Skip to content
Prev 248729 / 398506 Next

survreg 3-way interaction

The fire devoured my laptop, but the PC is still working... :)
Sorry for the lazyness, you are totally right.

Here goes a reproducible example, which resembles the main features of our
dataset:

# three experimental factors, 30 subjects
expgrid <- expand.grid(F1=0:1,F2=0:1,F3=0:1,id=1:30)

# set up the dataset (make it unbalanced)
inpdata <- expgrid[sample(nrow(expgrid),nrow(expgrid)*3,replace=T),]
inpdata <- inpdata[order(inpdata$id),]

# define coefficients
# fixed effects
coefs <- c(5,rnorm(7))
# random effects
subj.coefs <- rnorm(30)

# compute fixed effects
modmat <- model.matrix(logtimes~F1*F2*F3,
	data=data.frame(logtimes=0,inpdata))
fix.eff <- rep(modmat%*%coefs)

# compute logtimes (fixed effects + random effects + noise)
logtimes <- fix.eff+rep(subj.coefs,table(inpdata$id))+
	rnorm(nrow(inpdata),0,0.5)

# random censoring
cens <- runif(nrow(inpdata))<0.2
logtimes[cens] <- logtimes[cens]-rnorm(sum(cens))

# finalize dataset
inpdata$status <- !cens
inpdata$times <- exp(logtimes)
inpdata$id <- factor(inpdata$id)

# run survreg with fixed df in frailty
# (method="aic" leads to error, don't know why)
library(survival)
survreg3 <- survreg(Surv(times,status) ~
	F1*F2*F3 + frailty.gamma(id,df=29), data=inpdata,
	dist="lognormal")
# Error in survpenal.fit(X, Y, weights, offset, init = init, controlvals =
control,  :
#  Invalid pcols or pattr arg

# note, that without the frailty term or without the tree-way interaction
term, it runs smoothly


# the same with gamlss
library(gamlss.cens)
gamlss3 <- gamlss(Surv(times,status) ~
	F1*F2*F3 + random(id,df=29), data=inpdata,
	family=cens(LOGNO))