ERROR while using medsens for Sensitivity Analysis: subscript out of bounds

[Alaiti]

I'm having trouble with mediation:: medsens()

med.fit <- glm(ADM_ROTSUP_0_4S ~ fgrupo, data = dados, family="gaussian")

out.fit <- lm(VAS.0.8 ~ ADM_ROTSUP_0_4S +fgrupo + VAS.0 + ADM_ROTSUP_0_INICIAL + Idade..anos. + 
Tempo.de.dor..meses. + fsexo, data = dados)

med.out <- mediate(med.fit, out.fit, boot = TRUE, treat = "fgrupo", mediator = "ADM_ROTSUP_0_4S", robustSE = TRUE, sims = 1000)

The code works fine until med.out, but when I try to run a sensitivity analysis using medsens I get:

sens.out <- medsens(med.out, rho.by = 0.1, sims = 1000, effect.type = "indirect")
summary(sens.out)
ERROR: subscript out of bounds 

Someone might help me?

Thanks!

[kosukeimai]

It's difficult to know the reason, but I would not use the variable names like Idade..anos. which may cause some unexpeted behavior.

[liebscher]

I'm experiencing a similar error. Running mediation v4.5.0.

Case 1:

Replicating the framing data example from the mediation vignette, such that:

med.fit <- lm(emo ~ treat + age + educ + gender + income, data = framing)
out.fit <- glm(cong_mesg ~ emo + treat + age + educ + gender + income, data = framing, family = binomial("probit"))
med.out <- mediate(med.fit, out.fit, treat = "treat", mediator = "emo", robustSE = TRUE, sims = 100)

and then subsequently:

sens.out <- medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100)
summary(sens.out)

This works fine and as expected.

Case 2:

Simple mediation test, beginning with:

N <- 100

X <- c(rep("Treat", N/2), rep("Control", N/2))

sd <- 0.1
M <- c(rnorm(N/4, 4, sd), rnorm(N/4, 3, sd), rnorm(N/4, 2, sd), rnorm(N/4, 1, sd))

Y <- c(rnorm(N/4, 4, 0.5), rnorm(N/4, 3, 0.5), rnorm(N/4, 2, 0.5), rnorm(N/4, 1, 0.5))

then constructing some models:

mod_m_x <- lm(M ~ X)
mod_y_mx <- lm(Y ~ X*M)

and finally, building a simple mediation model:

s <- mediate(model.m = mod_m_x,
                model.y = mod_y_mx,
                treat = "X",
                mediator = "M",
                boot = T,
                boot.ci.type = "bca",
                sims = 1000)

Now, we find:

ms <- medsens(s)
summary(ms)

leads to the error: Error in m.coefs[T.out, ] : subscript out of bounds

While it could be due to variable names, I find that unlikely considering "X" and "M", with trivial levels names on the "X" factor, work fine for the rest of the package's functions and are as simple as one could hope.

I have tested not using bootstrap (I think one of the main differences between the two examples), which leads to no change. Another difference I see is the probit model in Case 1, compared to a gaussian model in the second. Advice would be appreciated, hopefully my code reproduces.

[kuriwaki]

(deleted my previous post) -- Correction: I can reproduce the case 2 error via reprex

library(mediation)
#> Loading required package: MASS
#> Loading required package: Matrix
#> Loading required package: mvtnorm
#> Loading required package: sandwich
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0

# Simple mediation test, beginning with:
N <- 100
X <- c(rep("Treat", N/2), rep("Control", N/2))
sd <- 0.1
M <- c(rnorm(N/4, 4, sd), rnorm(N/4, 3, sd), rnorm(N/4, 2, sd), rnorm(N/4, 1, sd))
Y <- c(rnorm(N/4, 4, 0.5), rnorm(N/4, 3, 0.5), rnorm(N/4, 2, 0.5), rnorm(N/4, 1, 0.5))

# then constructing some models:
mod_m_x <- lm(M ~ X)
mod_y_mx <- lm(Y ~ X*M)

# and finally, building a simple mediation model:

s <- mediate(model.m = mod_m_x,
             model.y = mod_y_mx,
             treat = "X",
             mediator = "M",
             boot = TRUE,
             boot.ci.type = "bca",
             sims = 5)
#> Warning in mediate(model.m = mod_m_x, model.y = mod_y_mx, treat = "X", mediator
#> = "M", : treatment and control values do not match factor levels; using Control
#> and Treat as control and treatment, respectively
#> Running nonparametric bootstrap

# medsens
medsens(s)
#> Error in m.coefs[T.out, ]: subscript out of bounds

^{Created on 2020-07-11 by the reprex package (v0.3.0)}

Session info

``` r devtools::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.0.0 (2020-04-24) #> os macOS Catalina 10.15.5 #> system x86_64, darwin17.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2020-07-11 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> acepack 1.4.1 2016-10-29 [1] CRAN (R 4.0.0) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.0) #> backports 1.1.8 2020-06-17 [1] CRAN (R 4.0.0) #> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.0.0) #> boot 1.3-24 2019-12-20 [1] CRAN (R 4.0.0) #> callr 3.4.3 2020-03-28 [1] CRAN (R 4.0.0) #> checkmate 2.0.0 2020-02-06 [1] CRAN (R 4.0.0) #> cli 2.0.2 2020-02-28 [1] CRAN (R 4.0.0) #> cluster 2.1.0 2019-06-19 [1] CRAN (R 4.0.0) #> colorspace 1.4-1 2019-03-18 [1] CRAN (R 4.0.0) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.0) #> data.table 1.12.8 2019-12-09 [1] CRAN (R 4.0.0) #> desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.0) #> devtools 2.3.0 2020-04-10 [1] CRAN (R 4.0.0) #> digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.0) #> dplyr 1.0.0.9000 2020-07-10 [1] Github (tidyverse/dplyr@f53e9ce) #> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.0) #> evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.0) #> fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.0) #> foreign 0.8-78 2020-04-13 [1] CRAN (R 4.0.0) #> Formula 1.2-3 2018-05-03 [1] CRAN (R 4.0.0) #> fs 1.4.1 2020-04-04 [1] CRAN (R 4.0.0) #> generics 0.0.2 2018-11-29 [1] CRAN (R 4.0.0) #> ggplot2 3.3.2 2020-06-19 [1] CRAN (R 4.0.0) #> glue 1.4.1 2020-05-13 [1] CRAN (R 4.0.0) #> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.0.0) #> gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.0) #> highr 0.8 2019-03-20 [1] CRAN (R 4.0.0) #> Hmisc 4.4-0 2020-03-23 [1] CRAN (R 4.0.0) #> htmlTable 2.0.1 2020-07-05 [1] CRAN (R 4.0.0) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 4.0.0) #> htmlwidgets 1.5.1 2019-10-08 [1] CRAN (R 4.0.0) #> jpeg 0.1-8.1 2019-10-24 [1] CRAN (R 4.0.0) #> knitr 1.28 2020-02-06 [1] CRAN (R 4.0.0) #> lattice 0.20-41 2020-04-02 [1] CRAN (R 4.0.0) #> latticeExtra 0.6-29 2019-12-19 [1] CRAN (R 4.0.0) #> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.0) #> lme4 1.1-23 2020-04-07 [1] CRAN (R 4.0.0) #> lpSolve 5.6.15 2020-01-24 [1] CRAN (R 4.0.0) #> magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.0) #> MASS * 7.3-51.5 2019-12-20 [1] CRAN (R 4.0.0) #> Matrix * 1.2-18 2019-11-27 [1] CRAN (R 4.0.0) #> mediation * 4.5.0 2019-10-08 [1] CRAN (R 4.0.0) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.0) #> minqa 1.2.4 2014-10-09 [1] CRAN (R 4.0.0) #> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.0) #> mvtnorm * 1.1-0 2020-02-24 [1] CRAN (R 4.0.0) #> nlme 3.1-147 2020-04-13 [1] CRAN (R 4.0.0) #> nloptr 1.2.2.1 2020-03-11 [1] CRAN (R 4.0.0) #> nnet 7.3-14 2020-04-26 [1] CRAN (R 4.0.0) #> pillar 1.4.5 2020-07-09 [1] CRAN (R 4.0.0) #> pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 4.0.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.0) #> pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.0) #> png 0.1-7 2013-12-03 [1] CRAN (R 4.0.0) #> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.0) #> processx 3.4.2 2020-02-09 [1] CRAN (R 4.0.0) #> ps 1.3.3 2020-05-08 [1] CRAN (R 4.0.0) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.0) #> R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.0) #> RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 4.0.0) #> Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.0) #> remotes 2.1.1 2020-02-15 [1] CRAN (R 4.0.0) #> rlang 0.4.7 2020-07-09 [1] CRAN (R 4.0.0) #> rmarkdown 2.1 2020-01-20 [1] CRAN (R 4.0.0) #> rpart 4.1-15 2019-04-12 [1] CRAN (R 4.0.0) #> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.0) #> rstudioapi 0.11 2020-02-07 [1] CRAN (R 4.0.0) #> sandwich * 2.5-1 2019-04-06 [1] CRAN (R 4.0.0) #> scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.0) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.0) #> statmod 1.4.34 2020-02-17 [1] CRAN (R 4.0.0) #> stringi 1.4.6 2020-02-17 [1] CRAN (R 4.0.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.0) #> survival 3.1-12 2020-04-10 [1] CRAN (R 4.0.0) #> testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.0) #> tibble 3.0.3.9000 2020-07-10 [1] Github (tidyverse/tibble@87bdf41) #> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.0) #> usethis 1.6.1 2020-04-29 [1] CRAN (R 4.0.0) #> vctrs 0.3.1 2020-06-05 [1] CRAN (R 4.0.0) #> withr 2.2.0 2020-04-20 [1] CRAN (R 4.0.0) #> xfun 0.13 2020-04-13 [1] CRAN (R 4.0.0) #> yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0) #> zoo 1.8-8 2020-05-02 [1] CRAN (R 4.0.0) #> #> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library ```

[kosukeimai]

Can you code the treatment as a factor variable and see if that solves the problem?

[kuriwaki]

That fixes it for this example! It can probably be documented in the help page.

library(mediation)
#> Loading required package: MASS
#> Loading required package: Matrix
#> Loading required package: mvtnorm
#> Loading required package: sandwich
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0

# Simple mediation test, beginning with:
N <- 100
X <- factor(c(rep("Treat", N/2), rep("Control", N/2)))
sd <- 0.1
M <- c(rnorm(N/4, 4, sd), rnorm(N/4, 3, sd), rnorm(N/4, 2, sd), rnorm(N/4, 1, sd))
Y <- c(rnorm(N/4, 4, 0.5), rnorm(N/4, 3, 0.5), rnorm(N/4, 2, 0.5), rnorm(N/4, 1, 0.5))

# then constructing some models:
mod_m_x <- lm(M ~ X)
mod_y_mx <- lm(Y ~ X*M)

# and finally, building a simple mediation model:

s <- mediate(model.m = mod_m_x,
             model.y = mod_y_mx,
             treat = "X",
             mediator = "M",
             boot = TRUE,
             boot.ci.type = "bca",
             sims = 5)
#> Warning in mediate(model.m = mod_m_x, model.y = mod_y_mx, treat = "X", mediator
#> = "M", : treatment and control values do not match factor levels; using Control
#> and Treat as control and treatment, respectively
#> Running nonparametric bootstrap

# medsens
summary(medsens(s))
#> 
#> Mediation Sensitivity Analysis: Average Mediation Effect
#> 
#> Sensitivity Region: ACME for Control Group
#> 
#>      Rho ACME(control) 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
#> [1,] 0.5        0.4905      -0.1423       1.1233         0.25       0.0113
#> [2,] 0.6        0.0922      -0.5388       0.7232         0.36       0.0162
#> [3,] 0.7       -0.4389      -1.0713       0.1935         0.49       0.0221
#> 
#> Rho at which ACME for Control Group = 0: 0.6
#> R^2_M*R^2_Y* at which ACME for Control Group = 0: 0.36
#> R^2_M~R^2_Y~ at which ACME for Control Group = 0: 0.0162 
#> 
#> 
#> Sensitivity Region: ACME for Treatment Group
#> 
#>      Rho ACME(treated) 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
#> [1,] 0.4        0.5516      -0.1058       1.2091         0.16       0.0072
#> [2,] 0.5        0.2265      -0.4290       0.8821         0.25       0.0113
#> [3,] 0.6       -0.1718      -0.8272       0.4836         0.36       0.0162
#> 
#> Rho at which ACME for Treatment Group = 0: 0.6
#> R^2_M*R^2_Y* at which ACME for Treatment Group = 0: 0.36
#> R^2_M~R^2_Y~ at which ACME for Treatment Group = 0: 0.0162

^{Created on 2020-07-11 by the reprex package (v0.3.0)}

[liebscher]

Good find, thank you!

[teppeiyamamoto]

Thanks for catching this. I'm guessing this is probably because stringsAsFactors now defaults to FALSE as of R 4.0. We'll look into it and release a fix.