library(SemiParBIVProbit)
############
## EXAMPLE 1
## Generate data
## Correlation between the two equations 0.5 - Sample size 400
set.seed(0)
n <- 400
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
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)
y1 <- -1.55 + 2*x1 + f1(x2) + u[,1]
y2 <- -0.25 - 1.25*x1 + f2(x2) + u[,2]
dataSim <- data.frame(y1, y2, x1, x2, x3)
resp.check(y1, "N")
resp.check(y2, "N")
eq.mu.1 <- y1 ~ x1 + s(x2) + s(x3)
eq.mu.2 <- y2 ~ x1 + s(x2) + s(x3)
eq.sigma2.1 <- ~ 1
eq.sigma2.2 <- ~ 1
eq.theta <- ~ x1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2.1, eq.sigma2.2, eq.theta)
# the order above is the one to follow when
# using more than two equations
out <- copulaReg(fl, data = dataSim)
conv.check(out)
post.check(out)
summary(out)
AIC(out)
BIC(out)
jc.probs(out, 1.4, 2.3, intervals = TRUE)[1:4,]
############
## EXAMPLE 2
############
## Generate data with one endogenous binary variable
## and continuous outcome
set.seed(0)
n <- 1000
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
cov <- rMVN(n, rep(0,2), Sigma)
cov <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
f2 <- function(x) x+exp(-30*(x-0.5)^2)
y1 <- ifelse(-1.55 + 2*x1 + f1(x2) + u[,1] > 0, 1, 0)
y2 <- -0.25 - 1.25*y1 + f2(x2) + u[,2]
dataSim <- data.frame(y1, y2, x1, x2)
## RECURSIVE Model
rc <- resp.check(y2, margin = "N", print.par = TRUE, loglik = TRUE)
AIC(rc); BIC(rc)
out <- copulaReg(list(y1 ~ x1 + x2,
y2 ~ y1 + x2),
data = dataSim, margins = c("probit","N"))
conv.check(out)
summary(out)
post.check(out)
## SEMIPARAMETRIC RECURSIVE Model
eq.mu.1 <- y1 ~ x1 + s(x2)
eq.mu.2 <- y2 ~ y1 + s(x2)
eq.sigma2 <- ~ 1
eq.theta <- ~ 1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)
out <- copulaReg(fl, data = dataSim,
margins = c("probit","N"), gamlssfit = TRUE)
conv.check(out)
summary(out)
post.check(out)
jc.probs(out, 1, 1.5, intervals = TRUE)[1:4,]
AT(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")
#
#
############
## EXAMPLE 3
############
## Generate data with one endogenous continuous exposure
## and binary outcome
set.seed(0)
n <- 1000
Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u <- rMVN(n, rep(0,2), Sigma)
cov <- rMVN(n, rep(0,2), Sigma)
cov <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
f2 <- function(x) x+exp(-30*(x-0.5)^2)
y1 <- -0.25 - 2*x1 + f2(x2) + u[,2]
y2 <- ifelse(-0.25 - 0.25*y1 + f1(x2) + u[,1] > 0, 1, 0)
dataSim <- data.frame(y1, y2, x1, x2)
eq.mu.1 <- y2 ~ y1 + s(x2)
eq.mu.2 <- y1 ~ x1 + s(x2)
eq.sigma2 <- ~ 1
eq.theta <- ~ 1
fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)
out <- copulaReg(fl, data = dataSim,
margins = c("probit","N"))
conv.check(out)
summary(out)
post.check(out)
AT(out, nm.end = "y1")
AT(out, nm.end = "y1", type = "univariate")
RR(out, nm.end = "y1", rr.plot = TRUE)
RR(out, nm.end = "y1", type = "univariate")
OR(out, nm.end = "y1", or.plot = TRUE)
OR(out, nm.end = "y1", type = "univariate")
#
#
############
## EXAMPLE 4
############
## Survival model
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
}
delta1 <- ifelse(t < c, 1, 0)
u1 <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u1, delta1, z1, z2)
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
}
delta2 <- ifelse(t < c,1, 0)
u2 <- apply(cbind(t, c), 1, min)
dataSim$delta2 <- delta2
dataSim$u2 <- u2
eq1 <- u1 ~ z1 + s(z2) + s(u1, bs = "mpi")
eq2 <- u2 ~ z + s(u2, bs = "mpi")
eq3 <- ~ s(z2)
out <- copulaReg(list(eq1, eq2), data = dataSim, surv = TRUE,
margins = c("PH", "PO"),
cens1 = delta1, cens2 = delta2)
# PH margin fit can also be compared with cox.ph from mgcv
conv.check(out)
post.check(out)
summary(out)
AIC(out); BIC(out)
plot(out, eq = 1, scale = 0, pages = 1)
plot(out, eq = 2, scale = 0, pages = 1)
hazsurv.plot(out, eq = 1, newdata = data.frame(z1 = 0, z2 = 0),
shade = TRUE, n.sim = 1000)
hazsurv.plot(out, eq = 1, newdata = data.frame(z1 = 0, z2 = 0),
shade = TRUE, n.sim = 1000, type = "hazard")
hazsurv.plot(out, eq = 2, newdata = data.frame(z = 0),
shade = TRUE, n.sim = 1000)
hazsurv.plot(out, eq = 2, newdata = data.frame(z = 0),
shade = TRUE, n.sim = 1000, type = "hazard")
out1 <- copulaReg(list(eq1, eq2, eq3), data = dataSim, surv = TRUE,
margins = c("PH", "PO"),
cens1 = delta1, cens2 = delta2, gamlssfit = TRUE)
eq1 <- u1 ~ z1 + s(z2)
eq2 <- u2 ~ z
eq3 <- ~ s(z2)
# note that Weibull is implemented as AFT model
out2 <- copulaReg(list(eq1, eq2, ~ 1, ~ 1, eq3), data = dataSim, surv = TRUE,
margins = c("WEI", "WEI"),
cens1 = delta1, cens2 = delta2)
Run the code above in your browser using DataLab