Dear R users:
I'm new to R and am trying to fit a mixed model
Cox regression model with coxme function.
I have one two-level factor (treat) and one
covariate (covar) and 32 different groups
(centers). I'd like to fit a random coefficients model, with treat and covar
as fixed factors and a random intercept, random
treat effect and random covar slope per center.
I haver a couple of doubts on how to use coxme function for this task:
* Following Therneau's indications (On
mixed-effect Cox models, sparse matrices, and
modeling data from large pedigrees, December 31,
2007) I've modelled a model with a random
intercept and a random treat effect, as follows:
y<-
rep(c(2,4,1,5,6,7,6,5,2,5,6,5,7,3,4,5,2,4,1,5,6,7,6,5,2,5,6,5,7,3,4,5,
6,5,7,8,5,6,8,6,5,6,8,6,5,8,6,9,6,5,7,8,5,6,8,6,5,6,8,6,5,8,6,9),2)
treat<-rep(c(0,1),each=64)
covar<-
rep(c(1,3,2,5,6,4,7,8,9,6,4,3,4,3,2,1,4,3,5,4,6,3,4,5,6,4,6,3,6,7,4,5,3,4,5,6,4,6,6,4,7,6,5,7,5,4,3,2,8,6,4,9,7,9,5,4,6,8,6,4,6,2,5,8),2)
uncens<-
rep(c(0,1,0,0,0,1,0,0,1,1,1,1,0,1,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,1,0
,0,1,1,0,0,1,0,0,1,0,1,0,0,0,1,1,0,1,1,0,0,1,0,0,1,0,1,0,0,0,1,1),2)
centers<- rep(rep(1:32, 2),2)
data1<-list(y, treat, covar, uncens, centers)
names(data1)<- c("y", "treat", "covar", "uncens", "centers")
data1<-as.data.frame(data1)[order(centers,treat),]
library(kinship)
ugroup<-paste(rep(1:32, each=2), rep(0:1, 32), sep='/') #unique groups
mat1<-bdsmatrix(rep(c(1,1,1), 32),
blocksize=rep(2,32), dimnames=list(ugroup, ugroup))
mat2<-bdsmatrix(rep(c(0,0,0,1), 32),
blocksize=rep(2,32), dimnames=list(ugroup, ugroup))
group1<-paste(data1$centers, data1$treat, sep="/")
fit2<- coxme(Surv(y, uncens)~ treat + covar,
data1,random=~1|group1, varlist=list(mat1, mat2),rescale=F, pdcheck=F)
However, I'm not sure about how to proceed to
introduce the random slope for covar. Should I
introduce a third variance matrix per center and
redefine the groups as follows?
ugroup<-paste(rep(1:32, each=3), rep(c(0:1,"z"), 32), sep='/') #unique groups
mat1_1<-bdsmatrix(rep(c(1,1,1,1,1,1), 32),
blocksize=rep(3,32), dimnames=list(ugroup, ugroup))
mat2_1<-bdsmatrix(rep(c(0,0,0,1,0,0), 32),
blocksize=rep(3,32), dimnames=list(ugroup, ugroup))
mat3_1<-bdsmatrix(rep(c(0,0,0,0,0,1), 32),
blocksize=rep(3,32), dimnames=list(ugroup, ugroup))
group<-paste(data1$centers, data1$treat, data1$covar, sep="/")
fit3<- coxme(Surv(y, uncens)~ treat + covar,
data1, random=~1|group1, varlist=list(mat1_1,
mat2_1, mat3_1), rescale=F, pdcheck=F)
This seems to work, but I don't think I'm getting what I want.
* What if the treatment factor has more than two
levels. Should I follow the same procedure, with just bigger block sizes?
* Coxme returns a variance per each of the
variance matrices I defined, but no residual
variance estimate. Is there a way to get it?
Thanks a lot in advance,
Best regards,
Asier
----------
Asier R. Larrinaga
Institut Mediterrani d'Estudis Avançats (IMEDEA)
Carrer Miquel Marqués, 21
07.190-Esporles
Mallorca- Illes Balears
Spain
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
