# NOT RUN {
#####
# Fit model with lung data set from survival
# Warning: this has a longer computation time
# }
# NOT RUN {
library(dynamichazard)
.lung <- lung[!is.na(lung$ph.ecog), ]
set.seed(43588155)
pf_fit <- PF_EM(
Surv(time, status == 2) ~ ph.ecog + age,
data = .lung, by = 50, id = 1:nrow(.lung),
Q_0 = diag(1, 3), Q = diag(1, 3),
max_T = 800,
control = PF_control(
N_fw_n_bw = 500,
N_first = 2500,
N_smooth = 2500,
n_max = 50,
n_threads = max(parallel::detectCores(logical = FALSE), 1)),
trace = 1)
# Plot state vector estimates
plot(pf_fit, cov_index = 1)
plot(pf_fit, cov_index = 2)
plot(pf_fit, cov_index = 3)
# Plot log-likelihood
plot(pf_fit$log_likes)
# }
# NOT RUN {
#####
# Can be compared with this example from ?coxph in R 3.4.1. Though, the above
# only has a linear effect for age
# }
# NOT RUN {
cox <- coxph(
Surv(time, status) ~ ph.ecog + tt(age), data= .lung,
tt=function(x,t,...) pspline(x + t/365.25))
cox
# }
# NOT RUN {
# }
# NOT RUN {
######
# example with fixed effects
# prepare data
temp <- subset(pbc, id <= 312, select=c(id, sex, time, status, edema, age))
pbc2 <- tmerge(temp, temp, id=id, death = event(time, status))
pbc2 <- tmerge(pbc2, pbcseq, id=id, albumin = tdc(day, albumin),
protime = tdc(day, protime), bili = tdc(day, bili))
pbc2 <- pbc2[, c("id", "tstart", "tstop", "death", "sex", "edema",
"age", "albumin", "protime", "bili")]
pbc2 <- within(pbc2, {
log_albumin <- log(albumin)
log_protime <- log(protime)
log_bili <- log(bili)
})
# standardize
for(c. in c("age", "log_albumin", "log_protime", "log_bili"))
pbc2[[c.]] <- drop(scale(pbc2[[c.]]))
# fit model with extended Kalman filter
ddfit <- ddhazard(
Surv(tstart, tstop, death == 2) ~ ddFixed_intercept() + ddFixed(age) +
ddFixed(edema) + ddFixed(log_albumin) + ddFixed(log_protime) + log_bili,
pbc2, Q_0 = 100, Q = 1e-2, by = 100, id = pbc2$id,
model = "exponential", max_T = 3600,
control = list(eps = 1e-5, NR_eps = 1e-4, n_max = 1e4))
summary(ddfit)
# fit model with particle filter
set.seed(88235076)
ppfit <- PF_EM(
Surv(tstart, tstop, death == 2) ~ ddFixed_intercept() + ddFixed(age) +
ddFixed(edema) + ddFixed(log_albumin) + ddFixed(log_protime) + log_bili,
pbc2, Q_0 = 100, Q = ddfit$Q * 100, # use estimate from before
by = 100, id = pbc2$id,
model = "exponential", max_T = 3600,
control = PF_control(
N_fw_n_bw = 250, N_smooth = 500, N_first = 1000, eps = 1e-3,
method = "AUX_normal_approx_w_cloud_mean",
n_max = 25, # just take a few iterations as an example
n_threads = max(parallel::detectCores(logical = FALSE), 1), trace = TRUE)
# compare results
plot(ddfit)
plot(ppfit)
sqrt(ddfit$Q * 100)
sqrt(ppfit$Q)
rbind(ddfit$fixed_effects, ppfit$fixed_effects)
# }
Run the code above in your browser using DataLab