Merge pull request #1 from nt-williams/devel

log and exp operators

Author: nt-williams

DESCRIPTION

@@ -1,23 +1,24 @@
 Package: ife
 Type: Package
-Title: Influence Function Based Estimate Objects
-Version: 0.1.12
+Title: Autodiff for Influence Function Based Estimates
+Version: 0.1.2
 Authors@R: 
     c(person(given = "Nicholas",
              family = "Williams",
              role = c("aut", "cre", "cph"),
              email = "ntwilliams.personal@gmail.com", 
              comment = c(ORCID = "0000-0002-1378-4831")))
 Maintainer: Nicholas Williams <ntwilliams.personal@gmail.com>
-Description: Implements an S7 class for estimates based on influence functions. 
-  Standard arithmetic operations are defined for the class.
+Description: Implements an S7 class for estimates based on influence functions, 
+  with forward mode automatic differentiation defined for standard arithmetic 
+  operations.
 License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true
 Imports: 
     cli,
     generics,
-    S7
+    S7 (>= 0.2.0)
 RoxygenNote: 7.3.2
 Suggests: 
     testthat (>= 3.0.0)

NAMESPACE

@@ -3,7 +3,6 @@
 export(ife)
 export(influence_func_estimate)
 export(tidy)
-if (getRversion() < "4.3.0") importFrom(S7,`@`)
 importFrom(S7,S7_object)
 importFrom(S7,`method<-`)
 importFrom(S7,class_character)
@@ -13,8 +12,10 @@ importFrom(S7,new_class)
 importFrom(S7,new_generic)
 importFrom(S7,new_object)
 importFrom(S7,new_property)
+importFrom(cli,cli)
 importFrom(cli,cli_div)
 importFrom(cli,cli_end)
 importFrom(cli,cli_text)
+importFrom(cli,format_inline)
 importFrom(generics,tidy)
 importFrom(stats,qnorm)

NEWS.md

@@ -1,3 +1,8 @@
+# ife 0.1.2
+
+* Adding `log` and `exp` methods
+* Printing now returns up to 3 digits
+
 # ife 0.1.12
 
 * Conditional `importFrom` statement to fix bug when R version is less than 4.3.0.

R/influence_func_estimand.R

@@ -1,6 +1,5 @@
-#' @importFrom cli cli_div cli_text cli_end
-#' @importFrom S7 `@` new_class new_generic new_property new_object S7_object class_double class_character `method<-` class_numeric
-#' @rawNamespace if (getRversion() < "4.3.0") importFrom(S7,`@`)
+#' @importFrom cli cli cli_div format_inline cli_text cli_end
+#' @importFrom S7 new_class new_generic new_property new_object S7_object class_double class_character `method<-` class_numeric
 #' @importFrom stats qnorm
 NULL
 
@@ -97,10 +96,14 @@ influence_func_estimate <- new_class("influence_func_estimate",
 
 # print
 method(print, influence_func_estimate) <- function(x, ...) {
-  div <- cli_div(theme = list(.val = list(digits = 2)))
-  cli_text(cat("      "), "Estimate: {.val {x@x}}")
-  cli_text(cat("    "), "Std. error: {.val {x@std_error}}")
-  cli_text("95% Conf. int.: {.val {x@conf_int[1]}}, {.val {x@conf_int[2]}}")
+  div <- cli_div(theme = list(.val = list(digits = 3)))
+  cli({
+    cat(format_inline("      Estimate: {.val {x@x}}\n"))
+    # cli_text(cat("      "), "Estimate: {.val {x@x}}")
+    cat(format_inline("    Std. error: {.val {x@std_error}}\n"))
+    # cli_text(cat("    "), "Std. error: {.val {x@std_error}}")
+    cli_text("95% Conf. int.: {.val {x@conf_int[1]}}, {.val {x@conf_int[2]}}")
+  })
   cli_end(div)
 }
 
@@ -156,3 +159,13 @@ method(`*`, list(class_numeric, influence_func_estimate)) <- function(e1, e2) {
 }
 
 method(`*`, list(influence_func_estimate, class_numeric)) <- function(e1, e2) e2 * e1
+
+# log(x)
+method(log, influence_func_estimate) <- function(x, base) {
+  influence_func_estimate(log(x@x), x@eif / x@x, x@weights, x@id)
+}
+
+# exp(x)
+method(exp, influence_func_estimate) <- function(x) {
+  influence_func_estimate(exp(x@x), exp(x@x) * x@eif, x@weights, x@id)
+}

README.Rmd

@@ -5,6 +5,8 @@ output: github_document
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
 ```{r, include = FALSE}
+options(digits = 3)
+
 knitr::opts_chunk$set(
   collapse = TRUE,
   comment = "#>",
@@ -16,45 +18,72 @@ knitr::opts_chunk$set(
 # ife
 
 <!-- badges: start -->
+[![CRAN status](https://www.r-pkg.org/badges/version/ife)](https://CRAN.R-project.org/package=ife)
 <!-- badges: end -->
 
-S7 class (with Ops) for influence function based estimands
+S7 class for influence function estimands with forward mode automatic differentiation for variance estimation.
 
 ## Installation
 
-You can install the development version of ife from [GitHub](https://github.com/) with:
+You can install the development version of *ife* from [GitHub](https://github.com/) with:
 
 ``` r
 # install.packages("pak")
 pak::pak("nt-williams/ife")
 ```
 
+## Example
+
+Consider estimating the population mean outcome under treatment $E[E[Y \mid A=1,W]]$ and control $E[E[Y \mid A=0,W]]$ using augmented inverse probability weighting (AIPW).
+
 ```{r example}
 library(ife)
 
+# Generate simulated data
 n <- 500
-w <- runif(n)
-a <- rbinom(n, 1, 0.5)
-y <- rbinom(n, 1, plogis(-0.75 + a + w))
+w <- runif(n)                                   # confounder
+a <- rbinom(n, 1, 0.5)                          # treatment (randomized)
+y <- rbinom(n, 1, plogis(-0.75 + a + w))        # outcome
 
+# Create data-frames for counterfactual predictions
 foo <- data.frame(w, a, y)
 foo1 <- foo0 <- foo
-foo1$a <- 1
-foo0$a <- 0
+foo1$a <- 1  # everyone treated
+foo0$a <- 0  # everyone untreated
 
-pi <- 0.5
+# Fit outcome model and generate predictions
+pi <- 0.5  # known propensity score
 m <- glm(y ~ a + w, data = foo, family = binomial())
+Qa <- predict(m, type = "response")                    # predicted outcomes
+Q1 <- predict(m, newdata = foo1, type = "response")    # under treatment
+Q0 <- predict(m, newdata = foo0, type = "response")    # under control
 
-Qa <- predict(m, type = "response")
-Q1 <- predict(m, newdata = foo1, type = "response")
-Q0 <- predict(m, newdata = foo0, type = "response")
+# Calculate un-centered influence functions
+if1 <- a / pi * (y - Qa) + Q1              
+if0 <- (1 - a) / (1 - pi) * (y - Qa) + Q0  
+```
 
-if1 <- a / pi * (y - Qa) + Q1
-if0 <- (1 - a) / pi * (y - Qa) + Q0
+Create *ife* objects for these estimates using `influence_func_estimate()` or `ife()`:
 
+```{r}
 ife1 <- influence_func_estimate(mean(if1), if1)
 ife0 <- ife(mean(if0), if0)
+```
+
+*ife* then allows you to estimate contrasts between estimates, with variance estimated using automatic differentiation. The additive effect (risk difference) can be calculated as:
 
+```{r riskdiff}
 ife1 - ife0
+```
+
+The multiplicative effect (risk ratio) can be estimated as:
+
+```{r}
 ife1 / ife0
 ```
+
+For the risk ratio, which is strictly positive, you can estimate the effect on the log scale and exponentiate the confidence intervals to ensure the lower bound is always positive:
+
+```{r}
+exp(log(ife1 / ife0)@conf_int)
+```

README.md

@@ -4,56 +4,90 @@
 # ife
 
 <!-- badges: start -->
+
+[![CRAN
+status](https://www.r-pkg.org/badges/version/ife)](https://CRAN.R-project.org/package=ife)
 <!-- badges: end -->
 
-S7 class (with Ops) for influence function based estimands
+S7 class for influence function estimands with forward mode automatic
+differentiation for variance estimation.
 
 ## Installation
 
-You can install the development version of ife from
+You can install the development version of *ife* from
 [GitHub](https://github.com/) with:
 
 ``` r
 # install.packages("pak")
 pak::pak("nt-williams/ife")
 ```
 
+## Example
+
+Consider estimating the population mean outcome under treatment
+$E[E[Y \mid A=1,W]]$ and control $E[E[Y \mid A=0,W]]$ using augmented
+inverse probability weighting (AIPW).
+
 ``` r
 library(ife)
 
+# Generate simulated data
 n <- 500
-w <- runif(n)
-a <- rbinom(n, 1, 0.5)
-y <- rbinom(n, 1, plogis(-0.75 + a + w))
+w <- runif(n)                                   # confounder
+a <- rbinom(n, 1, 0.5)                          # treatment (randomized)
+y <- rbinom(n, 1, plogis(-0.75 + a + w))        # outcome
 
+# Create data-frames for counterfactual predictions
 foo <- data.frame(w, a, y)
 foo1 <- foo0 <- foo
-foo1$a <- 1
-foo0$a <- 0
+foo1$a <- 1  # everyone treated
+foo0$a <- 0  # everyone untreated
 
-pi <- 0.5
+# Fit outcome model and generate predictions
+pi <- 0.5  # known propensity score
 m <- glm(y ~ a + w, data = foo, family = binomial())
+Qa <- predict(m, type = "response")                    # predicted outcomes
+Q1 <- predict(m, newdata = foo1, type = "response")    # under treatment
+Q0 <- predict(m, newdata = foo0, type = "response")    # under control
 
-Qa <- predict(m, type = "response")
-Q1 <- predict(m, newdata = foo1, type = "response")
-Q0 <- predict(m, newdata = foo0, type = "response")
+# Calculate un-centered influence functions
+if1 <- a / pi * (y - Qa) + Q1              
+if0 <- (1 - a) / (1 - pi) * (y - Qa) + Q0  
+```
 
-if1 <- a / pi * (y - Qa) + Q1
-if0 <- (1 - a) / pi * (y - Qa) + Q0
+Create *ife* objects for these estimates using
+`influence_func_estimate()` or `ife()`:
 
+``` r
 ife1 <- influence_func_estimate(mean(if1), if1)
 ife0 <- ife(mean(if0), if0)
+```
 
+*ife* then allows you to estimate contrasts between estimates, with
+variance estimated using automatic differentiation. The additive effect
+(risk difference) can be calculated as:
+
+``` r
 ife1 - ife0
+#>       Estimate: 0.254
+#>     Std. error: 0.042
+#> 95% Conf. int.: 0.172, 0.336
 ```
 
-    #> Estimate: 0.15
+The multiplicative effect (risk ratio) can be estimated as:
 
-    #> Std. error: 0.04
-    #> 95% Conf. int.: 0.07, 0.24
-    ife1 / ife0
+``` r
+ife1 / ife0
+#>       Estimate: 1.583
+#>     Std. error: 0.129
+#> 95% Conf. int.: 1.33, 1.837
+```
 
-    #> Estimate: 1.32
+For the risk ratio, which is strictly positive, you can estimate the
+effect on the log scale and exponentiate the confidence intervals to
+ensure the lower bound is always positive:
 
-    #> Std. error: 0.11
-    #> 95% Conf. int.: 1.11, 1.52
+``` r
+exp(log(ife1 / ife0)@conf_int)
+#> [1] 1.35 1.86
+```

tests/testthat/test-Ops.R

@@ -84,4 +84,12 @@ describe("operators", {
       expect_equal(z@std_error, 0.6581774, tolerance = tol)
     })
   })
+
+  describe("log", {
+    it("log(ife / ife)", {
+      z <- log(x / y)
+      log_eif <- (x@eif / x@x) - (y@eif / y@x)
+      expect_equal(z@eif, log_eif)
+    })
+  })
 })