[R-sig-ME] Toeplitz covariance structure

[Ben Bolker]

   I will admit that I am having trouble getting the connection between 
the phi parameters and the correlations.  However, the fitted models 
(and their implied correlation matrices) are the same:

all.equal(c(logLik(mmrm_toe)), logLik(mmrm_new_toe))

all.equal(unclass(v1 <- getVarCov(mmrm_toe)),
           unname(v2 <- VarCorr(mmrm_new_toe)),
           tolerance = 1e-5)

cs <- mmrm_toe$modelStruct$corStruct
cc


Pinheiro and Bates (2009, via Google Books) say

## for AR models of order greater than 1, the correlation function
## does not admit a simple representation ... being defined recursively
## through the difference equation (Box et al., 1994, Sec 3.2.2)

## h(k, Phi) = phi_1*h(|k-1|,Phi) + ... + phi_p*h(|k-p|,Phi), k = 1,2,...

I don't understand this yet ... would have expected correlations phi1, 
phi1^2 + phi2, ....

   glmmTMB *should* be able to do this as well but I haven't succeeded 
in making the results match ...

On 2023-05-26 5:33 a.m., Jackson, Daniel wrote:
> Dear r-sig-mixed-models
> 
> I am interested in fitting a model with the Toeplitz covariance structure, using the gls function from the nlme package. I took an example from here
> 
> Mixed model repeated measures (MMRM) in Stata, SAS and R - The Stats Geek<https://thestatsgeek.com/2020/12/30/mixed-model-repeated-measures-mmrm-in-stata-sas-and-r/>
> 
> And
> 
> Mixed models repeated measures (mmrm) package for R - The Stats Geek<https://thestatsgeek.com/2022/10/31/mixed-models-repeated-measures-mmrm-package-for-r/>
> 
> I have tried fitting the model using both the new mmrm and the gls function, codes are below. But I have some questions
> 
> 
>    1.  Am I correct in specifying the covariance structure using a corARMA structure with p = "number of visits" - 1 and q=0?
>    2.  All estimates using my proposed gls implementation and the mmrm package seem to agree well except the Phi reported by gls do not seem to be in agreement with the correlation coefficients implied by the mmrm package? Are the Phi a transformation of the correlations?
> 
> Many thanks, code below:
> 
> # Code taken from/inspired by The stats geek website
> # simulate data
> ###############################################################################################
> expit <- function(x) {
>    exp(x)/(1+exp(x))
> }
> 
> set.seed(1234)
> n <- 500
> library(MASS)
> #we will make correlation with baseline the same to visit 1 and visit 2
> corr <- matrix(1, nrow=4, ncol=4) + diag(0.5, nrow=4)
> corr
> data <- mvrnorm(n, mu=c(0,0,0,0), Sigma=corr)
> 
> trt <- 1*(runif(n)<0.5)
> y0 <- data[,1]
> y1 <- data[,2]
> y2 <- data[,3]
> y3 <- data[,4]
> 
> #add in effect of treatment
> y1 <- y1+trt*0.5
> y2 <- y2+trt*1
> y3 <- y3+trt*1.5
> 
> #now make some patients dropout before visit 1
> #r1=1 indicates visit 1 observed
> r1 <- 1*(runif(n)<expit(3-y0))
> 
> #dropout before visit 2, based on change from y0 to y1
> r2 <- 1*(runif(n)<expit(3-(y1-y0)))
> r2[r1==0] <- 0
> 
> #dropout before visit 3, based on change from y1 to y2
> r3 <- 1*(runif(n)<expit(3-(y2-y1)))
> r3[r2==0] <- 0
> 
> y1[r1==0] <- NA
> y2[r2==0] <- NA
> y3[r3==0] <- NA
> 
> wideData <- data.frame(id=1:n, trt=trt, y0=y0, y1=y1, y2=y2, y3=y3)
> summary(wideData)
> 
> longData <- reshape(wideData, varying=c("y1", "y2", "y3"),
>                      direction="long", sep="", idvar="id")
> longData <- longData[order(longData$trt,longData$id,longData$time),]
> 
> 
> ## use gls - but is this correct? - p=2 because 3 visits and q=0 is correct?
> require(nlme)
> mmrm_toe <- gls(y~y0*factor(time)+trt*factor(time),
>              na.action=na.omit, data=longData,
>              correlation=nlme::corARMA(p=2, q=0, form=~time | id),
>              weights=nlme::varIdent(form=~1|time))
> summary(mmrm_toe)
> ## use mmrm package
> library(mmrm)
> #To fit the same model using the new mmrm package, I discovered I first needed to convert the id, time and trt variables to factors. I then fitted the MMRM model using:
> longData_factor<-longData
> longData_factor$time<-factor(longData_factor$time)
> longData_factor$id<-factor(longData_factor$id)
> longData_factor$trt<-factor(longData_factor$trt)
> 
> 
> 
> mmrm_new_toe <- mmrm(y~y0*time+trt*time+toeph(time|id), data=longData_factor)
> summary(mmrm_new_toe)
> 
> ## Questions: is gls implementation correct and why do Phi reported not seem to agree with correlations from mmrm?
> 
> Dan Jackson
> 
> ________________________________
> 
> AstraZeneca UK Limited is a company incorporated in Engl...{{dropped:16}}
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models