phmm (version 0.4)

phmm: Proportional Hazards Model with Mixed Effects

Description

Fits a proportional hazards regression model incorporating random effects. The function implements an EM algorithm using Markov Chain Monte Carlo (MCMC) at the E-step as described in Vaida and Xu (2001).

Usage

phmm(formula, random, data, subset, na.action, Sigma, varcov, NINIT,
	VARSTART, MAXSTEP, CONVERG, emstep, Gbs, Gbsvar, verbose)

Arguments

formula

model formula for the fixed part of the model, i.e. the part that specifies $\beta'x_{ij}$. See survreg for further details. Intercept is implicitely included in the model by estimation

random

formula for the random part of the model, i.e. the part that specifies $w_{ij}'b_{i}$. If omitted, no random part is included in the model. E.g. to specify the model with a random intercept, say random=~1.

data

optional data frame in which to interpret the variables occuring in the formulas.

subset

subset of the observations to be used in the fit.

na.action

function to be used to handle any NAs in the data. The user is discouraged to change a default value na.fail.

Sigma

initial covariance matrix for the random effects. Defaults to "identity".

varcov

constraint on Sigma. Currently only "diagonal" is supported.

NINIT

number of starting values supplied to Adaptive Rejection Metropolis Sampling (ARMS) algorithm.

VARSTART

starting value of the variances of the random effects.

MAXSTEP

number of EM iterations.

CONVERG

iteration after which Gibbs sampling size changes from Gbs to Gbsvar.

emstep

initial step value.

Gbs

initial Gibbs sampling size (until CONVERG iterations).

Gbsvar

Gibbs sampling size after CONVERG iterations.

verbose

Set to TRUE to print EM steps.

Value

• The function produces an object of class "phmm" consisting of:
• steps
  a matrix of estimates at each EM step;
• bhat
  empirical Bayes estimates of expectation of random effects;
• sdbhat
  empirical Bayes estimates of standard deviation of random effects;
• coef
  the final parameter estimates for the fixed effects;
• var
  the estimated variance-covariance matrix;
• loglik
  a vector of length four with the conditional log-likelihood and marginal log-likelihood estimated by Laplace approximation, reciprocal importance sampling, and bridge sampling (only implemented for nreff < 3);
• lambda
  the estimated baseline hazard;
• Lambda
  the estimated cumulative baseline hazard.

Details

The proportional hazards model with mixed effects is equipped to handle clustered survival data. The model generalizes the usual frailty model by allowing log-linearl multivariate random effects. The software can only handle random effects from a multivariate normal distribution. Maximum likelihood estimates of the regression parameters and variance components is gotten by EM algorithm, with Markov chain Monte Carlo (MCMC) used in the E-step.

References

Cox, D.R. 1972. "Regression models and life tables" (with discussion). Journal of the Royal Statistical Society, Series B; 34:187-220.

Gilks, W. R. and Wild, P. (1992) Adaptive rejection sampling for Gibbs sampling. Applied Statistics 41, pp 337-348.

Vaida, F. and Xu, R. 2000. "Proportional hazards model with random effects", Statistics in Medicine, 19:3309-3324.

Xu, R, Gamst, A, Donohue, M, Vaida, F, & Harrington, DP. 2006. Using Profile Likelihood for Semiparametric Model Selection with Application to Proportional Hazards Mixed Models. Harvard University Biostatistics Working Paper Series, Working Paper 43.

See Also

survfit, Surv.

Examples

N <- 100
B <- 100
n <- 50
nclust <- 5
clusters <- rep(1:nclust,each=n/nclust)
beta0 <- c(1,2)
set.seed(13)
#generate phmm data set
Z <- cbind(Z1=sample(0:1,n,replace=TRUE),
           Z2=sample(0:1,n,replace=TRUE),
           Z3=sample(0:1,n,replace=TRUE))
b <- cbind(rep(rnorm(nclust),each=n/nclust),rep(rnorm(nclust),each=n/nclust))
Wb <- matrix(0,n,2)
for( j in 1:2) Wb[,j] <- Z[,j]*b[,j]
Wb <- apply(Wb,1,sum)
T <- -log(runif(n,0,1))*exp(-Z[,c('Z1','Z2')]%*%beta0-Wb)
C <- runif(n,0,1)
time <- ifelse(T<C,T,C)
event <- ifelse(T<=C,1,0)
mean(event)
phmmdata <- data.frame(Z)
phmmdata$cluster <- clusters
phmmdata$time <- time
phmmdata$event <- event

fit.phmm <- phmm(Surv(time, event)~Z1+Z2+cluster(cluster),
   ~-1+Z1+Z2, phmmdata, Gbs = 100, Gbsvar = 1000, VARSTART = 1,
   NINIT = 10, MAXSTEP = 100, CONVERG=90)
summary(fit.phmm)