mlfm: Semiparametric multilevel frailty models, estimated via EMPL...

Description Usage Arguments Note Author(s) References See Also Examples

View source: R/mlfm.R

Description

The mlfm function fits proportional hazard models with K>0 random effects. Inference is based on the EM-PPL algorithm proposed by Horny (2009).

Usage

1
2
3
4
5
6
  mlfm(formula, data, 
       eps = 1e-10, frailty.eps = 1e-11, maxit = 50,
       showtime = TRUE, verbose = FALSE, sparse = TRUE)
  
  ## S3 method for class 'mlfm'
print(x, digits = max(3, getOption("digits") - 3) - 1, ...)

Arguments

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 data.frame containing all the variables

eps

the threshold used to state convergence

frailty.eps

the threshold used to state convergence within PPL fits (see coxph.control)

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 frailty)

x

the object of class mlfm to print

digits

the number of digits to display in results

...

other arguments (see print)

Note

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.

Author(s)

Authors: Federico Rotolo and Guillaume Horny

Maintainer: Federico Rotolo <federico.rotolo@stat.unipd.it>

References

Horny G (2009) Inference in mixed proportional hazard models with K random effects Stat Papers 50:481–499

See Also

coxph, parfm, frailtypack

Examples

 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

mlfm documentation built on May 2, 2019, 4:16 p.m.

Related to mlfm in mlfm...