Published on :
Statistics CREATION_INTERNE

Advanced Statistical Analysis: GLIMMIX and NLMIXED

This code is also available in: Deutsch Español Français
Awaiting validation
This script first generates a simulated dataset ('hair'). It then performs a series of analyses using PROC GLIMMIX to compare linear and logistic models with different types of random effects (G-side vs R-side). Finally, it aggregates the data to fit a complex model via PROC NLMIXED. Technical note: A variable inconsistency ('id' vs 'subj_id') is present in the sorting and aggregation section.
Data Analysis

Type : CREATION_INTERNE


All data is generated by the first DATA STEP using random functions (normal, uniform).

1 Code Block
DATA STEP Data
Explanation :
Generation of the 'hair' working dataset containing simulated repeated measures (1000 subjects, 8 weeks).
Copied!
1DATA hair;
2 seed = 1999;
3 beta = 1.386294; *--- P is approximately 0.8;
4 DO Subj_id=1 to 1000;
5 s1 = 2*normal(seed);
6 s2 = sqrt(1.5)*normal(seed);
7 DO week=1 to 8;
8 score = round(s1 + normal(seed), 0.01);
9 p = 1/(1+exp(-beta - s2));
10 y = 0;
11 u = uniform(seed);
12 IF u < p THEN y = 1;
13 OUTPUT;
14 END;
15 END;
16keep subj_id score y;
17RUN;
2 Code Block
PROC GLIMMIX
Explanation :
Linear mixed model with random intercept per subject (G-side).
Copied!
1ods html;
2title "*** Linear Mixed Model with G-side Random Effects ***";
3 PROC GLIMMIX DATA=hair;
4 class subj_id;
5 model score = / s;
6 random int / subject=subj_id;
7RUN;
3 Code Block
PROC GLIMMIX
Explanation :
Linear mixed model modeling the residual covariance structure (R-side, Compound Symmetry type).
Copied!
1title "*** Linear Mixed Model with R-side Random Effects ***";
2PROC GLIMMIX DATA=hair;
3 class subj_id;
4 model score = / s;
5 random _residual_ / subject=subj_id type=cs;
6RUN;
4 Code Block
PROC GLIMMIX
Explanation :
Logistic mixed model with G-side random effects.
Copied!
1title "*** Logistic Model with G-side Random Effects ***";
2 PROC GLIMMIX DATA=hair;
3 class subj_id;
4 model y (ref='0') = / dist=binary link=logit s;
5 random int / subject=subj_id;
6 estimate 'P' int 1 / ilink;
7RUN;
5 Code Block
PROC GLIMMIX
Explanation :
Logistic mixed model with R-side random effects.
Copied!
1title "*** Logistic Model with R-side Random Effects ***";
2PROC GLIMMIX DATA=hair;
3 class subj_id;
4 model y (ref='0')= / dist=binary link=logit s;
5 random _residual_ / subject=subj_id type=cs;
6 estimate 'P' int 1 / ilink;
7RUN;
8ods html close;
6 Code Block
PROC SORT
Explanation :
Sorting data for aggregation. (Warning: the 'id' variable seems incorrect compared to the created 'subj_id').
Copied!
1PROC SORT DATA=hair;
2 BY id;
3RUN;
7 Code Block
PROC MEANS Data
Explanation :
Data aggregation: sum of the response variable 'y' by identifier.
Copied!
1PROC MEANS DATA=hair sum noprint;
2 BY id;
3 var y;
4 OUTPUT out=sum_hair sum=t;
5RUN;
8 Code Block
DATA STEP Data
Explanation :
Preparation of aggregated data for PROC NLMIXED (frequency calculation and renaming).
Copied!
1DATA sum_hair;
2 SET sum_hair;
3 m_t = _freq_ - t;
4 rename _freq_ = m;
5 drop _type_;
6RUN;
9 Code Block
PROC NLMIXED
Explanation :
Non-linear mixed modeling defining a custom log-likelihood function (ll).
Copied!
1PROC NLMIXED DATA=sum_hair;
2 parms beta=0 rho=0.1;
3 pi = 1/(1+exp(-beta));
4 pic = 1 - pi;
5 p1 = ( 1 - rho )*pi + rho;
6 p1c = 1 - p1;
7 p2 = p1 - rho;
8 p2c = 1 - p2;
9 z = lgamma(m+1) - lgamma(t+1) - lgamma(m_t+1);
10 ll = z + log( pi*p1**t*p1c**m_t + pic*p2**t*p2c**m_t );
11 model t ~ general( ll );
12 estimate 'P' 1/(1+exp(-beta));
13 estimate 'Rho*Rho' rho*Rho;
14RUN;
This material is provided "as is" by We Are Cas. There are no warranties, expressed or implied, as to merchantability or fitness for a particular purpose regarding the materials or code contained herein. We Are Cas is not responsible for errors in this material as it now exists or will exist, nor does We Are Cas provide technical support for it.