## Final update today [R for BE/BA]

 rm(list=ls(all=TRUE))   ###### ######   Section 1: Household things ###### CreateX=function(D) {  ##first the two treatments  TrtT=as.numeric(as.character(D$Trt)=="T") TrtR=as.numeric(as.character(D$Trt)=="R")  X=cbind(TrtT, TrtR)  R= qr(X)$rank ## subjects: in a mixed model with subj as random ## we do not do subjects also as fixed, therefore they are #'ed away here ## for (s in unique(D$Subj))  ## {  ## v=as.numeric(D$Subj==s) ## #print(v) ## XX=data.frame(X, v) ## names(XX)[ncol(XX)] = paste("S", s, sep="") ## rnk=qr(XX)$rank  ## if (rnk>R)   ##  {   ##     X=XX   ##     R=rnk   ##  }   ## }  ##now the Pers  for (p in unique(D$Per)) { v=as.numeric(D$Per==p)   #print(v)   XX=data.frame(X, v)   names(XX)[ncol(XX)] = paste("P", p, sep="")   rnk=qr(XX)$rank if (rnk>R) { X=XX R=rnk } } for (q in unique(D$Seq))  {   v=as.numeric(D$Seq==q) #print(v) XX=data.frame(X, v) names(XX)[ncol(XX)] = paste("Q", q, sep="") rnk=qr(XX)$rank   if (rnk>R)     {        X=XX        R=rnk     }  }  return(as.matrix(X)) } Create.CovM=function(Params) ##block diagonal covariance matrix {   varT=Params   varBR=Params   varWR=Params   varRT=Params   #cat("Vars:", varT, varBR, varWR,"\n")   Nobs=length(D$Y) V=matrix(0,ncol=Nobs, nrow=Nobs) for (iRow in 1:Nobs) for (iCol in 1:Nobs) { ## the diagonal if (iCol==iRow) { if (D$Trt[iRow]=="T") V[iRow,iCol]=V[iRow,iCol]+varT       if (D$Trt[iRow]=="R") V[iRow,iCol]=V[iRow,iCol]+varWR +varBR } ## off diagonal if (iCol!=iRow) if (D$Subj[iRow]==D$Subj[iCol]) { if (D$Trt[iCol]==D$Trt[iRow]) V[iRow,iCol]= V[iRow,iCol]+varBR if (D$Trt[iCol]!=D$Trt[iRow]) V[iRow,iCol]= V[iRow,iCol]+varRT } } return(as.matrix(V)) } ###### ###### Section 2: Matrix things ###### Obj.F12=function(Pars) ##this is line 3 of page 10 of: ##http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf { CovM=Create.CovM(Pars) A= -0.5*log(det(CovM)) B= -0.5*log(det(t(X) %*% solve(CovM) %*% X)) est.b = solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y tmp= y - X %*% est.b C=-0.5 *(t(tmp) %*% solve(CovM) %*% tmp) return(A+B+C) } #### #### section 3: execution #### Some.Initial.Guesses=function(foo) { #guess varT D1=subset(D, (D$Trt=="T"))  m=lm(log(Y)~factor(Seq)+factor(Per)   , data=D1)  varT=anova(m)["Residuals", "Mean Sq"]  #guess VarWR  D1=subset(D, (D$Trt=="R")) m=lm(log(Y)~factor(Seq)+factor(Per)+factor(Subj) , data=D1) varWR=anova(m)["Residuals", "Mean Sq"] #guess varBR S=unique(D$Subj)  v=NULL  for (s in S)  {   Dx=subset(D, (D$Subj==s) & (D$Trt=="R"))   v=c(v, mean(log(Dx$Y))) } varBR=var(v)-varWR varRT=0.5*(varT+varBR) rslt=c(varT, varBR, varWR, varRT) return(rslt) } D=read.csv("EMAII.a.csv" , header=T, sep="\t") ##public = easier X=CreateX(D) ##public = easier y=log(D$Y) ##public = easier MyREML=function(foo.bar)   {   L=Some.Initial.Guesses(0)   print(L)   parIni=L   F=optim(par=parIni, fn=Obj.F12,       control=list(reltol=1e-9, trace=T, fnscale=-1))   Nms=c("varT","varBR", "varWR", "varRT")   X=data.frame(Nms, F\$par, parIni)   names(X)=c("Var.Component", "Value", "Ini.value")   return(X) }   MyREML(0) 

which gives
   Var.Component      Value  Ini.value 1          varT 0.07050210 0.06841615 2         varBR 0.03621403 0.02776728 3         varWR 0.01324510 0.01240137 4         varRT 0.04563784 0.04809171 

By the way: You can adjust how many decimals you want on e.g. varWR (~s2WR) through the reltol setting. Use e.g. 1e-12 to get 0.0.01324652 with just a few more iterations.

Job done?

Pass or fail!
ElMaestro Ing. Helmut Schütz 