-
Notifications
You must be signed in to change notification settings - Fork 0
Regressions and mediations
By Etienne Aumont
- Packages
- Regression model
- Preparing for a mediation
- Single mediation
- Multiple mediation
library(lme4) #for linear models
library(effectsize) #to calculate effect sizes
library(lm.beta) #to calculate standardized beta
library(cocor) #for correlation coefficient comparisons
library(mediation) #for single mediations
library(mma) #for multiple mediations
Your model is defined using lm(dep_variable ~ indep_variable + indep_variable2(or covariates), data=your_dataframe) The order of the independent variables doesn't change anything
Model1 <- lm(AB_Jack~Braak1+ApoE4+Age+sex , data=data3)
m1<-summary(Model1)
Output for m1:
Call:
lm(formula = AB_Jack ~ Braak1 + ApoE4 + Age + sex, data = data3)
Residuals:
Min 1Q Median 3Q Max
-0.90695 -0.19775 -0.05647 0.13575 0.91952
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.161241 0.333614 -0.483 0.62959
Braak1 0.556425 0.041041 13.558 < 2e-16 ***
ApoE41 0.026035 0.057214 0.455 0.64975
Age 0.014212 0.004512 3.150 0.00198 **
sexM 0.051209 0.051036 1.003 0.31732
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3038 on 147 degrees of freedom
(21 observations effacées parce que manquantes)
Multiple R-squared: 0.6051, Adjusted R-squared: 0.5944
F-statistic: 56.32 on 4 and 147 DF, p-value: < 2.2e-16
The lines you should pay attention to are the ones starting with your independent variable names (here: Braak1).
You can consider that they have been corrected for all other independent variables included in the model.
You can perform interaction analyses by joining 2 or more variables with "*"
(lm(var1~var2*var3+var2+var3, data=data)
parameter confidence intervals can be obtained using confint()
Storing LM parameters
Parameters such as p-values can be stored in a matrix to be summarized later. Just use the coordinates of the coefficient you want. Here, we want the parameter on line #2 and column #4
pvalue_Amy<-c(m1$coefficients[2,4],m2$coefficients[2,4],m3$coefficients[2,4],m4$coefficients[2,4],m5$coefficients[2,4])
#For example, this is useful to perform FDR correction on your p-values when performing many analyses
pvalue_Amy_FDR <- p.adjust(pvalue_Amy, method = "fdr")
pvalue_Amy
output:
[1] 0.68134461 0.34078715 0.01992650 0.50145910 0.08192144
pvalue_Amy_FDR
output:
[1] 0.68134461 0.56797858 0.09963251 0.62682387 0.20480359
Calculate effect sizes
The basic package for regression effect sizes is this: lm.beta(Model1)
A <- effectsize(Model1) #This function is useful when the effect size you need is not just a beta, such as for ANOVA/ANCOVA
Output:
> A
# Standardization method: refit
Parameter | Std. Coef. | 95% CI
-----------------------------------------
(Intercept) | 0.19 | [-0.07, 0.44]
AB_Jack | -0.04 | [-0.23, 0.15]
ApoE41 | -0.61 | [-1.01, -0.21]
Age | -0.17 | [-0.36, 0.01]
sexM | 0.07 | [-0.30, 0.45]
To store the effect sizes like we did with the p-values (Here, we want line #2, column #2). Some trial and error might be needed.
beta_Amy<-c(A[2,2],B[2,2],C[2,2],D[2,2],E[2,2])
> beta_Amy
[1] -0.03904966 -0.09407310 -0.21374795 -0.06225126 -0.15826004
Compare correlations to one another
for this, you first need to remove all covariates. I recommend you obtain the data residualized for your covariates using something like this:
data_tPa$AB_Jack_resid2 <- residuals(lm(AB_Jack~ Age, data=data_tPa))
data_tPa$AB_Jack_resid2 <- residuals(lm(AB_Jack_resid2~ sex, data=data_tPa))
data_tPa$AB_Jack_resid2 <- residuals(lm(AB_Jack_resid2~ ApoE4, data=data_tPa))
If there are missing data, you may add ", na.action=na.exclude" at the end of each line. However, missing data has to be removed before running step 2
cocor.result <- cocor(~Braak2_DP_resid2 + AB_Jack_resid2 | Braak6_DP_resid2 + AB_Jack_resid2, data = data_tPa)
as.htest(cocor.result)
Here, you must define a predictor (the independent variable), a dependent variable (affected by the predictor) and a mediator (the variable that you believe might explain the link between the predictor and the dependent variable)
Step 1: remove all missing data
Data for mediation need to be devoid of missing data. For example, I removed here all data without RAVLT scores or both left and right hippocampal segmentations
data_RtR <- data_Rt[!is.na(data_Rt$RAVLT7),] #Removing all lines with NA in the RAVLT7 column
data_BtR <- data_RtR[!is.na(data_RtR$L_CA1),] #Removing all lines with missing left CA1 volume
Step 2: make sure that there is a main effect (independent variable effect on the dependent variable)
summary(lm(RAVLT7~Braak1+ApoE4+Age+edu+sex , data=data3))
If the Braak1 effect is significant, then you have an effect, and you can proceed
Step 3: make sure that the mediator makes sense
The mediator must be affected by the predictor to have any chance at mediating the effect. The mediator should also predict the dependent variable.
summary(lm(RAVLT7~L_CA1+Age+edu+sex , data=data_L2))
summary(lm(L_CA1~Braak1+ApoE4+Age+sex , data=data_L))
For this, your criteria can be more lenient than p < 0.05. p < 0.1, or even 0.2 is sufficient.
You will need to generate 2 linear models and plug them into the mediation model. Both LM need to be made using the same dataframe. LM #1 is the predictor predicting the mediator. LM #2 is the mediator and the predictor predicting the dependent variable
A<-lm(L_CA1~Braak1+ApoE4+Age+edu+sex , data=data_BtR)
A2<-lm(RAVLT7~L_CA1+Braak1+ApoE4+Age+edu+sex , data=data_BtR)
med.out <- mediate(A, A2, treat = "Braak1", mediator = "L_CA1")
summary(med.out)
The output will look like this:
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.0937 -0.4055 0.11 0.41
ADE -2.6932 -3.8018 -1.63 <2e-16 ***
Total Effect -2.7869 -3.8932 -1.71 <2e-16 ***
Prop. Mediated 0.0258 -0.0430 0.15 0.41
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 103
Simulations: 1000
Here, the total effect is the total of the direct effect (ADE) and the mediation effect (ACME) The ADE, is significant, which means that adding the mediator doesn't remove the effect of Braak1 MK on RAVLT but not the mediation effect (ACME). The total effect considers both the mediation and the direct effect
The prerequisites are the same as a single regression, but you include multiple mediators in the same model. I could only get it working by creating a new dataframe that includes the dependent variable, the predictor, the mediators and the covariates
new_data <- subset(data_BtR, select = c(RAVLT7, Braak2, L_DG.CA4_resid1,
R_DG.CA4_resid1, R_Sub_resid1, R_SRLM_resid1, Age, ApoE4, edu, sex))
x=new_data[,3:10] #This includes the covariates and mediators
y<-data.frame(new_data[,1]) #defining the dependent variable
pred <- new_data[, 2] #defining the predictor
cova <- data.frame(
Age = new_data$Age,
ApoE4 = new_data$ApoE4,
edu = new_data$edu,
sex = new_data$sex) #Defining the covariates
Then, you can plug this new dataframe in an MMA function
med_results <- MMA(
x = x,
y = y,
pred = pred,
mediator=1:4, #Based on X
jointm = list(n=1,j1=1:4), #list of mediators to include in a joint mediator (j1).
#Here, all 4 mediators are included, but you can select just 2 or 3.
cova = cova,
n2=500 #This is the number of iterations of the model. 500 iterations will require at least 15 minutes.
#You can start with 50 the first time to get an idea of the results.
)
summary(med_results, bymed = FALSE)
The output looks like this:
MMA Analysis: Estimated Mediation Effects Using GLM
For Predictor/Moderator at pred
$total.effect
est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
-6.839 -7.249 1.701 -3.916 -10.582 -4.403 -11.008 -4.133 -10.244 -3.059 -12.574 -4.474 -11.357 0.000 0.000
$direct.effect
est mean sd upbd lwbd upbd_q lwbd_q upbd_bcbi lwbd_bcbi upbd_b lwbd_b upbd_win lwbd_win p_norm p_quan
-6.625 -7.031 1.745 -3.611 -10.450 -4.165 -11.018 -3.700 -10.201 -2.797 -12.228 -4.202 -11.393 0.000 0.000
$indirect.effect
y1.all y1.L_DG.CA4_resid1 y1.R_DG.CA4_resid1 y1.R_Sub_resid1 y1.R_SRLM_resid1 y1.j1
est -0.207 -0.586 0.533 0.291 -0.444 -0.207
mean -0.219 -0.591 0.465 0.260 -0.353 -0.219
sd 0.613 0.559 0.745 0.403 0.536 0.613
upbd 0.982 0.503 1.926 1.049 0.698 0.982
lwbd -1.420 -1.686 -0.995 -0.530 -1.403 -1.420
upbd_q 0.872 0.172 2.075 1.258 0.500 0.872
lwbd_q -1.498 -2.008 -0.804 -0.367 -1.645 -1.498
upbd_bcbi 0.809 0.054 3.321 1.602 0.156 0.809
lwbd_bcbi -1.770 -2.501 -0.283 -0.118 -2.333 -1.770
upbd_b 1.445 0.383 2.226 1.405 0.888 1.445
lwbd_b -2.165 -2.381 -1.720 -0.946 -1.780 -2.165
upbd_win 0.845 0.120 2.288 1.332 0.445 0.845
lwbd_win -1.990 -2.747 -0.408 -0.147 -2.137 -1.990
p_norm 0.721 0.290 0.532 0.519 0.511 0.721
p_quan 0.752 0.180 0.500 0.460 0.496 0.752
The output is more complex, but the idea is exactly the same: ADE is the direct.effect and ACME is the indirect.effect Each mediator has a column and y1.all is the sum of the effect of all mediators together. j1 is the sum of the mediators you've included in the jointm list. The "est" is equivalent to the estimate in the single mediation. I recommend to use the p_norm instead of the p_quan for more details and ideas for personalization, see this guide