Update III [🇷 for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-12 23:46 (2152 d 12:55 ago) – Posting: # 21678
Views: 39,927

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

Complete thread:

UA Flag
Activity
 Admin contact
23,653 posts in 4,991 threads, 1,571 registered users;
250 visitors (0 registered, 250 guests [including 25 identified bots]).
Forum time: 12:41 CEST (Europe/Vienna)

I’m all in favor of the democratic principle
that one idiot is as good as one genius, but I draw the line
when someone takes the next step and concludes
that two idiots are better than one genius.    Leo Szilard

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