R help - Dec 2008 - Random coefficients model with a covariate: coxme function

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

> 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:
 example deleted
 
> * 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?
   The coxme function does not support random slopes.  It's been on my
"to do"
list for a long time.  I am supposed to teach an American Stat Assoc course at 
the end of March, however, which has escalated the urgency.

  If the covariate has only 2 levels, such as a random treatment effect when 
there are only 2 treatments, then by coding the treatment as 0/1 and creating 
just the right covariates you could "trick" coxme into fitting the
model.  This
is what is described in the report.  You essentially make treatment a nested 
effect.
  	fit1 <-  coxme(Surv(y, uncens) ~ treat + covar, data1,
  		    random= ~1 | centers)
        fit2 <-  coxme(Surv(y, uncens) ~ treat + covar, data1,
        	    random= ~1 | centers/treat)
  
  There is no residual variance for a Cox model.  
  
  Your example was very hard to read.  Consider using spaces, indentation, etc 
to make it easier for old eyes.
  
  	Terry T.