Skip to content

Instantly share code, notes, and snippets.

@peterdalle

peterdalle/medyad.r

Created February 4, 2020 22:31
Show Gist options
  • Save peterdalle/3f305194c2482351afd309e6deb25f76 to your computer and use it in GitHub Desktop.
Save peterdalle/3f305194c2482351afd309e6deb25f76 to your computer and use it in GitHub Desktop.
Download ZIP
MEDYAD - Easy statistical mediation analysis with distinguishable dyadic data
Raw
medyad.r
# Downloaded from http://www.jjcoutts.com/medyad
# Coutts, J., Hayes, A. F., & Jiang, T. (2019). Easy statistical mediation analysis with distinguishable dyadic data. Journal of Communication
# https://doi.org/10.1093/joc/jqz034
# -----------------------------------------------------------------
#MEDYAD for SAS beta version 1.1;
#current as of 18 December 2019;
#Copyright 2019 by Jacob J. Coutts and Andrew F. Hayes;
#Distribution of this code in any form, except through
#jjcouts.com or afhayes.com is prohibited;
#without the written permission of the copyright holder;
#THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND;
#EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF;
#MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT;
#IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, ;
#DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT;
#OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE ;
#SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE ;
#USE OF THIS SOFTWARE IMPLIES AGREEMENT WITH THESE TERMS ;
medyad <- function(data,y="xxxxx",x="xxxxx",m="xxxxx",cov="xxxxx",cmatrix=-999,mb=0,conf=95,hc=5,decimals=4,boot=5000,seed=0,contrast=0,save=0,bootmax=10000,describe=1,total=1)
{
widthex <- getOption("width")
options(width=160)
contrast<-floor(contrast)
covYES<-0;if(seed!=0){seed<-abs(floor(seed));set.seed(seed)}
errs<-0;notes<-1;criterr<-0;errcode<-matrix(0,100,1);notecode<-matrix(0,100,1);
adjust<-0;bad<-0;conseq<-" ";antec<-" ";goodboot<-0;badboot<-0;go<-0;total<-as.numeric(total==1);
covNUM<-0;yesxym<-0
hc<-trunc(hc);if((hc>5) | (hc<0)){hc<-5}
if((hc>=0) & (hc<5)){notecode[notes,1]<-2;notes<-notes+1}
describe<-as.numeric(describe==1); saveb<-(save==1);
specify_decimal <- function(x, k)
{trimws(format(round(x, k), nsmall=k))}
bootconf <- function(boottmp)
{
sbootind <- as.matrix(boottmp)
bootvar <- var(boottmp);
if(ncol(sbootind)>1)
{bootot <- colMeans(boottmp);seboot <- sqrt(diag(bootvar))}
if(ncol(sbootind)==1){bootot<-mean(boottmp);seboot<-sd(boottmp)}
temp<-sbootind;
for(bi in 1:ncol(sbootind))
{
temp[rank(sbootind[,bi]),bi] <- sbootind[,bi]
sbootind[,bi] <- temp[,bi]
llci <- as.matrix(sbootind[cilow,])
ulci <- as.matrix(sbootind[cihigh,])
}
bootmat<-cbind(bootot,seboot,llci,ulci);
return(bootmat)
}
hcestmd <- function(x,resid,hc,mse)
{
n1<-nrow(x);
invXtX <- solve(t(x)%*%x);varb<-as.numeric(mse)*invXtX;k3 <- ncol(x);xhc <- 0
if (hc!=5)
{
xhc <- x;hat <- xhc[,1];
for(i3 in 1:nrow(xhc))
{hat[i3]<-t(xhc[i3,])%*%invXtX%*%xhc[i3,]};
if((hc==0) | (hc==1))
{for(i3 in 1:k3){xhc[,i3] <- xhc[,i3]*resid}}
if ((hc==3) | (hc==2))
{
for(i3 in 1:k3)
{xhc[,i3] <- (resid/(1-hat)^(1/(4-hc)))*xhc[,i3];}
}
if (hc==4)
{
hcmn <- matrix(4,n,2);hcmn[,2] <- (n1*hat)/k3;
rminhcmn <- apply(hcmn,1,min)
for(i3 in 1:k3)
{xhc[,i3] <- (resid/(1-hat)^(rminhcmn/2))*xhc[,i3]}
}
varb<-invXtX%*%t(xhc)%*%xhc%*%invXtX;
if (hc==1)
{varb <- (n1/(n1-ncol(x)))*varb}
}
hclab <- c("se(HC0)","se(HC1)","se(HC2)","se(HC3)","se(HC4)","se")
hclab <- hclab[(hc+1)]
hcflab <- c("F(HC0)","F(HC1)","F(HC2)","F(HC3)","F(HC4)","F")
hcflab <- hcflab[(hc+1)]
hclist <- list("varb"=varb,"hclab"=hclab,"hcflab"=hcflab)
return(hclist)
}
modest <- function(y,x,type,full=0,contest=999)
{
if (type==1)
{
b <- solve(t(x)%*%x)%*%t(x)%*%y;
modres <- b
if (full>0)
{
n1 <- nrow(x);dfres <- n1-(ncol(x))
sstotal <- t(y-(sum(y)/n1))%*%(y-(sum(y)/n1))
resid <- y-x%*%b
ssresid <- t(resid)%*%resid
r2 <- (sstotal-ssresid)/sstotal
adjr2 <- 1-((1-r2)*(n1-1)/(dfres))
mse <- ssresid/(n1-ncol(x));
varb2 <- hcestmd(x=x,resid=resid,hc=hc,mse=mse);
seb <- sqrt(diag(varb2$varb))
trat <- b[,1]/seb
p <- 2*(1-pt(abs(trat),(dfres)))
tval <- sqrt(dfres*(exp((dfres-(5/6))*((xp2/(dfres-(2/3)+(.11/dfres)))*(xp2/(dfres-(2/3)+(.11/dfres)))))-1))
if(length(contest)>1)
{
contestv <- contest*b
secont <- sqrt(contest%*%varb%*%t(contest))
print(contestv);print(secont)
}
if (full>1)
{
modres <- cbind(modres,seb,trat,p)
modres <- cbind(modres,(b-tval*seb),(b+tval*seb))
modresl <- c("coeff",varb2$hclab,"t","p","LLCI","ULCI")
lmat <- diag(ncol(x))
lmat <- lmat[,2:ncol(lmat)]
fratio <- (t(t(lmat)%*%b)%*%solve(t(lmat)%*%varb2$varb%*%lmat)%*%((t(lmat)%*%b)))/(ncol(x)-1)
pfr <- 1-pf(fratio,(ncol(x)-1),dfres);
modsum <<- cbind(sqrt(r2),r2,mse,fratio,(ncol(x)-1),dfres,pfr)
modsuml <- c("R","R-sq","MSE",varb2$hcflab,"df1","df2","p")
modfin<-list("modres"=modres,"modsum"=modsum,"modresl"=modresl,"modsuml"=modsuml,"resid"=resid)
return(modfin)
}
}
if(full==0){return(modres)}
}
}
#bootstrapping checks
if(trunc(conf)>=100 | trunc(conf)<=50){conf<-95;notecode[notes,1]<-3;notes<-notes+1}
if ((abs(boot)>boot) | ((boot-trunc(boot)) != 0))
{boot<-abs(floor(boot));notecode[notes,1]<-4;notes<-notes+1;}
if(boot < 1000){boot<-1000;notecode[notes,1]<-5;notes<-notes+1;}
if (boot!=0)
{
bootsz<-boot;
cilow<-trunc(bootsz*(1-(conf/100))/2)
cihigh<-trunc((bootsz*(conf/100)+(bootsz*(1-(conf/100))/2)))+1
if(cilow < 0 | cihigh >= bootsz)
{
bad<-1;errs<-errs+1;errcode[errs,1]<-16;criterr<-1;
#while(bad==1)
#{
# bootsz<-floor((bootsz+1000)/1000)*1000
# cilowi<-trunc(bootsz*(1-(conf/100))/2)
# cihighi<-trunc((bootsz*(conf/100)+(bootsz*(1-(conf/100))/2)))+1
# if(cilowi > 0 & cihighi <= bootsz){bad<-0}
#}
}
}
if(boot > 0){boot<-bootsz}
if((adjust==1) & (boot > 0)){notecode[notes,1]<-8;notes<-notes+1;}
if((saveb==1) & (boot==0)){saveb<-0;notecode[notes,1]<-9;notes<-notes+1;}
bootmax<-trunc(bootmax)
if(bootmax<boot){bootmax<-2*boot;notecode[notes,1]<-12;notes<-notes+1;}
#end boot checks
p0=-.322232431088;p1=-1;p2=-.342242088547;p3=-.0204231210245;p4=-.0000453642210148;
q0=.0993484626060;q1=.588581570495;q2=.531103462366;q3=.103537752850;q4=.0038560700634;
alpha2<-(1-(conf/100))/2;cilm<-alpha2*2;
y5<-sqrt(-2*log(alpha2));
xp2<-(y5+((((y5*p4+p3)*y5+p2)*y5+p1)*y5+p0)/((((y5*q4+q3)*y5+q2)*y5+q1)*y5+q0));
if(missing(x)=="TRUE"|missing(y)=="TRUE"|missing(m)=="TRUE"){errs<-errs+1;errcode[errs,1]<-1;criterr<-1;}
if(criterr==0) #start bb
{
ytemp<-data[y];xtemp<-data[x];mtemp<-data[m]
ynames<-as.matrix(names(ytemp));xnames<-as.matrix(names(xtemp));mnames<-as.matrix(names(mtemp));
yNUM<-ncol(ytemp);xNUM<-ncol(xtemp);mNUM<-ncol(mtemp);
varnames<-as.matrix(rbind(ynames,xnames,mnames))
if(missing(cov) != "TRUE")
{
covtemp<-data[cov];covNUM<-ncol(covtemp);covnames<-as.matrix(names(covtemp))
varnames<-as.matrix(rbind(varnames,covnames));covYES=1
}
#mchecks
if(mNUM==1){mb=1}
if(mb>mNUM){errs<-errs+1;errcode[errs,1]<-6;criterr<-1}
mlabs <- rbind(" M1 :"," M2 :"," M3 :"," M4 :"," M5 :"," M6 :"," M7 :"," M8 :"," M9 :"," M10 :"," M11 :"," M12 :")
mblabs <- rbind(" MB1 :"," MB2 :"," MB3 :"," MB4 :"," MB5 :"," MB6 :"," MB7 :"," MB8 :"," MB9 :"," MB10 :"," MB11 :"," MB12 :")
xlabs <- rbind(" X1 :"," X2 :")
ylabs <- rbind(" Y1 :"," Y2 :")
pairNUM <- mNUM-mb
varlabs <- c(xlabs[1:xNUM], ylabs[1:yNUM])
if(pairNUM>0)
{pm <- mlabs[1:pairNUM];pmnames <- mnames[1:length(pm),]}
if(mb>0 & pairNUM>0)
{
npm <- mblabs[1:mb]
npmnames <- mnames[(length(pm)+1):mNUM,]
}
if(mb>0 & pairNUM==0)
{
npm <- mblabs[1:mb]
npmnames <- mnames[1:mNUM,]
}
varnl <- c(xnames,ynames);same<-0
for (k in 1:(length(varnames)-1))
{
for (p in (k+1):length(varnames))
{if(varnames[k,]==varnames[p,]){same<-1}}
}
if(same==1){errs<-errs+1;errcode[errs,1]<-12;criterr<-1}
} #end bb
if(criterr==0)#start aa
{
if(xNUM > 2){errs<-errs+1;errcode[errs,1]<-2;criterr<-1}
if(yNUM > 2){errs<-errs+1;errcode[errs,1]<-3;criterr<-1}
nonint<-mb-trunc(mb);medeven<-((mNUM-mb)/2)
if(nonint != 0){errs<-errs+1;errcode[errs,1]<-4;criterr<-1}
if(mb<0){errs<-errs+1;errcode[errs,1]<-5;criterr<-1}
if(mNUM>1)
{
if(medeven-trunc(medeven)!=0)
{errs<-errs+1;errcode[errs,1]<-7;criterr<-1}
}
if(mNUM>12){errs<-errs+1;errcode[errs,1]<-8;criterr<-1}
temp<-cbind(ytemp,xtemp,mtemp)
if(covYES==1){temp<-cbind(temp,covtemp)}
nrowC<-nrow(temp);datmat<-temp[complete.cases(temp),];nCHECK<-nrow(datmat);n<-nrow(datmat);
constant<-rep(1,n)
nMISSING <- nrowC-nCHECK
if(nMISSING>0){notes<-notes+1;notecode[notes,1]<-1}
ymat<-as.matrix(datmat[,1:yNUM]);xmat<-as.matrix(datmat[,(yNUM+1):(yNUM+xNUM)]);mmat<-as.matrix(datmat[,(yNUM+xNUM+1):(yNUM+xNUM+mNUM)])
if(covYES==1){covmat<-as.matrix(datmat[,(yNUM+xNUM+mNUM+1):(yNUM+xNUM+mNUM+covNUM)])}
vcheck<-sum(as.numeric((apply(datmat,2,FUN=min)-apply(datmat,2,FUN=max)==0)))
if(vcheck>0){errs<-errs+1;errcode[errs,1]<-11;criterr<-1}
} #end aa
if(criterr==0) #start cc
{
mdchk<-0;ydchk<-0
#check for dichomoous consequent
for(dc in 1:mNUM)
{
tempy <- as.matrix(mmat)
dicm<-length(table(tempy[,dc]));
if(dicm==2&mdchk==0){errs<-errs+1;errcode[errs,1]<-13;criterr<-1;mdchk<-1}
}
for(dc in 1:yNUM)
{
tempy <- as.matrix(ymat)
dicy<-length(table(tempy[,dc]))
if(dicy==2&ydchk==0){errs<-errs+1;errcode[errs,1]<-14;criterr<-1;ydchk<-1}
}
#covariates and cmatrix
if (covNUM > 0)
{
ccmat <- matrix(1,(mNUM+yNUM),covNUM)
if (cmatrix[1] != -999)
{
if (length(cmatrix) != ((mNUM+yNUM)*covNUM))
{errs<-errs+1;errcode[errs,1]<-9;criterr<-1} #cmatrix size error
if (criterr==0)
{
tmp<-1
for(i in 1:(mNUM+yNUM))
{
for(j in 1:covNUM)
{ccmat[i,j]<- 1-(as.numeric((cmatrix[tmp] == 0)));tmp<-tmp+1} #safeguarding >1 values
}
if(as.numeric(sum(ccmat)==0) != 0)
{errs<-errs+1;errcode[errs,1]<-10;criterr<-1} #must assign covariate
}
}
}
xtom<-matrix(-999, xNUM, mNUM);mtoy<-matrix(-999,yNUM,mNUM);xtoy<-matrix(-999, xNUM, yNUM)
seqa=c(1:12);seqa2<-c(2,3)
colxtom=seqa2[1:xNUM];colmtoy=seqa[1:mNUM];colxtoy<-colxtom
if(covNUM > 0)
{
if((sum(ccmat))<(nrow(ccmat)*ncol(ccmat)))
{total<-0;notecode[notes]<-11;notes=notes+1}
}
} #end cc
###note: type 1 is ols, full 0 is coeff, full 1 is everything
###note: contrast 0=no, 1 = dif con, 2 = abs con, 3 = dif all, 4 = abs all
cat("********************** MEDYAD Procedure for R Beta Version 1.1 **********************\n")
cat(" Written by Jacob Coutts, Andrew Hayes, and Tao Jiang \n")
cat("*************************************************************************************\n")
if(criterr==0)
{
obsmat<-matrix(-999,1,6);residmat<-matrix(0,n,1);totresid<-matrix(0,n,1);totmat<-matrix(-999,1,6)
ptemp=(varnl);cat("Model Variables:\n");ptemp <- as.matrix(ptemp);rownames(ptemp)<-varlabs;
cat(ptemp,fill=(xNUM+yNUM), labels=varlabs)
if(pairNUM>0){cat("\nPaired Mediators:\n");ptemp<-pmnames;cat(ptemp,fill=mNUM,labels=pm)}
if(mb>0){cat("\nUnpaired Mediators:\n");ptemp<-npmnames;cat(ptemp,fill=mb,labels=npm)}
if(covNUM>0){cat("\nCovariates:\n");cat(covnames,fill=covNUM)}
sn<-n;desclab<-"N:";
if(seed!=0){sn<-c(sn,seed);desclab=c(desclab,"Seed:")}
cat("\n");cat(fill=1,sn,labels=desclab);cat("\n")
if(describe!=0)
{
cat("************************ DESCRIPTIVES FOR MODEL VARIABLES ***************************\n")
datmatt<-as.matrix(cbind(xmat,mmat,ymat));varnamest<-c(xnames,mnames,ynames)
if(covNUM>0){
datmatt<-as.matrix(cbind(datmatt,covmat));varnamest<-c(varnamest,covnames)
}
sigmatal<-cov(datmatt)
SD<-sqrt(diag(sigmatal));sdall<-diag(1/SD);corall<-sdall%*%sigmatal%*%t(sdall)
Mean <- colMeans(datmatt)
Min<-apply(datmatt,2,FUN=min);Max<-apply(datmatt,2,FUN=max)
descr<-c("Mean","SD","Min","Max");desctot<-cbind(Mean,SD,Min,Max);rownames(desctot) <- varnamest
cat("\nDescriptive statistics of model variables\n");
(print(specify_decimal(desctot,decimals),quote=FALSE, right=TRUE));
cat("\n----------------\n")
mainnames<-c(xnames,mnames,ynames)
cormv<-corall[1:(xNUM+mNUM+yNUM),1:(xNUM+mNUM+yNUM)]
cat("\nCorrelations between antecendents and consequents\n");prmatrix(specify_decimal(cormv,decimals),rowlab=mainnames, collab=mainnames,quote=FALSE,right=TRUE);
if(covNUM>0)
{
coreall<-corall[1:(xNUM+mNUM+yNUM),(xNUM+mNUM+yNUM+1):ncol(corall)]
cat("\nCorrelations for model covariates\n");prmatrix(specify_decimal(coreall,decimals),rowlab=mainnames,collab=covnames,quote=FALSE, right=TRUE);cat("\n")
}
}
cat("*************************************************************************************\n")
rcoeff <- matrix("constant",1,1)
#pulling relevant model data
outmat <- cbind(mmat,ymat);bootcnum<-0;
for(i in 1:(mNUM+yNUM))
{
outname<-as.matrix(rbind(mnames,ynames));outname<-outname[i,];outv<-outmat[,i]
if(covNUM==0)
{
predname<-as.matrix(rbind(rcoeff,xnames))
predv<-as.matrix(cbind(constant,xmat))
if(i>mNUM)
{
predname<-as.matrix(rbind(predname,mnames));predv<-as.matrix(cbind(predv,mmat))
}
antect<-predname;antec<-rbind(antec,antect)
}
if(covNUM>0)
{
nc<-ncol(ccmat);cdata<-matrix(999,n,1);newcn<-matrix("999",1,1)
for(j2 in 1:nc)
{
if(ccmat[i,j2]==1){cdata<-cbind(cdata,covmat[,j2]);newcn<-rbind(newcn,covnames[j2,])}
}
ccmatmp<-ccmat[i,];
ccmatmp<-sum(ccmatmp)
if(ccmatmp>0){cdata<-cdata[,2:ncol(cdata)];newcn<-as.matrix(newcn[2:nrow(newcn),])}
predname<-rbind(rcoeff, xnames);predv<-as.matrix(cbind(constant,xmat));
ccmatmp<-ccmat[i,];
ccmatmp<-sum(ccmatmp)
if(ccmatmp>0)
{predname<-rbind(rcoeff,xnames,newcn);predv<-cbind(constant,xmat,cdata)}
if(i > mNUM)
{
predname<-rbind(rcoeff,xnames,mnames);predv<-as.matrix(cbind(constant,xmat,mmat));
ccmatmp<-ccmat[i,]
ccmatmp<-sum(ccmatmp)
if(ccmatmp>0)
{predname<-rbind(rcoeff,xnames,mnames,newcn);predv<-cbind(constant,xmat,mmat,cdata)}
}
antect<-predname
antec<-rbind(antec,antect)
}
blab<-t(predname)
modestmo <- modest(y=outv,x=predv,type=1,full=2);residmat<-cbind(residmat,modestmo$resid);
cat("Outcome:");cat("\n")
cat(outname)
cat("\n");cat("\n");cat("Model Summary");cat("\n")
rrlab <- "";
prmatrix(specify_decimal(modestmo$modsum,decimals), collab=modestmo$modsuml, rowlab=rrlab, quote=FALSE, right=TRUE)
cat("\n");cat("Model");cat("\n")
prmatrix(specify_decimal(modestmo$modres,decimals), collab=modestmo$modresl,rowlab=blab,quote=FALSE, right=TRUE)
if(i<(mNUM+yNUM))
{cat("*************************************************************************************\n")}
conseqt<-matrix(outname,nrow(modestmo$modres),1);conseq<-rbind(conseq,conseqt)
if(i <= mNUM)
{xtom[1:xNUM,i]<-colxtom;colxtom<-colxtom+nrow(modestmo$modres)}
if (i>mNUM)
{mtoy[(i-mNUM),]<-colmtoy+bootcnum+xNUM+1;xtoy[,(i-mNUM)]<-colxtoy+bootcnum}
bootcnum<-bootcnum+nrow(modestmo$modres);obsmat<-rbind(obsmat,modestmo$modres)
}
obsmat<-obsmat[2:nrow(obsmat),];obscoeff<-t(obsmat[,1])
cat("\n*************************** RESIDUAL CORRELATION MATRIX *****************************\n")
residmat<-residmat[,2:ncol(residmat)]
sigmatr<-(t(residmat)%*%(diag(n)-(1/n)*constant%*%t(constant))%*%residmat)/(n-1)
SDr<-sqrt(diag(sigmatr));sdr<-diag(1/SDr);corr<-sdr%*%sigmatr%*%t(sdr)
rclb<-c(mnames,ynames);rrlb<-t(rclb);
cat("\n");prmatrix(specify_decimal(corr,decimals),rowlab=rrlb, collab=rclb,quote=FALSE, right=TRUE);cat("\n");
if(total!=0)
{
totx<-as.matrix(cbind(constant,xmat))
totlab<-c("constant",xnames)
if(covNUM>0)
{totx<-as.matrix(cbind(totx,cdata));totlab<-c(totlab,covnames)}
cat("***************************** TOTAL EFFECT(S) MODEL(S) ******************************\n")
for(i in 1:yNUM)
{
totname<-ynames[i,]
ymat <- as.matrix(ymat);tymat<-ymat[,i]
modestmo2<-modest(x=totx,y=tymat,type=1,full=2)
totresid<-cbind(totresid,modestmo2$resid);
modres2<-modestmo2$modres[2:(1+xNUM),]
totmat<-rbind(totmat,modres2)
cat("Outcome:");cat("\n")
cat(totname)
cat("\n");cat("\n");cat("Model Summary");cat("\n")
rrlab <- ""
prmatrix(specify_decimal(modestmo2$modsum,decimals), collab=modestmo2$modsuml, rowlab = rrlab,quote=FALSE, right=TRUE)
cat("\n");cat("Model");cat("\n")
prmatrix(specify_decimal(modestmo2$modres,decimals), collab=modestmo2$modresl,rowlab=totlab,quote=FALSE, right=TRUE)
if(i < yNUM)
{cat("-------------------------------\n")}
}
totmat<-as.matrix(totmat[2:nrow(totmat),])
if(xNUM*yNUM==4)
{
totfix<-c(1,3,2,4)
totmat<-totmat[totfix,]
}
if (yNUM==2)
{
cat("-------------------------------\n")
totresid<-totresid[,2:ncol(totresid)]
sigmatot<-(t(totresid)%*%(diag(n)-(1/n)*constant%*%t(constant))%*%totresid)/(n-1)
sdtotr<-sqrt(diag(sigmatot));sdtotr<-diag(1/sdtotr);totcorr<-sdtotr%*%sigmatot%*%t(sdtotr)
cat("Correlation between residuals:\n");cat(specify_decimal(totcorr[1,2],decimals));cat("\n");
}
}
bootres <- matrix(-99, boot, bootcnum)
bs<-1
for(bs in 1:boot) #start aaa
{
nogo<-0;bootstmp<-0;s<-floor(runif(matrix(seed,n,1))*n)+1
bootm<-mmat[s,];bootx<-xmat[s,];
if(covNUM>0){bootcov<-as.matrix(covmat[s,])}
for(g in 1:(mNUM+yNUM))#start bbb
{
boutv <- outmat[s,g]
if(covNUM==0)
{bpredv<-as.matrix(cbind(constant,bootx));if(g>mNUM){bpredv<-as.matrix(cbind(bpredv,bootm))}}
if(covNUM>0)
{
bcdata<-matrix(999,n,1)
for(jc in 1:nc)
{if(ccmat[g,jc]==1){bcdata<-cbind(bcdata,bootcov[,jc])}}
if(sum(ccmat[g,])>0){bcdata<-bcdata[,2:ncol(bcdata)]}
bpredv<-as.matrix(cbind(constant,bootx))
if(sum(ccmat[g,])>0){bpredv<-cbind(constant,bootx,bcdata)}
if(g>mNUM)
{
bpredv<-as.matrix(cbind(constant,bootx,bootm))
ccmatmp<-ccmat[g,]
ccmatmp<-sum(ccmatmp)
if(ccmatmp>0){bpredv<-as.matrix(cbind(constant,bootx,bootm,bcdata))}
}
}
tmp<-cbind(bpredv,boutv)
totcheck<-sum(as.numeric((apply(tmp,2,FUN=min)-apply(tmp,2,FUN=max)>1)))
if(totcheck>0){nogo<-1}
if(nogo==0)
{
modestmo3 <- modest(y=boutv,x=bpredv,type=1,full=0)
bootstmp <- cbind(bootstmp,t(modestmo3))
}
}#end bbb
if(nogo==1){badboot<-badboot+1}
if(nogo==0){go<-go+1;goodboot<-goodboot+1;bootres[go,]<-bootstmp[1,2:ncol(bootstmp)]}
}#end aaa
if(badboot>0){notecode[notes,1]<-10;notes<-notes+1}
if(goodboot<0){errs<-errs+1;errcode[errs,1]<-15;criterr<-1}
}
if(criterr==0) #start a
{
if(total==0)
{cat("**************************** DIRECT AND INDIRECT EFFECTS ****************************\n")}
if(total==1)
{cat("\n************************ TOTAL, DIRECT, AND INDIRECT EFFECTS ************************\n")}
temp<-nrow(xtom)*nrow(mtoy)*ncol(mtoy);
indpaths<-matrix(" --> ",temp,7)
indcnt<-1;dircnt<-1
hclab <- c("se(HC0)","se(HC1)","se(HC2)","se(HC3)","se(HC4)","se")
hclab <- hclab[(hc+1)]
befflab<-c("effect","BootSE","BootLLCI","BootULCI");efflabs<-c("effect",hclab,"t","p","LLCI","ULCI")
obsind<-matrix(-999,1,(xNUM*mNUM*yNUM))
if(boot>0){booti<-matrix(-999,boot,(xNUM*mNUM*yNUM))}
dirmat<-matrix(-999,(xNUM*yNUM),6)
ind<-c("Ind");lab<-1:48;indlab<-paste(ind,lab,sep="")
indlab<-indlab[1:(xNUM*mNUM*yNUM)]
if(contrast<0 | contrast > 4){notes<-notes+1;notecode[notes,1]<-7;contrast<-0}
if((length(varnames)==3) & (contrast!=0)){notes<-notes+1;notecode[notes,1]<-6;contrast<-0}
if(contrast!=0)
{
cpairs<-(mNUM*yNUM);cpairsc<-((cpairs*(cpairs-1))/2)
nconb<-matrix(999,4,3)
left<-"(C";mid<-1:600;end<-")";clabt<-paste(left,mid,end,sep="")
npairmed<-pairNUM/2
pairloc<-c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10)
npairloc<-c(1:10)
pairmtch<-999
if(npairmed>0){pairmtch<-c(pairmtch,pairloc[1:(npairmed*2)])}
if(mb>0){pairmtch<-c(pairmtch,(npairmed+npairloc[1:mb]))}
pairmtch<-t(as.matrix(pairmtch[2:length(pairmtch)]))
if(yNUM>1){pairmtch<-cbind(pairmtch,pairmtch)}
}
cstart<-1
for(indlp1 in 1:nrow(xtom)) #start c
{
contot <-matrix(999,1,4);ckey<-matrix(" ",1,4)
for(indlp2 in 1:nrow(mtoy))
{
for(indlp3 in 1:ncol(mtoy))
{
obsind[1,indcnt]<-obscoeff[1,xtom[indlp1,indlp3]]%*%obscoeff[1,mtoy[indlp2,indlp3]]
booti[,indcnt]<-bootres[,xtom[indlp1,indlp3]]*bootres[,mtoy[indlp2,indlp3]]
indpaths[indcnt,1]<-indlab[indcnt]
indpaths[indcnt,2]<-":"
indpaths[indcnt,3]<-xnames[indlp1,]
indpaths[indcnt,5]<-mnames[indlp3,]
indpaths[indcnt,7]<-ynames[indlp2,]
indcnt<-indcnt+1
}
}
for(dirlp in 1:ncol(xtoy)){dirmat[dircnt,]<-obsmat[xtoy[indlp1,dirlp],];dircnt<-dircnt+1}
xntemp<-xnames[indlp1,]
if(total==0)
{cat("DIRECT AND INDIRECT EFFECTS OF\n");cat(xntemp);cat("\n")}
if(total==1)
{
cat("TOTAL, DIRECT, AND INDIRECT EFFECTS OF\n");cat(xntemp);cat("\n")
if(xNUM*yNUM*mNUM==1 | (xNUM==1 & yNUM==1)){
totmat2<-t(totmat)
} else{
totmat2<-as.matrix(totmat[(((indlp1-1)*yNUM)+1):(yNUM*indlp1),])}
if(ncol(totmat2)==1){totmat2<-t(totmat2)}
cat("\nTotal effect(s) on\n");prmatrix(specify_decimal(totmat2,decimals), rowlab = ynames, collab=efflabs, quote=FALSE, right=TRUE)
}
indres<-t(obsind)
indres2<-as.matrix(indres[(((indlp1-1)*(mNUM*yNUM))+1):(indlp1*(mNUM*yNUM)),])
indlab2<-indlab[(((indlp1-1)*(mNUM*yNUM))+1):(indlp1*(mNUM*yNUM))]
indpath2<-indpaths[(((indlp1-1)*(mNUM*yNUM))+1):(indlp1*(mNUM*yNUM)),]
if(xNUM*yNUM*mNUM==1){
dirmat2<-dirmat
} else{
dirmat2<-as.matrix((dirmat[(((indlp1-1)*yNUM)+1):(yNUM*indlp1),]))}
if(ncol(dirmat2)==1){dirmat2<-t(dirmat2)}
cat("\nDirect effect(s) on\n");prmatrix(specify_decimal(dirmat2,decimals), rowlab=ynames,collab=efflabs,quote=FALSE, right=TRUE)
bootind<-booti[,(((indlp1-1)*(mNUM*yNUM))+1):(indlp1*(mNUM*yNUM))];
bootmat<-as.matrix(bootconf(boottmp=bootind))
if((xNUM==1 & mNUM==1 & yNUM==1) | (xNUM==2 & mNUM==1 & yNUM==1)){
indres2<-as.matrix(t(c(indres2,bootmat[,2:ncol(bootmat)])));} else{
indres2<-cbind(indres2,bootmat[,2:ncol(bootmat)])
}
if((contrast!=0)&((mNUM*yNUM)>1))
{
counter<-1
for(o in 1:(cpairs-1))
{
for(b in (o+1):cpairs)
{
doit<-1
if((pairmtch[1,o] != pairmtch[1,b]) & (yNUM > 1) & (contrast <3))
{if((o<=(cpairs/2)) & (b>(cpairs/2))){doit<-0}}
if(doit==1)
{
conboot<-as.matrix(bootind[,o]-bootind[,b])
contindt<-as.matrix(indres2[o,1]-indres2[b,1])
if((contrast==2)|(contrast==4))
{
conboot<-as.matrix(abs(bootind[,o])-abs(bootind[,b]))
contindt<-as.matrix(abs(indres2[o,1])-abs(indres2[b,1]))
}
ckeyt<- cbind(" : ",indpath2[o,1], " minus ",indpath2[b,1])
ckey<-rbind(ckey,ckeyt)
bootmat2<-bootconf(boottmp=conboot)
contott<-cbind(contindt,t(as.matrix(bootmat2[,2:ncol(bootmat2)])))
contot<-rbind(contot,contott);counter<-counter+1
}
}
}
contot<-as.matrix(contot[2:nrow(contot),])
if(ncol(contot)==1){contot<-t(contot)}
conlab<-as.matrix(clabt[cstart:(cstart+nrow(contot)-1)])
cstart<-cstart+nrow(contot)
conlabt<-rbind(" ", conlab)
indres2<-rbind(indres2,contot)
indlab2<-rbind(as.matrix(indlab2),conlab)
}
cat("\nIndirect Effect(s):\n");prmatrix(specify_decimal(indres2,decimals), rowlab=indlab2,collab=befflab, quote=FALSE, right=TRUE)
if((xNUM==1 & mNUM==1 & yNUM==1) | (xNUM==2 & mNUM ==1 & yNUM==1)){indpath2<-as.matrix(t(indpath2))}
cat("\nIndirect Effect Key:");prmatrix(indpath2, quote=FALSE, rowlab=rep(" ", length(indlab2)), collab=rep(" ",7));
if((contrast!=0)&((mNUM*yNUM)>1))
{
ckey<-ckey[2:nrow(ckey),]
if(mNUM==1){ckey<-t(ckey)}
if(mNUM==2 & yNUM==1){ckey<-t(as.matrix(ckey))}
cat("\nContrast of Indirect Effects Key:");prmatrix(ckey, rowlab=conlab, collab=rep(" ",7), quote=FALSE)
}
if(indlp1 < nrow(xtom)){cat("\n--------------------\n")}
}#end c
if((xNUM==2) & (contrast!=0)) #start b
{
cat("\n--------------------\n")
contmp<-matrix(-99,1,4)
ckey<-matrix(" ",1,4)
for(i in 1:ncol(pairmtch))
{
for(j in 1:ncol(pairmtch))
{
doit<-1
if((pairmtch[1,i] != pairmtch[1,j]) & (contrast<3)){doit<-0}
if(doit==1)
{
if((contrast==1)|(contrast==3))
{
conboot<-as.matrix(booti[,i]-booti[,(j+ncol(pairmtch))])
bootmat3<-bootconf(boottmp=conboot)
bootmat3[1,1]<-as.matrix(obsind[1,i]-obsind[1,(j+ncol(pairmtch))])
}
if((contrast==2)|(contrast==4))
{
conboot<-as.matrix(abs(booti[,i])-abs(booti[,(j+ncol(pairmtch))]))
bootmat3<-bootconf(boottmp=conboot)
bootmat3[1,1]<-as.matrix(abs(obsind[1,i])-abs(obsind[1,(j+ncol(pairmtch))]))
}
ckeyt<-cbind(" : ",indlab[i]," minus ",indlab[((length(indlab)/2)+j)])
ckey<-rbind(ckey,ckeyt);contmp<-rbind(contmp,bootmat3)
}
}
}
contmp<-contmp[2:nrow(contmp),]
if((mNUM*yNUM)>1)
{
conlab2<-clabt[((2*nrow(conlab))+1):(2*nrow(conlab)+nrow(contmp))]
cat("CONTRASTS BETWEEN INDIRECT EFFECTS OF X1 and X2\n");prmatrix(specify_decimal(contmp,decimals),rowlab=conlab2,collab=befflab,quote=FALSE, right=TRUE)
}
if(yNUM==1){contmp<-t(contmp)}
if((mNUM*yNUM)==1)
{
conlab2<-clabt[1]
cat("CONTRASTS BETWEEN INDIRECT EFFECTS OF X1 and X2\n");prmatrix(specify_decimal(contmp,decimals),rowlab=conlab2,collab=befflab,quote=FALSE, right=TRUE)
}
if((xNUM==2 & mNUM==1 & yNUM==1)){
ckey<-(t(as.matrix(ckey[2,])))} else{
ckey<-ckey[2:nrow(ckey),]
}
cat("\nContrast of Indirect Effects Key:");prmatrix(ckey,rowlab=conlab2,collab=rep(" ",7),quote=FALSE)
}#end b
if(saveb==1 & boot > 0 )
{
cat("\n**************************************************************************\n")
write.csv(bootres, file="medyadres.csv", row.names=FALSE)
cat("Bootstrap estimates were saved to a file in your working directory (medyadres.csv)\n")
save<-c("V");lab<-1:400;savelab<-paste(save,lab,sep="")
savelab<-savelab[1:ncol(bootres)]
conseq<-conseq[2:nrow(conseq),]
antec<-antec[2:nrow(antec),]
savelab<-cbind(savelab,conseq,antec)
cat("\nColumn names key for saved model coefficients:\n");prmatrix(savelab, collab=c(" ","Conseqnt","Antecdnt"),quote=FALSE)
}
} #end a
cat("\n**************************** ANALYSIS NOTES AND ERRORS ******************************\n")
if(criterr==0)
{
cat(sprintf("The percentage for all confidence intervals in the output was: %s\n",specify_decimal(conf,decimals)))
if(boot!=0)
{cat(sprintf("\nThe total number of bootstrap samples used to construct percentile bootstrap confidence intervals was: %s\n",boot))}
if (contrast!=0)
{
if ((contrast==1) | (contrast==3))
{cat("\nContrast(s) computed through difference of indirect effects.\n")}
if ((contrast==2) | (contrast==4))
{cat("\nContrast(s) computed through difference of absolute values of indirect effects.\n")}
}
for(i in 1:notes)
{
if(notecode[i,1]==1)
{cat(sprintf("\nNOTE: Some cases were deleted due to missing data. The number of such cases was: %s\n",nMISSING))}
if(notecode[i,1]==2)
{cat("\nNOTE: A heteroscedasticity-consistent standard error and covariance matrix estimator was used.\n")}
if(notecode[i,1]==3)
{cat("\nNOTE: Confidence level restricted to between 50 and 99.9999%. 95% confidence is provided in the output.\n")}
if(notecode[i,1]==4)
{cat("\nNOTE: Number of bootstrap samples created was rounded down to the nearest integer.\n")}
if(notecode[i,1]==5)
{cat("\nNOTE: MEDYAD requires a minimum of 1000 bootstrap samples. One thousand was used.\n")}
if(notecode[i,1]==6)
{cat("\nNOTE: There are no constrasts to conduct.\n")}
if(notecode[i,1]==7)
{cat("\nNOTE: Invalid contrast option requested. No contrasts are conducted.\n")}
if(notecode[i,1]==8)
{cat("\nNOTE: The number of bootstrap samples was adjusted upward given your desired confidence.\n")}
if(notecode[i,1]==10)
{cat(sprintf("\nNOTE: Some bootstrap samples had to be replaced. The number of times this happened was: %s.\n",badboot))}
if(notecode[i,1]==11)
{cat("\nNOTE: Total effect(s) are generated only when all covariates are assigned to all model equations.\n")}
if(notecode[i,1]==12)
{cat("\nNOTE: BOOTMAX set to twice the number of bootstrap samples requested.\n")}
}
}
if (errs>0)
{
for(i in 1:errs)
{
if(errcode[i,1]==1)
{cat("\nERROR: You must specify a Y, X, and M variable.\n")}
if(errcode[i,1]==2)
{cat("\nERROR: Too many variables entered for X. No more than two allowed.\n")}
if(errcode[i,1]==3)
{cat("\nERROR: Too many variables entered for Y. No more than two allowed.\n")}
if(errcode[i,1]==4)
{cat("\nERROR: Number of mediators cannot contain fractions.")}
if(errcode[i,1]==5)
{cat("\nERROR: There cannot be a negative number of mediators.\n")}
if(errcode[i,1]==6)
{cat("\nERROR: MB cannot exceed the number of mediators.\n")}
if(errcode[i,1]==7)
{cat("\nERROR: Improper number of mediators entered.\n")}
if(errcode[i,1]==8)
{cat("\nERROR: The total number of mediators may not exceed 12.\n")}
if(errcode[i,1]==9)
{cat("\nERROR: Incorrect number of 1s or 0s entered for cmatrix.\n")}
if(errcode[i,1]==10)
{cat("\nERROR: You must assign a covariate to M or Y.\n")}
if(errcode[i,1]==11)
{cat("\nERROR: One of the variables in the model is a constant.\n")}
if(errcode[i,1]==12)
{cat("\nERROR: A variable cannot play more than one role (e.g., X and M).\n")}
if(errcode[i,1]==13)
{cat("\nERROR: M is a dichotomous variable. MEDYAD does not allow this.\n")}
if(errcode[i,1]==14)
{cat("\nERROR: Y is a dichotomous variable. MEDYAD does not allow this.\n")}
if(errcode[i,1]==15)
{cat("\nERROR: Bootstrapping did not finish. Results have been suppressed. Try increasing bootmax.\n")}
if(errcode[i,1]==16)
{cat("\nNOTE: Your requested level of confidence is too high given your number of bootstrap samples. Adjust conf or boot.\n")}
}
}
options(width=widthex)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment