/*--------------------------------------------------------------------*/ /* */ /* This macro is based on the article "Response Transformations */ /* in Repeated Measures and Growth Curve Models," by R. J. Boik, */ /* C. Todd, and S. K. Hyde (Communications in Statistics---Theory */ /* and Methods, 2000, 29, 669-773). */ /* */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /* This iml program computes the MLE of the Box-Cox or folded-power */ /* transformation parameter, lambda, in a repeated measures or in a */ /* growth curve model. The user must customize the program for her/ */ /* his particular application. A maximum of three within-subject */ /* factors are allowed. */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /*-------------------Required Arguments-------------------------------*/ /*--------------------------------------------------------------------*/ /* */ /* response = var, where var is a variable or a list of variables. */ /* */ /* data = dataset, where dataset is the name of the sas data set that */ /* contains the data. */ /* */ /*--------------------------------------------------------------------*/ /*-------------------Optional Arguments-------------------------------*/ /*--------------------------------------------------------------------*/ /* */ /* class = var, where var a list of classification variables. */ /* */ /* explanatory = var, where var is a variable or a list */ /* of variables. If explanatory is omitted, */ /* then explanatory = ; is assumed. */ /* */ /* family = value, where value is either boxcox or foldpow */ /* If family is omitted, then boxcox is assumed. */ /* */ /* plot_like = value, where value is either yes or no. If there are */ /* any repeated measures, then the plot option is not available.*/ /* */ /* plot_title = value, where value is the plot title that is desired. */ /* */ /* levels = values, where values is a row vector that specifies the */ /* number of levels of the within subject factors. The levels */ /* should be specified in the same order as on the REPEATED */ /* statement in PROC GLM. For example, if there are two */ /* within-subject factors with 5 and 4 levels, then levels */ /* is specified as "levels = 4 5" or as "levels = 5 4" */ /* depending on which index varies fastest. If the linear */ /* model has no repeated measures, then levels can be omitted. */ /* */ /* Polydeg: This is a matrix that specifies the polynomial degree of */ /* the growth curve that is desired for each within-subject */ /* factor. To obtain a repeated measures model, a saturated */ /* polynomial curve should be specified. That is, the */ /* polynomial degree for each factor should be one fewer that */ /* the number of levels of the factor. For example, if */ /* "levels = 5 4" then a repeated measures model is obtained */ /* by setting "Polydeg = 4 3." If the linear model has no */ /* repeated measures, then Polydeg can be omitted. /* */ /* To obtain a growth curve model, a lower order polynomial */ /* curve should be specified for one or more within-subject */ /* factor. Multiple models can be requested by specifying */ /* more than one set of polynomials. For example, if */ /* levels = 4 5, then "Polydeg = 3 4 3 3 2 3" specifies the */ /* polynomial degrees for three growth curve models. In the */ /* first model, Factor 1 is modeled as a third degree */ /* polynomial and Factor 2 is modeled as a fourth degree */ /* polynomial. Note that this is a saturated model and is */ /* equivalent to a repeated measures model. In the second */ /* model, both Factors are modeled as third degree polynomials.*/ /* A polynomial curve with degree zero contains an intercept */ /* term only. */ /* */ /* time1: This is a row vector that specifies the time points at */ /* which observations were made on Factor 1. The arguments */ /* time2 and time3 are defined analogously. If the levels of */ /* Factor 1 represent equal time intervals, then time1 can be */ /* omitted. If the levels of Factor 1 represent time points */ /* 1, 3, 4, 6, then time1 can be specified as */ /* "time1 = 1 3 4 6" If there is only one within-subject */ /* factor, then time2 and time3 can be omitted. */ /* */ /* alpha: This is a scalar that represents the desired confidence */ /* coefficient. Typically, alpha is set to 0.05 or 0.01. */ /* The default value is 0.05. The value of confid must be in */ /* the interval (0,1). */ /* */ /* guesslam: This is a scalar that represents an initial guess for */ /* lambda. The default value 1 appears to work well, but for */ /* some data sets, the initial guess might need to be */ /* modified. */ /* */ /* converge: This is a scalar convergence criterion used when */ /* computing the MLE of the covariance matrix under local */ /* sphericity. Iteration stops when improvement in the */ /* evaluated log likelihood function is smaller than "converge.*/ /* The default value 1e-6 appears to work well and might not */ /* need to be modified. */ /* */ /* small: This is a scalar that is used to determine the rank of a */ /* matrix. The rank of a matrix is computed as the number of */ /* singular values that are larger than the maximum singular */ /* value multiplied by small. The default value 1.e-10 appears*/ /* to work well and might not need to be modified. */ /* */ %macro transform(response = , explanatory = , class = , alpha = 0.05, guesslam = 1, converge = 1.e-6, small = 1.e-10, data = , levels = 1, Polydeg = 0, time1 = 1, time2 = 1, time3 = 1, family = boxcox, plot_like = no); * plot_title = Log Likelihood Function,; /*------------------------------------------------------------*/ /*----------- PROC GLMMOD: Create Design Matrix --------------*/ /*------------------------------------------------------------*/ proc glmmod outparm=outp outdesign=outX noprint data = &data; class &class; model &response = &explanatory; run; /*-----------------------------------------------------------------*/ /*-- PROC IML: Compute MLE of Lambda & Perform Lack of Fit Tests --*/ /*-----------------------------------------------------------------*/ proc iml; nlength = %length(&levels); levels=J(1,nlength,0); %let j1=1; %do %while (%scan(&levels,%eval(&j1)) > 0); jj=%eval(&j1); levels[jj]=%scan(&levels,%eval(&j1)); %let j1 = %eval(&j1+1); %end; nlevels = %eval(&j1)-1; levels=shape(levels[1:nlevels],1,nlevels); nlength=%length(&Polydeg); vecPolydeg = J(1,nlength,0); %let j2=1; %do %while (%scan(&Polydeg,%eval(&j2)) > 0); jj=%eval(&j2); vecPolydeg[jj]=%scan(&Polydeg,%eval(&j2)); %let j2 = %eval(&j2+1); %end; npoly = %eval(&j2)-1; if npoly > 0 then Polydeg=shape(vecPolydeg[1:npoly],npoly/nlevels,nlevels); else Polydeg = {0}; nlength=%length(&time1); time1=J(1,nlength,0); %let j4=1; %do %while (%scan(&time1,%eval(&j4)) > 0); jj=%eval(&j4); time1[jj]=%scan(&time1,%eval(&j4)); %let j4 = %eval(&j4+1); %end; time1=time1[1:jj]; if ncol(time1) < levels[1] then time1=1:levels[1]; time1=shape(time1,1,levels[1]); if nlevels > 1 then do; nlength=%length(&time2); time2=J(1,nlength,0); %let j4=1; %do %while (%scan(&time2,%eval(&j4)) > 0); jj=%eval(&j4); time2[jj]=%scan(&time2,%eval(&j4)); %let j4 = %eval(&j4+1); %end; time2=time2[1:jj]; if ncol(time2) < levels[2] then time2=1:levels[2]; time2=shape(time2,1,levels[2]); end; else time2 = {0}; if nlevels > 2 then do; nlength=%length(&time3); time3=J(1,nlength,0); %let j4=1; %do %while (%scan(&time3,%eval(&j4)) > 0); jj=%eval(&j4); time3[jj]=%scan(&time3,%eval(&j4)); %let j4 = %eval(&j4+1); %end; time3=time3[1:jj]; if ncol(time3) < levels[3] then time3=1:levels[3]; time3=shape(time3,1,levels[3]); end; else time3 = {0}; confid = 1-%sysevalf(&alpha); boxcox = 'yes'; foldpow = 'no'; %if &family = foldpow %then %do; boxcox='no'; foldpow='yes'; %end; guesslam = %sysevalf(&guesslam); converge = %sysevalf(&converge); small = %sysevalf(&small); /*------------------------------------------------------------*/ /*------------End of User Defined Variables ------------------*/ /*------------------------------------------------------------*/ /*-------------------------------------------------*/ /*-------------- Define Sub-Modules ---------------*/ /*-------------------------------------------------*/ /*-------------------------------------------------*/ /* Module to stack columns of matrix into a vector */ /*-------------------------------------------------*/ start vec(mat); n = nrow(mat); p = ncol(mat); tot = n*p; *; * Caution shape(M,pq,1) = vec(M`).; *; Y = shape(mat`,tot,1); return(Y); finish vec; /* ----------------------------------------------------------*/ /* Module to unstack vector into a matrix with the specified */ /* number of rows. No checking is made if elements are lost */ /* or duplicated at the end */ /* ----------------------------------------------------------*/ start dvec(mat,nr); tot = nrow(mat); nc = tot/nr; *; * Caution shape(vecM,p,q) = dvec(vecM,p,q)`; *; Y = shape(mat,nc,nr)`; return(Y); finish dvec; /*-----------------------------------------*/ /* Module to compute the rank of a matrix */ /*-----------------------------------------*/ start rankM(A) global(small); m = nrow(A); n = ncol(A); call svd(U,S,V,A); if m = 1 then S = S[1]; else if m = 0 then S = {0}; tol = max(S) * small; r = max(loc(S>tol)); /*rank of A = number of nonzero singular values */ return(r); finish rankM; /*-------------------------------------------*/ /* Module to compute the commutation matrix. */ /*-------------------------------------------*/ start commute(p,q); Ip = I(p); Iq = I(q); k = shape({0},p*q,p*q); if p <= q then do i = 1 to p; k = k + (Ip[,i]`@Iq)@Ip[,i]; end; else do i = 1 to q; k = k + (Iq[,i]@Ip)@Iq[,i]`; end; return(k); finish commute; /*-------------------------------------------*/ /* Module to compute elementary matrices */ /*-------------------------------------------*/ start makeE(m,i); nr = max(ncol(m),nrow(m)); if nr = 1 then do; E = I(m[1]); return(E); end; if i = 1 then do; s2 = m[2:nr][+,+]; E = I(m[1]) // J(s2,m[1],0); return(E); end; if i = nr then do; s1 = m[1:nr-1][+,+]; E = J(s1,m[nr],0) // I(m[nr]); return(E); end; s1 = m[1:i-1][+,+]; s2 = m[i+1:nr][+,+]; E = J(s1,m[i],0) // I(m[i]) // J(s2,m[i],0); return(E); finish makeE; /*-------------------------------------------*/ /* Module to compute the C and G matrices */ /*-------------------------------------------*/ start makeCG(C,G,q,d,levels,w) global(time1,time2,time3); nfact = max(nrow(levels),ncol(levels)); if nfact = 1 then do; t1 = J(levels[1],1,1)/sqrt(levels[1]); if w[1] = 0 then do; G=t1; d=0; q=0; return; end; G = orpol(time1,w[1]); C = G[,2:w[1]+1]; q = w[1]; d = 1; return; end; else if nfact = 2 then do; t1 = J(levels[1],1,1)/sqrt(levels[1]); t2 = J(levels[2],1,1)/sqrt(levels[2]); if w[+] = 0 then goto AB00; if w[1]*w[2] > 0 then goto AB11; if w[1] > 0 then goto AB10; if w[2] > 0 then goto AB01; AB00: G=t1@t2; d=0; q=0; return; AB11: G1 = orpol(time1,w[1]); G2 = orpol(time2,w[2]); C1 = G1[,2:w[1]+1]; C2 = G2[,2:w[2]+1]; G = G1@G2; C= (C1@t2) || (t1@C2) || (C1@C2); q = w[1] || w[2] || w[1]*w[2]; d = 3; return; AB10: G1 = orpol(time1,w[1]); C1 = G1[,2:w[1]+1]; G=G1@t2; C=C1@t2; q = w[1]; d = 1; return; AB01: G2 = orpol(time2,w[2]); C2 = G2[,2:w[2]+1]; G=t1@G2; C=t1@C2; q = w[2]; d = 1; return; end; else if nfact = 3 then do; t1 = J(levels[1],1,1)/sqrt(levels[1]); t2 = J(levels[2],1,1)/sqrt(levels[2]); t3 = J(levels[3],1,1)/sqrt(levels[3]); if w[+] = 0 then goto ABC000; if w[1]*w[2]*w[3] > 0 then goto ABC111; if w[1]*w[2] > 0 then goto ABC110; if w[1]*w[3] > 0 then goto ABC101; if w[2]*w[3] > 0 then goto ABC011; if w[1] > 0 then goto ABC100; if w[2] > 0 then goto ABC010; if w[3] > 0 then goto ABC001; ABC000: G=(t1@t2)@t3; q=0; d=0; return; ABC111: G1 = orpol(time1,w[1]); G2 = orpol(time2,w[2]); G3 = orpol(time3,w[3]); C1 = G1[,2:w[1]+1]; C2 = G2[,2:w[2]+1]; C3 = G3[,2:w[3]+1]; G = (G1@G2)@G3; C = (C1@t2)@t3 || (t1@C2)@t3 || (t1@t2)@C3 || (C1@C2)@t3 || (C1@t2)@C3 || (t1@C2)@C3 || (C1@C2)@C3; q = w[1] || w[2] || w[3] || w[1]*w[2] || w[1]*w[3] || w[2]*w[3] || w[1]*w[2]*w[3]; d = 7; return; ABC110: G1 = orpol(time1,w[1]); G2 = orpol(time2,w[2]); C1 = G1[,2:w[1]+1]; C2 = G2[,2:w[2]+1]; G = (G1@G2)@t3; C = (C1@t2)@t3 || (t1@C2)@t3 || (C1@C2)@t3; q = w[1] || w[2] || w[1]*w[2]; d = 3; return; ABC101: G1 = orpol(time1,w[1]); G3 = orpol(time3,w[3]); C1 = G1[,2:w[1]+1]; C3 = G3[,2:w[3]+1]; G = (G1@t2)@G3; C = (C1@t2)@t3 || (t1@t2)@C3 || (C1@t2)@C3; q = w[1] || w[3] || w[1]*w[3]; d = 3; return; ABC011: G2 = orpol(time2,w[2]); G3 = orpol(time3,w[3]); C2 = G2[,2:w[2]+1]; C3 = G3[,2:w[3]+1]; G = (t1@G2)@G3; C = (t1@C2)@t3 || (t1@t2)@C3 || (t1@C2)@C3; q = w[2] || w[3] || w[2]*w[3]; d = 3; return; ABC100: G1 = orpol(time1,w[1]); C1 = G1[,2:w[1]+1]; G = (G1@t2)@t3; C = (C1@t2)@t3; q = w[1]; d = 1; return; ABC010: G2 = orpol(time2,w[2]); C2 = G2[,2:w[2]+1]; G = (t1@G2)@t3; C = (t1@C2)@t3; q = w[2]; d = 1; return; ABC001: G3 = orpol(time3,w[3]); C3 = G3[,2:w[3]+1]; G = (t1@t2)@G3; C = (t1@t2)@C3; q = w[3]; d = 1; return; end; finish makeCG; /*-----------------------------------------------------------------*/ /* Module to compute the mle of Omega under local sphericity */ /* The algorithm is based on Fisher scoring */ /*-----------------------------------------------------------------*/ start mlesigma(A,n,m,converge) global(sig_iter,GG_save, G2_save); mt = m[+,+]; Vhat = A/n; d = max(ncol(m),nrow(m)); if d = 1 then do; v = trace(Vhat)/mt; Sigmahat = v*I(mt); return(Sigmahat); end; if m[##] = d then return(Vhat); if sig_iter = 0 then do; GG = J(1,mt*mt,0); /* Initialize with a row of zeros */ do i =1 to d; Ei = makeE(m,i); do j = i+1 to d; Ej = makeE(m,j); GG = GG // (Ej` @ Ei`); end; end; nr = nrow(GG); GG = GG[2:nr,]; /* Remove row of zeros */ mmt=mt; G2 = (I(mt*mt)+commute(mt,mmt))*GG`; GG_save=GG; G2_save=G2; end; old = {0}; do i = 1 to d; if i = 1 then s = 0; else s = m[1:i-1][+,+]; mi = m[i]; Temp=A[s+1:s+mi,s+1:s+mi]; Temp=(Temp+Temp`)/2; call eigen(eigs,dumb,Temp); vv = eigs[<>]/n; old = old // vv; Vhat[s+1:s+mi,s+1:s+mi] = I(mi)*vv; end; old = old[2:d+1]; if sig_iter > 0 then do; GG=GG_save; G2=G2_save; end; sig_iter=sig_iter+1; old2 = GG*vec(Vhat); old = old // old2; del = 1; guess = Vhat; Vi = inv(Vhat); logl = -trace(A*Vi)/2 - n*log(det(Vhat))/2; d1 = (mt*mt-m[##])/2; iter = 0; do while(del > converge); H11 = J(d,d,0); TT = Vi*A*Vi/n; mi1 = m[1]; Vi1 = Vi[,1:mi1]; Vii = Vi[1:mi1,1:mi1]; Tii = TT[1:mi1,1:mi1]; h1 = trace(Tii-Vii)/2; U = Vi1*Vi1`; H21 = GG*vec(U); do j1 = 1 to d; if j1 = 1 then sj1 = 0; else sj1 = m[1:j1-1][+,+]; mj1 = m[j1]; Vi1j1 = Vi[1:mi1,sj1+1:sj1+mj1]; v = vec(Vi1j1); H11[1,j1] = v`*v/2; H11[j1,1] = H11[1,j1]; end; do i1 = 2 to d; si1 = m[1:i1-1][+,+]; mi1 = m[i1]; Vi1 = Vi[,si1+1:si1+mi1]; Vii = Vi[si1+1:si1+mi1,si1+1:si1+mi1]; Tii = TT[si1+1:si1+mi1,si1+1:si1+mi1]; h1 = h1 // trace(Tii-Vii)/2; U = Vi1*Vi1`; H21 = H21 || GG*vec(U); do j1 = i1 to d; sj1 = m[1:j1-1][+,+]; mj1 = m[j1]; Vi1j1 = Vi[si1+1:si1+mi1,sj1+1:sj1+mj1]; v = vec(Vi1j1); H11[i1,j1] = v`*v/2; H11[j1,i1] = H11[i1,j1]; end; end; H22 = GG*(Vi@Vi)*G2; H = - ((H11 || H21`) // (H21 || H22)); h2 = GG*vec(TT-Vi); hh = h1 // h2; change = -inv(H)*hh; c = 2; dd = -1; do while(dd < -1.e-10); c = c/2; new = old+c*change; new1 = new[1:d]; new2 = new[d+1:d+d1]; Vhat = dvec(GG`*new2,mt); Vhat = Vhat + Vhat`; do i = 1 to d; if i = 1 then s = 0; else s = m[1:i-1][+,+]; mi = m[i]; Vhat[s+1:s+mi,s+1:s+mi] = I(mi)*new1[i]; end; call eigen(eigs,dumb,Vhat); dd = eigs[><]; if dd > 0 then do; dd1=dd; vi=inv(Vhat); newl = -trace(A*Vi)/2-log(det(Vhat))*n/2; dd = newl - logl; end; end; iter = iter + 1; old = new; del = abs(logl-newl); logl = newl; end; return(Vhat); finish mlesigma; /*---------------------------------------------------------------*/ /* Module for computing the log likelihood for the following */ /* four covariance structures. In the growth curve model, these */ /* are conditional covariance structures. */ /* 1. Unconstrained. */ /* 2. Local Sphericity */ /* 3. Local Sphericity with Independence. */ /* 4. Omnibus Sphericity. */ /*---------------------------------------------------------------*/ start lnlike(lam) global(Y,X,XPXiX,C,G,logg,n,q,t,wd,ci, boxcox,foldpow,converge,del_max,maxlnlik,lnlik); if boxcox = 'yes' then do; lnJ=(lam-1)*logg; if abs(lam) > .0001 then Z = (Y##lam-1)/lam; else Z = log(Y); end; if foldpow = 'yes' then do; lnJ=log(Y##(lam-1) + (J(n,t,1)-Y)##(lam-1)); lnJ=lnJ[+,+]; if abs(lam) > .0001 then Z = (Y##lam -(J(n,t,1)-Y)##lam)/lam; else Z= log(Y/(J(n,t,1)-Y)); end; B = XPXiX*Z; S = (Z`*Z - Z`*X*B)/n; grow = 'no'; if ncol(G) < t then grow = 'yes'; SigGQ = S; if grow = 'yes' then do; Sinv = inv(S); temp1 = G*ginv(G`*Sinv*G); Bg = B*Sinv*temp1; SigGQ = S + (B - Bg*G`)`*X`*X*(B - Bg*G`)/n; end; if q[+] = 0 then do; lnlik = n*t/2 + n/2*log(det(SigGQ)) - lnJ; /*Confidence interval calculation*/ if ci=1 then lnlik = 10000*(lnlik - maxlnlik)**2; else lnlik=lnlik-del_max; return(lnlik); end; /* Compute Omegahat and Sigmahat for various covariance structures */ if wd = 1 then Sigmag = SigGQ; /* Unconstrained */ else do; Sig2=S; if grow='yes' then Sig2 = temp1*G`; Omegat = C`*Sig2*C; Omegati = inv(Omegat); /* Local Sphericity*/ if wd = 2 then Omegah = mlesigma(n*Omegat,n,q,converge); /* Local with independence*/ if wd = 3 then do; mt = q[+,+]; Omegah = J(mt,mt,0); s2 = 0; do i = 1 to max(nrow(q),ncol(q)); s1 = s2 + 1; s2 = s2 + q[i]; thetai = trace(Omegat[s1:s2,s1:s2])/q[i]; Omegah[s1:s2,s1:s2] = thetai*I(q[i]); end; end; /* Omnibus */ if wd = 4 then do; sizeO = ncol(C); Omegah = trace(Omegat)*I(sizeO)/sizeO; end; Sigmag = SigGQ + Sig2*C*Omegati*(Omegah-Omegat)*Omegati*C`*Sig2; end; lnlik = n*t/2 + n/2*log(det(Sigmag)) - lnJ; /*Confidence interval calculation*/ if ci=1 then lnlik = 10000*(lnlik - maxlnlik)**2; else lnlik=lnlik-del_max; return(lnlik); finish lnlike; /*-------------------------------------------------*/ /*------------------ Main Module ----------------- */ /*-------------------------------------------------*/ do; use outx; read all var _num_ into X; t = levels[#,#]; /*Product of levels of within-subject factors*/ n = nrow(X); p2 = ncol(X); Y = X[,1:t]; tmp=vec(Y); tmp2=dvec(tmp,n); X = X[,t+1:p2]; r = rankM(X); maxY=1; error='no'; if boxcox = 'yes' then do; maxY=max(Y); /* Normalization of Data for numerical stability*/ Y = Y/maxY; /* for Box-Cox family only */ if foldpow = 'yes' then do; print 'Input error: Both transformation families have been selected.'; print 'Select either Box-Cox or folded-power,'; error = 'yes'; end; end; minY=min(Y); if minY <= 0 then do; print 'Data Error: One or more data values are non-positive'; error='yes'; end; if foldpow='yes' then do; maxYY = max(Y); if maxYY > 1 then do; print 'Data Error: One or more data values are greater than 1'; print 'For the folded power family, all data values must be '; print 'in (0, 1).'; error = 'yes'; end; end; checkdeg=J(nrow(Polydeg),1,1)*levels-Polydeg; if min(checkdeg) < 1 then do; print 'Input Error: Polynomial degree cannot exceed'; print 'the number of levels minus one.'; print levels; print Polydeg; error = 'yes'; end; if error = 'yes' then goto exit; logg = sum(log(Y)); XPXiX = ginv(X`*X)*X`; /*----------------------------------------------------------------*/ /* Main module to compute mle of lambda (transformation parameter)*/ /* and a confidence interval for lambda under various covariance */ /* structures. In the growth curve model, these are conditional */ /* covariance structures. */ /* wd = 1 --> Unconstrained */ /* wd = 2 --> Local Sphericity */ /* wd = 3 --> Local Sphericity with Independence */ /* wd = 4 --> Omnibus Sphericity */ /*----------------------------------------------------------------*/ sizew = nrow(Polydeg); crit = 2*gaminv(confid,.5); /*percentile of Chi-square 1 distribution*/ critn = sqrt(crit); /*percentile of standard normal dist*/ do wj = 0 to sizew; if wj = 0 then do; /*Find the mle for unconstrained full model */ wd = 1; /*Unconstrained*/ w = levels - J(1,ncol(levels),1); run makeCG(C,G,q,d,levels,w); f = ncol(G); guess = guesslam; /*Initial guess for lambda */ ci = 0; /* Flag for finding the mle */ maxlnlik = 0; sig_iter=0; del_max=0; call nlpnrr(rc,lambda,"lnlike",guess); lam_full=lambda; guess=lambda; bdim = r*t + t*(t+1)/2 + 1; blnlik = -lnlik -n*t*log(maxY); lnlik_m=-lnlik; end; else do; lofv = { 0 0 0 0 0 }; estm = { 0 0 0 0 0 }; w = Polydeg[wj,]; run makeCG(C,G,q,d,levels,w); f = ncol(G); nwd=4; if max(nrow(levels),ncol(levels)) = 1 then nwd=2; if levels[+] = 1 then nwd=1; do wdtemp = 1 to nwd; wd=wdtemp; if nwd=2 & wdtemp=2 then wd=4; /* Find the mle estimate of the transformation parameter */ ci = 0; /* Flag for finding the mle */ maxlnlik = 0; sig_iter=0; full='no'; if wd=1 & f=t then do; lambda=lam_full; logl=lnlik_m; full = 'yes'; end; if full = 'no' then do; del_max=-lnlik_m; call nlpnrr(rc,lambda,"lnlike",guess); guess=lambda; logl = -lnlik-del_max; end; /* Estimate the standard error of lambdahat */ sig_iter=0; del_max=0; call nlpfdd(maxlik,gradient,se,"lnlike",lambda); se = 1/sqrt(se); /* Find lower confidence interval limit for lambda */ ci = 1; /* Flag for finding an endpoint of the ci */ maxlnlik = maxlik + crit/2; guess1 = lambda - critn*se; sig_iter=0; low=guess1; tc=J(13,1,0); tc[1]=200; tc[2]=500; tc[3]=1.e-6; tc[4]=1.e-8; tc[6]=1.e-5; tc[7]=2.2204e-16; tc[10]=1; call nlpnrr(rc,low,"lnlike",guess1) tc=tc; /* Find upper confidence interval limit for lambda */ guess2 = lambda + critn*se; sig_iter=0; high=guess2; call nlpnrr(rc,high,"lnlike",guess2) tc=tc; /* Calculate Lack of Fit Statistic and corresponding p-value*/ /* plus the dimension of the model and AIC */ /* Calculate dimensions*/ qq = q[+,+]; if wd = 1 then dim = r*f + t*(t+1)/2 + 1; if wd = 2 then dim = r*f + (t*(t+1) - qq - ssq(q))/2 + d + 1; if wd = 3 then dim = r*f + (t-qq)*(t+qq+1)/2 + d + 1; if wd = 4 then dim = r*f + (t-qq)*(t+qq+1)/2 + d + 1; logl=logl-n*t*log(maxY); lof = -2*(logl - blnlik); df = bdim - dim; if lof > 0 then pv = 1 - probchi(lof,df); else do; df=.; pv = .; lof = .; end; aic = logl - dim; lofv = lofv // ( lambda || logl || lof || df || pv ); estm = estm // ( lambda || low || high || dim || aic ); end; lofv = lofv[2:nwd+1,]; estm = estm[2:nwd+1,]; rows1 = {'Unconstrained','Local Sphericity','Local w Indep','Omnibus'}; rows2 = {'Unconstrained','Omnibus Spher.'}; rows=rows1; if nwd=2 then rows=rows2; if w = levels - J(1,ncol(levels),1) then do; if levels[+] > 1 then do; print 'Repeated Measures Model'; rows3=rows; end; else do; print 'Linear Model with No Repeated Measures'; rows3={' '}; end; mattrib lofv rowname = rows3 colname = {'Estimate','LogLik','LoF','df','p-value'} label = ' ' format = F9.4; end; else do; print 'Growth Curve Model'; mattrib lofv rowname = rows colname = {'Estimate','LogLik','LoF','df','p-value'} label = 'Cond Cov Struct' format = F9.4; mattrib w label = ' ' colname = ' ' rowname={'Polynomial Degree ='} format=F2.0; print w; end; print lofv; clevel=confid*100; dum2 = char(clevel,3); lowlim = concat(dum2,'% Low'); highlim = concat(dum2,'% High'); cols = {'Estimate'} // lowlim // highlim // {'Dimen', 'AIC'}; rows3=rows; if levels[+] = 1 then rows3 = {' '}; mattrib estm rowname = rows3 colname = cols label = ' ' format = F9.4; print estm; if wj < sizew then do; print /; end; end; end; if levels[+] = 1 then do; /* Compute loglikelihood values for plot */ if "&plot_like" = 'yes' then do; wd = 1; /*Unconstrained*/ w = levels - J(1,ncol(levels),1); run makeCG(C,G,q,d,levels,w); f = ncol(G); ci = 0; /* Flag for finding the mle */ maxlnlik = 0; sig_iter=0; del_max=0; lnlik_mle=lnlik_m-n*t*log(maxY); lnlik_ci = lnlik_mle - crit/2; lb = low - (high-low)/2; ub = high+ (high-low)/2; ss=(ub-lb)/100; do lam=lb to ub by ss; loglike = lnlike(lam); loglike = -loglike -n*t*log(maxY); tab = tab // (lam || loglike || 1); end; minlik = min(tab[,2]); tab = tab // ( low || lnlik_ci || 2 ) // (low || minlik || 2); tab = tab // ( high || lnlik_ci || 3 ) // (high || minlik || 3); tab = tab // ( lam_full || lnlik_mle || 4) // (lam_full || minlik || 4); create plotdata from tab [ colname = { Lambda LogLikel graph } ]; append from tab; end; end; /*------------------------------------------------------------*/ exit: if error = 'yes' then abort; end; quit; /*-----------------------------------------------------------------*/ /*-- PROC gplot: Plot likelihood function if there are no----------*/ /*------------------repeated measures--------------------------- --*/ /*-----------------------------------------------------------------*/ %if &plot_like = yes %then %do; /*This is the main SAS code for doing the Plot of the */ /*Log Likelihood Function*/ symbol1 color=black interpol=join width=2 value=none height=3; symbol2 color=green interpol=join width=2 value=none line=2 height=3; symbol3 color=red interpol=join width=2 value=none line=3 height=3; symbol4 color=blue interpol=join width=2 value=none line=1 height=3; axis1 label=(font=greek 'l'); axis2 label=(font=centx 'ln L' font=greek '(l)'); legend across=1 position=(top inside left) label=none cborder=black value=(font=zapf height=1 tick=1 'Loglikelihood' tick=2 'Lower Bound' tick=3 'Upper Bound' tick=4 'MLE of ' font=greek 'l') mode=protect shape=symbol(3,1); * proc print data = plotdata; proc gplot data=plotdata; plot LogLikel*Lambda=graph /legend=legend1 haxis=axis1 vaxis=axis2; /*===============================================================*/ run; %end; %mend transform; /* title font=zapf height=2 "&plot_title"; */