# NOT RUN {
# }
# NOT RUN {
library(JRM)
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.s2 <- ~ s(x3)
fl <- list(eq.mu, eq.s2)
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)
##
eq.s2 <- ~ x3
fl <- list(eq.mu, eq.s2)
out <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE)
conv.check(out)
summary(out)
##
eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s2 <- ~ s(x3)
fl <- list(eq.mu, eq.s2)
out1 <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE)
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 example
##################################
library(JRM)
########################################
## 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
out <- gamlss(list(u ~ z1 + s(z2) + s(u, bs = "mpi") ), 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.plot(out, newdata = data.frame(z1 = 0, z2 = 0), shade = TRUE, n.sim = 1000)
hazsurv.plot(out, type = "hazard", newdata = data.frame(z1 = 0, z2 = 0),
shade = TRUE, n.sim = 1000)
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 ~ z + s(u, bs = "mpi") ), data = dataSim, surv = TRUE,
margin = "PO", cens = delta)
post.check(out)
summary(out)
AIC(out)
BIC(out)
plot(out, eq = 1, scale = 0)
hazsurv.plot(out, newdata = data.frame(z =0), shade = TRUE, n.sim = 1000)
hazsurv.plot(out, type = "hazard", newdata = data.frame(z = 0),
shade = TRUE, n.sim = 1000)
# note that the Fisk is implemented as AFT
# as using the PH parametrisation makes
# computation unstable
out1 <- gamlss(list(u ~ z), data = dataSim, surv = TRUE,
margin = "FISK", cens = delta)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab