Nothing
tests/testthat/test-msmm.R
In OneSampleMR: One Sample Mendelian Randomization and Instrumental Variable Analyses
# tests for MSMM
# Test equivalence to Stata output ----
test_that("Stata output check", {
skip_on_cran()
skip_if_not_installed("haven")
library(haven)
dat <- read_dta("https://www.stata-press.com/data/r17/trip.dta")
fit1 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
cbd + ptn + worker + weekend + pt,
data = dat, estmethod = "gmmalt")
expect_equal(log(fit1$crrci["tcost",1]), log(1.036), tolerance = 0.05)
# fit2 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat)
# fit3 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat, estmethod = "gmmalt",
# t0 = c(.265, -.008, -.011, .662, .301, .035))
# fit4 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat,
# t0 = c(exp(.265), -.008, -.011, .662, .301, .035))
})
# . * Setup
# . webuse trip, clear
# (Household trips)
#
# .
# . * Generalized method of moments: multiplicative errors
# . ivpoisson gmm trips cbd ptn worker weekend (tcost=pt), multiplicative nolog
#
# Final GMM criterion Q(b) = 1.27e-31
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 6 Number of obs = 5,000
# Number of moments = 6
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# _cons | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * Display incidence-rate ratios
# . ivpoisson, irr
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 6 Number of obs = 5,000
# Number of moments = 6
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | IRR Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# cbd | .9916371 .0020003 -4.16 0.000 .9877243 .9955655
# ptn | .9887491 .0021573 -5.19 0.000 .9845299 .9929864
# worker | 1.939251 .1008234 12.74 0.000 1.751376 2.14728
# weekend | 1.351118 .0490026 8.30 0.000 1.258409 1.450657
# _cons | 1.304008 .2021377 1.71 0.087 .9623501 1.766962
# ------------------------------------------------------------------------------*
# Note: _cons estimates baseline incidence rate.
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * Control-function method
# . ivpoisson cfunction trips cbd ptn worker weekend (tcost=pt), nolog
#
# Final GMM criterion Q(b) = 9.78e-27
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 13 Number of obs = 5,000
# Number of moments = 13
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | -.0082567 .0020005 -4.13 0.000 -.0121777 -.0043357
# ptn | -.0113719 .0021625 -5.26 0.000 -.0156102 -.0071335
# worker | .6903044 .0521642 13.23 0.000 .5880645 .7925444
# weekend | .2978149 .0356474 8.35 0.000 .2279472 .3676825
# tcost | .0320718 .0092738 3.46 0.001 .0138955 .0502481
# _cons | .2145986 .1359327 1.58 0.114 -.0518246 .4810218
# -------------+----------------------------------------------------------------*
# tcost |
# cbd | .0165466 .0043693 3.79 0.000 .0079829 .0251102
# ptn | -.040652 .0045946 -8.85 0.000 -.0496573 -.0316467
# worker | 1.550985 .0996496 15.56 0.000 1.355675 1.746294
# weekend | .0423009 .0779101 0.54 0.587 -.1104002 .1950019
# pt | .7739176 .0150072 51.57 0.000 .7445041 .8033312
# _cons | 12.13934 .1123471 108.05 0.000 11.91915 12.35954
# -------------+----------------------------------------------------------------*
# /c_tcost | .1599984 .0111752 14.32 0.000 .1380954 .1819014
# ------------------------------------------------------------------------------*
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# . ivpoisson, irr
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 13 Number of obs = 5,000
# Number of moments = 13
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | IRR Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | .9917773 .0019841 -4.13 0.000 .9878961 .9956737
# ptn | .9886926 .002138 -5.26 0.000 .984511 .9928919
# worker | 1.994323 .1040322 13.23 0.000 1.8005 2.20901
# weekend | 1.346912 .0480139 8.35 0.000 1.256019 1.444383
# tcost | 1.032592 .009576 3.46 0.001 1.013992 1.051532
# _cons | 1.239364 .1684701 1.58 0.114 .9494954 1.617727
# -------------+----------------------------------------------------------------*
# tcost |
# cbd | .0165466 .0043693 3.79 0.000 .0079829 .0251102
# ptn | -.040652 .0045946 -8.85 0.000 -.0496573 -.0316467
# worker | 1.550985 .0996496 15.56 0.000 1.355675 1.746294
# weekend | .0423009 .0779101 0.54 0.587 -.1104002 .1950019
# pt | .7739176 .0150072 51.57 0.000 .7445041 .8033312
# _cons | 12.13934 .1123471 108.05 0.000 11.91915 12.35954
# -------------+----------------------------------------------------------------*
# /c_tcost | .1599984 .0111752 14.32 0.000 .1380954 .1819014
# ------------------------------------------------------------------------------*
# Note: Estimates are transformed only in the first equation.
# Note: _cons estimates baseline incidence rate.
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * ivpois for comparison
# . ivpois trips cbd ptn worker weekend, endog(tcost) exog(pt)
# ------------------------------------------------------------------------------*
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# _cons | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
#
# . * coding msmm using the gmm command
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*`x'*{psi} - {xb:`covars'}) - {ey0}), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 9.64e-34
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# psi |
# _cons | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# -------------+----------------------------------------------------------------*
# xb |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /ey0 | 1.304008 .2021377 6.45 0.000 .907825 1.70019
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . lincom [psi]_cons, eform // causal risk ratio
#
# ( 1) [psi]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# ------------------------------------------------------------------------------*
#
# . // estat overid
# .
# . * coding msmm using the gmm command with first derivatives
# . local y trips
#
# . local xlist tcost cbd ptn worker weekend
#
# . local zlist pt cbd ptn worker weekend
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*{xb:`xlist'}) - {ey0}), ///
# > instruments(`zlist') tracelevel("none") ///
# > deriv(/xb = -1*`y'*exp(-1*{xb:})) /// remembering that deriv automatically multiplies by x in xb
# > deriv(/ey0 = -1)
#
# Final GMM criterion Q(b) = 5.31e-33
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /ey0 | 1.304008 .2021377 6.45 0.000 .907825 1.70019
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# .
# . * coding msmm using the alternative moment condition using the gmm command
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*`x'*{psi} - {xb:`covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 1.09e-28
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# psi |
# _cons | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# -------------+----------------------------------------------------------------*
# xb |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /logey0 | .2654423 .1550127 1.71 0.087 -.038377 .5692616
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . lincom [psi]_cons, eform // causal risk ratio
#
# ( 1) [psi]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# ------------------------------------------------------------------------------*
#
# . lincom [logey0]_cons, eform // causal baseline risk
#
# ( 1) [logey0]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.304008 .2021377 1.71 0.087 .9623501 1.766962
# ------------------------------------------------------------------------------*
#
# . // estat overid
# .
# . // with first derivatives
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none") ///
# > deriv(/xb = -1*`y'*exp(-1*{xb:} - {logey0})) ///
# > deriv(/logey0 = -1*`y'*exp(-1*{xb:} - {logey0}))
#
# Final GMM criterion Q(b) = 1.09e-28
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /logey0 | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . * multiple phenotype examples
# .
# . webuse trip, clear
# (Household trips)
#
# .
# . ** without covariates
# .
# . ** ivpoisson fit
# . ivpoisson gmm trips (tcost cbd = ptn worker), multiplicative nolog
#
# Final GMM criterion Q(b) = 3.12e-29
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 3 Number of obs = 5,000
# Number of moments = 3
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4214074 .0700091 6.02 0.000 .2841921 .5586228
# cbd | .075126 .1304184 0.58 0.565 -.1804893 .3307413
# _cons | -5.359227 1.47909 -3.62 0.000 -8.258191 -2.460263
# ------------------------------------------------------------------------------*
# Instrumented: tcost cbd
# Instruments: ptn worker
#
# . lincom _cons, eform
#
# ( 1) [trips]_cons = 0
#
# ------------------------------------------------------------------------------*
# trips | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | .0047045 .0069584 -3.62 0.000 .0002591 .0854124
# ------------------------------------------------------------------------------*
#
# .
# . *** moment conditions 1
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . gmm (`y'*exp(-1*{xb:`x'}) - {ey0}), ///
# > instruments(`z') tracelevel("none")
#
# Final GMM criterion Q(b) = 5.27e-09
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 3
# Number of moments = 3
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4213734 .0699533 6.02 0.000 .2842675 .5584793
# cbd | .0751167 .1303905 0.58 0.565 -.1804441 .3306774
# -------------+----------------------------------------------------------------*
# /ey0 | .0047057 .0069581 0.68 0.499 -.0089319 .0183433
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker _cons
#
# .
# . *** moment conditions 2
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . gmm (`y'*exp(-1*{xb:`x'} - {logey0}) - 1), ///
# > instruments(`z') tracelevel("none")
#
# Final GMM criterion Q(b) = 1.21e-33
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 3
# Number of moments = 3
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4214074 .0700091 6.02 0.000 .2841921 .5586228
# cbd | .075126 .1304184 0.58 0.565 -.1804893 .3307413
# -------------+----------------------------------------------------------------*
# /logey0 | -5.359227 1.47909 -3.62 0.000 -8.258191 -2.460264
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker _cons
#
# .
# . ** including covariates
# . ivpoisson gmm trips weekend (tcost cbd = ptn worker), multiplicative nolog
#
# Final GMM criterion Q(b) = 2.33e-32
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 4 Number of obs = 5,000
# Number of moments = 4
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252243 .0750159 5.67 0.000 .2781958 .5722528
# cbd | .0773626 .1414159 0.55 0.584 -.1998075 .3545328
# weekend | .2157507 .0804855 2.68 0.007 .058002 .3734994
# _cons | -5.484431 1.57664 -3.48 0.001 -8.574588 -2.394274
# ------------------------------------------------------------------------------*
# Instrumented: tcost cbd
# Instruments: weekend ptn worker
#
# . lincom _cons, eform
#
# ( 1) [trips]_cons = 0
#
# ------------------------------------------------------------------------------*
# trips | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | .0041509 .0065445 -3.48 0.001 .0001888 .0912389
# ------------------------------------------------------------------------------*
#
# .
# . *** moment conditions 1
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . local covars weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'}) - {ey0}), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 2.73e-12
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 4
# Number of moments = 4
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252238 .0750177 5.67 0.000 .2781917 .5722559
# cbd | .0773625 .1414325 0.55 0.584 -.1998401 .354565
# weekend | .2157505 .0804888 2.68 0.007 .0579953 .3735057
# -------------+----------------------------------------------------------------*
# /ey0 | .0041509 .006545 0.63 0.526 -.0086771 .0169789
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker weekend _cons
#
# .
# . *** moment conditions 2
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . local covars weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 2.43e-32
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 4
# Number of moments = 4
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252243 .0750159 5.67 0.000 .2781958 .5722528
# cbd | .0773626 .1414159 0.55 0.584 -.1998075 .3545328
# weekend | .2157507 .0804855 2.68 0.007 .058002 .3734994
# -------------+----------------------------------------------------------------*
# /logey0 | -5.484431 1.57664 -3.48 0.001 -8.574588 -2.394274
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker weekend _cons
# Single instrument example ----
test_that("Single instrument example", {
skip_on_cran()
# skip_if_not_installed("ivtools")
# Data generation from the example in the ivtools ivglm() helpfile
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
# ivtools for comparison fit
# library(ivtools)
# fitZ.L <- glm(Z ~ 1, data = dat)
# fitY.LZX <- glm(Y ~ X + Z + X*Z, family = "binomial", data = dat)
# fitLogGest <- ivglm(estmethod = "g",
# X = "X",
# fitZ.L = fitZ.L,
# fitY.LZX = fitY.LZX,
# data = dat,
# link = "log",
# Y = "Y")
logcrr <- 0.1314654 # fitLogGest$est["X"]
logcrrse <- 0.06035374 # sqrt(fitLogGest$vcov)
fit01 <- msmm(Y ~ X | Z, data = dat)
expect_equal(log(fit01$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit02 <- msmm(Y ~ X | Z, data = dat, estmethod = "gmm")
expect_equal(log(fit02$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit03 <- msmm(Y ~ X | Z, data = dat, estmethod = "gmmalt")
expect_equal(log(fit03$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit04 <- msmm(Y ~ X | Z, data = dat, estmethod = "tsls")
expect_equal(log(fit04$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit05 <- msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt")
expect_equal(log(fit05$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
expect_s3_class(fit01, "msmm")
smy <- summary(fit01)
expect_s3_class(smy, "summary.msmm")
expect_output(print(fit01))
expect_output(print(smy))
})
# check subset argument ----
test_that("Check subset argument", {
skip_on_cran()
# Data generation from the example in the ivtools ivglm() helpfile
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
datfifty <- dat[1:50,]
fitcompare <- msmm(Y ~ X | Z, dat = datfifty)
fitcheck <- msmm(Y ~ X | Z, dat = dat, subset = 1:50)
expect_equal(fitcheck$crrci, fitcompare$crrci)
})
# check variables with different names ----
test_that("Check using different variable names", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$E <- dat$X
dat$R <- dat$Y
dat$W <- dat$Z
fit06 <- msmm(R ~ E | W, data = dat)
expect_equal(log(fit06$crrci[1]), 0.1314654, tolerance = 0.05, ignore_attr = "names")
})
# Multiple instrument example ----
test_that("Multiple instrument example", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX <- plogis(0.7*G1 + G2 - G3 + U)
X <- rbinom(n, 1, pX)
pY <- plogis(-2 + psi0*X + U)
Y <- rbinom(n, 1, pY)
dat <- data.frame(G1, G2, G3, X, Y)
fit11 <- msmm(Y ~ X | G1 + G2 + G3, data = dat)
expect_equal(log(fit11$crrci[1]), 0.45, tolerance = 0.05)
fit12 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "gmm")
expect_equal(log(fit12$crrci[1]), 0.45, tolerance = 0.05)
fit13 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "gmmalt")
expect_equal(log(fit13$crrci[1]), 0.45, tolerance = 0.1)
fit14 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "tsls")
expect_equal(log(fit14$crrci[1]), 0.45, tolerance = 0.02)
fit15 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "tslsalt")
expect_equal(log(fit15$crrci[1]), 0.45, tolerance = 0.15)
})
# non-binary x with tsls, tslsalt methods fail ----
test_that("Non-binary x with tsls and tslsalt methods should produce an error", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$X[1] <- 2
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# non-binary y fail for tsls methods ----
test_that("Y needs to be binary for the TSLS methods", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$Y[1] <- 2
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# All methods fail for a non-integer in Y ----
test_that("Methods fail for non-integer Y", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$Y[1] <- 1.5
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "gmm"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "gmmalt"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# Multiple exposure example ----
test_that("Multiple exposure example", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U)
Y <- rbinom(n, 1, pY)
dat <- data.frame(G1, G2, G3, X1, X2, Y)
fit21 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat)
expect_equal(log(fit21$crrci[1,1]), 0.08, tolerance = 0.02)
fit22 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "gmm")
expect_equal(log(fit22$crrci[1,1]), 0.08, tolerance = 0.02)
fit23 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "gmmalt")
expect_equal(log(fit23$crrci[1,1]), 0.08, tolerance = 0.1)
expect_error(msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "tslsalt"))
})
# Multiple exposure example with different variable names ----
test_that("Multiple exposure example with different variable names", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U)
Y <- rbinom(n, 1, pY)
E1 <- X1
E2 <- X2
R <- Y
dat <- data.frame(G1, G2, G3, E1, E2, R)
fit24 <- msmm(Y ~ E1 + E2 | G1 + G2 + G3, data = dat)
expect_equal(log(fit24$crrci[1,1]), 0.08, tolerance = 0.01)
})
# Example adjusting for a covariate ----
test_that("Adjusting for covariate", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
C <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U + C)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U + C)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U + C)
Y <- rbinom(n, 1, pY)
E1 <- X1
E2 <- X2
R <- Y
dat <- data.frame(G1, G2, G3, E1, E2, R, C)
fit25 <- msmm(Y ~ E1 + C | G1 + G2 + G3 + C, data = dat)
expect_equal(log(fit25$crrci[1,1]), .173, tolerance = .01)
fit26 <- msmm(Y ~ E1 + E2 + C | G1 + G2 + G3 + C, data = dat)
expect_equal(log(fit26$crrci[1,1]), .184, tolerance = .01)
})
Try the OneSampleMR package in your browser
Any scripts or data that you put into this service are public.
OneSampleMR documentation built on June 27, 2024, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# tests for MSMM
# Test equivalence to Stata output ----
test_that("Stata output check", {
skip_on_cran()
skip_if_not_installed("haven")
library(haven)
dat <- read_dta("https://www.stata-press.com/data/r17/trip.dta")
fit1 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
cbd + ptn + worker + weekend + pt,
data = dat, estmethod = "gmmalt")
expect_equal(log(fit1$crrci["tcost",1]), log(1.036), tolerance = 0.05)
# fit2 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat)
# fit3 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat, estmethod = "gmmalt",
# t0 = c(.265, -.008, -.011, .662, .301, .035))
# fit4 <- msmm(trips ~ cbd + ptn + worker + weekend + tcost |
# cbd + ptn + worker + weekend + pt,
# data = dat,
# t0 = c(exp(.265), -.008, -.011, .662, .301, .035))
})
# . * Setup
# . webuse trip, clear
# (Household trips)
#
# .
# . * Generalized method of moments: multiplicative errors
# . ivpoisson gmm trips cbd ptn worker weekend (tcost=pt), multiplicative nolog
#
# Final GMM criterion Q(b) = 1.27e-31
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 6 Number of obs = 5,000
# Number of moments = 6
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# _cons | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * Display incidence-rate ratios
# . ivpoisson, irr
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 6 Number of obs = 5,000
# Number of moments = 6
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | IRR Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# cbd | .9916371 .0020003 -4.16 0.000 .9877243 .9955655
# ptn | .9887491 .0021573 -5.19 0.000 .9845299 .9929864
# worker | 1.939251 .1008234 12.74 0.000 1.751376 2.14728
# weekend | 1.351118 .0490026 8.30 0.000 1.258409 1.450657
# _cons | 1.304008 .2021377 1.71 0.087 .9623501 1.766962
# ------------------------------------------------------------------------------*
# Note: _cons estimates baseline incidence rate.
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * Control-function method
# . ivpoisson cfunction trips cbd ptn worker weekend (tcost=pt), nolog
#
# Final GMM criterion Q(b) = 9.78e-27
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 13 Number of obs = 5,000
# Number of moments = 13
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | -.0082567 .0020005 -4.13 0.000 -.0121777 -.0043357
# ptn | -.0113719 .0021625 -5.26 0.000 -.0156102 -.0071335
# worker | .6903044 .0521642 13.23 0.000 .5880645 .7925444
# weekend | .2978149 .0356474 8.35 0.000 .2279472 .3676825
# tcost | .0320718 .0092738 3.46 0.001 .0138955 .0502481
# _cons | .2145986 .1359327 1.58 0.114 -.0518246 .4810218
# -------------+----------------------------------------------------------------*
# tcost |
# cbd | .0165466 .0043693 3.79 0.000 .0079829 .0251102
# ptn | -.040652 .0045946 -8.85 0.000 -.0496573 -.0316467
# worker | 1.550985 .0996496 15.56 0.000 1.355675 1.746294
# weekend | .0423009 .0779101 0.54 0.587 -.1104002 .1950019
# pt | .7739176 .0150072 51.57 0.000 .7445041 .8033312
# _cons | 12.13934 .1123471 108.05 0.000 11.91915 12.35954
# -------------+----------------------------------------------------------------*
# /c_tcost | .1599984 .0111752 14.32 0.000 .1380954 .1819014
# ------------------------------------------------------------------------------*
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# . ivpoisson, irr
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 13 Number of obs = 5,000
# Number of moments = 13
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | IRR Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | .9917773 .0019841 -4.13 0.000 .9878961 .9956737
# ptn | .9886926 .002138 -5.26 0.000 .984511 .9928919
# worker | 1.994323 .1040322 13.23 0.000 1.8005 2.20901
# weekend | 1.346912 .0480139 8.35 0.000 1.256019 1.444383
# tcost | 1.032592 .009576 3.46 0.001 1.013992 1.051532
# _cons | 1.239364 .1684701 1.58 0.114 .9494954 1.617727
# -------------+----------------------------------------------------------------*
# tcost |
# cbd | .0165466 .0043693 3.79 0.000 .0079829 .0251102
# ptn | -.040652 .0045946 -8.85 0.000 -.0496573 -.0316467
# worker | 1.550985 .0996496 15.56 0.000 1.355675 1.746294
# weekend | .0423009 .0779101 0.54 0.587 -.1104002 .1950019
# pt | .7739176 .0150072 51.57 0.000 .7445041 .8033312
# _cons | 12.13934 .1123471 108.05 0.000 11.91915 12.35954
# -------------+----------------------------------------------------------------*
# /c_tcost | .1599984 .0111752 14.32 0.000 .1380954 .1819014
# ------------------------------------------------------------------------------*
# Note: Estimates are transformed only in the first equation.
# Note: _cons estimates baseline incidence rate.
# Instrumented: tcost
# Instruments: cbd ptn worker weekend pt
#
# .
# . * ivpois for comparison
# . ivpois trips cbd ptn worker weekend, endog(tcost) exog(pt)
# ------------------------------------------------------------------------------*
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# trips |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# _cons | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
#
# . * coding msmm using the gmm command
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*`x'*{psi} - {xb:`covars'}) - {ey0}), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 9.64e-34
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# psi |
# _cons | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# -------------+----------------------------------------------------------------*
# xb |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /ey0 | 1.304008 .2021377 6.45 0.000 .907825 1.70019
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . lincom [psi]_cons, eform // causal risk ratio
#
# ( 1) [psi]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# ------------------------------------------------------------------------------*
#
# . // estat overid
# .
# . * coding msmm using the gmm command with first derivatives
# . local y trips
#
# . local xlist tcost cbd ptn worker weekend
#
# . local zlist pt cbd ptn worker weekend
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*{xb:`xlist'}) - {ey0}), ///
# > instruments(`zlist') tracelevel("none") ///
# > deriv(/xb = -1*`y'*exp(-1*{xb:})) /// remembering that deriv automatically multiplies by x in xb
# > deriv(/ey0 = -1)
#
# Final GMM criterion Q(b) = 5.31e-33
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /ey0 | 1.304008 .2021377 6.45 0.000 .907825 1.70019
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# .
# . * coding msmm using the alternative moment condition using the gmm command
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*`x'*{psi} - {xb:`covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 1.09e-28
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# psi |
# _cons | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# -------------+----------------------------------------------------------------*
# xb |
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /logey0 | .2654423 .1550127 1.71 0.087 -.038377 .5692616
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . lincom [psi]_cons, eform // causal risk ratio
#
# ( 1) [psi]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.035846 .0101701 3.59 0.000 1.016103 1.055972
# ------------------------------------------------------------------------------*
#
# . lincom [logey0]_cons, eform // causal baseline risk
#
# ( 1) [logey0]_cons = 0
#
# ------------------------------------------------------------------------------*
# | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | 1.304008 .2021377 1.71 0.087 .9623501 1.766962
# ------------------------------------------------------------------------------*
#
# . // estat overid
# .
# . // with first derivatives
# . local y trips
#
# . local x tcost
#
# . local z pt
#
# . local covars cbd ptn worker weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none") ///
# > deriv(/xb = -1*`y'*exp(-1*{xb:} - {logey0})) ///
# > deriv(/logey0 = -1*`y'*exp(-1*{xb:} - {logey0}))
#
# Final GMM criterion Q(b) = 1.09e-28
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 6
# Number of moments = 6
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .0352185 .0098182 3.59 0.000 .0159752 .0544617
# cbd | -.008398 .0020172 -4.16 0.000 -.0123517 -.0044444
# ptn | -.0113146 .0021819 -5.19 0.000 -.015591 -.0070383
# worker | .6623018 .0519909 12.74 0.000 .5604015 .764202
# weekend | .3009323 .0362682 8.30 0.000 .2298479 .3720167
# -------------+----------------------------------------------------------------*
# /logey0 | .2654423 .1550127 1.71 0.087 -.0383769 .5692616
# ------------------------------------------------------------------------------*
# Instruments for equation 1: pt cbd ptn worker weekend _cons
#
# . * multiple phenotype examples
# .
# . webuse trip, clear
# (Household trips)
#
# .
# . ** without covariates
# .
# . ** ivpoisson fit
# . ivpoisson gmm trips (tcost cbd = ptn worker), multiplicative nolog
#
# Final GMM criterion Q(b) = 3.12e-29
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 3 Number of obs = 5,000
# Number of moments = 3
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4214074 .0700091 6.02 0.000 .2841921 .5586228
# cbd | .075126 .1304184 0.58 0.565 -.1804893 .3307413
# _cons | -5.359227 1.47909 -3.62 0.000 -8.258191 -2.460263
# ------------------------------------------------------------------------------*
# Instrumented: tcost cbd
# Instruments: ptn worker
#
# . lincom _cons, eform
#
# ( 1) [trips]_cons = 0
#
# ------------------------------------------------------------------------------*
# trips | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | .0047045 .0069584 -3.62 0.000 .0002591 .0854124
# ------------------------------------------------------------------------------*
#
# .
# . *** moment conditions 1
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . gmm (`y'*exp(-1*{xb:`x'}) - {ey0}), ///
# > instruments(`z') tracelevel("none")
#
# Final GMM criterion Q(b) = 5.27e-09
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 3
# Number of moments = 3
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4213734 .0699533 6.02 0.000 .2842675 .5584793
# cbd | .0751167 .1303905 0.58 0.565 -.1804441 .3306774
# -------------+----------------------------------------------------------------*
# /ey0 | .0047057 .0069581 0.68 0.499 -.0089319 .0183433
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker _cons
#
# .
# . *** moment conditions 2
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . gmm (`y'*exp(-1*{xb:`x'} - {logey0}) - 1), ///
# > instruments(`z') tracelevel("none")
#
# Final GMM criterion Q(b) = 1.21e-33
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 3
# Number of moments = 3
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4214074 .0700091 6.02 0.000 .2841921 .5586228
# cbd | .075126 .1304184 0.58 0.565 -.1804893 .3307413
# -------------+----------------------------------------------------------------*
# /logey0 | -5.359227 1.47909 -3.62 0.000 -8.258191 -2.460264
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker _cons
#
# .
# . ** including covariates
# . ivpoisson gmm trips weekend (tcost cbd = ptn worker), multiplicative nolog
#
# Final GMM criterion Q(b) = 2.33e-32
#
# note: model is exactly identified
#
# Exponential mean model with endogenous regressors
#
# Number of parameters = 4 Number of obs = 5,000
# Number of moments = 4
# Initial weight matrix: Unadjusted
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# trips | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252243 .0750159 5.67 0.000 .2781958 .5722528
# cbd | .0773626 .1414159 0.55 0.584 -.1998075 .3545328
# weekend | .2157507 .0804855 2.68 0.007 .058002 .3734994
# _cons | -5.484431 1.57664 -3.48 0.001 -8.574588 -2.394274
# ------------------------------------------------------------------------------*
# Instrumented: tcost cbd
# Instruments: weekend ptn worker
#
# . lincom _cons, eform
#
# ( 1) [trips]_cons = 0
#
# ------------------------------------------------------------------------------*
# trips | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# (1) | .0041509 .0065445 -3.48 0.001 .0001888 .0912389
# ------------------------------------------------------------------------------*
#
# .
# . *** moment conditions 1
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . local covars weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'}) - {ey0}), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 2.73e-12
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 4
# Number of moments = 4
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252238 .0750177 5.67 0.000 .2781917 .5722559
# cbd | .0773625 .1414325 0.55 0.584 -.1998401 .354565
# weekend | .2157505 .0804888 2.68 0.007 .0579953 .3735057
# -------------+----------------------------------------------------------------*
# /ey0 | .0041509 .006545 0.63 0.526 -.0086771 .0169789
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker weekend _cons
#
# .
# . *** moment conditions 2
# . local y trips
#
# . local x tcost cbd
#
# . local z ptn worker
#
# . local covars weekend
#
# . gmm (`y'*exp(-1*{xb:`x' `covars'} - {logey0}) - 1), ///
# > instruments(`z' `covars') tracelevel("none")
#
# Final GMM criterion Q(b) = 2.43e-32
#
# note: model is exactly identified
#
# GMM estimation
#
# Number of parameters = 4
# Number of moments = 4
# Initial weight matrix: Unadjusted Number of obs = 5,000
# GMM weight matrix: Robust
#
# ------------------------------------------------------------------------------*
# | Robust
# | Coef. Std. Err. z P>|z| [95% Conf. Interval]
# -------------+----------------------------------------------------------------*
# tcost | .4252243 .0750159 5.67 0.000 .2781958 .5722528
# cbd | .0773626 .1414159 0.55 0.584 -.1998075 .3545328
# weekend | .2157507 .0804855 2.68 0.007 .058002 .3734994
# -------------+----------------------------------------------------------------*
# /logey0 | -5.484431 1.57664 -3.48 0.001 -8.574588 -2.394274
# ------------------------------------------------------------------------------*
# Instruments for equation 1: ptn worker weekend _cons
# Single instrument example ----
test_that("Single instrument example", {
skip_on_cran()
# skip_if_not_installed("ivtools")
# Data generation from the example in the ivtools ivglm() helpfile
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
# ivtools for comparison fit
# library(ivtools)
# fitZ.L <- glm(Z ~ 1, data = dat)
# fitY.LZX <- glm(Y ~ X + Z + X*Z, family = "binomial", data = dat)
# fitLogGest <- ivglm(estmethod = "g",
# X = "X",
# fitZ.L = fitZ.L,
# fitY.LZX = fitY.LZX,
# data = dat,
# link = "log",
# Y = "Y")
logcrr <- 0.1314654 # fitLogGest$est["X"]
logcrrse <- 0.06035374 # sqrt(fitLogGest$vcov)
fit01 <- msmm(Y ~ X | Z, data = dat)
expect_equal(log(fit01$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit02 <- msmm(Y ~ X | Z, data = dat, estmethod = "gmm")
expect_equal(log(fit02$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit03 <- msmm(Y ~ X | Z, data = dat, estmethod = "gmmalt")
expect_equal(log(fit03$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit04 <- msmm(Y ~ X | Z, data = dat, estmethod = "tsls")
expect_equal(log(fit04$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
fit05 <- msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt")
expect_equal(log(fit05$crrci[1]), logcrr, tolerance = 0.05, ignore_attr = "names")
expect_s3_class(fit01, "msmm")
smy <- summary(fit01)
expect_s3_class(smy, "summary.msmm")
expect_output(print(fit01))
expect_output(print(smy))
})
# check subset argument ----
test_that("Check subset argument", {
skip_on_cran()
# Data generation from the example in the ivtools ivglm() helpfile
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
datfifty <- dat[1:50,]
fitcompare <- msmm(Y ~ X | Z, dat = datfifty)
fitcheck <- msmm(Y ~ X | Z, dat = dat, subset = 1:50)
expect_equal(fitcheck$crrci, fitcompare$crrci)
})
# check variables with different names ----
test_that("Check using different variable names", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$E <- dat$X
dat$R <- dat$Y
dat$W <- dat$Z
fit06 <- msmm(R ~ E | W, data = dat)
expect_equal(log(fit06$crrci[1]), 0.1314654, tolerance = 0.05, ignore_attr = "names")
})
# Multiple instrument example ----
test_that("Multiple instrument example", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX <- plogis(0.7*G1 + G2 - G3 + U)
X <- rbinom(n, 1, pX)
pY <- plogis(-2 + psi0*X + U)
Y <- rbinom(n, 1, pY)
dat <- data.frame(G1, G2, G3, X, Y)
fit11 <- msmm(Y ~ X | G1 + G2 + G3, data = dat)
expect_equal(log(fit11$crrci[1]), 0.45, tolerance = 0.05)
fit12 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "gmm")
expect_equal(log(fit12$crrci[1]), 0.45, tolerance = 0.05)
fit13 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "gmmalt")
expect_equal(log(fit13$crrci[1]), 0.45, tolerance = 0.1)
fit14 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "tsls")
expect_equal(log(fit14$crrci[1]), 0.45, tolerance = 0.02)
fit15 <- msmm(Y ~ X | G1 + G2 + G3, data = dat, estmethod = "tslsalt")
expect_equal(log(fit15$crrci[1]), 0.45, tolerance = 0.15)
})
# non-binary x with tsls, tslsalt methods fail ----
test_that("Non-binary x with tsls and tslsalt methods should produce an error", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$X[1] <- 2
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# non-binary y fail for tsls methods ----
test_that("Y needs to be binary for the TSLS methods", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$Y[1] <- 2
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# All methods fail for a non-integer in Y ----
test_that("Methods fail for non-integer Y", {
skip_on_cran()
set.seed(9)
n <- 1000
psi0 <- 0.5
Z <- rbinom(n, 1, 0.5)
X <- rbinom(n, 1, 0.7*Z + 0.2*(1 - Z))
m0 <- plogis(1 + 0.8*X - 0.39*Z)
Y <- rbinom(n, 1, plogis(psi0*X + log(m0/(1 - m0))))
dat <- data.frame(Z, X, Y)
dat$Y[1] <- 1.5
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "gmm"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "gmmalt"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X | Z, data = dat, estmethod = "tslsalt"))
})
# Multiple exposure example ----
test_that("Multiple exposure example", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U)
Y <- rbinom(n, 1, pY)
dat <- data.frame(G1, G2, G3, X1, X2, Y)
fit21 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat)
expect_equal(log(fit21$crrci[1,1]), 0.08, tolerance = 0.02)
fit22 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "gmm")
expect_equal(log(fit22$crrci[1,1]), 0.08, tolerance = 0.02)
fit23 <- msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "gmmalt")
expect_equal(log(fit23$crrci[1,1]), 0.08, tolerance = 0.1)
expect_error(msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "tsls"))
expect_error(msmm(Y ~ X1 + X2 | G1 + G2 + G3, data = dat, estmethod = "tslsalt"))
})
# Multiple exposure example with different variable names ----
test_that("Multiple exposure example with different variable names", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U)
Y <- rbinom(n, 1, pY)
E1 <- X1
E2 <- X2
R <- Y
dat <- data.frame(G1, G2, G3, E1, E2, R)
fit24 <- msmm(Y ~ E1 + E2 | G1 + G2 + G3, data = dat)
expect_equal(log(fit24$crrci[1,1]), 0.08, tolerance = 0.01)
})
# Example adjusting for a covariate ----
test_that("Adjusting for covariate", {
skip_on_cran()
set.seed(123456)
n <- 1000
psi0 <- 0.5
psi1 <- 0.4
G1 <- rbinom(n, 2, 0.5)
G2 <- rbinom(n, 2, 0.3)
G3 <- rbinom(n, 2, 0.4)
U <- runif(n)
C <- runif(n)
pX1 <- plogis(0.7*G1 + G2 - G3 + U + C)
X1 <- rbinom(n, 1, pX1)
pX2 <- plogis(-1 + 0.2*G1 - 0.2*G2 + 0.4*G3 + U + C)
X2 <- rbinom(n, 1, pX2)
pY <- plogis(-2 + psi0*X1 + psi1*X2 + U + C)
Y <- rbinom(n, 1, pY)
E1 <- X1
E2 <- X2
R <- Y
dat <- data.frame(G1, G2, G3, E1, E2, R, C)
fit25 <- msmm(Y ~ E1 + C | G1 + G2 + G3 + C, data = dat)
expect_equal(log(fit25$crrci[1,1]), .173, tolerance = .01)
fit26 <- msmm(Y ~ E1 + E2 + C | G1 + G2 + G3 + C, data = dat)
expect_equal(log(fit26$crrci[1,1]), .184, tolerance = .01)
})
Try the OneSampleMR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.