[R-sig-ME] Toeplitz covariance structure

[Jackson, Daniel]

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}}