Methods for targeted and semiparametric inference.

Author

Klaus K. Holst (Maintainer) klaus@holst.it

Examples

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+ }