kosukeimai / mediation Goto Github PK
View Code? Open in Web Editor NEWR package mediation
R package mediation
I have what seems a simple case that cannot properly fit.
I have a dependent variable that is measured 10 times per subject. However, the mediator is only once after the intervention (there are two groups: Control and Interventions)
I thought that the following would work:
model_m <- lm(mediator ~ group, data=unique (df) )
model_dv <- lmer(dv ~ group + mediator + (1 | subj), data = df)
mediate(model_m, model_dv, treat='group', mediator='mediator', boot=F)
To be sure: model_m is fit using lm
because the mediator is only measured once (note that data=unique(df)
), whereas model_dv is fitted using lmer
because it is measured 10 times (within subject).
It triggers this error, which makes sense: one group has 10 times more rows.
Instead, I tried the following:
model_m <- lmer(mediator ~ group + (1|subj), data=df )
model_dv <- lmer(dv ~ group + mediator + (1 | subj), data = df)
mediate(model_m, model_dv, treat='group', mediator='mediator', boot=F)
This works except that model.m
is incorrect because it is fitted in 10 times the measurements that actually were taken. As a consequence, there are several warnings displayed (e.g., model failed to converge which makes sense), *the estimates are the same as the ones obtained with lm
but the the standard errors are smaller.
Is it OK to still the model_m
(fitted with lmer
) even though some of the statistics are wrong? In other words, which statistics of the lmerMod object are used inside the mediate
function? (I can then check if they are all the same and it is just about the format of the data)
Hi there,
thanks for the package. It is insightful and helpful.
My outcome variable has different characteristics for the bottom and upper quantiles. Unfortunately, multimed function doesn't have tau.y argument unlike mediate function. Do you have any idea how to deal with this issue?
Sincerely,
At present, summary(mediation(...)) produces output like this:
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.1477 0.0237 0.28 0.018 *
ADE 0.2998 0.1166 0.50 <2e-16 ***
Total Effect 0.4475 0.2293 0.66 <2e-16 ***
Prop. Mediated 0.3295 0.0741 0.60 0.018 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 200
Simulations: 1000
It would be nice if this could be supplemented by a tidy data frame, such that this would work:
as_data_frame(summary(mediation(...)))
To produce a data frame with these columns:
(these chosen to match those typically produced by broom::tidy)
Thank you for this great package.
The current mediate() function may be extended further by adding a function that computes odds ratios and/or risk ratios (with confidence intervals) for a dichotomous outcome. In the medical literature, we are often asked to report effect measures in these scales in addition the difference. Thank you very much for your consideration.
This is most probably my own misunderstanding, but I'm wondering whether bootstrap CIs are correctly implemented. Two questions:
z.inv <- length(theta[theta < mean(theta)])/sims
However, eq. 2.8 in DiCiccio & Efron, 1996 is different: theta is compared to the original, non-bootstrapped statistic. The current implementation assumes that the mean of the bootstrap distribution equals to the original statistic. Is this necessarily true?
U <- (sims - 1) * (mean(theta) - theta)
However, eq. 6.7 uses jackknife estimator, not the bootstrap estimator. Is the usage of the bootstrap estimator here justified?
This has caused some problems from some users.
Can be replicated from the package example file (below). I suggest using the same rounding for upper as we do for lower.
summary(contcont)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.0183 -0.0397 0.00 0.04 *
ADE -0.0450 -0.1417 0.03 0.40
Total Effect -0.0634 -0.1636 0.01 0.16
Prop. Mediated 0.2686 -0.9015 2.54 0.20
Hello,
I'm running a mediation analysis and trying to increase the number of simulations, but I run into an error when my number of simulations reaches n ~= 700.
M2 <- lm( mediator~ temperature, dsite )
M3 <- glmer( rate~ temperature+ mediator + (1|Site), d, family="binomial" )
med <- mediation::mediate( M2,M3, sims=700, treat="temperature", mediator="mediator" )
summary(med)
This produces the following error: Error in if (xhat == 0) out <- 1 else { : missing value where TRUE/FALSE needed
Note that the mediator and treatment here are both continuous variables, so the medation model is a linear model. The mediator model also comes from data at a higher level than the outcome model, so the datasets are slightly different. Interesting though that no error is produced when the number of simulations is lower (e.g., 500). Can you help me understand why I might be getting these errors?
Here is the output of a working mediation model with a low number of simulations:
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediator Groups:
Outcome Groups: Site
Output Based on Overall Averages Across Groups
Estimate 95% CI Lower 95% CI Upper p-value
ACME (control) 4.49e-03 2.28e-06 0.03 <2e-16 ***
ACME (treated) 3.80e-03 1.77e-06 0.02 <2e-16 ***
ADE (control) -2.42e-03 -1.65e-02 0.00 0.02 *
ADE (treated) -3.11e-03 -2.06e-02 0.00 0.02 *
Total Effect 1.38e-03 1.27e-06 0.01 0.02 *
Prop. Mediated (control) 1.85e+00 1.05e+00 4.93 0.02 *
Prop. Mediated (treated) 1.51e+00 1.03e+00 4.23 0.02 *
ACME (average) 4.15e-03 2.03e-06 0.02 <2e-16 ***
ADE (average) -2.76e-03 -1.85e-02 0.00 0.02 *
Prop. Mediated (average) 1.68e+00 1.04e+00 4.56 0.02 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 234
Simulations: 100
I'm also interested to know what estimates of "Prop. Mediated" is greater than 1. Should I be worried about that? I see some examples in the vignette where at least the upper confidence limit exceeds 1, but I'm not sure what this reveals about the modeling exercise or the underlying data.
Thank you!
I am conducting a mediation analysis with multiply imputed 30 data sets using "amelidiate {mediation}".
variables: mediators= m,outcome=y, treatment=rep(e,30),covariates=z, control_val=a0, treat_val=a1,
medout <- mediations(datasets, treatment, mediators, outcome, covariates,
families=c("binomial","gaussian"), control.value = control_val,
treat.value = treat_val, interaction=TRUE, conf.level=.95, sims=1000, set.seed(1))
med_pool <- amelidiate(medout)
summary(med_pool)
The package provides pooled CI for all estimates but does not provide pooled p-values. Would there any way to calculate pooled p-values?
Thank you so much!
Hi!
I love the package, thanks for it (and the theoretical work underlying it)
Is it possible to extract a standard path diagram from the output of mediation()
?
I believe reviewers will require these standard graphs because the readers are used to those.
Thanks for the help!
Kind regards
Jana
Hi,
I just re-ran an old analysis which functioned without problems in November 2020 but returns this error now Error in eval(extras, data, env) : object 'Call.M' not found
.
I did not change the code at all.
I already tried reinstalling the package, but that didn't help.
Could you help me with this issue?
Best regards
Leo
This works fine with boot=F, but throws an error when boot=T.
Seems related to #29
library(mediation)
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0
test <- function(f1, f2, boot){
mtcars$mpg <- mtcars$mpg/max(mtcars$mpg)
m_med <- lm(f1, mtcars)
m_out <- lm(f2, mtcars)
mediate(m_med, m_out, treat = "wt", mediator = 'am', boot = boot)
}
x <- test('am ~ wt', 'mpg ~ wt + am', F)
#success!
x <- test('am ~ wt', 'mpg ~ wt + am', T)
#> Running nonparametric bootstrap
#> Error in stats::model.frame(formula = f1, data = structure(list(am = c(1, : object 'f1' not found
Hello,
Thank you for developing such a great package. It is highly useful in implementing many different kinds of mediation analyses, and I'm getting a lot of help for my projects.
But recently, I ran into trouble while using the "ivmediate" function.
The codes I ran was like this:
`############################################## IV Mediation with mediation package ##############################
library(mediation)
library(car)
library(fastDummies)
treatment = 'rh_adi_bi_0y'
instrument = 'section8_0y'
mediator = 'discount_rate_1y_z'
outcome = 'totalscore_pps_1y_z'
data<-dummy_cols(data, select_columns = c('sex_0y'))
data$sex_0y<-as.factor(data$sex_0y)
data$married_0y<-as.factor(data$married_0y)
data$race_g<-as.factor(data$race_g)
data$parent_identity_0y<-as.factor(data$parent_identity_0y)
data$gender_identity_0y<-as.factor(data$gender_identity_0y)
data$foreign_born_0y<-as.factor(data$foreign_born_0y)
data$foreign_born_family_0y<-as.factor(data$foreign_born_family_0y)
data$gay_parent_0y<-as.factor(data$gay_parent_0y)
data$gay_youth_0y<-as.factor(data$gay_youth_0y)
data$race_ethnicity_0y<-as.factor(data$race_ethnicity_0y)
fmla1 <- as.formula("rh_adi_bi_0y~section8_0y+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")
fmla2 <- as.formula("discount_rate_1y_z~rh_adi_bi_0y+section8_0y+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")
fmla3 <-as.formula("totalscore_pps_1y_z ~ discount_rate_1y_z*(rh_adi_bi_0y+section8_0y)+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")
model_treatment <- glm(fmla1, data = data, family = binomial)
model_mediator <- lm(fmla2, data = data)
model_outcome <- lm(fmla3, data = data)
ivmediation_result <- ivmediate(model_treatment, model_mediator, model_outcome, sims = 500,
enc = 'section8_0y', treat='rh_adi_bi_0y', mediator='discount_rate_1y_z')
summary(ivmediation_result)
plot(ivmediation_result)`
and the error messages that I got were like this:
Error in integrate(LACME1.integrand, -Inf, Inf, i = i, stop.on.error = F) : evaluation of function gave a result of wrong type In addition: Warning messages: 1: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 2: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 3: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 4: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 5: glm.fit: algorithm did not converge 6: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 7: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead.
I think I've done the same thing as shown in the example below, but this just doesn't seem to work for some reason.
https://www.rdocumentation.org/packages/mediation/versions/4.5.0/topics/ivmediate
`a <- lm(comply ~ treat + sex + age + marital + nonwhite + educ + income,
data = jobs)
b <- glm(job_dich ~ comply + treat + sex + age + marital + nonwhite + educ + income,
data = jobs, family = binomial)
c <- lm(depress2 ~ job_dich * (comply + treat) + sex + age + marital + nonwhite + educ + income,
data = jobs)
out <- ivmediate(a, b, c, sims = 50, boot = FALSE,
enc = "treat", treat = "comply", mediator = "job_dich")
summary(out)
plot(out)`
Have I missed something, or is there some problem with the ivmediate function that makes ends up with this error?
I'll really appreciate it if you can help.
Thanks in advance!
I believe there must be an error which is causing the sign of the reported ACME
effect to change depending on the contrast coding matrix for X
, when both X->M (a
) and M->Y (b
) are negative, regardless of how control.value
and treat.value
are specified in the call to mediate
.
In the experimental data during the analysis of which I encountered this issue, X
represents a two-contrast factor with english
and dutch
as the levels. I am interested in setting english
as the control and testing for the effect of dutch
compared with english
. Note that this implies changing the factor contrast coding from the default, because E comes after D, and would usually default as the treatment level.
I simulate these data below:
simulate_data = function(x1){
# x2, x3, and x4 in a matrix, these will be modified to meet the criteria
x2 <- scale(matrix( rnorm(length(x1)), ncol=1))
# put all into 1 matrix for simplicity
x12 <- cbind(scale(x1),x2)
# find the current correlation matrix
c1 <- var(x12)
# cholesky decomposition to get independence
chol1 <- solve(chol(c1))
newx <- x12 %*% chol1
# create new correlation structure (zeros can be replaced with other rvals)
newc <- matrix(c(1, -0.5,
-0.5, 1),
ncol=2)
chol2 <- chol(newc)
finalx <- newx %*% chol2 * sd(x1) + mean(x1)
finalx[,2]
}
data = data.frame("X"=c(rep('B_english',300),rep('A_dutch',300)),
"M"=scale(c(rnorm(300,mean=150, sd = 50),
rnorm(300,mean=100,sd=50))
))
data$Y = simulate_data(data[,'M'])
#set contrasts
data$X_rev = factor(data$X)
contrasts(data$X_rev) = contr.treatment(2,base=2)
colnames(contrasts(data$X_rev)) = 'A_dutch'
data$X = factor(ifelse(data$X == 'B_english','A_english','B_dutch'))
> contrasts(data$X)
B_dutch
A_english 0
B_dutch 1
> contrasts(data$X_rev)
A_dutch
A_dutch 1
B_english 0
The only difference between data$X
and data$X_rev
is the letter they start with. In both of the contrast matrices, the level representing dutch
is coded as the treatment contrast, so I would expect that this difference should not result in any substantive changes to model outputs. However this does not appear to be true.
When X = data$X
, the sign for ACME
is correct (i.e. positive since a
and b
are both negative):
### X = data$X
med.fit = lm(M ~ X,data = data)
out.fit = lm(Y ~ M + X,data = data)
med.out = mediate(med.fit,out.fit,treat='X',mediator='M',sims=100,
control.value = 'A_english', treat.value = 'B_dutch')
>summary(med.out)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.4568 0.3632 0.55 <2e-16 ***
ADE -0.0764 -0.2247 0.07 0.38
Total Effect 0.3804 0.2186 0.54 <2e-16 ***
Prop. Mediated 1.2186 0.8517 1.98 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 600
Simulations: 100
However, in the case of X = data$X_rev
, the ACME
is reported as negative when it should be positive:
med.fit.rev = lm(M ~ X_rev,data = data)
out.fit.rev = lm(Y ~ M + X_rev,data = data)
med.out.rev = mediate(med.fit.rev,out.fit.rev,treat='X_rev',mediator='M',sims=100,
control.value = 'B_english', treat.value = 'A_dutch')
> summary(med.out.rev)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.4670 -0.5757 -0.35 <2e-16 ***
ADE 0.0787 -0.0637 0.19 0.24
Total Effect -0.3883 -0.5335 -0.22 <2e-16 ***
Prop. Mediated 1.2161 0.8620 1.83 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 600
Simulations: 100
This would seem to be undesirable and/or confusing behaviour. I'm sorry if I missed somewhere in the docs where it specifies that factor contrasts must be coded to the R
default (i.e. in alphabetical order). I would suggest that this point either be made more salient in the vignette, or that this behaviour be corrected.
Mediate returns NaNs for the CIs if I run a model with ordered outcomes and the bca option.
Here a reproducible example:
library(mediate)
library(MASS)
dv <- sample(1:4, size=200, replace=TRUE)
m <- sample(1:4, size=200, replace=TRUE)
t <- sample(0:1, size=200, replace=TRUE)
ex.dat <- data.frame(dv,m,t)
ex.dat$dv <- as.factor(ex.dat$dv)
ex.dat$m <- as.factor(ex.dat$m)
med.m <- polr(m ~ t, data = ex.dat)
summary(med.m)
out.m <- polr(dv ~ t + m, data = ex.dat)
summary(out.m)
out.ord <- mediate(med.m, out.m, treat = "t",
mediator = "m", sims = 100,
boot = T, boot.ci.type = "bca")
summary(out.ord)
It produces an output, but the CIs are missing.
It returns the following warning:
Warning messages:
1: In qnorm(z.inv) : NaNs generated
2: In qnorm(z.inv) : NaNs generated
3: In qnorm(z.inv) : NaNs generated
4: In qnorm(z.inv) : NaNs generated
5: In qnorm(z.inv) : NaNs generated
Could Package "mediate" support multi-categorical treatment?
I could find arguments 'control.value' & 'treat.value' in the function mediate, it seems that the function could support binary treatment. I am not sure whether mediate can support multi-categorical treatment (# groups >= 3) as well. Thank you for your help in advance!
David
How can I perform a mediation analysis with a continuous treatment. I want to esimate indirect effect when exposure increase 1 unit instead of specifying control.value and treat.value.
Hi!
Thank you very much for providing the mediation package to R, I really enjoy working with it! I was wondering, however, if it is somehow possible to extend the mediation package to also accept GAMM models, in one way or another.
Being aware that "mediation" does not handle "gamm" models directly, I tried to circumvent the problem by implementing my design with the mgcv "gam" function by using bs="re" flag as specification (as suggested by S. Wood , https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/random.effects.html). Unfortunately, it appears that using this approach the "predict.gam" (as currently used in "mediation") fails, and I get the following error message:
"Error in predict.gam(new.fit.M, type = "response", newdata = pred.data.t) :
403, 407, 418, 421,..."
In advance, thank you for your time!
Best regards
René
I was trying to fit a multilevel mediation analysis with crossed random factors using the mediate function from the mediation package in R. For the med.fit and out.fit models, I used lmer. However, I received the following error message:
"Error in mediate(med.fit, out.fit, treat = "condition_numfct", mediator = "attr_rating", :
mediate does not support more than two levels per model."
I am not sure, but I suspect that this error may be due to the mediate function not supporting crossed random factors. Can anyone confirm this? If so, does anyone know of any alternative packages or functions that can handle multilevel mediation analyses with crossed random factors?
Thank you for any help or insights you can provide!
I'm running two mediations on the same data, once using prior (sampling) weights when fitting the models (through the survey
package), and once without those prior weights (just using the glm()
function). My outcome is binary. My mediator is a proportion out of 3, so I also pass 3 as the weights argument (or rather, the name of a column that is equal to 3 for every subject). When using the models fit with survey:svyglm()
inside mediate()
, there is no issue; however, when using the glm()
models, mediate()
complains that "weights on outcome and mediator models not identical".
After checking the documentation, I see that you're not supposed to use a multiple-trial binomial mediator, as that is not implemented yet. My question is, why does it work when the models are fit through survey::svyglm()? And can I trust those results at all?
I am having some trouble using your mediate() function with a Generalized Additive Model(GAM) and see that your documentation specifies it can be done, but does not explain the appropriate syntax. Specifically, I am trying to run a causal mediation analysis with a GAM with an interaction between a continuous variable (time) and two factor categorical variable (illness group). The model aims to asses the relationship of how gene expression (continuous variable) differs over time (also a continuous variable) by a two factor categorical variable (Illness group) and whether this is mediated by viral load (a continuous variable), which differs overtime between the two illness groups. The model syntax for the GAMs is here:
model.mediate <- gam(viral_load~illnessgroup+ s(time) +s(time, by = illnessgroup), data=data, method = "REML")
model.outcome <- gam( gene_exp~illnessgroup + viral_load + s(time) +s(time, by = illnessgroup) ,data=data, method = "REML")
mediation_result <- mediate(model.m = model.mediate, model.y = model.outcome, mediator = "viral_load", treat = "illnessgroup", treat.value = "Severe", control.value = "Not_Severe", boot = TRUE, control = "Not_Severe")
This gives me an output with the average ACME, ADE, Prop. mediated and Total Effects but I was expecting to see separate statistics for the treatment (Severe) and control (Not Severe) conditions. This is what we would expect if using a linear model rather than a GAM. Any advice on how to appropriately run the mediate() function with a GAM model using an interaction for 2 comparison groups to get the output for both the treatment and control groups would be much appreciated.
Thank you,
Basilin Benson
When calling summary
on a mediate
object, the underlying code uses the printCoefmat
method. By default, this assumes the results of the second-last column (i.e. the upper-confidence interval here) are a test statistic, which is reported with less precision than other values such the lower confidence interval. That is inconsistent between the bounds, but can also lead to rounding to 0, which can lead to inconsistencies between the confidence intervals and the p-values. For instance
> dat <- read.csv("https://mspeekenbrink.github.io/sdam-jasp-companion/data/redist2015.csv")
> mod2 <- lm(SC_mean ~ income, data=dat)
> mod3 <- lm(redist ~ SC_mean + income, data=dat)
> set.seed(20201027)
> med <- mediate(model.m = mod2, model.y = mod3, sims=1000, boot = TRUE, boot.ci.type = "bca", treat = "income", mediator = "SC_mean")
Running nonparametric bootstrap
> summary(med)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.00229 -0.00463 0.00 0.014 *
ADE -0.00292 -0.00650 0.00 0.078 .
Total Effect -0.00521 -0.00801 0.00 <2e-16 ***
Prop. Mediated 0.43953 0.12895 1.29 0.014 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 301
Simulations: 1000
It would be better if the upper confidence interval had the same precision as the lower CI. An easy fix (I believe) is to indicate that there is no test statistic, replacing all calls to printCoefmat
to
printCoefmat(smat, tst.ind=NULL)
Removing the digits=3
argument that is used at the moment also allows the digits to be controlled by the general options, rather than hardcoded. E.g.
> x <- med
> smat <- c(x$d1, x$d1.ci, x$d1.p)
> smat <- rbind(smat, c(x$z0, x$z0.ci, x$z0.p))
> smat <- rbind(smat, c(x$tau.coef, x$tau.ci, x$tau.p))
> smat <- rbind(smat, c(x$n0, x$n0.ci, x$n0.p))
> rownames(smat) <- c("ACME", "ADE",
+ "Total Effect", "Prop. Mediated")
> clp <- "95"
> colnames(smat) <- c("Estimate", paste(clp, "% CI Lower", sep=""),
+ paste(clp, "% CI Upper", sep=""), "p-value")
> printCoefmat(smat, tst.ind=NULL)
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.00228812 -0.00462947 -0.00046455 0.014 *
ADE -0.00291772 -0.00649948 0.00024205 0.078 .
Total Effect -0.00520584 -0.00801499 -0.00287072 <2e-16 ***
Prop. Mediated 0.43953008 0.12894917 1.29021831 0.014 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Hi there,
thanks for the package. I managed to use it for a single mediator analysis, but I can't seem to get it working for two mediators, using the mediators function I constantly get this error - Error in rep(1, nrow(dataarg)) : invalid 'times' argument
This is my code :
mediations(datasets = positive, treatment = c("factor_scores"), mediators = c("downward_comparison", "upward_comparison"), outcome = c("self_esteem"), families=c("gaussian","gaussian"), LowerY = NULL, UpperY = NULL, interaction=FALSE, conf.level=.90, sims=100, boot = FALSE, weights=NULL)
There are no NA values in the columns. What am I missing?
Thanks,
Ivaylo
I reduced the code as far as I could and inserted a browser() statement into med.fun.
med.fun has 3 arguments: y.data, m.data and index.
boot requires the function to have an argument "data", but med.fun doesn't.
Therefore, y.data and m.data are not resampled.
Also, new.fit.M <- eval.parent(Call.M)
is problematic, as is the equivalent for Y.
In particular, Call.M and Call.Y reference the original data in the stacked environments. So both models are refit to the original data, not to y.data and m.data.
One approach to fixing this is to add
y.data <- y.data[index, ]
m.data <- m.data[index, ]
near the top of med.fun and to do
model.m <- update(model.m, data = m.data)
model.y <- update(model.y, data = y.data)
instead of re-evaluating Call.M and Call.Y.
Finally, because we're dealing with a treatment and control, I think the observed group sizes should be conditioned on, so the function should default to passing strata = y.data[, treat] into the call to boot (though you can argue I'm wrong about that).
Harry
When trying to run the following code for non-parametric mediation the function "mediate.binary" no longer works. Has the name of this function been updated or changed? Below is the replication code for 'figure 3' that I am having issues with
Imai, Keele, & Tingley (2010):
#Initialize set of models
a<-lm(depress2 ~ treat + depress1 +age+educ+sex+nonwhite+econ_hard+marital+occp+income, data=job)
b <- lm(job_seek ~ treat + depress1 +age+educ+sex+nonwhite+econ_hard+marital+occp+income, data=job)
c <- gam(depress2 ~ treat +s(job_seek, bs="cr")+ depress1+age+educ+sex+nonwhite+econ_hard+marital+occp+income , data=job)
#No interaction assumption
mod.1 <- mediate.binary(b, c, sims=sims, boot=TRUE, T="treat", M="job_seek")
summary(mod.1)
Error in mediate.binary(b, c, sims = sims, boot = TRUE, T = "treat", M = "job_seek") :
could not find function "mediate.binary"
I recently updated R and my packages and the below code, which was previously working, is no longer working:
fit.dv <- lme4::lmer(opti3 ~ no_recalled_notes + N + condition + log(attempt) + log(attempt):N + (1|unique_melody_name), data = na.omit(main2))
fit.mediator <- lme4::lmer(no_recalled_notes ~ N + condition + log(attempt) + log(attempt):N + (1|unique_melody_name), data = na.omit(main2))
results <- mediation::mediate(fit.mediator, fit.dv, treat = c('N', 'condition'), mediator = 'no_recalled_notes')
The error is:
Error in if(!INT & isGam.y) { : the condition has length > 1
As title.
Hi, I know there is a function mediations() to perform "multiple" mediation models, and each of them is with a single mediator.
However, if I want to perform a mediaiton model with multiple mediators (like multiple mechanisms in a model, and they are uncorrelated with each other), how could I do that in this package? Thanks!
The formula
object in R allows variables to be specified flexibly e.g. y ~ factor(X)
or y ~ log(X)
. Doing this when fitting models for mediate
, however, results in error.
This example from the package documentation works as expected:
library("mediation")
data("jobs")
set.seed(4646)
model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
mediator = "job_seek")
out.2 <- mediate(model.m, model.y, sims = 10, treat = "treat",
mediator = "job_seek")
But would break if we replace sex
with factor(sex)
:
> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + factor(sex) + age
+ + occp + marital + nonwhite + educ + income, data = jobs)
> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ + occp + marital + nonwhite + educ + income, data = jobs)
>
> out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
+ mediator = "job_seek")
Running nonparametric bootstrap
Error in factor(sex) : object 'sex' not found
Or, if we try replacing age with log(age):
> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + log(age)
+ + occp + marital + nonwhite + educ + income, data = jobs)
> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ + occp + marital + nonwhite + educ + income, data = jobs)
>
> out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
+ mediator = "job_seek")
Running nonparametric bootstrap
Error in eval(predvars, data, env) : object 'age' not found
The root cause has to do with how model.frame()
handles column names. For example:
> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + factor(sex) + age
+ + occp + marital + nonwhite + educ + income, data = jobs)
> test <- predict(model.frame())
Error in terms.formula(formula, data = data) :
argument is not a valid model
> test <- predict(model.m)
> test <- predict(model.m, newdata = model.frame(model.m))
Error in factor(sex) : object 'sex' not found
I am using the mediations() and amelidiate() functions to estimate causal mediated effects with 5 imputed datasets. After running the amelidiate() function on the 5 mediation analyses from the mediations() function, the estimates are equivalent to the mediation based on the first imputed dataset. I did some digging in the source code and realized that the amelidiate() function uses the names from mediations() to organize the data but the names of the 5 mediation analyses are identical so it is only using the first dataset. I am going to rename the resulting list from the mediations() function so that I can use the amelidiate() function. The main application for these functions seems to be a scenario with multiple outcome or mediators and so that resulting names would be different in that case but multiple imputations is also mentioned as an application so I thought this might be a useful comment for future releases of the package or for anyone who is having a similar problem.
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!
Lines 78 to 84 in a9a5c7b
Hello,
I am performing a mediation analysis with mixed models in the presence of a (potential) treatment * mediator interaction. I would like to test if this interaction is significant. Is there anything wrong with using the same logic in this code to preform this test? My model.m and model.y that I provide to the mediate function are both lmer objects, and test.TMint says that it is not implemented for these models quite yet.
So my question is, is there anything theoretically wrong with this:
mod.m <- lmer(M ~ X + (1|id), data = data)
mod.y <- lmer(Y ~ M * X + (1|id), data = data)
med.results <- mediate(mod.m, mod.y, mediator = "M", treat = "X" )
acme_diff = med.results$d1.sims - med.results$d0.sims
quantile(acme_diff, probs = c(0.025, 0.975))
[Error in Surv(time, event) : cannot find object "time"] in survreg()
Here I am doing mediation analysis among X, Y, and m.
Here model.m is linear model and model.y is survival model.
The code is as follow:
model.m <- lm(m ~ X, data = data)
model.y <- survreg(Surv(time, event) ~ X+m, data = data)
mediation <- mediate(model.m, model.y, treat = "X", mediator = "m" , boot = T, outcome = "event")
However, error occurred as Error in Surv(time, event) : cannot find object "time
.
Running the following model, with MapChange as a binary factor
medModel <-
lme4::glmer( MapChange ~ TallCondition + ( 1 | V1.ResponseID ), data=bin.tallSonaK, family=binomial(link="logit") )
outModel <-
lme4::lmer( SGAIN ~ MapChange + TallCondition + ( 1 | V1.ResponseID ), data=bin.tallSonaK )
med <-
mediation::mediate(model.m=medModel, model.y=outModel,treat='TallCondition',mediator='MapChange',data=bin.tallSonaK)
summary(med)
yields the error message Error in if (xhat == 0) out <- 1 else { : missing value where TRUE/FALSE needed
I believe this is an error in function pval
:
pval <- function(x, xhat){
## Compute p-values
if (xhat == 0) out <- 1
else {
out <- 2 * min(sum(x > 0), sum(x < 0)) / length(x)
}
return(min(out, 1))
}
caused by xhat
having value NA.
Interestingly I can remove this error by keeping MapChange
as num
instead of factor
Is it possible to use the mediate function for multi-level outcome variables?
I want to do a mediation analysis, the outcome model is a survival model, but the document of this packages shows the outcome model can be survreg object instead of coxph object. Cox ph model can be used as outcome model?
Hi,
I was wondering if there was a potential bug with how headers are reported out for the summary of a mediation object?
I have a 4-level ordinal outcome (No (1), Mild (2), Moderate (3), Severe (4)) analyzed as a factor, a continuous mediator & a binary treatment.
I ran a glm for my mediator ~ treatment & a polr for my outcome ~ treatment + mediator.
When I run the mediation analysis using the 4-level factor outcome where the levels are coded as text, and run summary(), it reports out this way:
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Pr(Y=Mild) Pr(Y=Moderate) Pr(Y=No) Pr(Y=Severe)
ACME (control) -0.01801 0.00874 0.00447 0.00480
2.5% -0.02936 0.00262 0.00134 0.00139
97.5% -0.00561 0.01455 0.00773 0.00842
p-value 0.00200 0.00200 0.00200 0.00200
But when I change my outcome to a 4 level factor outcome where the levels are coded as numbers (1 (No), 2 (Mild), 3 (Moderate), 4 (Severe), it reports out this way:
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Pr(Y=1) Pr(Y=2) Pr(Y=3) Pr(Y=4)
ACME (control) -0.01801 0.00874 0.00447 0.00480
2.5% -0.02936 0.00262 0.00134 0.00139
97.5% -0.00561 0.01455 0.00773 0.00842
p-value 0.00200 0.00200 0.00200 0.00200
Is this expected behavior or have I made a mistake?
Thanks!
-Stephanie
I am working on a package that takes user input and within its functions, uses lme4 to fit and return a "merMod" object.
I tried to pass these models to mediate and was getting a "mediator model is not yet implemented" error. Poking through the source code, I see why. There are a series of if statements like
if(isMer.m && getCall(model.m)[[1]]=="lmer")
The models I've passed it, though, have a different 1st element, such that getCall(model.m)[[1]]
would be "lme4::lmer"
rather than "lmer"
. Of course, mediate has the technical capacity to fit these models but just doesn't recognize them for what they are.
I think it may be as simple of a fix as changing the if statements to something like
if(isMer.m && class(model.m)[[1]]=="lmerMod")
I am running into an issue when using mediations command to run the mediation analysis on multiple imputation data. I have weights in each of my datasets, but the packages seems to not be location them, as the error reads:
Error in eval(extras, data, env) : object 'weight' not found
The variable ipw_response_t are the truncated weights and are a numeric column in each of my imputed datasets. I specify it as below:
olsols <- mediations(datasets, treatment, mediators, outcome, covariates,
families=c("gaussian","gaussian")
, interaction=FALSE
, conf.level=.95
, sims=100,weights=c("ipw_response_t"))
I'm having trouble with mediation::mediate()
nb.terror.y.anx <- glm(data = terror_removed2,
certain~cond + event + WC + anx,
control = glm.control(maxit = 50000),
family = "poisson")
nb.terror.m.anx <- glm(data = terror_removed2,
anx~cond + event + WC + certain,
control = glm.control(maxit = 50000),
family = "poisson")
med.anx <- mediation::mediate(nb.terror.m.anx,nb.terror.y.anx, treat = "cond",
mediator = "anx", sim = 2)
Works fine, but whenever sim>2, I get:
Error in if (xhat == 0) out <- 1 else { :
missing value where TRUE/FALSE needed
I saw there's still an open topic in a realted issue
Hope it's ok to flag it up again.
Thanks!
Can you please create a conda wheel so I can install this R package using anaconda? Thanks!
Hello,
I am trying to run the mediate function using a mediator model of class "lmer" (from lme4) and an outcome model of class "gamm", due to nested random effects. However, only a class "gam" is compatible. Is it possible to add "gamm" as an allowable model class?
Thanks!
Ashley
Hi, does medsens function support the sensitivity analysis of causal mediation analysis with multilevel data? Thank you so much!
Due to the essence of multilevel data, I use lmer function to obtain med.fit and out.fit. Then, use these two outputs as inputs of the mediate function to get med.out. Finally, I try to use medsens function to perform the sensitivity analysis.
However, I get the error message below:
Error in medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100) :
mediate object fitted with non-supported model combinations
thx!
Hi! Thank you so much for this wonderful package. I was wondering if it would be possible to add support for models that require cluster-robust estimates, specifically lm_robust
models from the estimatr
package?
Using the mediate function (please see my code below), I receive the following error message. I would like to estimate a multilevel mediation for the following 2x5 mixed design. While the IV is manipulated between-subjects, the DV is a repeated measure (binary). In particular, each participant made 5 binary choices for 5 different choice sets (within-subject). The mediator (i.e., perceived expertise) is a repeatedly measured continuous variable, which was measured after each choice.
Do you have any idea what could be wrong with my code?
Error in if (xhat == 0) out <- 1 else { :
missing value where TRUE/FALSE needed
med.fit <- lmer(Expertise ~ Condition + (1|SubjectID), data = data)
out.fit <- glmer(Choice ~ Expertise + Condition + (1 + Expertise|SubjectID),
family = binomial(link = "logit"), data = data)
med.out <- mediate(med.fit, out.fit, treat = "Condition", mediator = "Expertise",
sims = 5000)
summary(med.out)
The data looks like the following:
SubjectID I ChoiceSet I Condition I Expertise I Choice
1 I 1 I 0 I 5 I 1
1 I 2 I 0 I 4 I 1
1 I 3 I 0 I 4 I 0
1 I 4 I 0 I 4 I 0
1 I 5 I 0 I 3 I 0
2 I 1 I 1 I 3 I 1
2 I 2 I 1 I 3 I 1
2 I 3 I 1 I 4 I 1
2 I 4 I 1 I 4 I 1
2 I 5 I 1 I 5 I 0
The mediate function does not work when package lmerTest is loaded, because lmerTest changes the class of the model from lmerMod to lmerModLmerTest.
The function returns the error: mediator model is not yet implemented
Can the function be upgraded to be compatible with lmerTest package?
Thanks
Hello,
thank you for making this great package!
I have a dataset with missing values so I am using the mediations() and amelidiate() functions. At the same time, I would also like to test whether the mediated effect differs by demographic so I would like to use the test.modmed() function. Can they be somehow used together?
Hello,
I am currently running a causal mediation analysis for a clinical trial where my mediator model uses a negative binomial regression model and my outcome model uses a logistic regression model.
model.x <- glm.nb(total.malaria.incident.events ~ Txarm + offset(years.in.study), data = delivery.analysis)
model.y = glm(compositeBO ~ Txarm + total.malaria.incident.events + offset(years.in.study), family = "binomial", data = delivery.analysis)
m.out = mediate(malaria.glm.nb, model.3.malaria.compositeBO, treat = "Txarm", mediator = "total.malaria.incident.events", boot=T)
However, I get the following error: Error in offset(years.in.study) : object 'years.in.study' not found
Does the mediate package allow for an "offset" function to be in the model because the time each participant contributes varies? If so, how can I correct this error? If not, is there another way to handle these models?
Thanks in advance!
I have a question regarding the functionality of the medsens function.
I am conducting a causal mediation analysis with a binary outcome and a continous mediator. My data is subject to missing data and I have opted to use multiple imputation in my workflow (through MICE). The mediations package workflow allows me to conduct the mediation analysis on imputed data in a very straightforward manner (through the mediations and amelidiate functions). The summary and plot functions work appropiately; I have obtained the results of the main analysis.
I am now in the stage where I would like to pass the mediate object (that results from the amelidiate function, which I use to pool results from the imputed datasets) through the medsens function in order to conduct the sensitivity analysis. Unfortunately, this is where I run into issues. I receive an error message which mentions that the model combination used in my analysis is not supported. The model combination should however be supported if I am not mistaken (continous mediator, binary outcome. Families in the mediations call are correctly specified as "gaussian" and "binomial", respectively; the mediate object confirms this. No treatment-mediator interaction was specified).
The question: should the medsens function in theory support the workflow described in the vignet accompanying the 'mediation' package on multiple-imputed data? I have no issues running the workflow on the non-imputed data, nor when I run the analysis on the individual imputed datasets (I have not gone through pooling the results from the individual analyses so far).
Looking forward as to any reply. Thanks in advance!
The demo for code here for gam
interactions no longer runs and gives error "Error in predict.gam(new.fit.Y, type = "response", newdata = pred.data.t) : 0 not in original fit". This appears to be an error that happens when predict.gam is given unknown factor levels.
Strangely if I comment out control
as below it runs or if I set "control=treat" it runs. However I get different results in with these two methods. The commenting out method gives a summary result that is closer to the old vignette:
model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp + marital
+ nonwhite + educ + income, data = jobs)
model.y <- gam(depress2 ~ treat + s(job_seek, by = treat) + s(job_seek, by = control)
+ depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income,
data = jobs)
out.5 <- mediation::mediate(model.m, model.y, sims = 1000, boot = TRUE,
treat = "treat", mediator = "job_seek",
#control = "control"
)
summary(out.5)
However the result isn't close enough for me to feel confident that this is correct (e.g. within .01).
Both the paper referenced in mediate
's documentation for gam
and the old vignette imply that the by
variables should be indicator variables. The current jobs
data does not have this form: treat
is binary but control
is a factor. If I create a new variable jobs$control.1 <- 1 - jobs$treat
I can then set the by
options sensibly, but again I don't get the values from the old vignette. I've also tried just using the spline by factor control
since the gam
documentation suggests that creates an interaction by itself anyways. Altogether I've explored a number of alternatives, and the closest to the old vignette seems to be the "comment out control" listing above.
This works fine with boot=F, but throws an error when boot=T. A real corner case, I admit:
library(mediation)
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0
test <- function(f1, f2, boot){
mtcars$mpg <- mtcars$mpg/max(mtcars$mpg)
m_med <- glm(f1, family=binomial, mtcars)
m_out <- glm(f2, family=binomial, mtcars)
mediate(m_med, m_out, treat = "wt", mediator = 'am', boot = boot)
}
x <- test('am ~ wt', 'mpg ~ wt + am', F)
#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
x <- test('am ~ wt', 'mpg ~ wt + am', T)
#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
#> Running nonparametric bootstrap
#> Error in stats::model.frame(formula = f1, data = structure(list(am = c(1, : object 'f1' not found
Created on 2021-03-17 by the reprex package (v0.3.0)
Hi,
Im trying to use the multimed function.
I can run the example code for multimed, but when I run it on my dataset, or this 'fake' dataset, it errors (see below).
treat is a 0/1 numeric variable.
What am I doing wrong? I assume its user error?
thanks!
library(palmerpenguins)
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
library(tidyverse)
library(reprex)
glimpse(penguins)
#> Rows: 344
#> Columns: 8
#> $ species Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel~
#> $ island Torgersen, Torgersen, Torgersen, Torgersen, Torgerse~
#> $ bill_length_mm 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, ~
#> $ bill_depth_mm 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, ~
#> $ flipper_length_mm 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186~
#> $ body_mass_g 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, ~
#> $ sex male, female, female, NA, female, male, female, male~
#> $ year 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007~
penguins %>%
filter(species!="Gentoo") %>%
mutate(treat = as.numeric(rbernoulli(n=220, 0.5))) %>%
drop_na()->df
summary(df)
#> species island bill_length_mm bill_depth_mm
#> Adelie :146 Biscoe : 44 Min. :32.1 Min. :15.50
#> Chinstrap: 68 Dream :123 1st Qu.:37.8 1st Qu.:17.50
#> Gentoo : 0 Torgersen: 47 Median :40.6 Median :18.40
#> Mean :42.0 Mean :18.37
#> 3rd Qu.:46.0 3rd Qu.:19.10
#> Max. :58.0 Max. :21.50
#> flipper_length_mm body_mass_g sex year treat
#> Min. :172.0 Min. :2700 female:107 Min. :2007 Min. :0.000
#> 1st Qu.:187.0 1st Qu.:3400 male :107 1st Qu.:2007 1st Qu.:0.000
#> Median :191.0 Median :3700 Median :2008 Median :1.000
#> Mean :191.9 Mean :3715 Mean :2008 Mean :0.514
#> 3rd Qu.:196.0 3rd Qu.:3994 3rd Qu.:2009 3rd Qu.:1.000
#> Max. :212.0 Max. :4800 Max. :2009 Max. :1.000
test<-multimed(outcome = "bill_length_mm",
med.main = "flipper_length_mm",
med.alt = "bill_depth_mm",
treat = "treat",
covariates = "year",
data=df, sims=100)
#> Warning in mean.default(data.1[, treat] * data.1[, med.main]^2): argument is not
#> numeric or logical: returning NA
#> Warning in ETM2 * sigma^2/VY: Recycling array of length 1 in vector-array arithmetic is deprecated.
#> Use c() or as.vector() instead.
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
Created on 2022-05-18 by the reprex package (v2.0.1)
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.