if (FALSE) {
library(GJRM)
set.seed(0)
n <- 400
x1 <- round(runif(n))
x2 <- runif(n)
x3 <- runif(n)
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
y1 <- -1.55 + 2*x1 + f1(x2) + rnorm(n)
dataSim <- data.frame(y1, x1, x2, x3)
resp.check(y1, "N")
eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s <- ~ s(x3)
fl <- list(eq.mu, eq.s)
out <- gamlss(fl, data = dataSim)
conv.check(out)
post.check(out)
plot(out, eq = 1, scale = 0, pages = 1, seWithMean = TRUE)
plot(out, eq = 2, seWithMean = TRUE)
summary(out)
AIC(out)
BIC(out)
################
# Robust example
################
eq.mu <- y1 ~ x1 + x2 + x3
fl <- list(eq.mu)
out <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE,
rc = 3, lB = -Inf, uB = Inf)
conv.check(out)
summary(out)
rob.const(out, 100)
##
eq.s <- ~ x3
fl <- list(eq.mu, eq.s)
out <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE)
conv.check(out)
summary(out)
##
eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s <- ~ s(x3)
fl <- list(eq.mu, eq.s)
out1 <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE,
sp.method = "efs")
conv.check(out1)
summary(out1)
AIC(out, out1)
plot(out1, eq = 1, all.terms = TRUE, pages = 1, seWithMean = TRUE)
plot(out1, eq = 2, seWithMean = TRUE)
##########################
## GEV link binary example
##########################
# this incorporates the bgeva
# model implemented in the bgeva package
# however this implementation is more general
# stable and efficient
set.seed(0)
n <- 400
x1 <- round(runif(n)); x2 <- runif(n); x3 <- runif(n)
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
f2 <- function(x) x+exp(-30*(x-0.5)^2)
y <- ifelse(-3.55 + 2*x1 + f1(x2) + rnorm(n) > 0, 1, 0)
dataSim <- data.frame(y, x1, x2, x3)
out1 <- gamlss(list(y ~ x1 + x2 + x3), margin = "GEVlink", data = dataSim)
out2 <- gamlss(list(y ~ x1 + s(x2) + s(x3)), margin = "GEVlink", data = dataSim)
conv.check(out1)
conv.check(out2)
summary(out1)
summary(out2)
AIC(out1, out2)
BIC(out1, out2)
plot(out2, eq = 1, all.terms = TRUE, pages = 1, seWithMean = TRUE)
##################
# prediction of Pr
##################
# Calculate eta (that is, X*model.coef)
# For a new data set the argument newdata should be used
eta <- predict(out2, eq = 1, type = "link")
# extract gev tail parameter
gev.par <- out2$gev.par
# multiply gev tail parameter by eta
gevpeta <- gev.par*eta
# establish for which values the model is defined
gevpetaIND <- ifelse(gevpeta < -1, FALSE, TRUE)
gevpeta <- gevpeta[gevpetaIND]
# estimate probabilities
pr <- exp(-(1 + gevpeta)^(-1/gev.par))
###################################
## Flexible survival model examples
###################################
library(GJRM)
########################################
## Simulate proportional hazards data ##
########################################
set.seed(0)
n <- 2000
c <- runif(n, 3, 8)
u <- runif(n, 0, 1)
z1 <- rbinom(n, 1, 0.5)
z2 <- runif(n, 0, 1)
t <- rep(NA, n)
beta_0 <- -0.2357
beta_1 <- 1
f <- function(t, beta_0, beta_1, u, z1, z2){
S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
exp(-exp(log(-log(S_0))+beta_0*z1 + beta_1*z2))-u
}
for (i in 1:n){
t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5,
beta_0 = beta_0, beta_1 = beta_1, u = u[i],
z1 = z1[i], z2 = z2[i], extendInt = "yes" )$root
}
delta <- ifelse(t < c, 1, 0)
u <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u, delta, z1, z2)
1-mean(delta) # average censoring rate
# log(u) helps obtaining smoother hazards
out <- gamlss(list(u ~ s(log(u), bs = "mpi") + z1 + s(z2) ), data = dataSim,
surv = TRUE, margin = "PH", cens = delta)
post.check(out)
summary(out)
AIC(out)
BIC(out)
plot(out, eq = 1, scale = 0, pages = 1)
hazsurv(out, newdata = data.frame(z1 = 0, z2 = 0), shade = TRUE,
n.sim = 1000, baseline = TRUE)
hazsurv(out, type = "hazard", newdata = data.frame(z1 = 0, z2 = 0),
shade = TRUE, n.sim = 1000, baseline = TRUE)
out1 <- gam(u ~ z1 + s(z2), family = cox.ph(),
data = dataSim, weights = delta)
summary(out1)
# estimates of z1 and s(z2) are
# nearly identical between out and out1
# note that the Weibull is implemented as AFT
# as using the PH parametrisation makes
# computation unstable
out2 <- gamlss(list(u ~ z1 + s(z2) ), data = dataSim, surv = TRUE,
margin = "WEI", cens = delta)
#####################################
## Simulate proportional odds data ##
#####################################
set.seed(0)
n <- 2000
c <- runif(n, 4, 8)
u <- runif(n, 0, 1)
z <- rbinom(n, 1, 0.5)
beta_0 <- -1.05
t <- rep(NA, n)
f <- function(t, beta_0, u, z){
S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
1/(1 + exp(log((1-S_0)/S_0)+beta_0*z))-u
}
for (i in 1:n){
t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5,
beta_0 = beta_0, u = u[i], z = z[i],
extendInt="yes" )$root
}
delta <- ifelse(t < c,1, 0)
u <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u, delta, z)
1-mean(delta) # average censoring rate
out <- gamlss(list(u ~ s(log(u), bs = "mpi") + z ), data = dataSim, surv = TRUE,
margin = "PO", cens = delta)
post.check(out)
summary(out)
AIC(out)
BIC(out)
plot(out, eq = 1, scale = 0)
hazsurv(out, newdata = data.frame(z = 0), shade = TRUE, n.sim = 1000,
baseline = TRUE)
hazsurv(out, type = "hazard", newdata = data.frame(z = 0),
shade = TRUE, n.sim = 1000)
#############################
## Mixed censoring example ##
#############################
f1 <- function(t, u, z1, z2, z3, z4, s1, s2){
S_0 <- 0.7 * exp(-0.03*t^1.8) + 0.3*exp(-0.3*t^2.5)
exp( -exp(log(-log(S_0)) + 1.3*z1 + 0.5*z2 + s1(z3) + s2(z4) ) ) - u
}
datagen <- function(n, z1, z2, z3, z4, s1, s2, f1){
u <- runif(n, 0, 1)
t <- rep(NA, n)
for (i in 1:n) t[i] <- uniroot(f1, c(0, 100), tol = .Machine$double.eps^0.5,
u = u[i], s1 = s1, s2 = s2, z1 = z1[i], z2 = z2[i],
z3 = z3[i], z4 = z4[i], extendInt = "yes")$root
c1 <- runif(n, 0, 2)
c2 <- c1 + runif(n, 0, 6)
df <- data.frame(u1 = t, u2 = t, cens = character(n), stringsAsFactors = FALSE)
for (i in 1:n){
if(t[i] <= c1[i]) {
df[i, 1] <- c1[i]
df[i, 2] <- NA
df[i, 3] <- "L"
}else if(c1[i] < t[i] && t[i] <= c2[i]){
df[i, 1] <- c1[i]
df[i, 2] <- c2[i]
df[i, 3] <- "I"
}else if(t[i] > c2[i]){
df[i, 1] <- c2[i]
df[i, 2] <- NA
df[i, 3] <- "R"}
}
uncens <- (df[, 3] %in% c("L", "I")) + (rbinom(n, 1, 0.2) == 1) == 2
df[uncens, 1] <- t[uncens]
df[uncens, 2] <- NA
df[uncens, 3] <- "U"
dataSim <- data.frame(u1 = df$u1, u2 = df$u2, cens = as.factor(df$cens), z1, z2, z3, z4, t)
dataSim
}
set.seed(0)
n <- 1000
SigmaC <- matrix(0.5, 4, 4); diag(SigmaC) <- 1
cov <- rMVN(n, rep(0,4), SigmaC)
cov <- pnorm(cov)
z1 <- round(cov[, 1])
z2 <- round(cov[, 2])
z3 <- cov[, 3]
z4 <- cov[, 4]
s1 <- function(x) -0.075*exp(3.2 * x)
s2 <- function(x) sin(2*pi*x)
eq1 <- u1 ~ s(log(u1), bs = "mpi") + z1 + z2 + s(z3) + s(z4)
dataSim <- datagen(n, z1, z2, z3, z4, s1, s2, f1)
out <- gamlss(list(eq1), data = dataSim, surv = TRUE, margin = "PH",
cens = cens, type.cen = "mixed", upperB = "u2")
conv.check(out)
summary(out)
plot(out, eq = 1, scale = 0, pages = 1)
ndf <- data.frame(z1 = 1, z2 = 0, z3 = 0.2, z4 = 0.5)
hazsurv(out, eq = 1, newdata = ndf, type = "surv")
hazsurv(out, eq = 1, newdata = ndf, type = "hazard", n.sim = 1000)
}
Run the code above in your browser using DataLab