Methods for targeted and semiparametric inference.
example(riskreg)
#>
#> riskrg> m <- lvm(a[-2] ~ x,
#> riskrg+ z ~ 1,
#> riskrg+ lp.target[1] ~ 1,
#> riskrg+ lp.nuisance[-1] ~ 2*x)
#>
#> riskrg> distribution(m,~a) <- binomial.lvm("logit")
#>
#> riskrg> m <- binomial.rr(m, "y","a","lp.target","lp.nuisance")
#>
#> riskrg> d <- sim(m,5e2,seed=1)
#>
#> riskrg> I <- model.matrix(~1, d)
#>
#> riskrg> X <- model.matrix(~1+x, d)
#>
#> riskrg> with(d, riskreg_mle(y, a, I, X, type="rr"))
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9453 0.1150 0.7198 1.1707 2.092e-16
#> odds-product:(Intercept) -0.8068 0.3031 -1.4009 -0.2126 7.783e-03
#> odds-product:x 2.3296 0.3914 1.5624 3.0969 2.659e-09
#>
#> riskrg> with(d, riskreg_fit(y, a, nuisance=X, propensity=I, type="rr"))
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9469 0.113 0.7254 1.168 5.418e-17
#>
#> riskrg> riskreg(y ~ a | 1, nuisance=~x , data=d, type="rr")
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9469 0.113 0.7254 1.168 5.418e-17
#>
#> riskrg> ## Model with same design matrix for nuisance and propensity model:
#> riskrg> with(d, riskreg_fit(y, a, nuisance=X, type="rr"))
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9402 0.1173 0.7104 1.17 1.082e-15
#>
#> riskrg> ## a <- riskreg(y ~ a, target=~z, nuisance=~x, propensity=~x, data=d, type="rr")
#> riskrg> a <- riskreg(y ~ a | z, nuisance=~x, propensity=~x, data=d, type="rr")
#>
#> riskrg> a
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9797 0.1212 0.7421 1.217291 6.365e-16
#> z -0.2218 0.1115 -0.4404 -0.003204 4.674e-02
#>
#> riskrg> predict(a, d[1:5,])
#> [1] 0.1491919 0.6897799 0.1040211 0.3736667 0.3304784
#>
#> riskrg> riskreg(y ~ a, nuisance=~x, data=d, type="rr", mle=TRUE)
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9453 0.1150 0.7198 1.1707 2.092e-16
#> odds-product:(Intercept) -0.8068 0.3031 -1.4009 -0.2126 7.783e-03
#> odds-product:x 2.3296 0.3914 1.5624 3.0969 2.659e-09
example(cate)
#>
#> cate> sim1 <- function(n=1e4,
#> cate+ seed=NULL,
#> cate+ return_model=FALSE, ...) {
#> cate+ suppressPackageStartupMessages(require("lava"))
#> cate+ if (!is.null(seed)) set.seed(seed)
#> cate+ m <- lava::lvm()
#> cate+ regression(m, ~a) <- function(z1,z2,z3,z4,z5)
#> cate+ cos(z1)+sin(z1*z2)+z3+z4+z5^2
#> cate+ regression(m, ~u) <- function(a,z1,z2,z3,z4,z5)
#> cate+ (z1+z2+z3)*a + z1+z2+z3 + a
#> cate+ distribution(m, ~a) <- binomial.lvm()
#> cate+ if (return_model) return(m)
#> cate+ lava::sim(m, n, p=par)
#> cate+ }
#>
#> cate> d <- sim1(200)
#>
#> cate> e <- cate(a ~ z1+z2+z3, response=u~., data=d)
#>
#> cate> e
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.02724 0.4805 -0.9144 0.9689 9.548e-01
#> z1 2.08298 0.5099 1.0836 3.0824 4.407e-05
#> z2 0.89715 0.3777 0.1568 1.6375 1.754e-02
#> z3 0.89901 0.6210 -0.3181 2.1161 1.477e-01
example(ate)
#>
#> ate> m <- lvm(y ~ a+x, a~x)
#>
#> ate> distribution(m,~ a+y) <- binomial.lvm()
#>
#> ate> d <- sim(m,1e4,seed=1)
#>
#> ate> a <- ate(y ~ a, nuisance=~x, data=d)
#>
#> ate> summary(a)
#>
#> Augmented Inverse Probability Weighting estimator
#> Response y (Outcome model: gaussian):
#> y ~ x + a
#> Exposure a (Propensity model: logistic regression):
#> a ~ x
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> a=0 0.496843 0.007478 0.48219 0.51150 0.000e+00
#> a=1 0.706396 0.007364 0.69196 0.72083 0.000e+00
#> Outcome model:
#> (Intercept) 0.485158 0.006898 0.47164 0.49868 0.000e+00
#> x 0.175278 0.004231 0.16699 0.18357 0.000e+00
#> a 0.224335 0.010050 0.20464 0.24403 2.287e-110
#> Propensity model:
#> (Intercept) -0.001888 0.022109 -0.04522 0.04144 9.320e-01
#> x -1.019511 0.026347 -1.07115 -0.96787 0.000e+00
#>
#> Average Treatment Effect (constrast: 'a=1' - 'a=0'):
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> ATE 0.2096 0.01013 0.1897 0.2294 4.494e-95
#>
#>
#> ate> b <- cate(a ~ 1, y ~ a*x, a ~ x, data=d)
#>
#> ate> # Multiple treatments
#> ate> m <- lvm(y ~ a+x, a~x)
#>
#> ate> distribution(m,~ y) <- binomial.lvm()
#>
#> ate> m <- ordinal(m, K=4, ~a)
#>
#> ate> transform(m, ~a) <- factor
#>
#> ate> d <- sim(m, 1e4, seed=1)
#>
#> ate> (a <- ate(y~a|a*x|x, data=d))
#> Estimate Std.Err 2.5% 97.5% P-value
#> a=0 0.2121 0.01695 0.1788 0.2453 6.738e-36
#> a=1 0.3422 0.01831 0.3063 0.3781 6.005e-78
#> a=2 0.3922 0.01617 0.3605 0.4239 6.820e-130
#> a=3 0.6052 0.00771 0.5901 0.6203 0.000e+00
#>
#> ate> # Comparison with randomized experiment
#> ate> m0 <- cancel(m, a~x)
#>
#> ate> d0 <- sim(m0,2e5)
#>
#> ate> lm(y~a-1,d0)
#>
#> Call:
#> lm(formula = y ~ a - 1, data = d0)
#>
#> Coefficients:
#> a0 a1 a2 a3
#> 0.2245 0.3405 0.4111 0.6206
#>
#>
#> ate> # Choosing a different contrast for the association measures
#> ate> summary(a, contrast=c(2,4))
#>
#> Augmented Inverse Probability Weighting estimator
#> Response y (Outcome model: gaussian):
#> y ~ a * x
#> Exposure a (Propensity model: logistic regression):
#> a ~ x
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> a=0 0.2121 0.01695 0.1788 0.2453 6.738e-36
#> a=1 0.3422 0.01831 0.3063 0.3781 6.005e-78
#> a=2 0.3922 0.01617 0.3605 0.4239 6.820e-130
#> a=3 0.6052 0.00771 0.5901 0.6203 0.000e+00
#>
#> Average Treatment Effect (constrast: 'a=1' - 'a=3'):
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> ATE -0.263 0.01955 -0.3013 -0.2247 3.016e-41
#>
example(calibration)
#>
#> clbrtn> sim1 <- function(n, beta=c(-3, rep(.5,10)), rho=.5) {
#> clbrtn+ p <- length(beta)-1
#> clbrtn+ xx <- lava::rmvn0(n,sigma=diag(nrow=p)*(1-rho)+rho)
#> clbrtn+ y <- rbinom(n, 1, lava::expit(cbind(1,xx)%*%beta))
#> clbrtn+ d <- data.frame(y=y, xx)
#> clbrtn+ names(d) <- c("y",paste0("x",1:p))
#> clbrtn+ return(d)
#> clbrtn+ }
#>
#> clbrtn> set.seed(1)
#>
#> clbrtn> beta <- c(-2,rep(1,10))
#>
#> clbrtn> d <- sim1(1e4, beta=beta)
#>
#> clbrtn> a1 <- NB(y ~ ., data=d)
#>
#> clbrtn> a2 <- glm(y ~ ., data=d, family=binomial)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#>
#> clbrtn> ## a3 <- randomForest(factor(y) ~ ., data=d, family=binomial)
#> clbrtn>
#> clbrtn> d0 <- sim1(1e5, beta=beta)
#>
#> clbrtn> p1 <- predict(a1, newdata=d0)
#>
#> clbrtn> p2 <- predict(a2, newdata=d0, type="response")
#>
#> clbrtn> ## p3 <- predict(a3, newdata=d0, type="prob")
#> clbrtn>
#> clbrtn> c2 <- calibration(p2, d0$y, method="isotonic")
#>
#> clbrtn> c1 <- calibration(p1, d0$y, breaks=100)
#>
#> clbrtn> if (interactive()) {
#> clbrtn+ plot(c1)
#> clbrtn+ plot(c2,col="red",add=TRUE)
#> clbrtn+ abline(a=0,b=1)##'
#> clbrtn+ with(c1$xy[[1]], points(pred,freq,type="b", col="red"))
#> clbrtn+ }
#>
#> clbrtn> set.seed(1)
#>
#> clbrtn> beta <- c(-2,rep(1,10))
#>
#> clbrtn> dd <- lava::csplit(sim1(6e4, beta=beta), k=3)
#>
#> clbrtn> mod <- NB(y ~ ., data=dd[[1]])
#>
#> clbrtn> p1 <- predict(mod, newdata=dd[[2]])
#>
#> clbrtn> cal <- calibration(p1, dd[[2]]$y)
#>
#> clbrtn> p2 <- predict(mod, newdata=dd[[3]])
#>
#> clbrtn> pp <- predict(c1, p2)
#>
#> clbrtn> cc <- calibration(pp, dd[[3]]$y)
#>
#> clbrtn> if (interactive()) {##'
#> clbrtn+ plot(cal)
#> clbrtn+ plot(cc, add=TRUE, col="blue")
#> clbrtn+ }