Description Usage Arguments Note Author(s) References See Also Examples
The mlfm
function fits proportional hazard models
with K>0 random effects.
Inference is based on the EM-PPL algorithm proposed by Horny (2009).
1 2 3 4 5 6 |
formula |
a formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the Surv() function |
data |
a |
eps |
the threshold used to state convergence |
frailty.eps |
the threshold used to state convergence within PPL fits
(see |
maxit |
the maximum number of iteration for the outer loop (EM) |
showtime |
display the execution time |
verbose |
display the progression status over EM steps |
sparse |
wether to use sparse computation
(see |
x |
the object of class |
digits |
the number of digits to display in results |
... |
other arguments (see |
The standard errors of regression coefficients
as given by the EMPL algorithm are underestimated, so are not returned.
A correct estimation is, for example, the one based on Louis (1982) approach,
still not implemented in mlfm
.
Authors: Federico Rotolo and Guillaume Horny
Maintainer: Federico Rotolo <federico.rotolo@stat.unipd.it>
Horny G (2009) Inference in mixed proportional hazard models with K random effects Stat Papers 50:481–499
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | ############################## - Data simulation - #############################
set.seed(1)
nsim <- 2000
data <- data.frame(id = 1:nsim)
### - Clustering - ###
grVar <- cbind(nclust = c(50, 20, 80),
Fpar = c(1, .5, 1.5))
for (gr in 1:nrow(grVar)) {
data[[paste("gr", gr, sep="")]] <- apply(
rmultinom(n= nsim, 1, rep(1, grVar[gr, 'nclust'])), 2, which.max)
data[[paste("u", gr, sep="")]] <- rgamma(grVar[gr, 'nclust'],
shape=1/grVar[gr, 'Fpar'],
scale=grVar[gr, 'Fpar'])[
data[[paste("gr", gr, sep="")]]]
}
head(data)
### - Covariates - ###
data$Treat <- rbinom(nsim, 1, .5)
data$Age <- rnorm(nsim, 60, 4)
### - Times - ###
lambda <- 2
rho <- 1.2
tcens <- 1.5
simBetas <- c(.5, log(1.2) / 10)
data$Time <-pmin(tcens, (- log(runif(nsim, 0, 1)) / (lambda *
apply(data[, paste("u", 1:nrow(grVar), sep="")], 1, prod) *
as.vector(exp(simBetas %*% t(data[, c("Treat", "Age")])))
))^(1 / rho))
data$Status <- 0 + (data$Time < tcens)
summary(data)
head(data)
############################## - Model fitting - #############################
mlfmod <- mlfm(Surv(Time, Status) ~ Treat + Age +
frailty(gr1) +frailty(gr2) +frailty(gr3) , data=data, verbose=TRUE)
mlfmod
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.