## Update III [🇷 for BE/BA]

Look at this baby:

 rm(list=ls(all=TRUE)) ###### par vector is c(varT, varBR, varWR) ###### ######   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[1]   varBR=Params[2]   varWR=Params[3]   #varRT=Params[4]   #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) } Some.Initial.Guesses=function(foo) { #guess variance within R D1=subset(D, D$Trt=="R")   m=lm(log(Y)~factor(Subj)+factor(Seq)+factor(Per) , data=D1)   #print(anova(m))   varWR= anova(m)["Residuals", "Mean Sq"]   m=lm(log(Y)~factor(Subj)+factor(Seq)+factor(Per)+factor(Trt) , data=D)   varBR= anova(m)["factor(Subj)", "Mean Sq"]   ##total var T, cheaply   ##actually, this may be much better: varT=varBR+varWR   D2=subset(D, D$Trt=="T") varT=var(log(D2$Y))   varT=varBR+varWR   L=list(varT=varT, varBR=varBR, varWR=varWR)   return(L) } #### #### section 3: execution #### 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=c(0.06, 0.032, 0.012) F=optim(par=parIni, fn=Obj.F12, #method="BFGS", #lower=c(0.003,0.003,0.003), upper=c(1,1,1), control=list(reltol=1e-11, trace=T, fnscale=-1)) print(F) } MyREML(0)  On my machine it gives: $par [1] 0.06776119 0.03636979 0.01237020

Pass or fail!
ElMaestro