minr bug fixes to npOR

R/Param_npCATE.R

@@ -1,4 +1,8 @@
-#' Average Treatment Effect
+#' Nonparametric inference for user-specified parametric working models for the conditional treatment effect.
+#' The true conditional average treatment effect is projected onto a parametric working model using least-squares regression.
+#' Unlike \code{Param_npCATT}, this function uses all observations to compute the projection.
+#' This can be used to assess heterogeneity of the average treatment effect.
+#' We note that `formula_CATE = ~ 1` gives an estimator of the nonparametric average treatment effect (ATE).
 #'
 #' Parameter definition for the Average Treatment Effect (ATE).
 #' @importFrom R6 R6Class

R/Param_npCATT.R

@@ -1,6 +1,7 @@
-#' Average Treatment Effect
-#'
-#' Parameter definition for the Average Treatment Effect (ATE).
+#' Nonparametric inference for user-specified parametric working models for the conditional treatment effect.
+#' The true conditional average treatment effect is projected onto a parametric working model using only individuals with `A=1` (among the treated).
+#' This can be used to assess heterogeneity of the average treatment effect and avoids positivity issues by focusing on best approximating the conditional average treatment effect ampng the treated.
+#' We note that `formula_CATT = ~ 1` gives an estimator of the nonparametric average treatment effect among the treated (ATT).
 #' @importFrom R6 R6Class
 #' @importFrom uuid UUIDgenerate
 #' @importFrom methods is

R/Param_npOR.R

@@ -1,5 +1,8 @@
-#' Average Treatment Effect
-#'
+#' Nonparametric inference for user-specified parametric working models for the conditional odds ratio between two binary variables
+#' The true conditional odds ratio is projected onto a parametric working model using logistic-regression.
+#' This can be used to assess heterogeneity of the odds ratio.
+#' We note that `formula_logOR = ~ 1` gives an estimator of the nonparametric marginal odds ratio among the treated.
+#' The parametric model is at the log-scale and therefore the coefficients returned code the linear predictor for the `log`-conditional odds ratio.
 #' Parameter definition for the Average Treatment Effect (ATE).
 #' @importFrom R6 R6Class
 #' @importFrom uuid UUIDgenerate
@@ -17,7 +20,7 @@
 #'   \describe{
 #'     \item{\code{observed_likelihood}}{A \code{\link{Likelihood}} corresponding to the observed likelihood
 #'     }
-#'     \item{\code{formula_OR}}{...
+#'     \item{\code{formula_logOR}}{...
 #'     }
 #'     \item{\code{intervention_list_treatment}}{A list of objects inheriting from \code{\link{LF_base}}, representing the treatment intervention.
 #'     }
@@ -48,7 +51,7 @@ Param_npOR <- R6Class(
   class = TRUE,
   inherit = Param_base,
   public = list(
-    initialize = function(observed_likelihood,  formula_OR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
+    initialize = function(observed_likelihood,  formula_logOR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
       super$initialize(observed_likelihood, list(), outcome_node)
       if (!is.null(observed_likelihood$censoring_nodes[[outcome_node]])) {
         # add delta_Y=0 to intervention lists
@@ -57,7 +60,7 @@ Param_npOR <- R6Class(
         intervention_list_treatment <- c(intervention_list_treatment, censoring_intervention)
         intervention_list_control <- c(intervention_list_control, censoring_intervention)
       }
-      private$.formula_OR <- formula_OR
+      private$.formula_logOR <- formula_logOR
       private$.cf_likelihood_treatment <- CF_Likelihood$new(observed_likelihood, intervention_list_treatment)
       private$.cf_likelihood_control <- CF_Likelihood$new(observed_likelihood, intervention_list_control)
     },
@@ -75,11 +78,11 @@ Param_npOR <- R6Class(
       intervention_nodes <- union(names(self$intervention_list_treatment), names(self$intervention_list_control))
 
       W <- tmle_task$get_tmle_node("W")
-      V <- model.matrix(self$formula_OR, as.data.frame(W))
+      V <- model.matrix(self$formula_logOR, as.data.frame(W))
       A <- tmle_task$get_tmle_node("A", format = TRUE)[[1]]
       Y <- tmle_task$get_tmle_node("Y", format = TRUE)[[1]]
       W_train <- training_task$get_tmle_node("W")
-      V_train <- model.matrix(self$formula_OR, as.data.frame(W_train))
+      V_train <- model.matrix(self$formula_logOR, as.data.frame(W_train))
       A_train <- training_task$get_tmle_node("A", format = TRUE)[[1]]
       Y_train <- training_task$get_tmle_node("Y", format = TRUE)[[1]]
 
@@ -93,29 +96,33 @@ Param_npOR <- R6Class(
       Qorig <- Q
       Q0 <- bound(Q0, 0.005)
       Q1 <- bound(Q1, 0.005)
-      beta <- get_beta(W_train, A_train, self$formula_OR, Q1, Q0, family = binomial(), weights = self$weights)
+      beta <- get_beta(W_train, A_train, self$formula_logOR, Q1, Q0, family = binomial(), weights = self$weights)
+
       Q1beta <- plogis(qlogis(Q0) + V%*%beta)
-      ORbeta <- Q1beta*(1-Q1beta) / (Q0*(1-Q0))
-      omega <- (g0 + g1*ORbeta) / (g0)
+
+      sigma_rel <- Q1beta*(1-Q1beta) / (Q0*(1-Q0))
 
 
-      h_star <-  -1*as.vector((g1*ORbeta) / (g1*ORbeta + (1-g1)))
+      omega <- (g0 + g1*sigma_rel) / (g0)
+
+      h_star <-  -1*as.vector((g1*sigma_rel) / (g1*sigma_rel + (1-g1)))
       H <- as.matrix(omega*V*(A  + h_star))
 
       # Store EIF component
       EIF_Y <- NULL
       EIFWA <- NULL
       if(is_training_task) {
         tryCatch({
-          scale <- apply(V,2, function(v){apply(self$weights*as.vector( Q1beta*(1-Q1beta) * Q0*(1-Q0) * g1 * (1-g1) / (g1 * Q1beta*(1-Q1beta) + (1-g1) *Q0*(1-Q0) )) * v*V,2,mean)})
+          scale <- apply(V,2, function(v){apply(self$weights*(A * Q1beta*(1-Q1beta) * v*V),2,mean)})
           scaleinv <- solve(scale)
           EIF_Y <- self$weights * (H%*% scaleinv) * (Y-Q)
           EIFWA <- apply(V, 2, function(v) {
-            (weights*(A*v*(Q1 - Q1beta)) - mean( self$weights*(A*v*(Q1 - Q1beta))))
-          }) %*% scale_inv
-        }, error = function(...){
+            (self$weights*(A*v*(Q1 - Q1beta)) - mean( self$weights*(A*v*(Q1 - Q1beta))))
+          })
+
+          EIFWA <- EIFWA %*% scaleinv
 
-        })
+        } ,error = function(...) {} )
       }
 
       return(list(Y = H, EIF = list(Y = EIF_Y, WA = EIFWA)))
@@ -147,8 +154,8 @@ Param_npOR <- R6Class(
       # Q <- Q_packed[[3]]
       Q0 <- bound(Q0, 0.0005)
       Q1 <- bound(Q1, 0.0005)
-      beta <- get_beta(W, A, self$formula_OR, Q1, Q0, family = binomial(), weights = weights)
-      V <- model.matrix(self$formula_OR, as.data.frame(W))
+      beta <- get_beta(W, A, self$formula_logOR, Q1, Q0, family = binomial(), weights = weights)
+      V <- model.matrix(self$formula_logOR, as.data.frame(W))
       OR <- exp(V%*%beta)
 
       IC <- EIF
@@ -177,16 +184,16 @@ Param_npOR <- R6Class(
     update_nodes = function() {
       return(c(self$outcome_node))
     },
-    formula_OR = function(){
-      return(private$.formula_OR)
+    formula_logOR = function(){
+      return(private$.formula_logOR)
     }
   ),
   private = list(
     .type = "OR",
     .cf_likelihood_treatment = NULL,
     .cf_likelihood_control = NULL,
     .supports_outcome_censoring = TRUE,
-    .formula_OR = NULL,
+    .formula_logOR = NULL,
     .submodel = list(Y = "gaussian_identity")
   )
 )

R/Param_spCATE.R

@@ -1,6 +1,7 @@
-#' Average Treatment Effect
-#'
-#' Parameter definition for the Average Treatment Effect (ATE).
+#' Semiparametric estimation of the conditonal average treatment effect for arbitrary partially-linear least-squares regression models.
+#' This is a semiparametric version of \code{Param_npCATT} and \code{Param_npCATE} where the parametric model for the CATE is assumed correct.
+#' Assuming the semiparametric model to be true allows for some efficiency gain (when true) but may lead to less robust estimates due to misspecification.
+#' Note a linear-link is used for the CATE.
 #' @importFrom R6 R6Class
 #' @importFrom uuid UUIDgenerate
 #' @importFrom methods is

R/Param_spOR.R

@@ -1,6 +1,7 @@
-#' Average Treatment Effect
-#'
-#' Parameter definition for the Average Treatment Effect (ATE).
+#' Semiparametric estimation of the conditonal odds ratio for arbitrary partially-linear logistic regression models.
+#' This is a semiparametric version of \code{Param_npOR} where the parametric model for the OR is assumed correct.
+#' Assuming the semiparametric model to be true allows for some efficiency gain (when true) but may lead to less robust estimates due to misspecification.
+#' The parametric model is at the log-scale and therefore the coefficients returned code the linear predictor for the `log`-conditional odds ratio.
 #' @importFrom R6 R6Class
 #' @importFrom uuid UUIDgenerate
 #' @importFrom methods is
@@ -17,7 +18,7 @@
 #'   \describe{
 #'     \item{\code{observed_likelihood}}{A \code{\link{Likelihood}} corresponding to the observed likelihood
 #'     }
-#'     \item{\code{formula_OR}}{...
+#'     \item{\code{formula_logOR}}{...
 #'     }
 #'     \item{\code{intervention_list_treatment}}{A list of objects inheriting from \code{\link{LF_base}}, representing the treatment intervention.
 #'     }
@@ -48,7 +49,7 @@ Param_spOR <- R6Class(
   class = TRUE,
   inherit = Param_base,
   public = list(
-    initialize = function(observed_likelihood,  formula_OR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
+    initialize = function(observed_likelihood,  formula_logOR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
       super$initialize(observed_likelihood, list(), outcome_node)
       if (!is.null(observed_likelihood$censoring_nodes[[outcome_node]])) {
         # add delta_Y=0 to intervention lists
@@ -57,7 +58,7 @@ Param_spOR <- R6Class(
         intervention_list_treatment <- c(intervention_list_treatment, censoring_intervention)
         intervention_list_control <- c(intervention_list_control, censoring_intervention)
       }
-      private$.formula_OR <- formula_OR
+      private$.formula_logOR <- formula_logOR
       private$.cf_likelihood_treatment <- CF_Likelihood$new(observed_likelihood, intervention_list_treatment)
       private$.cf_likelihood_control <- CF_Likelihood$new(observed_likelihood, intervention_list_control)
     },
@@ -75,7 +76,7 @@ Param_spOR <- R6Class(
       intervention_nodes <- union(names(self$intervention_list_treatment), names(self$intervention_list_control))
 
       W <- tmle_task$get_tmle_node("W")
-      V <- model.matrix(self$formula_OR, as.data.frame(W))
+      V <- model.matrix(self$formula_logOR, as.data.frame(W))
       A <- tmle_task$get_tmle_node("A", format = TRUE)[[1]]
       Y <- tmle_task$get_tmle_node("Y", format = TRUE)[[1]]
       g <- self$observed_likelihood$get_likelihoods(tmle_task, "A", fold_number)
@@ -91,10 +92,10 @@ Param_spOR <- R6Class(
       Qorig <- Q
       Q0 <- bound(Q0, 0.005)
       Q1 <- bound(Q1, 0.005)
-      OR <- Q1*(1-Q1) / (Q0*(1-Q0))
+      sigma_rel <- Q1*(1-Q1) / (Q0*(1-Q0))
 
 
-      h_star <-  -1*as.vector((g1*OR) / (g1*OR + (1-g1)))
+      h_star <-  -1*as.vector((g1*sigma_rel) / (g1*sigma_rel + (1-g1)))
       H <- as.matrix(V*(A  + h_star))
 
       # Store EIF component
@@ -137,8 +138,8 @@ Param_spOR <- R6Class(
       # Q <- Q_packed[[3]]
       Q0 <- bound(Q0, 0.0005)
       Q1 <- bound(Q1, 0.0005)
-      beta <- get_beta(W, A, self$formula_OR, Q1, Q0, family = binomial(), weights = weights)
-      V <- model.matrix(self$formula_OR, as.data.frame(W))
+      beta <- get_beta(W, A, self$formula_logOR, Q1, Q0, family = binomial(), weights = weights)
+      V <- model.matrix(self$formula_logOR, as.data.frame(W))
       OR <- exp(V%*%beta)
 
       IC <- EIF
@@ -167,16 +168,16 @@ Param_spOR <- R6Class(
     update_nodes = function() {
       return(c(self$outcome_node))
     },
-    formula_OR = function(){
-      return(private$.formula_OR)
+    formula_logOR = function(){
+      return(private$.formula_logOR)
     }
   ),
   private = list(
     .type = "OR",
     .cf_likelihood_treatment = NULL,
     .cf_likelihood_control = NULL,
     .supports_outcome_censoring = TRUE,
-    .formula_OR = NULL,
+    .formula_logOR = NULL,
     .submodel = list(Y = "gaussian_identity")
   )
 )

R/Param_spRR.R

@@ -1,6 +1,8 @@
-#' Average Treatment Effect
-#'
-#' Parameter definition for the Average Treatment Effect (ATE).
+#' Semiparametric estimation of the conditonal relative risk/treatment-effect for arbitrary partially-linear log-linear/link regression models.
+#' Arbitrary user-specified parametric models for the conditional relative-risk are supported.
+#` This method implements semiparametric efficient relative-risk regression for nonnegative outcomes.
+#' Assuming the semiparametric model to be true allows for some efficiency gain (when true) but may lead to less robust estimates due to misspecification.
+#' The parametric model is at the log-scale and therefore the coefficients returned code the linear predictor for the `log`-relative-risk.
 #' @importFrom R6 R6Class
 #' @importFrom uuid UUIDgenerate
 #' @importFrom methods is
@@ -17,7 +19,7 @@
 #'   \describe{
 #'     \item{\code{observed_likelihood}}{A \code{\link{Likelihood}} corresponding to the observed likelihood
 #'     }
-#'     \item{\code{formula_RR}}{...
+#'     \item{\code{formula_logRR}}{...
 #'     }
 #'     \item{\code{intervention_list_treatment}}{A list of objects inheriting from \code{\link{LF_base}}, representing the treatment intervention.
 #'     }
@@ -48,7 +50,7 @@ Param_spRR <- R6Class(
   class = TRUE,
   inherit = Param_base,
   public = list(
-    initialize = function(observed_likelihood,  formula_RR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
+    initialize = function(observed_likelihood,  formula_logRR =~ 1, intervention_list_treatment, intervention_list_control, outcome_node = "Y") {
       super$initialize(observed_likelihood, list(), outcome_node)
       if (!is.null(observed_likelihood$censoring_nodes[[outcome_node]])) {
         # add delta_Y=0 to intervention lists
@@ -57,7 +59,7 @@ Param_spRR <- R6Class(
         intervention_list_treatment <- c(intervention_list_treatment, censoring_intervention)
         intervention_list_control <- c(intervention_list_control, censoring_intervention)
       }
-      private$.formula_RR <- formula_RR
+      private$.formula_logRR <- formula_logRR
       private$.cf_likelihood_treatment <- CF_Likelihood$new(observed_likelihood, intervention_list_treatment)
       private$.cf_likelihood_control <- CF_Likelihood$new(observed_likelihood, intervention_list_control)
     },
@@ -75,7 +77,7 @@ Param_spRR <- R6Class(
       intervention_nodes <- union(names(self$intervention_list_treatment), names(self$intervention_list_control))
 
       W <- tmle_task$get_tmle_node("W")
-      V <- model.matrix(self$formula_RR, as.data.frame(W))
+      V <- model.matrix(self$formula_logRR, as.data.frame(W))
       A <- tmle_task$get_tmle_node("A", format = TRUE)[[1]]
       Y <- tmle_task$get_tmle_node("Y", format = TRUE)[[1]]
 
@@ -144,8 +146,8 @@ Param_spRR <- R6Class(
 
       Q0 <- pmax(Q0, 0.0005)
       Q1 <- pmax(Q1, 0.0005)
-      beta <- get_beta(W, A, self$formula_RR, Q1, Q0, family = poisson(), weights = weights)
-      V <- model.matrix(self$formula_RR, as.data.frame(W))
+      beta <- get_beta(W, A, self$formula_logRR, Q1, Q0, family = poisson(), weights = weights)
+      V <- model.matrix(self$formula_logRR, as.data.frame(W))
       RR <- exp(V%*%beta)
 
       IC <- as.matrix(EIF)
@@ -174,16 +176,16 @@ Param_spRR <- R6Class(
     update_nodes = function() {
       return(c(self$outcome_node))
     },
-    formula_RR = function(){
-      return(private$.formula_RR)
+    formula_logRR = function(){
+      return(private$.formula_logRR)
     }
   ),
   private = list(
     .type = "RR",
     .cf_likelihood_treatment = NULL,
     .cf_likelihood_control = NULL,
     .supports_outcome_censoring = TRUE,
-    .formula_RR = NULL,
+    .formula_logRR = NULL,
     .submodel = list(Y = "poisson_log")
   )
 )

R/tmle3_Spec_npCausalGLM.R

@@ -41,13 +41,13 @@ tmle3_Spec_npCausalGLM <- R6Class(
 
       return(likelihood)
     },
-    make_updater = function(convergence_type = "sample_size", verbose = TRUE,...) {
+    make_updater = function(convergence_type = "sample_size", verbose = F,...) {
       if(self$options$estimand == "CATE" || self$options$estimand == "CATT"){
         updater <- tmle3_Update$new(maxit=100,one_dimensional = FALSE,   verbose = verbose, constrain_step = FALSE, bounds = c(-Inf, Inf), ...)
       } else if (self$options$estimand == "OR"){
-        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE, convergence_type = convergence_type, verbose = verbose,delta_epsilon = 0.001, constrain_step = TRUE, bounds = 0.0025, ...)
+        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE, convergence_type = convergence_type, verbose = verbose,delta_epsilon = 0.01, constrain_step = TRUE, bounds = 0.0025, ...)
       } else if (self$options$estimand == "RR"){
-        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE,  convergence_type = convergence_type, verbose = verbose, delta_epsilon = 0.001, constrain_step = TRUE, bounds = c(0.0025, Inf), ...)
+        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE,  convergence_type = convergence_type, verbose = verbose, delta_epsilon = 0.01, constrain_step = TRUE, bounds = c(0.0025, Inf), ...)
       }
       return(updater)
     },

R/tmle3_spec_spCausalGLM.R

@@ -43,13 +43,13 @@ tmle3_Spec_spCausalGLM <- R6Class(
 
       return(likelihood)
     },
-    make_updater = function(convergence_type = "sample_size", verbose = TRUE,...) {
+    make_updater = function(convergence_type = "sample_size", verbose = F,...) {
       if(self$options$estimand == "CATE"){
         updater <- tmle3_Update$new(maxit=100,one_dimensional = FALSE,   verbose = verbose, constrain_step = FALSE, bounds = c(-Inf, Inf), ...)
       } else if (self$options$estimand == "OR"){
-        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE, convergence_type = convergence_type, verbose = verbose,delta_epsilon = 0.001, constrain_step = TRUE, bounds = 0.0025, ...)
+        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE, convergence_type = convergence_type, verbose = verbose,delta_epsilon = 0.01, constrain_step = TRUE, bounds = 0.0025, ...)
       } else if (self$options$estimand == "RR"){
-        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE,  convergence_type = convergence_type, verbose = verbose, delta_epsilon = 0.001, constrain_step = TRUE, bounds = c(0.0025, Inf), ...)
+        updater <- tmle3_Update$new(maxit = 200, one_dimensional = TRUE,  convergence_type = convergence_type, verbose = verbose, delta_epsilon = 0.01, constrain_step = TRUE, bounds = c(0.0025, Inf), ...)
       }
       return(updater)
     },

vignettes/testing.Rmd

@@ -11,8 +11,14 @@ knitr::opts_chunk$set(echo = TRUE)
 
 
 ```{r}
-library(sl3)
-n <- 200
+passes <- c()
+passes1 <- c()
+passes2 <- c()
+
+for(i in 1:100){
+  print(i)
+
+n <- 500
 W <- runif(n, -1, 1)
 A <- rbinom(n, size = 1, prob = plogis(W))
 Y <- rnorm(n, mean =  A+W, sd = 0.3)
@@ -25,37 +31,62 @@ learner_list <- list(A  = lrnr_A, Y = lrnr_Y0W, var_Y = Lrnr_mean$new())
 # spec_spCATE <- tmle3_Spec_spCausalGLM$new(~1, "CATE")
 # out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) 
 spec_spCATE <- tmle3_Spec_npCausalGLM$new(~1, "CATE")
-out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) 
-out
+suppressWarnings(out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) )
+ out <- out$summary
+passes <- c(passes , out$lower <= 1 & out$upper >= 1)
+
 
 spec_spCATE <- tmle3_Spec_npCausalGLM$new(~1, "CATT")
-out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) 
-out
+suppressWarnings(out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) )
+ out <- out$summary
+passes1 <- c(passes1 , out$lower <= 1 & out$upper >= 1)
+
 
 spec_spCATE <- tmle3_Spec_spCausalGLM$new(~1, "CATE")
-out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) 
-out
+suppressWarnings(out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) )
+ out <- out$summary
+passes2 <- c(passes2 , out$lower <= 1 & out$upper >= 1)
+
+print(mean(passes))
+print(mean(passes1))
+print(mean(passes2))
+}
 ```
 
 
 
 
-```{r}
+```{r, include = F}
+
+passes <- c()
+passes1 <- c()
+for(i in 1:100){
+  print(i)
 library(sl3)
-n <- 200
+n <- 500
 W <- runif(n, -1, 1)
-A <- rbinom(n, size = 1, prob = plogis(W))
+A <- rbinom(n, size = 1, prob = plogis(0))
 Y <-  rbinom(n, size = 1, prob = plogis(A + W))
+ quantile(plogis(1 + W) * (1-plogis(1 + W)) / ( plogis( W) * (1-plogis( W))))
 data <- data.table(W,A,Y)
 lrnr_Y0W <- Lrnr_glmnet$new()
 lrnr_A <- Lrnr_glm$new()
 node_list <- list (W = "W", A = "A", Y= "Y")
 learner_list <- list(A  = lrnr_A, Y = lrnr_Y0W)
 spec_spCATE <- tmle3_Spec_spCausalGLM$new(~1, "OR")
-out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) 
-out
+ suppressWarnings(out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list))
+ out <- out$summary
+passes <- c(passes , out$lower <= 1 & out$upper >= 1)
+
 
+spec_spCATE <- tmle3_Spec_npCausalGLM$new(~1, "OR")
+suppressWarnings(out <- tmle3(spec_spCATE, data, node_list, learner_list = learner_list) )
+ out <- out$summary
+passes1 <- c(passes1 , out$lower <= 1 & out$upper >= 1)
 
+print(mean(passes))
+print(mean(passes1))
+}
 ```