-
Notifications
You must be signed in to change notification settings - Fork 6
/
README.Rmd
148 lines (112 loc) · 5.58 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
---
output:
md_document:
variant: markdown_github
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.align = "center",
fig.path = "man/figures/README-"
)
```
<!-- badges: start -->
[](https://CRAN.R-project.org/package=investr)
[](https://cranlogs.r-pkg.org/badges/grand-total/investr)
<!-- badges: end -->
# investr: Inverse Estimation in R
Inverse estimation, also referred to as the calibration problem, is a classical and well-known problem in regression. In simple terms, it involves the use of an observed value of the response (or specified value of the mean response) to make inference on the corresponding unknown value of an explanatory variable.
A detailed introduction to investr has been published in The R Journal: ["investr: An R Package for Inverse Estimation"](https://journal.r-project.org/archive/2014/RJ-2014-009/index.html). You can track development at https://github.com/bgreenwell/investr. To report bugs or issues, contact the main author directly or submit them to https://github.com/bgreenwell/investr/issues.
As of right now, `investr` supports (univariate) inverse estimation with objects of class:
* `"lm"` - linear models (multiple predictor variables allowed)
* `"glm"` - generalized linear models (multiple predictor variables allowed)
* `"nls"` - nonlinear least-squares models
* `"lme"` - linear mixed-effects models (fit using the `nlme` package)
## Installation
The package is [currently listed on CRAN](https://cran.r-project.org/package=investr) and can easily be installed:
```{r, eval=FALSE}
# Install from CRAN
install.packages("investr", dep = TRUE)
# Alternatively, install the development version from GitHub
devtools::install_github("bgreenwell/investr")
```
The package is also part of the [ChemPhys task view](https://cran.r-project.org/view=ChemPhys), a collection of R packages useful for analyzing data from chemistry and physics experiments. These packages can all be installed at once (including `investr`) using the `ctv` package (Zeileis, 2005):
```{r, eval=FALSE}
# Install the ChemPhys task view
install.packages("ctv")
ctv::install.views("ChemPhys")
```
## Examples
### Dobson's Beetle Data
In binomial regression, the estimated lethal dose corresponding to a specific probability _p_ of death is often referred to as _LDp_. `invest` obtains an estimate of _LDp_ by inverting the fitted mean response on the link scale. Similarly, a confidence interval for _LDp_ can be obtained by inverting a confidence interval for the mean response on the link scale.
```{r glm-beetle-invest}
# Load required packages
library(investr)
# Binomial regression
beetle.glm <- glm(cbind(y, n-y) ~ ldose, data = beetle,
family = binomial(link = "cloglog"))
plotFit(beetle.glm, lwd.fit = 2, cex = 1.2, pch = 21, bg = "lightskyblue",
lwd = 2, xlab = "Log dose", ylab = "Probability")
# Median lethal dose
invest(beetle.glm, y0 = 0.5)
# 90% lethal dose
invest(beetle.glm, y0 = 0.9)
# 99% lethal dose
invest(beetle.glm, y0 = 0.99)
```
To obtain an estimate of the standard error, we can use the Wald method:
```{r glm-beetle-invest-se}
invest(beetle.glm, y0 = 0.5, interval = "Wald")
# The MASS package function dose.p can be used too
MASS::dose.p(beetle.glm, p = 0.5)
```
### Including a factor variable
Multiple predictor variables are allowed for objects of class `lm` and `glm`.
For instance, the example from `?MASS::dose.p` can be re-created as follows:
```{r glm-beetle-investr-vs-MASS}
# Load required packages
library(MASS)
# Data
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
sex <- factor(rep(c("M", "F"), c(6, 6)))
SF <- cbind(numdead, numalive = 20 - numdead)
budworm <- data.frame(ldose, numdead, sex, SF)
# Logistic regression
budworm.glm <- glm(SF ~ sex + ldose - 1, family = binomial, data = budworm)
# Using dose.p function from package MASS
dose.p(budworm.glm, cf = c(1, 3), p = 1/4)
# Using invest function from package investr
invest(budworm.glm, y0 = 1/4,
interval = "Wald",
x0.name = "ldose",
newdata = data.frame(sex = "F"))
```
### Bioassay on Nasturtium
The data here contain the actual concentrations of an agrochemical present in soil samples versus the weight of the plant after three weeks of growth. These data are stored in the data frame `nasturtium` and are loaded with the package. A simple
log-logistic model describes the data well:
```{r nls-nasturtium-fit}
# Log-logistic model for the nasturtium data
nas.nls <- nls(weight ~ theta1/(1 + exp(theta2 + theta3 * log(conc))),
start = list(theta1 = 1000, theta2 = -1, theta3 = 1),
data = nasturtium)
# Plot the fitted model
plotFit(nas.nls, lwd.fit = 2)
```
Three new replicates of the response (309, 296, 419) at an unknown concentration of interest ($x_0$) are measured. It is desired to estimate $x_0$.
```{r nls-nasturtium-invest}
# Inversion method
invest(nas.nls, y0 = c(309, 296, 419), interval = "inversion")
# Wald method
invest(nas.nls, y0 = c(309, 296, 419), interval = "Wald")
```
The intervals both rely on large sample results and normality. In practice, the bootstrap may be more reliable:
```{r nls-nasturtium-bootstrap}
# Bootstrap calibration intervals (may take a few seconds)
boo <- invest(nas.nls, y0 = c(309, 296, 419), interval = "percentile",
nsim = 9999, seed = 101)
boo # print bootstrap summary
plot(boo) # plot results
```