data one ; input y group @@ ; cards ; 3 1 5 1 6 1 8 2 9.5 2 11 2 10 3 11 3 10.3 3 14 4 16 4 15.7 4 ; proc glm data = one ; class group ; model y = group / solution ; lsmeans group ; run ; proc mixed data = one ; class group ; model y = / solution ; random group / solution ; estimate 'rand1' | group 1 0 0 0 ; estimate 'rand2' | group 0 1 0 0 ; estimate 'rand3' | group 0 0 1 0 ; estimate 'rand4' | group 0 0 0 1 ; estimate 'prand1' intercept 1 | group 1 0 0 0 ; estimate 'prand2' intercept 1 | group 0 1 0 0 ; estimate 'prand3' intercept 1 | group 0 0 1 0 ; estimate 'prand4' intercept 1 | group 0 0 0 1 ; run ; proc iml ; y = {3, 5, 6, 8, 9.5, 11, 10, 11, 10.3, 14, 16, 15.7} ; xfixed = { 1 0 0 0, 1 0 0 0, 1 0 0 0, 0 1 0 0, 0 1 0 0, 0 1 0 0, 0 0 1 0, 0 0 1 0, 0 0 1 0, 0 0 0 1, 0 0 0 1, 0 0 0 1 } ; n = nrow(xfixed) ; pfixed = ncol(xfixed) ; betafixed = inv(xfixed`*xfixed)*xfixed`*y ; sstotal = y` * (i(n) - (1/n)*j(n,n) ) * y ; ssmodel = y` * (xfixed*inv(xfixed`*xfixed)*xfixed` - (1/n)*j(n,n) ) * y ; ssresid = y` * (i(n) - xfixed*inv(xfixed`*xfixed)*xfixed` ) * y ; print xfixed y ; print n pfixed ; print betafixed ; print sstotal ssmodel ssresid ; * defining the random part of the model ; xrandom = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} ; z = { 1 0 0 0, 1 0 0 0, 1 0 0 0, 0 1 0 0, 0 1 0 0, 0 1 0 0, 0 0 1 0, 0 0 1 0, 0 0 1 0, 0 0 0 1, 0 0 0 1, 0 0 0 1 } ; * defining the REML function ; start remlcr(theta) global(xrandom,y,z) ; prandom = ncol(z) ; n = nrow(xrandom) ; G = i(prandom) * theta[1] ; R = i(n) * theta[2] ; V = z * G * z` + R ; bhat = inv(xrandom` * inv(V) * xrandom) * xrandom` * inv(V) * y ; lreml = -(n-prandom)*log(2*3.14159)/2 -log(det(V)) -log(det(xrandom` * inv(V) * xrandom )) -(y - xrandom*bhat)` * inv(V) * (y - xrandom*bhat) ; return(lreml) ; finish remlcr ; * checking the log likelihood function ; par1 = {16, 1} ; reml1 = remlcr(par1) ; par2 = {18, 1.5} ; reml2 = remlcr(par2) ; print reml1 reml2 ; * this optimization uses a derivative-free Nelder-Mead simplex method, which is not very efficient and for demonstration purposes only ; start = {10, 3} ; optn = {1 4 } ; call nlpnms(rreml,thetahat,"remlcr",start,optn,,,,,); print thetahat ; print rreml ; prandom = ncol(z) ; G = i(prandom) * thetahat[1] ; R = i(n) * thetahat[2] ; V = z * G * z` + R ; betahat = inv(xrandom` * inv(V) * xrandom) * xrandom` * inv(V) * y ; predhat = G * z` * inv(V) * (y - xrandom*betahat) ; predmeanhat = betahat*j(4,1) + predhat ; print betahat ; print predhat ; print predmeanhat ; quit ;