Bioequivalence and Bioavailability Forum

Hi all,

long post, I apologise.

Some time ago a weirdo published a paper about fitting mixed models for semi-replicated designs, where he showed how it is possible to work directly in the observational covariance matrix V, rather than using V=ZGZ^t+R and working in G and R.
This allows us to avoid fitting a model that involves a within-T variance component, since it does not exist with the semi-replicated designs. Great concept, eh?

The author gave some working R code for it. The code did not, however, include derivation of a 90% CI with Satterthwaite's df approximation (or any other approximation for that matter). I have been through literature incl. SAS's helpfiles for Proc Mixed and read into papers on mixed models. As I am not a mathematician, much of it is quite inaccessible to me.


So, I asked https://deepai.org for help with deriving the 90% CI for T-R based on Satt dfs, and so far I failed miserably.

At my end Deepai readily takes on the challenge, and produces all sorts of complicated code which ends up with extremely wrong estimates of the df. For example anywhere from -37165346958870.86 to 1e10.
It ends up in anobscure debate with itself about Hessians, ill-conditioned matrices, ridge-stabilisation, regularisation constants, Fishbutt adjustment, Cojones-Krasilnikoff correction, L-BFGS-B, and much other incomprehensible stuff. It then oscillates between dysfunctional solutions. Whatever I do I seem to end up with severely ill-conditioned Hessians which I guess may be the (or a) reason for badly estmated df's.

Can someone find a way to do it, with or without AI help?

Below is some code to get you started, the code is a slightly adjusted version of the weirdo's spew.
All we need is "the rest" which derives the 90% CI with Satt df.

rm(list=ls(all=TRUE))



##the following lines correspond to EMA dataset II.
Subj=c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4,
       5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8,
       9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12,
       13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16,
       17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20,
       21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24)
Trt=c("R", "T", "R", "R", "T", "R", "R", "R", "T", "T",
      "R", "R", "T", "R", "R", "R", "R", "T", "R", "T",
      "R", "T", "R", "R", "R", "T", "R", "T", "R", "R",
      "R", "R", "T", "R", "R", "T", "R", "R", "T", "R",
      "T", "R", "R", "T", "R", "T", "R", "R", "T", "R",
      "R", "R", "R", "T", "R", "T", "R", "R", "R", "T",
      "R", "R", "T", "T", "R", "R", "T", "R", "R", "R",
      "T", "R")
Per=c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,
      1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,
      1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,
      1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,
      1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3)
Seq=c("RTR", "RTR", "RTR", "RTR", "RTR", "RTR", "RRT",
      "RRT", "RRT", "TRR", "TRR", "TRR", "TRR", "TRR",
      "TRR", "RRT", "RRT", "RRT", "RTR", "RTR", "RTR",
      "TRR", "TRR", "TRR", "RTR", "RTR", "RTR", "TRR",
      "TRR", "TRR", "RRT", "RRT", "RRT", "RRT", "RRT",
      "RRT", "RRT", "RRT", "RRT", "RTR", "RTR", "RTR",
      "RTR", "RTR", "RTR", "TRR", "TRR", "TRR", "TRR",
      "TRR", "TRR", "RRT", "RRT", "RRT", "RTR", "RTR",
      "RTR", "RRT", "RRT", "RRT", "RRT", "RRT", "RRT",
      "TRR", "TRR", "TRR", "TRR", "TRR", "TRR", "RTR",
      "RTR", "RTR")
Y=c(4053.6, 3970.4, 3748.8, 2986.2, 2378.8, 2804.6,
    3464.4, 3340.2, 4028.8, 4105, 3191.2, 3803.6,
    4767.8, 4542.6, 3940, 3050.8, 3027.2, 2419.6,
    2530.2, 3072, 2962.6, 2205, 2041.4, 2018, 4647.6,
    4159.6, 3400, 2228.2, 2360.4, 2221.2, 1863.8,
    2212.4, 2394.4, 2278.4, 3170.4, 3927.2, 2640.4,
    2430.4, 2869.2, 3030.8, 2459.8, 2970.4, 2254.4,
    1994.8, 2724.4, 2959.6, 3442, 3342.6, 2396.8,
    2659.4, 2172, 2725, 2805.6, 3146.6, 2418.8,
    2749.8, 2504, 2662.4, 2929.8, 3037.2, 2869.6,
    2666.4, 3069, 2949, 2926.8, 2855.4, 3154.8,
    3185.6, 3548.6, 1874.8, 1808.8, 2730.8)

D=data.frame(Subj, Seq, Trt, Per, Y)

Mod=lm(log(Y) ~ factor(Subj)+factor(Per)+factor(Seq)+factor(Trt), data=D)
anova(Mod)

Create.Model.Matrix.Fixed.Effects=function(D, Want.Intercept=F)
##creates a full column-rank model matrix
##fixed: Treatment, Period, Sequence
{
 ##first the two treatments
 Int=rep(1, nrow(D))
 TrtT=as.numeric(as.character(D$Trt)=="T")
 TrtR=as.numeric(as.character(D$Trt)=="R")
 X=cbind(TrtT, TrtR)
 if (Want.Intercept) X=cbind(Int, TrtR)
 R= qr(X)$rank

 ##now the Periods
 for (p in unique(D$Per))
 {
  v=as.numeric(D$Per==p)
  XX=data.frame(X, v)
  names(XX)[ncol(XX)] = paste("Per", p, sep="")
  rnk=qr(XX)$rank
  if (rnk>R)
    {
       X=XX
       R=rnk
    }
 }
 
 ##and sequence
 for (q in unique(D$Seq))
 {
  v=as.numeric(D$Seq==q)
  XX=data.frame(X, v)
  names(XX)[ncol(XX)] = paste("Seq", 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]   Nobs=length(D$Y)
  V=matrix(0,ncol=Nobs, nrow=Nobs) ##creates the square matrix of zeros
  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))
}

Log.REML=function(Pars)
{
  CovM=Create.CovM(Pars) ##this is V
  k=  -((length(y)-ncol(X))/2)*log(2*pi)
  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(k+A+B+C)
}


Some.Initial.Guesses=function()
{
 #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
 D1=subset(D, (D$Trt=="R"))
 m=lm(log(Y)~factor(Seq)+factor(Per)   , data=D1)
 varBR=anova(m)["Residuals", "Mean Sq"] - varWR

 #guess varRT this is also a between-variance;
 #might be perhaps close to varBR?
 varRT=varBR     

 rslt=c(varT, varBR, varWR, varRT)
 return(rslt)
}

Run.REML=function(parIni=NULL, Want.Intercept=F)
## call this with a vector of four custom guesses
## of varT, varBR, varWR, and  varRT (in this sequence).
## Or leave it out to allow the script to use a set of qualified guesses
{
  X<<-Create.Model.Matrix.Fixed.Effects(D, Want.Intercept) ##public = easier
  y<<-log(D$Y)                                             ##public = easier
  if (is.null(parIni)) parIni=Some.Initial.Guesses()
 
  F=optim(par=parIni, fn=Log.REML, 
          control=list(reltol=1e-9, maxit=2000,
          trace=F, fnscale= -1))
  if (F$convergence ==0) s="The model converged." else s="The model did not converge."
  Nms=c("varT","varBR", "varWR", "varBRT")
  X=data.frame(Nms, parIni, F$par)
  names(X)=c("Var.Component", "Ini.guess", "REML.Estimate")
  rslt=list(Outcome=s, LogLikelihood= F$value,
            Fixed.eff=est.b, Random.eff=X, V=Create.CovM(F$par))
  return(rslt)
}

L=Run.REML() ; print(L)