Exercise 6 instructions.

Download the data ‘PSE_MOL_Doors.sav’, these are the data from an experiment by Kamila Irvine and Piers Cornelissen. This file contains data on 95 women performing various scales and body image-related tasks. doors_front is the score from a gap estimation task, w_dn is the actual gap a participant can pass through. The (estimated) Point of subjective equality or PSE (the BMI they believe themselves to be) when viewing an imageset varying in BMI. Participants used the method of adjustment to estimate their body size with the same stimulus set as for the yes-no task (MOL). BMI is the participant’s actual BMI.

Test the mediation model: doors_front –> PSE –> BMI via the causal steps method by Baron & Kenny. Report as you would do in a paper.

Make a diagram (use ‘mediate’).

Calculate a Sobel z test and report.

Test the mediation via Preacher & Hayes method.

Now test a mediation model with 2 mediators (PSE and MOL) but with the same independent and dependent variables.

Export a figure for that mediation model.

Test the mediation via Imai et al.’s method.

BONUS: perform the sensitivity analysis via Imai et al.’s method.

Load the data and describe.

Load the data and get some descriptives.

setwd("~/Dropbox/Teaching_MRes_Northumbria/Lecture6")
require(haven)
## Loading required package: haven
Data<-read_spss("PSE_MOL_DOORS.sav")
require(skimr)
## Loading required package: skimr
skim(Data)
Data summary
Name Data
Number of rows 98
Number of columns 23
_______________________
Column type frequency:
character 2
numeric 21
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
widest 0 1 0 8 2 3 0
an_history 0 1 1 1 0 3 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
subject 0 1.00 49.50 28.43 1.00 25.25 49.50 73.75 98.00 ▇▇▇▇▇
age 2 0.98 22.44 4.69 18.00 19.00 20.50 24.00 35.00 ▇▂▁▁▁
height 2 0.98 1.65 0.06 1.46 1.61 1.65 1.69 1.78 ▁▁▇▆▂
weight 2 0.98 62.58 11.11 44.20 55.88 60.00 68.45 100.10 ▅▇▃▁▁
BMI 2 0.98 22.98 4.15 15.83 20.37 22.32 24.56 38.62 ▅▇▂▁▁
PSE 2 0.98 23.77 4.43 16.02 20.58 23.42 26.19 41.50 ▆▇▃▁▁
MOL 2 0.98 24.34 5.20 17.03 20.32 23.22 25.49 41.86 ▆▇▂▁▁
w_dn 2 0.98 45.83 3.52 39.20 43.45 45.55 47.50 59.80 ▅▇▃▁▁
w_ax 2 0.98 40.00 2.89 34.00 38.27 40.00 41.00 57.20 ▅▇▁▁▁
depth 2 0.98 27.19 4.31 21.60 24.67 26.20 28.25 48.40 ▇▃▁▁▁
RSE 2 0.98 18.26 4.29 8.00 16.00 18.50 21.00 28.00 ▂▅▇▅▁
BSQ 2 0.98 44.43 17.05 19.00 30.75 41.00 56.50 87.00 ▇▇▆▃▂
EDEQ_res 2 0.98 1.43 1.31 0.00 0.35 1.10 2.40 5.00 ▇▃▂▂▁
EDEQ_eat 2 0.98 0.93 0.93 0.00 0.20 0.80 1.40 4.20 ▇▅▁▁▁
EDEQ_sc 2 0.98 2.56 1.61 0.00 1.25 2.44 3.87 8.00 ▇▇▆▂▁
EDEQ_wc 2 0.98 1.99 1.49 0.00 0.75 1.80 2.80 5.40 ▇▆▅▃▂
EDEQ_global 2 0.98 1.73 1.20 0.06 0.78 1.41 2.74 4.58 ▇▇▃▅▂
BDI 2 0.98 12.03 9.02 0.00 5.75 10.00 16.00 47.00 ▇▆▃▁▁
doors_front 3 0.97 52.77 8.55 31.55 47.53 52.65 56.53 87.60 ▁▇▆▁▁
doors_side 3 0.97 36.48 7.48 24.30 31.55 35.10 39.65 67.70 ▆▇▂▁▁
doors_ball 3 0.97 69.02 8.78 34.80 64.05 67.90 73.50 96.10 ▁▁▇▃▁

Baron & Kenny.

Bear in mind that here I have not checked the assumptions.

Step 1.

There is support for Step 1.

step1<-lm(BMI~doors_front, data=Data)
summary(step1)
## 
## Call:
## lm(formula = BMI ~ doors_front, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0480 -2.3237 -0.7612  1.2384 14.4915 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  14.7787     2.5658   5.760 1.08e-07 ***
## doors_front   0.1555     0.0480   3.239  0.00167 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.979 on 93 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.1014, Adjusted R-squared:  0.0917 
## F-statistic: 10.49 on 1 and 93 DF,  p-value: 0.001665

Step 2.

We also find support for step 2.

step2<-lm(PSE~doors_front, data= Data)
summary(step2)
## 
## Call:
## lm(formula = PSE ~ doors_front, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4181 -2.5253 -0.0993  2.2904 16.0526 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.59476    2.66410   5.103 1.76e-06 ***
## doors_front  0.19367    0.04984   3.886 0.000191 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.131 on 93 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.1397, Adjusted R-squared:  0.1304 
## F-statistic:  15.1 on 1 and 93 DF,  p-value: 0.0001909

Step 3.

Step 3 shows that the effect of doors_front was reduced from .156 to .020 and no longer significant.

step3<-lm(BMI~PSE+doors_front, data= Data)
summary(step3)
## 
## Call:
## lm(formula = BMI ~ PSE + doors_front, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4484 -2.0667 -0.0722  1.2345  7.8752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.23679    1.99847   2.620   0.0103 *  
## PSE          0.70188    0.06875  10.209   <2e-16 ***
## doors_front  0.01952    0.03563   0.548   0.5850    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.739 on 92 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.5787, Adjusted R-squared:  0.5695 
## F-statistic: 63.17 on 2 and 92 DF,  p-value: < 2.2e-16

Step 4

We want this so we can make a diagram.

step4<-lm(BMI~PSE, data= Data)
summary(step4)
## 
## Call:
## lm(formula = BMI ~ PSE, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6002 -2.0903 -0.0908  1.3892  7.6851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.13169    1.52850   4.012 0.000121 ***
## PSE          0.70885    0.06323  11.210  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.731 on 94 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.5721, Adjusted R-squared:  0.5675 
## F-statistic: 125.7 on 1 and 94 DF,  p-value: < 2.2e-16

Make a diagram

Here is an example. That’s admittedly not great.

require(DiagrammeR)
## Loading required package: DiagrammeR
mermaid(" graph LR
            D(Doors)-->|.19***|P(PSE)
            P-->|.71***|B(BMI)
          D-->|.16** / .02|B
            ")

This is perhaps better. Notice that these are B’s.

library(diagram)
## Loading required package: shape
data <- c(0, "'.19***'", 0,
          0, 0, 0, 
          "'.71***'", "'.16** (.02)'", 0)
M<- matrix (nrow=3, ncol=3, byrow = TRUE, data=data)
plot<- plotmat (M, pos=c(1,2), 
                name= c( "PSE","Doors", "BMI"), 
                box.type = "rect", box.size = 0.12, box.prop=0.5,  curve=0)

Note that ‘mediate’ will standardize coefficients, which is why they are different from the above graph. (You can replace the values from the above as it is a better figure)

require(psych)
## Loading required package: psych
mediationmodel1<-mediate(BMI~doors_front+(PSE), std=TRUE, data=Data, n.iter=10000, plot=F)
mediate.diagram(mediationmodel1)

Summary: A mediation analysis via Baron & Kenny’s (1986) method showed that there was mediation (Figure X).

Sobel z test.

The mediator contains missing values! Therefore we remove those… .

# create reduced dataframe. Remove missings.
Data_red<-na.omit(cbind.data.frame(Data$PSE,Data$doors_front,Data$BMI))
# rename the columns.
colnames(Data_red)<-c("PSE","doors","BMI")
require(bda)
## Loading required package: bda
## Loading required package: boot
## 
## Attaching package: 'boot'
## The following object is masked from 'package:psych':
## 
##     logit
## bda v17 (Bin Wang, 2023)
mediation.test(Data_red$PSE,Data_red$doors, Data_red$BMI)
##                Sobel       Aroian      Goodman
## z.value 3.6317053592 3.6165813558 3.6470207025
## p.value 0.0002815545 0.0002985195 0.0002652986

A Sobel z test showed that the mediation effect reported in Fig. X was significant (Sobel z= 3.63, p=.0003).

Preacher & Hayes method.

Standardized indirect effects were computed for 10,000 bootstrapped samples, and the 95% percentile confidence interval was examined via the ‘psych’ package. The bootstrapped standardized indirect effect was 0.27, and the 95% confidence interval ranged from 0.12 to 0.41. Thus, the indirect effect was statistically significant and mediation is supported.

require(psych)
set.seed(12345)
mediate(BMI~doors_front+(PSE), std=T, data= Data, n.iter=10000)
## 
## Mediation/Moderation Analysis 
## Call: mediate(y = BMI ~ doors_front + (PSE), data = Data, n.iter = 10000, 
##     std = T)
## 
## The DV (Y) was  BMI . The IV (X) was  doors_front . The mediating variable(s) =  PSE .
## 
## Total effect(c) of  doors_front  on  BMI  =  0.32   S.E. =  0.1  t  =  3.31  df=  96   with p =  0.0013
## Direct effect (c') of  doors_front  on  BMI  removing  PSE  =  0.04   S.E. =  0.07  t  =  0.6  df=  95   with p =  0.55
## Indirect effect (ab) of  doors_front  on  BMI  through  PSE   =  0.28 
## Mean bootstrapped indirect effect =  0.27  with standard error =  0.07  Lower CI =  0.12    Upper CI =  0.41
## R = 0.76 R2 = 0.57   F = 63.93 on 2 and 95 DF   p-value:  1.04e-22 
## 
##  To see the longer output, specify short = FALSE in the print statement or ask for the summary
# below calls a longer output
set.seed(12345)
summary(mediate(BMI~doors_front+(PSE), std=T, data= Data, n.iter=10000), short=F)

## Call: mediate(y = BMI ~ doors_front + (PSE), data = Data, n.iter = 10000, 
##     std = T)
## 
## Direct effect estimates (traditional regression)    (c') X + M on Y 
##              BMI   se     t df     Prob
## Intercept   0.00 0.07  0.00 95 1.00e+00
## doors_front 0.04 0.07  0.60 95 5.48e-01
## PSE         0.74 0.07 10.25 95 4.87e-17
## 
## R = 0.76 R2 = 0.57   F = 63.93 on 2 and 95 DF   p-value:  2.57e-18 
## 
##  Total effect estimates (c) (X on Y) 
##              BMI  se    t df    Prob
## Intercept   0.00 0.1 0.00 96 1.00000
## doors_front 0.32 0.1 3.31 96 0.00132
## 
##  'a'  effect estimates (X on M) 
##              PSE   se    t df     Prob
## Intercept   0.00 0.09 0.00 96 1.000000
## doors_front 0.37 0.09 3.95 96 0.000151
## 
##  'b'  effect estimates (M on Y controlling for X) 
##      BMI   se     t df     Prob
## PSE 0.74 0.07 10.25 95 4.87e-17
## 
##  'ab'  effect estimates (through all  mediators)
##              BMI boot   sd lower upper
## doors_front 0.28 0.27 0.07  0.12  0.41

Both PSE and MOL.

Note that this averages the effect of the mediators.

require(psych)
dual_path<-mediate(BMI~doors_front+(PSE)+(MOL), std=T, data= Data, n.iter=10000)

dual_path
## 
## Mediation/Moderation Analysis 
## Call: mediate(y = BMI ~ doors_front + (PSE) + (MOL), data = Data, n.iter = 10000, 
##     std = T)
## 
## The DV (Y) was  BMI . The IV (X) was  doors_front . The mediating variable(s) =  PSE MOL .
## 
## Total effect(c) of  doors_front  on  BMI  =  0.32   S.E. =  0.1  t  =  3.31  df=  96   with p =  0.0013
## Direct effect (c') of  doors_front  on  BMI  removing  PSE MOL  =  0.01   S.E. =  0.07  t  =  0.11  df=  94   with p =  0.91
## Indirect effect (ab) of  doors_front  on  BMI  through  PSE MOL   =  0.31 
## Mean bootstrapped indirect effect =  0.31  with standard error =  0.09  Lower CI =  0.13    Upper CI =  0.47
## R = 0.8 R2 = 0.64   F = 55.7 on 3 and 94 DF   p-value:  5.4e-24 
## 
##  To see the longer output, specify short = FALSE in the print statement or ask for the summary
setEPS()
postscript("dual_path.eps", horizontal = FALSE, onefile = FALSE, paper = "special")
par(mar=c(1,1,1,1))
mediate.diagram(dual_path)
dev.off()
## quartz_off_screen 
##                 2

Mediate. (+ Bonus: sensitivity analysis)

The ‘mediate’ function cannot handle multiple mediators, we therefore test them sequentially. ‘multimed’ can handle multiple mediators but cannot take continuous variables as ‘treat’.

Important: below will take time.

require(mediation)
## Loading required package: 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
## 
## Attaching package: 'mediation'
## The following object is masked from 'package:psych':
## 
##     mediate
med.fit<-lm(PSE~doors_front, data=Data)
out.fit<-lm(BMI~doors_front+PSE, data=Data)
# Robust SE is ignored for Bootstrap. Otherwise omit boot=TRUE.
set.seed(1984)
med.out <- mediate(med.fit, out.fit, treat = "doors_front", mediator = "PSE", boot=TRUE, sims = 10000)
## Running nonparametric bootstrap
summary(med.out)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME             0.1359       0.0570         0.21  0.0010 ***
## ADE              0.0195      -0.0629         0.10  0.6400    
## Total Effect     0.1555       0.0590         0.25  0.0028 ** 
## Prop. Mediated   0.8744       0.4641         1.75  0.0034 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 95 
## 
## 
## Simulations: 10000
plot(med.out) 

sensitivity_analysis<-medsens(med.out) 
summary(sensitivity_analysis)
## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##      Rho   ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] 0.7 0.0107      -0.0156       0.0371         0.49       0.1776
## 
## Rho at which ACME = 0: 0.7
## R^2_M*R^2_Y* at which ACME = 0: 0.49
## R^2_M~R^2_Y~ at which ACME = 0: 0.1776

Using a 10,000 bootstraps, the mediation analysis via Imai, Keele & Tingley’s method (2011) showed a significant average causal mediation effect (ACME): 0.14, 95%CI [0.06, 0.21]. However, the average direct effect was no longer significant (ADE):.02, 95%CI [-0.06, 0.10].

Now we do the same via MOL.

require(mediation)
med.fit2<-lm(MOL~doors_front, data=Data)
out.fit2<-lm(BMI~doors_front+MOL, data=Data)
# Robust SE is ignored for Bootstrap. Otherwise omit boot=TRUE.
set.seed(1984)
med.out2 <- mediate(med.fit2, out.fit2, treat = "doors_front", mediator = "MOL", boot=TRUE, sims = 10000)
## Running nonparametric bootstrap

So, we find support for mediation via both paths. The ‘hidden confounder’ has to be stronger for MOL than for PSE in order to explain away the mediation effect.

summary(med.out2)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            0.14665      0.05486         0.24  0.0032 **
## ADE             0.00881     -0.06330         0.08  0.7756   
## Total Effect    0.15546      0.05900         0.25  0.0028 **
## Prop. Mediated  0.94334      0.52073         1.64  0.0052 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 95 
## 
## 
## Simulations: 10000
plot(med.out2) 

sensitivity_analysis2<-medsens(med.out2) 
summary(sensitivity_analysis2)
## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##      Rho    ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] 0.7  0.0246      -0.0033       0.0525         0.49       0.1573
## [2,] 0.8 -0.0194      -0.0463       0.0075         0.64       0.2054
## 
## Rho at which ACME = 0: 0.8
## R^2_M*R^2_Y* at which ACME = 0: 0.64
## R^2_M~R^2_Y~ at which ACME = 0: 0.2054

MBESS

If mediation does not work, here is an alternative to achieve roughly the same thing.

MBESS needs complete data.

Data_no_miss<-na.omit(Data)

Below does the calculation with PSE as a mediator.

Note the warning points to the fact that there are ‘extreme order statistics as endpoints’. This points to potential issues with outliers in the bootstrap.

require(MBESS)
## Loading required package: MBESS
## 
## Attaching package: 'MBESS'
## The following object is masked from 'package:psych':
## 
##     cor2cov
set.seed(1234)
MBESS_1<-MBESS::mediation(x = Data_no_miss$doors_front, mediator = Data_no_miss$PSE, dv = Data_no_miss$BMI, bootstrap = TRUE, which.boot = "BCa", B = 10000)
## [1] "Bootstrap resampling has begun. This process may take a considerable amount of time if the number of replications is large, which is optimal for the bootstrap procedure."
## Warning in norm.inter(t, adj.alpha): extreme order statistics used as endpoints
MBESS_1
##                                          Estimate CI.Lower_BCa CI.Upper_BCa
## Indirect.Effect                        0.13593392   0.06784993    0.2237093
## Indirect.Effect.Partially.Standardized 0.03256002   0.01510387    0.0471279
## Index.of.Mediation                     0.27838375   0.13045123    0.4246567
## R2_4.5                                 0.09998371   0.00985396    0.2372505
## R2_4.6                                 0.07419206   0.01622114    0.1564360
## R2_4.7                                 0.12821577   0.02967435    0.2443907
## Ratio.of.Indirect.to.Total.Effect      0.87440417   0.48783460    1.8774763
## Ratio.of.Indirect.to.Direct.Effect     6.96204759   1.85296907 6242.2252560
## Success.of.Surrogate.Endpoint          0.80269625   0.40308496    1.3353193
## Residual.Based_Gamma                   0.08207922   0.01580217    0.1829540
## Residual.Based.Standardized_gamma      0.08172430   0.01561751    0.1817403
## SOS                                    0.98642918   0.86811368    0.9999995

Using a 10,000 bootstraps via Bias-Corrected Acceleration, the MBESS mediation method (Preacher & Kelley 2011) showed a significant indirect effect: 0.136, 95%CI [0.07, 0.22].

So, the indirect effect estimate (‘ACME in mediate’) is quite similar to the above (the confidence intervals differ somewhat due to differen methods of estimation). You can derive the total effect as follows:

# Indirect.Effect / Ratio.of.Indirect.to.Total.Effect
0.13593392 / 0.87440417 
## [1] 0.1554589

It does not print the average direct effect, like the Imai et al. method, when bootstrapping but you can still get an estimate (it is c.prime = .019). Note to the similarity to the estimate in the causal steps approach. These are unstandardised effects (hence the difference with the output from ‘psych’, for example). The 95% CI is not bootstrapped but is a normal approximation.

MBESS_2<-MBESS::mediation(x = Data_no_miss$doors_front, mediator = Data_no_miss$PSE, dv = Data_no_miss$BMI)
MBESS_2
## $Y.on.X
## $Y.on.X$Regression.Table
##                 Estimate Std. Error  t value      p(>|t|) Low Conf Limit
## Intercept.Y_X 14.7786954 2.56575643 5.759976 1.080008e-07     9.68361161
## c (Regressor)  0.1554589 0.04799937 3.238770 1.665163e-03     0.06014168
##               Up Conf Limit
## Intercept.Y_X    19.8737792
## c (Regressor)     0.2507762
## 
## $Y.on.X$Model.Fit
##        Residual standard error (RMSE) numerator df denomenator df F-Statistic
## Values                       3.978861            1             93    10.48963
##        p-value (F)       R^2    Adj R^2 Low Conf Limit Up Conf Limit
## Values 0.001665163 0.1013592 0.09169643     0.01538183     0.2367254
## 
## 
## $M.on.X
## $M.on.X$Regression.Table
##                 Estimate Std. Error  t value      p(>|t|) Low Conf Limit
## Intercept.M_X 13.5947632 2.66409947 5.102949 1.762926e-06     8.30438967
## a (Regressor)  0.1936709 0.04983914 3.885920 1.908838e-04     0.09470025
##               Up Conf Limit
## Intercept.M_X    18.8851368
## a (Regressor)     0.2926416
## 
## $M.on.X$Model.Fit
##        Residual standard error (RMSE) numerator df denomenator df F-Statistic
## Values                       4.131367            1             93    15.10037
##         p-value (F)       R^2   Adj R^2 Low Conf Limit Up Conf Limit
## Values 0.0001908838 0.1396884 0.1304378     0.03427217     0.2839665
## 
## 
## $Y.on.X.and.M
## $Y.on.X.and.M$Regression.Table
##                       Estimate Std. Error    t value      p(>|t|)
## Intercept.Y_XM      5.23679013 1.99846509  2.6204061 1.027445e-02
## c.prime (Regressor) 0.01952499 0.03562736  0.5480336 5.849961e-01
## b (Mediator)        0.70188094 0.06875424 10.2085471 8.334789e-17
##                     Low Conf Limit Up Conf Limit
## Intercept.Y_XM          1.26766595    9.20591431
## c.prime (Regressor)    -0.05123402    0.09028401
## b (Mediator)            0.56532908    0.83843281
## 
## $Y.on.X.and.M$Model.Fit
##        Residual standard error (RMSE) numerator df denomenator df F-Statistic
## Values                       2.739269            2             92     63.1729
##        p-value (F)     R^2   Adj R^2 Low Conf Limit Up Conf Limit
## Values           0 0.57865 0.5694902       0.427292     0.6903835
## 
## 
## $Effect.Sizes
##                                              [,1]
## Indirect.Effect                        0.13593392
## Indirect.Effect.Partially.Standardized 0.03256002
## Index.of.Mediation                     0.27838375
## R2_4.5                                 0.09998371
## R2_4.6                                 0.07419206
## R2_4.7                                 0.12821577
## Ratio.of.Indirect.to.Total.Effect      0.87440417
## Ratio.of.Indirect.to.Direct.Effect     6.96204759
## Success.of.Surrogate.Endpoint          0.80269625
## Residual.Based_Gamma                   0.08207922
## Residual.Based.Standardized_gamma      0.08172430
## SOS                                    0.98642918

Robust Mediation.

Here is yet another way to do a mediation, this time via the ‘robmed’ package. This method relaxes assumptions about normality which commonly apply to regression methods of mediation (Alfons, Ates, and Groenen, 2018). This package also allows controlling for a variable via specifying covariate =. Check the manual.

Using a 10,000 bootstraps via Bias-Corrected Acceleration, the robust mediation method (2011) showed a significant indirect effect: 0.11, 95%CI [0.04, 0.22]. The output also reports the bootstrapped direct effect estimate (akin to c’ or ADE) but not the 95% CI.

require(robmed)
## Loading required package: robmed
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## Loading required package: robustbase
## 
## Attaching package: 'robustbase'
## The following object is masked from 'package:boot':
## 
##     salinity
set.seed(1234)
robmed <- test_mediation(Data,
                       x = "doors_front",
                       y = "BMI",
                       m = "PSE", R=10000, type='bca')
summary(robmed)
## Robust bootstrap test for indirect effect via regression
## 
## x = doors_front
## y = BMI
## m = PSE
## 
## Sample size: 95
## ---
## Outcome variable: PSE
## 
## Coefficients:
##                 Data     Boot Std. Error z value Pr(>|z|)    
## (Intercept) 12.55650 12.70577    2.10689   6.031 1.63e-09 ***
## doors_front  0.20790  0.20478    0.04131   4.957 7.15e-07 ***
## 
## Robust residual standard error: 3.649 on 93 degrees of freedom
## Robust R-squared:  0.1976,   Adjusted robust R-squared:  0.1889
## Robust F-statistic: 16.27 on 1 and Inf DF,  p-value: 5.504e-05
## 
## Robustness weights:
## Observation 71 is a potential outlier with weight 0
## ---
## Outcome variable: BMI
## 
## Coefficients:
##                Data    Boot Std. Error z value Pr(>|z|)   
## (Intercept) 7.55362 7.32943    4.25312   1.723  0.08483 . 
## PSE         0.53157 0.53624    0.17782   3.016  0.00256 **
## doors_front 0.04143 0.04396    0.04453   0.987  0.32364   
## 
## Robust residual standard error: 2.382 on 92 degrees of freedom
## Robust R-squared:  0.4743,   Adjusted robust R-squared:  0.4629
## Robust F-statistic: 11.45 on 2 and Inf DF,  p-value: 1.066e-05
## 
## Robustness weights:
## Observation 68 is a potential outlier with weight 0
## ---
## Total effect of x on y:
##                Data    Boot Std. Error z value Pr(>|z|)   
## doors_front 0.15194 0.15440    0.05771   2.675  0.00746 **
## 
## Direct effect of x on y:
##                Data    Boot Std. Error z value Pr(>|z|)
## doors_front 0.04143 0.04396    0.04453   0.987    0.324
## 
## Indirect effect of x on y:
##       Data   Boot   Lower  Upper
## PSE 0.1105 0.1104 0.03861 0.2202
## ---
## Level of confidence: 95 %
## 
## Number of bootstrap replicates: 10000
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Do check the manual for further information.

A word of caution: Causal ordering!

This was a toy example. As you might have spotted the logical chain might not make sense at all: we are predicting participants’ actual BMI as an outcome measure. However, it might be more sensible that BMI is a predictor of people’s perceptions of whether they pass through a gap (rather than the other way round!). Have a look at the paper. It just goes to show, that we can build all sorts of models with these mediation packages and get “significant” results, but that they might make no theoretical sense!

THE END!