phmm (version 0.4)

AIC: Akaike Information Criterion for PHMM

Description

Function calculating the Akaike information criterion for PHMM fitted model objects, according to the formula $-2*log-likelihood + k*rho$, where $npar$ represents the number of parameters in the fitted model. The function returns a list of AIC calculations corresponding different likelihood estimations: conditional and marginal likelihoods calculated by Laplace approximation, reciprocal importance sampling, and bridge sampling (only implemented for nreff < 3). The default k = 2, is for the usual AIC.

Usage

## S3 method for class 'phmm':
AIC(object, ..., k = 2)

Arguments

object

a fitted PHMM model object of class phmm,

...

optionally more fitted model objects.

k

numeric, the penalty per parameter to be used; the default k = 2 is the classical AIC.

Value

• Returns a list of AIC values corresonding to all available log-likelihood values from the fit. See phmm for details of the log-likelihood values.

References

Whitehead, J. (1980). Fitting Cox's Regression Model to Survival Data using GLIM. Journal of the Royal Statistical Society. Series C, Applied statistics, 29(3), 268-.

See Also

phmm, AIC

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)

# Same data can be fit with lmer,
# though the correlation structures are different.
poisphmmdata <- pseudoPoisPHMM(fit.phmm)

library(lme4)
fit.lmer <- lmer(m~-1+as.factor(time)+z1+z2+
  (-1+w1+w2|cluster)+offset(log(N)), 
  poisphmmdata, family=poisson)

fixef(fit.lmer)[c("z1","z2")]
fit.phmm$coef

VarCorr(fit.lmer)$cluster
fit.phmm$Sigma

logLik(fit.lmer)
fit.phmm$loglik

traceHat(fit.phmm)

summary(fit.lmer)@AICtab
AIC(fit.phmm)