Nothing
## ----knitup, include = FALSE---------------------------------------------
library(knitr)
require(printr, quietly = TRUE)
set.seed(225452)
## ----eval = FALSE--------------------------------------------------------
# devtools::install_github("sachsmc/pseval")
## ----eval = FALSE--------------------------------------------------------
# install.packages("pseval")
## ----setup---------------------------------------------------------------
library(pseval)
library(survival)
## ------------------------------------------------------------------------
set.seed(1492)
fakedata <- generate_example_data(n = 800)
head(fakedata)
## ----psdesign------------------------------------------------------------
binary.ps <- psdesign(data = fakedata, Z = Z, Y = Y.obs, S = S.obs, BIP = BIP)
binary.ps
## ----casecont------------------------------------------------------------
fakedata.cc <- fakedata
missdex <- sample((1:nrow(fakedata.cc))[fakedata.cc$Y.obs == 0],
size = floor(sum(fakedata.cc$Y.obs == 0) * .8))
fakedata.cc[missdex, ]$S.obs <- NA
fakedata.cc$weights <- ifelse(fakedata.cc$Y.obs == 1, 1, .2)
## ----casecont2-----------------------------------------------------------
binary.cc <- psdesign(data = fakedata.cc, Z = Z, Y = Y.obs, S = S.obs,
BIP = BIP, weights = weights)
## ----survex--------------------------------------------------------------
surv.ps <- psdesign(data = fakedata, Z = Z, Y = Surv(time.obs, event.obs), S = S.obs,
BIP = BIP, CPV = CPV, BSM = BSM)
## ----inthelp-------------------------------------------------------------
?add_integration
## ----impp----------------------------------------------------------------
binary.ps <- binary.ps + integrate_parametric(S.1 ~ BIP)
binary.ps
## ----imp2----------------------------------------------------------------
binary.ps + integrate_parametric(S.0 ~ BIP)
## ----imp3, eval = FALSE--------------------------------------------------
# library(splines)
# binary.ps + integrate_parametric(S.1 ~ BIP^2)
# binary.ps + integrate_parametric(S.1 ~ bs(BIP, df = 3))
# binary.ps + integrate_parametric(S.1 ~ BIP + BSM + BSM^2)
## ----riskhelp------------------------------------------------------------
?add_riskmodel
## ----riskbin-------------------------------------------------------------
binary.ps <- binary.ps + risk_binary(model = Y ~ S.1 * Z, D = 5, risk = risk.logit)
binary.ps
## ----est-----------------------------------------------------------------
binary.est <- binary.ps + ps_estimate(method = "BFGS")
binary.boot <- binary.est + ps_bootstrap(n.boots = 10, progress.bar = FALSE,
start = binary.est$estimates$par, method = "BFGS")
binary.boot
## ----alltog, eval = FALSE------------------------------------------------
# binary.est <- psdesign(data = fakedata, Z = Z, Y = Y.obs, S = S.obs, BIP = BIP) +
# integrate_parametric(S.1 ~ BIP) +
# risk_binary(model = Y ~ S.1 * Z, D = 50, risk = risk.logit) +
# ps_estimate(method = "BFGS")
## ----stg, eval = FALSE---------------------------------------------------
# calc_STG(binary.boot, progress.bar = FALSE)
## ----summary-------------------------------------------------------------
smary <- summary(binary.boot)
## ----calcrisk------------------------------------------------------------
head(calc_risk(binary.boot, contrast = "TE", n.samps = 20), 3)
head(calc_risk(binary.boot, contrast = function(R0, R1) 1 - R1/R0, n.samps = 20), 3)
## ----plot1, fig.cap = "Plot showing the estimates using the example data, along with confidence bands (CB), and the true treatment efficacy (TE) curve."----
plot(binary.boot, contrast = "TE", lwd = 2)
abline(h = smary$TE.estimates[2], lty = 3)
expit <- function(x) exp(x)/(1 + exp(x))
trueTE <- function(s){
r0 <- expit(-1 - 0 * s)
r1 <- expit(-1 - .75 * s)
1 - r1/r0
}
rug(binary.boot$augdata$S.1)
curve(trueTE(x), add = TRUE, col = "red")
legend("bottomright", legend = c("estimated TE", "95\\% CB", "marginal TE", "true TE"),
col = c("black", "black", "black", "red"), lty = c(1, 2, 3, 1), lwd = c(2, 2, 1, 1))
## ----logrr, echo = -1, fig.cap = "Plot illustrating ways that different risk contrast functions can be plotted. "----
par(mfrow = c(2, 2), mar = c(4, 4, 1, .5))
plot(binary.boot, contrast = "logRR", lwd = 2, col = c("black", "grey75", "grey75"))
plot(binary.boot, contrast = "RR", log = "y", lwd = 2, col = c("black", "grey75", "grey75"))
plot(binary.boot, contrast = "RD", lwd = 2, col = c("black", "grey75", "grey75"))
plot(binary.boot, contrast = "risk", lwd = 2, lty = c(1, 0, 0, 2, 0, 0))
legend("topright", legend = c("R0", "R1"), lty = c(1, 2), lwd = 2)
## ----ve------------------------------------------------------------------
te.est <- calc_risk(binary.boot, CI.type = "pointwise", n.samps = 200)
head(te.est, 3)
## ----plotciag, eval = FALSE----------------------------------------------
# plot(binary.boot, contrast = "TE", lwd = 2, CI.type = "band")
# sbs <- calc_risk(binary.boot, CI.type = "pointwise", n.samps = 200)
# lines(Y.lower.CL.2.5 ~ S.1, data = sbs, lty = 3, lwd = 2)
# lines(Y.upper.CL.97.5 ~ S.1, data = sbs, lty = 3, lwd = 2)
# legend("bottomright", lwd = 2, lty = 1:3,
# legend = c("estimate", "simultaneous CI", "pointwise CI"))
## ----ggpt----------------------------------------------------------------
library(ggplot2)
TE.est <- calc_risk(binary.boot, n.samps = 200)
ggplot(TE.est,
aes(x = S.1, y = Y, ymin = Y.lower.CL.0.95, ymax = Y.upper.CL.0.95)) +
geom_line() + geom_ribbon(alpha = .2) + ylab(attr(TE.est, "Y.function"))
## ----ccest, eval = FALSE-------------------------------------------------
# cc.fit <- binary.cc + integrate_parametric(S.1 ~ BIP) +
# risk_binary(D = 10) + ps_estimate()
# cc.fit
## ----surv1, eval = FALSE-------------------------------------------------
# surv.fit <- psdesign(fakedata, Z = Z, Y = Surv(time.obs, event.obs),
# S = S.obs, BIP = BIP, CPV = CPV) +
# integrate_semiparametric(formula.location = S.1 ~ BIP, formula.scale = S.1 ~ 1) +
# risk_exponential(D = 10) + ps_estimate(method = "BFGS") + ps_bootstrap(n.boots = 20)
# surv.fit
# plot(surv.fit)
## ----cont1---------------------------------------------------------------
fakedata$Y.cont <- log(fakedata$time.obs + 0.01)
cont.fit <- psdesign(fakedata, Z = Z, Y = Y.cont,
S = S.obs, BIP = BIP, CPV = CPV) +
integrate_semiparametric(formula.location = S.1 ~ BIP, formula.scale = S.1 ~ 1) +
risk_continuous(D = 10) + ps_estimate(method = "BFGS") #+ ps_bootstrap(n.boots = 20)
cont.fit
plot(cont.fit, contrast = "risk")
## ------------------------------------------------------------------------
with(fakedata, table(S.obs.cat, BIP.cat))
## ----catfit--------------------------------------------------------------
cat.fit <- psdesign(fakedata, Z = Z, Y = Y.obs,
S = S.obs.cat, BIP = BIP.cat) +
integrate_nonparametric(formula = S.1 ~ BIP) +
risk_binary(Y ~ S.1 * Z, D = 10, risk = risk.probit) + ps_estimate(method = "BFGS")
cat.fit
plot(cat.fit)
## ----catfitps------------------------------------------------------------
cat.fit.ps <- psdesign(fakedata, Z = Z, Y = Y.obs,
S = S.obs, BIP = BIP.cat) +
integrate_nonparametric(formula = S.1 ~ BIP) +
risk_binary(Y ~ S.1 * Z, D = 10, risk = risk.logit) + ps_estimate(method = "pseudo-score") +
ps_bootstrap(n.boots = 20, method = "pseudo-score")
summary(cat.fit.ps)
plot(cat.fit.ps)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.