/* GLMM.PRG -- EM program for linear mixed models Author: Brian M. Steele Dept. of Math. Sci. University of Montana Missoula MT, 59812 Email: bsteele@selway.umt.edu The author gives permission for redistribution of this software for non- commercial purposes only. */ /* Instructions are in the file glmm.hlp */ library glmm; #include \gauss\src\glmm.dec; format /m1/rd 4,0; /**************** User-specified variables **********************/ output file = glmm.out reset; /* Control sequence for Schall's data - 2 normal random factors *********/ zflag = 2; {r,z,q,x,y,nvec,name,dname} = DATAS(zflag); respv = "B"; eff = "N"; s = {.2,.1}; s = {5,5}; test = 1E-5; lflag = 1; /************************************************************************/ /* Control sequence for Schall's data - 1 normal random factor ********* zflag = 1; {r,z,q,x,y,nvec,name,dname} = DATAS(zflag); respv = "B"; eff = "N"; s = {.2}; test = 1E-5; lflag = 1; ************************************************************************/ /**** Control sequence for Schall's data - 2 log-gamma random factors *** zflag = 2; {r,z,q,x,y,nvec,name,dname} = DATAS(zflag); respv = "B"; eff = "LG"; s = 5 | 5; sn = 1; test = 1E-5; lflag = 1; ************************************************************************/ /****** Control sequence for Salamander data ************************* zflag = 2; @ zflag = 1, 2, or 3 @ {r,z,q,x,y,nvec,name,dname} = DATASAL(zflag); respv = "B"; eff = "N"; s = {1,1}; test = 1E-5; lflag = 1; ************************************************************************/ /****** Control sequence for pooled Salamander data ************************* {r,z,q,x,y,nvec,name,dname} = JSAL; respv = "B"; eff = "N"; s = {1,1}; test = 1E-5; lflag = 1; ************************************************************************/ /* Control sequence for Thall and Vail's data - 1 normal factor ******** zflag = 1; respv = "P"; eff = "N"; {r,z,q,x,y,name,dname} = DATAT(zflag); s = {.8}; @ initial estimate @ test = 1E-5; lflag = 1; ***********************************************************************/ /* Control sequence for Thall and Vail's data - 2 normal factor ******** zflag = 2; respv = "P"; eff = "N"; {r,z,q,x,y,name,dname} = DATAT(zflag); s = {.8,.8}; @ initial estimate @ test = 1E-5; lflag = 1; ***********************************************************************/ /* Control sequence for Thall and Vail's data - 1 log-gamma factor ** zflag = 1; respv = "P"; eff = "LG"; {r,z,q,x,y,name,dname} = DATAT(zflag); s = 4; @ initial estimate @ sn = 1; @ sign parameter for log-gamma (controls direction of skew) @ test = 1E-5; lflag = 1; ***********************************************************************/ /* Control sequence for Thall and Vail's data - 2 log-gamma factors ** zflag = 2; respv = "P"; eff = "LG"; {r,z,q,x,y,name,dname} = DATAT(zflag); s = 4|4; @ initial estimate @ sn = 1; @ sign parameter for log-gamma (controls direction of skew) @ test = 1E-5; lflag = 1; ***********************************************************************/ /*********************************** Program ***************************/ a = MEMINIT(y,x,respv,sn); if eff $== "N"; {a,s,b} = MEMN(a,s,r,z,q,x,y,nvec,respv,eff,s,test,lflag); {se,ion} = IOMEMN(a,b,s,respv,x,z,q,names,dname); else; {a,s,b} = MEMLG(a,s,r,z,q,x,y,nvec,respv,eff,s,sn,test,lflag); {se,ion} = IOMEMLG(a,b,s,respv,x,z,q,names,dname); endif; /************************** Data procedures ***************************/ PROC(8) = DATAS(zflag); @ ref: Schall, R. 1991. Estimation in generalized linear models with random effects. Biometrika, 78, 719-727. @ dat = { 1 1 178, 1 2 193, 1 3 217, 2 4 109, 2 5 112, 2 6 115, 3 7 66, 3 8 75, 3 9 80, 4 10 118, 4 11 125, 4 12 137, 5 13 123, 5 14 146, 5 15 170, 6 16 115, 6 17 130, 6 18 133, 7 19 200, 7 20 189, 7 21 173, 8 22 88, 8 23 76, 8 24 90, 9 25 121, 9 26 124, 9 27 136}; dname = "Schall's data"; n = rows(dat); y = dat[.,3]; occ = dat[.,1]; dish = dat[.,2]; let name = const; nvec = 400*ones(n,1); x = ones(n,1); s = seqa(1,1,n/2); s = seqa(1,1,9); z = (occ.==s'); r = 1; q = 9; if zflag == 2; s = seqa(1,1,27); z = z~(dish.==s'); r = 2; q = q|27; endif; output on; $dname; format 8,8;" Variables = ";;$name'; format 4,0;" N random factors = ";;r; " N random effects = ";;q';?; output off; RETP(r,z,q,x,y,nvec,name,dname); ENDP; PROC(7) = DATAT(zflag); @ ref: Thall, P.F. and Vail, S.C. 1990. Some covariance models for longitudinal count data with overdispersion. Biometrics, 46, 657-671.@ dname = "Thall's data"; dat = { 104 5 3 3 3 0 11 31, 106 3 5 3 3 0 11 30, 107 2 4 0 5 0 6 25, 114 4 4 1 4 0 8 36, 116 7 18 9 21 0 66 22, 118 5 2 8 7 0 27 29, 123 6 4 0 2 0 12 31, 126 40 20 23 12 0 52 42, 130 5 6 6 5 0 23 37, 135 14 13 6 0 0 10 28, 141 26 12 6 22 0 52 36, 145 12 6 8 4 0 33 24, 201 4 4 6 2 0 18 23, 202 7 9 12 14 0 42 36, 205 16 24 10 9 0 87 26, 206 11 0 0 5 0 50 26, 210 0 0 3 3 0 18 28, 213 37 29 28 29 0 111 31, 215 3 5 2 5 0 18 32, 217 3 0 6 7 0 20 21, 219 3 4 3 4 0 12 29, 220 3 4 3 4 0 9 21, 222 2 3 3 5 0 17 32, 226 8 12 2 8 0 28 25, 227 18 24 76 25 0 55 30, 230 2 1 2 1 0 9 40, 234 3 1 4 2 0 10 19, 238 13 15 13 12 0 47 22, 101 11 14 9 8 1 76 18, 102 8 7 9 4 1 38 32, 103 0 4 3 0 1 19 20, 108 3 6 1 3 1 10 30, 110 2 6 7 4 1 19 18, 111 4 3 1 3 1 24 24, 112 22 17 19 16 1 31 30, 113 5 4 7 4 1 14 35, 117 2 4 0 4 1 11 27, 121 3 7 7 7 1 67 20, 122 4 18 2 5 1 41 22, 124 2 1 1 0 1 7 28, 128 0 2 4 0 1 22 23, 129 5 4 0 3 1 13 40, 137 11 14 25 15 1 46 33, 139 10 5 3 8 1 36 21, 143 19 7 6 7 1 38 35, 147 1 1 2 3 1 7 25, 203 6 10 8 8 1 36 26, 204 2 1 0 0 1 11 25, 207 102 65 72 63 1 151 22, 208 4 3 2 4 1 22 32, 209 8 6 5 7 1 41 25, 211 1 3 1 5 1 32 35, 214 18 11 28 13 1 56 21, 218 6 3 4 0 1 24 41, 221 3 5 4 3 1 16 32, 225 1 23 19 8 1 22 26, 228 2 3 0 1 1 25 21, 232 0 0 0 0 1 13 36, 236 1 4 3 2 1 12 37}; let names = ID Y1 Y2 Y3 Y4 TRT BASE AGE; makevars(dat,0,names); dat[.,7] = ln(dat[.,7]/4); dat[.,8] = ln(dat[.,8]); dat = dat~(dat[.,6].*dat[.,7]); x = {}; i = 2; do until i == 5; x = x|(dat[.,1 i 6:9]~zeros(rows(dat),1)); i = i + 1; endo; x = x|(dat[.,1 i 6:9]~(ones(rows(dat),1))); let namen = id count trt base age base_trt visit4; makevars(x,0,namen); n = rows(count); y = count; i = ones(4,1); ii = eye(59); z = i.*.ii; q = cols(z); if zflag == 1; x = ones(n,1)~base~trt~base_trt~age~visit4; let name = const base trt base_trt age visit4; endif; if zflag == 2; t = -0.3|-0.1|0.1|0.3; z = z~((i.*t).*.ii); q = q|cols(ii); i = eye(4); ii = ones(59,1); visit10 = (i.*.ii)*t; x = ones(n,1)~base~trt~base_trt~age~visit10; let name = const base trt base_trt age visit10; endif; r = rows(q); n = rows(y); output on; $dname; format 4,0;" N random factors = ";;r; " N random effects = ";;q';?; output off; RETP(r,z,q,x,y,name,dname); ENDP; PROC(8) = DATASAL(zflag); @ ref: McCullagh, P. and Nelder, J.A. 1989. Generalized linear models, 2nd Ed. Chapman and Hall. Tables 14.4 - 14.6. @ male = { 1 4 5 1 4 5, 5 5 3 3 1 2, 2 1 1 4 3 3, 4 2 2 5 5 4, 3 3 4 2 2 1, 9 9 10 7 6 8, 8 8 9 9 7 6, 6 6 7 10 10 9, 10 7 8 6 9 10, 7 10 6 8 8 7, 9 9 7 10 10 8, 7 6 9 7 6 10, 8 7 6 9 7 6, 10 10 8 8 9 9, 6 8 10 6 8 7, 5 2 3 4 2 1, 4 1 5 2 1 5, 1 4 2 5 5 3, 3 3 1 1 4 4, 2 5 4 3 3 2}; y1 = { 1 1 1 0 1 1, 1 1 1 1 1 1, 1 0 1 1 1 1, 1 1 1 0 1 1, 1 1 1 1 1 1, 1 1 1 0 1 1, 0 0 0 1 0 0, 0 1 0 0 1 1, 0 0 1 1 1 1, 0 0 1 0 1 0, 0 1 1 1 0 1, 0 0 0 1 0 0, 0 0 0 0 0 1, 0 1 1 1 0 1, 0 1 0 0 0 0, 0 0 1 0 0 0, 1 1 1 0 1 1, 1 0 1 0 1 0, 1 1 1 1 1 0, 1 0 0 1 1 0}; y2 = { 1 0 1 0 1 0, 1 1 1 0 0 1, 0 1 1 1 1 1, 1 0 0 0 0 0, 1 0 1 0 0 0, 1 1 1 0 0 0, 1 1 0 1 0 1, 1 0 0 0 1 0, 1 1 1 1 0 1, 0 0 1 1 1 0, 0 1 0 0 0 1, 0 1 0 0 0 1, 0 1 0 1 0 1, 0 0 1 1 0 1, 0 0 0 0 0 0, 1 0 1 1 1 1, 1 0 1 0 0 0, 0 0 0 0 0 0, 1 1 1 1 1 1, 1 0 1 1 1 0}; y3 = { 1 1 1 0 1 1, 0 0 0 1 0 0, 1 1 1 0 1 1, 1 0 1 1 0 0, 0 1 1 0 1 0, 0 1 1 1 1 1, 0 1 0 0 0 0, 1 1 1 1 1 0, 0 0 1 0 1 1, 0 1 1 1 1 1, 0 0 0 1 0 1, 0 1 0 1 0 1, 1 0 1 0 0 1, 1 1 1 0 0 0, 0 1 1 1 0 1, 0 0 1 0 0 0, 1 0 1 0 0 0, 0 0 1 0 1 0, 1 0 0 0 1 0, 0 0 1 0 1 0}; if zflag == 1; y = y1; dname = "Summer 86"; endif; if zflag == 2; y = y2; dname = "Fall 86 (re-runs)"; endif; if zflag == 3; y = y3; dname = "Fall 86 "; endif; z = zeros(5,1); o = ones(5,1); malesp = (z~o~z~o~z~o)|(o~z~o~z~o~z)|(z~o~z~o~z~o)|(o~z~o~z~o~z); s = seqa(1,1,10); female = (s|s).*ones(20,6); femalesp = zeros(10,6)|ones(10,6); y = vec(y); malesp = vec(malesp); femalesp = vec(femalesp); x = ones(120,1)~femalesp~malesp~femalesp.*malesp; s = seqa(1,1,20); z = ((vec(female)+10*femalesp).==s')~((vec(male)+10*malesp).==s'); q = rows(s)|rows(s); r = 2; nvec = ones(120,1); let name = const fws mws fws_mws; output on; ?;$dname; format 4,0;" N random factors = ";;r; " N random effects = ";;q';?; output off; RETP(r,z,q,x,y,nvec,name,dname); ENDP; PROC(8) = JSAL; {r,z,q,x,y,nvec,name,dname} = DATASAL(1); n = rows(y); ym = y; xm = x; zm = zeros(360,120); zm[1:120,1:20] = z[.,1:20]; zm[1:120,61:80] = z[.,21:40]; {r,z,q,x,y,nvec,name,dname} = DATASAL(2); ym = ym | y; xm = xm | x; zm[121:240,21:40] = z[.,1:20]; zm[121:240,81:100] = z[.,21:40]; {r,z,q,x,y,nvec,name,dname} = DATASAL(3); ym = ym | y; xm = xm | x; zm[241:360,41:60] = z[.,1:20]; zm[241:360,101:120] = z[.,21:40]; y = ym; x = xm; z = zm; nvec = nvec|nvec|nvec; q = 60|60; cls; dname = "Pooled salamander data"; RETP(r,z,q,x,y,nvec,name,dname); ENDP; END;