#' The Sobel mediation test
#'
#' To compute statistics and p-values for the Sobel test. Results for three versions of "Sobel test" are provided: Sobel test, Aroian test and Goodman test.
#'@param mv The mediator variable
#'@param iv The independent variable
#'@param dv The dependent variable
bda.mediation.test=function (mv, iv, dv)
{
if (any(is.na(mv)))
stop("Mediator contains missing value(s)")
if (any(is.na(iv)))
stop("Mediator contains missing value(s)")
if (any(is.na(dv)))
stop("Mediator contains missing value(s)")
nm = length(mv)
ni = length(iv)
nd = length(dv)
if (nm != ni | nm != nd | ni != nd)
stop("Variables have different lengths.")
tmp = summary(lm(mv ~ iv))
a = tmp$coef[2, 1]
sa = tmp$coef[2, 2]
tmp = summary(lm(dv ~ mv + iv))
b = tmp$coef[2, 1]
sb = tmp$coef[2, 2]
tmp1 = b^2 * sa^2 + a^2 * sb^2
tmp2 = sa^2 * sb^2
zsob = a * b/sqrt(tmp1)
psob = pnorm(-abs(zsob)) * 2
zaro = a * b/sqrt(tmp1 + tmp2)
paro = pnorm(-abs(zaro)) * 2
if (tmp1 > tmp2) {
zgm = a * b/sqrt(tmp1 - tmp2)
pgm = pnorm(-abs(zgm)) * 2
}
else {
zgm = NA
pgm = NA
}
p.value = c(psob, paro, pgm)
z.value = c(zsob, zaro, zgm)
out = data.frame(rbind(z.value, p.value))
names(out) = c("Sobel", "Aroian", "Goodman")
out
}
#' Perform mediation analysis by Baron and Kenny Method
#' @param X name of independent variable
#' @param M name of mediator variable
#' @param Y name of dependent variable
#' @param labels An optional list of variable names
#' @param data A data.frame
#' @param silent Logical. Whether or not show summary of regression tests
#' @param indirect.test Logical. Whether or not show results of bda.mediation.test
#' @param sig significant level. default value is 0.05
#' @importFrom stats pnorm qnorm
#' @export
#' @examples
#' labels=list(X="cond",M="pmi",Y="reaction")
#' result=mediationBK(labels=labels,data=pmi,silent=FALSE)
#' result
mediationBK=function(X=NULL,M=NULL,Y=NULL,labels=list(),data,silent=TRUE,indirect.test=TRUE,sig=0.05){
# X=NULL;M=NULL;Y=NULL;silent=FALSE;indirect.test=TRUE;data=pmi
if(is.null(X)) X=labels$X
if(is.null(M)) M=labels$M
if(is.null(Y)) Y=labels$Y
dataname=substitute(data)
Paths=c("Path C (Total Effect)","Path A (X on M)","Path B (M on Y, controlling for X)","Path C' (Direct Effect, X on Y, controlling for M)")
if(!silent) cat("Step 1: ",Paths[1],"\n-Estimate the relationship between X on Y\n")
temp1=paste0("lm(",Y,"~",X,",data=",dataname,")")
fit1=eval(parse(text=temp1))
if(!silent) print(summary(fit1))
if(!silent) cat("\nStep 2:",Paths[2],"\n-Estimate the relationship between X on M\n")
temp2=paste0("lm(",M,"~",X,",data=",dataname,")")
fit2=eval(parse(text=temp2))
if(!silent) print(summary(fit2))
a=fit2$coef[2]
sea=summary(fit2)$coef[2,2]
if(!silent) cat("\nStep 3:",Paths[3],"\n-Estimate the relationship between M on Y, controlling for X\n")
temp3=paste0("lm(",Y,"~",M,"+", X,",data=",dataname,")")
fit3=eval(parse(text=temp3))
if(!silent) print(summary(fit3))
b=fit3$coef[2]
seb=summary(fit3)$coef[2,2]
if(!silent) cat("\nStep 4:",Paths[4],"\n-Estimate the relationship between Y on X, controlling for M\n")
fit=list(fit1,fit2,fit3)
equations=list(temp1,temp2,temp3)
results=list()
coef=c()
pvalue=c()
for(i in 1:3){
coef=c(coef,summary(fit[[i]])$coef[2,1])
pvalue=c(pvalue,summary(fit[[i]])$coef[2,4])
}
coef=c(coef,summary(fit[[i]])$coef[3,1])
pvalue=c(pvalue,summary(fit[[i]])$coef[3,4])
results[[1]]=ifelse(pvalue[1]<0.05,"Acceptable","Not satisfied")
results[[2]]=ifelse(pvalue[2]<0.05,"Acceptable","Not satisfied")
results[[3]]=ifelse(pvalue[3]<0.05,"Acceptable","Not satisfied")
if(pvalue[4]>0.05) {
results[[4]]="Complete mediation"
} else if(coef[4]<coef[1]) {
results[[4]]="Partial mediation"
} else results[[4]]="Not satisfied"
# cat("\nResults of Baron and Kenny Method\n")
# name=c("c","a","b","c'")
# for(i in 1:4){
# temp=pvalue[i]
# if(temp<0.001) temp="< 0.001"
# else temp=paste0("= ",sprintf("%0.3f",temp))
# cat("Step",i,":", Paths[i],":",name[i],"=",sprintf("%0.3f",coef[i]),"( p",temp,")\n")
# }
# cat("Result :",results[[4]],"\n")
indirect=bda.mediation.test(data[[M]],data[[X]],data[[Y]])
# if(indirect.test){
# cat("\nResults of bda.mediation.test\n\n")
# print(indirect)
# }
seab=sqrt((a^2)*(seb^2)+(b^2)*(sea^2)+(sea^2)*(seb^2))
Z=(a*b)/seab
normalTheory=c(ab=a*b,seab=seab,z=Z,p=pnorm(Z,lower.tail = FALSE)*2,
ci.lower=a*b+qnorm(sig/2)*seab,ci.upper=a*b-qnorm(sig/2)*seab)
names(normalTheory)=c("ab","seab","z","p","ci.lower","ci.upper")
result=list(labels=labels,fit=fit,equations=equations,coef=coef,pvalue=pvalue,results=results,indirect=indirect,normalTheory=normalTheory)
class(result)="mediationBK"
invisible(result)
}
#' S3 method for class mediationBK
#' @param x An object of class mediationBK
#' @param ... Further arguments to be passed to print()
#' @export
print.mediationBK=function(x,...){
Paths=c("Path C (Total Effect)","Path A (X on M)",
"Path B (M on Y, controlling for X)","Path C' (Direct Effect, X on Y, controlling for M)")
cat("\nResults of Baron and Kenny Method\n")
name=c("c","a","b","c'")
for(i in 1:4){
temp=x$pvalue[i]
if(temp<0.001) temp="< 0.001"
else temp=paste0("= ",sprintf("%0.3f",temp))
cat("Step",i,":", Paths[i],":",name[i],"=",sprintf("%0.3f",x$coef[i]),"( p",temp,")\n")
}
cat("Result :",x$results[[4]],"\n")
cat("\nNormal theory test for indirect effect(s) :\n\n")
print(x$normalTheory)
cat("\nResults of bda.mediation.test\n\n")
print(x$indirect)
}
#' S3 method for class mediationBK
#' @param object An object of class mediationBK
#' @param ... Further arguments to be passed to summary()
#' @export
summary.mediationBK=function(object,...){
modelsSummary(labels=object$labels,object$fit)
}
#' S3 method for class mediationBK
#' @param x An object of class mediationBK
#' @param ... Further arguments to be passed to plot()
#' @importFrom rlang enexprs
#' @examples
#' labels=list(X="cond",M="pmi",Y="reaction")
#' result=mediationBK(labels=labels,data=pmi)
#' plot(result,type=1)
#' plot(result)
#' plot(result,type=1,whatLabel="label",arrowslabels="c",addprime=FALSE)
#' plot(result,whatLabel="label",arrowslabels=c("a","b","c"))
#' @export
plot.mediationBK=function(x,...){
vars=rlang::enexprs(...)
type=0
if(!is.null(vars$type)){
if(vars$type==1) type=1
}
whatLabel="est"
if(!is.null(vars$whatLabel)){
whatLabel=vars$whatLabel
}
if(type==1){
lty=ifelse(x$pvalue[1]<0.05,1,2)
if(whatLabel=="est"){
statisticalDiagram(0,labels=x$labels,arrowslabels=round(x$coef[1],3),arrowslty=lty,whatLabel = "label2",...)
} else{
statisticalDiagram(0,labels=x$labels,...)
}
} else{
lty=c()
for(i in 2:4){
lty=c(lty,ifelse(x$pvalue[i]<0.05,1,2))
}
if(whatLabel=="est"){
statisticalDiagram(4,labels=x$labels,arrowslabels=round(x$coef[2:4],3),arrowslty=lty,whatLabel = "label2",...)
} else{
statisticalDiagram(4,labels=x$labels,...)
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.