Lecture Notes: Multilevel Analysis/Hierarchical Linear Modeling

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Multilevel Analysis/Hierarchical Linear Modeling

Edpsych/Psych/Stat 587

C.J. Anderson

Fall 2013

http://courses.education.illinois.edu/edpsy587/

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Hierarchical Multilevel Modeling Analysis Analysi

谢谢分享！

1. proc sort;
   
2. by case diagnose treat time outcome;
   
3. 
   
4. data tmp1 tmp2 tmp3;
   
5. set depress;
   
6. if time=0 then output tmp1;
   
7. else if time=1 then output tmp2;
   
8. else if time=2 then output tmp3;
   
9. data tmp1a;
   
10. set tmp1;
    
11. if outcome=0 then Y1='A';
    
12. else if outcome=1 then Y1='N';
    
13. drop outcome;
    
14. run;
    
15. proc sort data=tmp1a; by case ;
    
16. data tmp2a;
    
17. set tmp2;
    
18. if outcome=0 then Y2='A';
    
19. else if outcome=1 then Y2='N';
    
20. drop outcome;
    
21. proc sort data=tmp1a; by case ;
    
22. data tmp3a;
    
23. set tmp3;
    
24. if outcome=0 then Y3='A';
    
25. else if outcome=1 then Y3='N';
    
26. drop outcome;
    
27. proc sort data=tmp1a; by case ;
    
28. 
    
29. data Table11_2;
    
30. merge tmp1a tmp2a tmp3a;
    
31. by case ;
    
32. /*
    
33. if Y1='N' then do;
    
34.      if Y2='N' then do;
    
35.        if Y3='N' then Y='NNN';
    
36.            else Y='NNA';
    
37.          end;
    
38.          else do;
    
39.            if Y3='N' then Y='NAN';
    
40.            else Y='NAA';
    
41.          end;
    
42. end;
    
43. else do;
    
44.       if Y2='N' then do;
    
45.        if Y3='N' then Y='ANN';
    
46.            else Y='ANA';
    
47.          end;
    
48.          else do;
    
49.            if Y3='N' then Y='AAN';
    
50.            else Y='AAA';
    
51.          end;
    
52. end;*/
    
53. run;
    
54. 
    
55. 
    
56. 
    
57. proc freq ;
    
58. tables diagnose*treat*Y1*y2*y3/ nopercent norow nocol list sparse;
    
59. run;
    
60. 
    
61. proc genmod descending data=depress;  
    
62.   model outcome = diagnose treat time treat*time / dist=bin link=logit type3;
    
63.   title 'MLE ignoring repeated aspect of the data';
    
64. run;
    
65. 
    
66. proc logistic descending data=depress;
    
67. model outcome = diagnose treat time treat*time / Lackfit ;
    
68. run;
    
69. 
    
70. proc genmod descending data=depress;  
    
71.   class case;
    
72.   model outcome = diagnose treat time treat*time / dist=bin link=logit type3;
    
73.   repeated subject=case / type=exch corrw corrb;
    
74.   title 'GEE with Exchangable';
    
75. run;
    
76. 
    
77. 
    
78. proc genmod descending data=depress;  
    
79.   class case;
    
80.   model outcome = diagnose treat time treat*time / dist=bin link=logit type3;
    
81.   repeated subject=case / type=AR(1) corrw;
    
82.   title 'GEE with AR(1)';
    
83. run;
    
84. 
    
85. 
    
86. proc genmod descending data=depress;  
    
87.   class case;
    
88.   model outcome = diagnose treat time treat*time / dist=bin link=logit type3;
    
89.   repeated subject=case / type=unstr corrw;
    
90.   title 'GEE with Unstructured';
    
91. run;
    
92. 
    
93. 
    
94. proc nlmixed  data=depress qpoints=20;
    
95.   parms alpha=-.03 beta1=-1.3 beta2=-.06 beta3=.48 beta4=1.02 sigma=.01;   
    
96.   eta = alpha + beta1*diagnose + beta2*treat + beta3*time + beta4*treat*time + u;
    
97.   p = exp(eta)/(1 + exp(eta));
    
98.   model outcome ~ binary(p);
    
99.   random u ~ normal(0, sigma*sigma) subject = case;
    
100.   estimate 'var(U_o)' sigma**2;
     
101.   title 'Results from nlmixed';
     
102. run;
     
103. 
     
104. proc nlmixed  data=depress qpoints=20;
     
105.   parms alpha=-.03 beta1=-1.3 beta2=-.06 beta3=.48 beta4=1.02 ;   
     
106.   eta = alpha + beta1*diagnose + beta2*treat + beta3*time + beta4*treat*time ;
     
107.   p = exp(eta)/(1 + exp(eta));
     
108.   model outcome ~ binary(p);
     
109.   title 'Results from nlmixed; w/o random';
     
110. run;
     
111. 
     
112. proc glimmix data=depress gradient method=RSPL;
     
113. class case ;
     
114. model outcome = diagnose treat time treat*time /dist=bin link=logit solution;
     
115. random intercept / subject=case;
     
116. title 'How about GLIMMIX?';
     
117. run;
     
118. 
     
119. proc glimmix data=depress;
     
120. class case ;
     
121. model outcome = diagnose treat time treat*time /dist=bin link=logit solution;
     
122. title 'How about GLIMMIX? without random';
     
123. run;

复制代码

Longitudinal_depression.zip (3.4 KB) 本附件包括：

• Longitudinal_depression.sas

2016-7-2 21:46:44 上传

1. /* Simple model */
   
2. 
   
3. /* Ignoring within child dependency */
   
4. proc genmod data=respire descending;
   
5. model resp = age /link=logit dist=bin type3 ;
   
6. output out=preds pred=fitted;
   
7. title 'Simple Model ignoring repeated measures';
   
8. run;
   
9. 
   
10. /* Simple Model:  GEE */
    
11. proc genmod data=respire descending;
    
12. class id ;
    
13. model resp = age /link=logit dist=bin type3 ;
    
14. repeated subject=id  / type=exch;
    
15. title 'Simple Model GEE (exchangable)';
    
16. run;
    
17. 
    
18. /* Compute fitted probabilities to graph*/
    
19. data preds;
    
20. lambda=-2.3436 ;
    
21. lamAge=-0.0248;
    
22. do age = -40 to 50 by .5;
    
23.   psimple = exp(lambda + lamAge*Age)/(1+ exp(lambda + lamAge));
    
24.   pnoinf = 1-psimple;
    
25.   id=1;
    
26.   output;
    
27. end;
    
28. run;
    
29. 
    
30. goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate
    
31.          hpos=0 vpos=0 htext=2.5 ftext=swiss ctext= target= gaccess= gsfmode= ;
    
32. goptions device=WIN  ctext=blue
    
33.          graphrc interpol=join;
    
34. axis1   color=blue   width=2.0  label=('Age in Months')  ;
    
35. axis2   color=blue   width=2.0  label=(angle=90 'Fitted Probability of Infection')
    
36.   order=0 to 1 by .1;
    
37. axis3   color=blue   width=2.0  label=(angle=90 'Fitted Probability of NO Infection')
    
38.   order=0 to 1 by .1;
    
39. proc gplot data=WORK.PREDS;
    
40.    plot psimple * age / haxis=axis1 vaxis=axis2 frame ;
    
41.    title 'Probability of Infection';
    
42. run;
    
43. 
    
44. proc gplot data=WORK.PREDS;
    
45.    plot pnoinf * age / haxis=axis1 vaxis=axis2 frame ;
    
46.    title 'Probability of NO Infection';
    
47. run;
    
48. 
    
49. /* Complex model */
    
50. /* Ignoring within child dependency */
    
51. proc genmod data=respire descending;
    
52. class xero(descending) female(descending) stunted(descending);
    
53. model resp = age xero female cosine sine height stunted /link=logit dist=bin type3 ;
    
54. title 'Complex Model ignoring repeated measures';
    
55. run;
    
56. 
    
57. 
    
58. 
    
59. 
    
60. 
    
61. 
    
62. /* Complex Model: GEE */
    
63. proc genmod data=respire descending;
    
64. class id xero(descending) female(descending) stunted(descending);
    
65. model resp = age xero female cosine sine height stunted /link=logit dist=bin type3 ;
    
66. repeated subject=id  / type=exch corrw;
    
67. title 'Complex Model GEE (exchangable)';
    
68. run;
    
69. 
    
70. /* Complex Model: GEE */
    
71. proc genmod data=respire descending;
    
72. class id xero(descending) female(descending) stunted(descending);
    
73. model resp = age xero female cosine sine height stunted /link=logit dist=bin type3 ;
    
74. repeated subject=id  / type=unstr corrw;
    
75. title 'Complex Model GEE (unstructured)';
    
76. run;
    
77. 
    
78. 
    
79. /* Subject specific model */
    
80. /* Simple */
    
81. 
    
82. 
    
83. proc nlmixed data=respire qpoints=10;   *<- I am starting in a good place so 10 should do it;
    
84.   parms lam=-2.3 bAge=0.02 tau0= 0.8;  *<- starting values from GEE;
    
85.   eta = lam + bAge*Age + u;             *<- linear predictor;
    
86.   p = exp(eta)/(1 + exp(eta));          *<- probabiliy;
    
87.   model resp ~ binary(p);
    
88.   random u ~ normal(0, tau0*tau0) subject=id out=ebuR;
    
89.   estimate 'Var(Uo)' tau0**2;
    
90.   estimate 'odds ratio' exp(bAge);
    
91.   title 'Random Intercept, Simple Model';
    
92. run;
    
93. 
    
94. proc glimmix data=respire method=mmpl gradient noclprint;
    
95. class id ;
    
96. model resp = age / solution link=logit dist=bin;
    
97. random intercept / subject=id;
    
98. title 'Random Intercept, Simple Model: Default estimation method';
    
99. run;
    
100. 
     
101. proc glimmix data=respire  method=rmpl gradient noclprint;
     
102. class id ;
     
103. model resp = age / solution link=logit dist=bin;
     
104. random intercept / subject=id;
     
105. title 'Random Intercept, Simple Model: RMPL';
     
106. run;
     
107. 
     
108. 
     
109. proc glimmix data=respire  method=mspl gradient noclprint;
     
110. class id ;
     
111. model resp = age / solution link=logit dist=bin;
     
112. random intercept / subject=id;
     
113. title 'Random Intercept, Simple Model: mSPL';
     
114. run;
     
115. 
     
116. 
     
117. 
     
118. proc glimmix data=respire  method=rspl gradient noclprint;
     
119. class id ;
     
120. model resp = age / solution link=logit dist=bin;
     
121. random intercept / subject=id;
     
122. title 'Random Intercept, Simple Model: mSPL';
     
123. run;
     
124. 
     
125. proc glimmix data=respire method=laplace gradient noclprint;
     
126. class id ;
     
127. model resp = age / solution link=logit dist=bin;
     
128. random intercept / subject=id;
     
129. title 'Random Intercept, Simple Model: laPlace';
     
130. run;
     
131. 
     
132. proc glimmix data=respire method=quad gradient noclprint;
     
133. class id ;
     
134. model resp = age / solution link=logit dist=bin;
     
135. random intercept / subject=id;
     
136. title 'Random Intercept, Simple Model:  QUAD';
     
137. run;
     
138. 
     
139. data random;
     
140. set ebuR;
     
141. u = estimate;
     
142. lam=-2.6130;
     
143. bAge=-0.02676;
     
144. do age=-40 to 50 by .5;
     
145.    pperson = exp(lam + bAge*age + u)/(1+exp(lam + bAge*age + u));
     
146.    output;
     
147. end;
     
148. run;
     
149. 
     
150. 
     
151. goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate
     
152.          hpos=0 vpos=0 htext=2.5 ftext=swiss ctext= target= gaccess= gsfmode= ;
     
153. goptions device=WIN  ctext=blue
     
154.          graphrc interpol=join;
     
155. axis1   color=blue   width=2.0  label=('Age in Months')  ;
     
156. axis2   color=blue   width=2.0  label=(angle=90 'Fitted Probability of Infection')
     
157.   order=0 to .5 by .1;
     
158. axis3   color=blue   width=2.0  label=(angle=90 'Fitted Probability of NO Infection')
     
159.   order=0 to .5 by .1;
     
160. legend1 position=(top left inside) label=none frame cshadow=pink
     
161.         value=('Ignoring dependency' 'GEE (exchangable)');
     
162. proc gplot data=random;
     
163.    plot pperson*age=id / overlay haxis=axis1 vaxis=axis2 frame nolegend;
     
164.    title 'Probability of Infection';
     
165. run;
     
166. 
     
167. 
     
168. data predtwo;
     
169. /* From nlmixed */
     
170. lam=-2.6130;
     
171. bAge=-0.02676;
     
172. /* From gendmode using GEE */
     
173. lambda= -2.3355;
     
174. lamAge= -0.0243;
     
175. do age= -40 to 50 by .5;
     
176.   pGEE = exp(lambda + lamAge*Age)/(1+exp(lambda+lamAge*age));
     
177.   pperson = exp(lam + bAge*age)/(1+exp(lam + bAge*age ));
     
178.   pnoGEE= 1-pGEE;
     
179.   pnoperson = 1 -pperson;
     
180.   id=1;
     
181.   output;
     
182. end;
     
183. run;
     
184. 
     
185. proc sort data=preds; by age;
     
186. proc sort data=predGEE; by age;
     
187. 
     
188. data tmp;
     
189. merge preds predtwo;
     
190. by age;
     
191. pnosimple=1-psimple;
     
192. run;
     
193. 
     
194. 
     
195. goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate
     
196.          hpos=0 vpos=0 htext=2.5 ftext=swiss ctext= target= gaccess= gsfmode= ;
     
197. goptions device=WIN  ctext=blue
     
198.          graphrc interpol=join;
     
199. axis1   color=blue   width=2.0  label=('Age in Months')  ;
     
200. axis2   color=blue   width=2.0  label=(angle=90 'Fitted Probability of Infection')
     
201.   order=0 to 1 by .1;
     
202. axis3   color=blue   width=2.0  label=(angle=90 'Fitted Probability of NO Infection')
     
203.   order=0 to 1 by .1;
     
204. legend1 position=(top left inside) label=none frame cshadow=pink
     
205.         value=('Ignoring dependency' 'GEE (exchangable)' 'Random effects');
     
206. legend2 position=(bottom left inside) label=none frame cshadow=pink
     
207.         value=('Ignoring dependency' 'GEE (exchangable)' 'Random effects');
     
208. 
     
209. proc gplot data=tmp;
     
210.    plot psimple*age pGEE*age pperson*age/ overlay haxis=axis1 vaxis=axis2 frame legend=legend1;
     
211.    title 'Probability of Infection';
     
212. run;
     
213. 
     
214. proc gplot data=tmp;
     
215.    plot pnosimple*age pnoGEE*age pnoperson*age/ overlay haxis=axis1 vaxis=axis3 frame legend=legend2;
     
216.    title 'Probability of No Infection';
     
217. run;
     
218. 
     
219. 
     
220. 
     
221.    /* Complex  */
     
222. 
     
223. 
     
224. proc nlmixed data=respire method=gauss qpoints=10;   
     
225.   parms lam=-2.42 bAge=-.03 bxero=0.6 bfemale=-.42 bcos=-.6 bsin=-.2 bheight=-.05
     
226.         bstunt=0.15 sigma= 0.1;  
     
227.   eta = lam + bAge*Age + bxero*xero + bfemale*female + bcos*cosine + bsin*sine
     
228.         + bheight*height + bstunt*stunted +  u;            
     
229.   p = exp(eta)/(1 + exp(eta));         
     
230.   model resp ~ binary(p);
     
231.   random u ~ normal(0, sigma*sigma) subject=id;
     
232.   title 'Random Intercept, Complex Model';
     
233.   estimate 'exp(bAge)' exp(bAge);
     
234.   estimate 'exp(bxero)' exp(bxero);
     
235.   estimate 'exp(bfemale)' exp(bfemale);
     
236.   estimate 'exp(bcosine)' exp(bcos);
     
237.   estimate 'exp(bheight)' exp(eheight);
     
238.   estimate 'exp(bstunt)'  exp(bstunt);
     
239.   estimate 'var(u)'  sigma**2;
     
240. 
     
241. run;
     
242. 
     
243. 
     
244. proc glimmix data=respire method=quad gradient noclprint;
     
245. class id;
     
246. model resp = age xero female cosine sine height stunted / link=logit dis=binary solution;
     
247. random intercept / subject=id;
     
248. run;

复制代码

respirtory.zip (12.62 KB) 本附件包括：

• respirtory.sas

2016-7-2 21:50:10 上传

谢谢分享！