Final update today [R for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-12 22:27 (298 d 22:48 ago) – Posting: # 21679
Views: 11,463
(edited by ElMaestro on 2020-07-12 23:51)


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:

Activity
 Admin contact
21,445 posts in 4,482 threads, 1,511 registered users;
online 18 (0 registered, 18 guests [including 3 identified bots]).
Forum time: Friday 21:16 CEST (Europe/Vienna)

The plural of anecdote is not data.    Roger Brinner

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