data(multcif)
multcif$cause[multcif$cause==0] <- 2
addm<-comp.risk(Event(time,cause)~const(X)+cluster(id),data=multcif,
cause=1,n.sim=0)
### making group indidcator
g.des<-data.frame(group2=rep(rbinom(200,1,0.5),rep(2,200)))
theta.des <- model.matrix(~-1+factor(group2),g.des)
out1m<-random.cif(addm,data=multcif,cause1=1,cause2=1,Nit=15,detail=0,
theta=2,theta.des=theta.des,step=1.0)
summary(out1m)
## this model can also be formulated as a random effects model
## but with different parameters
out2m<-Grandom.cif(addm,data=multcif,cause1=1,cause2=1,Nit=10,detail=0,
random.design=theta.des,step=1.0)
summary(out2m)
1/out2m$theta
out1m$theta
####################################################################
################### ACE modelling of twin data #####################
####################################################################
### assume that zygbin gives the zygosity of mono and dizygotic twins
### 0 for mono and 1 for dizygotic twins. We now formulate and AC model
zygbin <- g.des$group2 ## indicator of dizygotic twins
n <- nrow(multcif)
### random effects for each cluster
des.rv <- cbind(theta.des,(zygbin==1)*rep(c(1,0)),(zygbin==1)*rep(c(0,1)),1)
### design making parameters half the variance for dizygotic components
pardes <- rbind(c(1,0), c(0.5,0),c(0.5,0), c(0.5,0), c(0,1))
outacem <-Grandom.cif(addm,data=multcif,causeS=1,Nit=30,detail=0,
theta=c(-1.21,2.1),theta.des=pardes,step=1.0,random.design=des.rv)
summary(outacem)
### genetic variance is
exp(outacem$theta[1])/sum(exp(outacem$theta))^2Run the code above in your browser using DataLab