Final update today [🇷 for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-13 00:27 (1411 d 05:01 ago) – Posting: # 21679
Views: 17,516


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[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)
}






####
#### 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

Complete thread:

UA Flag
Activity
 Admin contact
23,035 posts in 4,835 threads, 1,645 registered users;
43 visitors (0 registered, 43 guests [including 5 identified bots]).
Forum time: 05:28 CEST (Europe/Vienna)

The mediocre teacher tells.
The good teacher explains.
The superior teacher demonstrates.
The great teacher inspires.    William Arthur Ward

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5