Who can help? [🇷 for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-12 15:26 (1377 d 03:47 ago) – Posting: # 21676
Views: 17,539

Important update, I forgot vBR on the diagonal. Now the singuality is less of a heacahe:

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


Check.Surf=function(i, x1,x2, N)
## Checks the log likelihood surface of the i'th variance component
## (given the other variance components) for the interval var = x1 to x2
## with a granularity of N points.
{
  L=Some.Initial.Guesses(0)
  parIni=c(L$varT, L$varBR, L$varWR )
  x=x1;
  vX=NULL
  vY=NULL
  dx=(x2-x1)/N
  for (j in 1:N)
  {
    vX=c(vX, x)
    parIni[i]=x
    vY=c(vY, Obj.F12(parIni))
    x=x+dx
  }
  plot(vX, vY)
}

MyREML=function(foo.bar)
## ideally this would work, but it doesn't right now :-(
## also, I have a feeling the surface may be quite rugged,
## so I may switch to a homecooked optimiser
## once I have ironed out the singularity issue, see below.
{
  L=Some.Initial.Guesses(0)
  print(L)
  parIni=c(L$varT, L$varBR, L$varWR)
  F=optim(par=parIni, fn=Obj.F12, #method="L-BFGS-B",lower=c(0,0,0),
      control=list(reltol=1e-10, trace=T))
  print(F)
}
# if ti worked it'd be called through
 MyREML(0)

# however:

Check.Surf(3, 0.01, 0.2, 50) ##checks the surface in one dimension


#
#
#


—
Pass or fail!
ElMaestro

Complete thread:

UA Flag
Activity
 Admin contact
22,988 posts in 4,825 threads, 1,661 registered users;
84 visitors (0 registered, 84 guests [including 3 identified bots]).
Forum time: 19:14 CEST (Europe/Vienna)

The only way to comprehend what mathematicians mean by Infinity
is to contemplate the extent of human stupidity.    Voltaire

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