#Libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Average DI models

No Covariates

#####Latinx ###
#  class x instrumental x moral 



#### White #### 
# Class x moral 

#class x instrumental x moral 

Pre-registered

### Asian ### 
#Instrumental  x moral 
# (0.279811631433574)^2  ### standard errors that need to be squared, taken from SPSS file
# (.018)^2 
# (.023)^2
# (.024)^2

z1=-2  #supply lower bound for w2 here     ###change z1 to -2 and z2 to 2
z2=2 #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.022+-0.049*z)+(1.9602*sqrt(0.000324+(2*z*0)+((z^2)*0.000576)))
y2 <- (-0.022+-0.049*z)-(1.9602*sqrt(0.000324+(2*z*0)+((z^2)*0.000576)))
fy <- c(y1,y2)
fline <- (-0.022+-0.049*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moral',ylab='Simple Slope',main='Confidence Bands for Asian Participants', 
     sub= "Model with pre-registered covariates")
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-11.8223,col=4,lty=2)
abline(v=0.3427,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.6371,1.9773)   #  <-- change to alter plot dims
leg <- c(-2,1.6796)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8425,1.9488)
y2 <- c(1.9099,1.8213)
y3 <- c(1.9773,1.6938)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.022167198640804+-0.048726847459291*z)+(1.9602*sqrt(0.000323154+(2*z*0.000134427123867)+((z^2)*0.000567663)))
y2 <- (-0.022167198640804+-0.048726847459291*z)-(1.9602*sqrt(0.000323154+(2*z*0.000134427123867)+((z^2)*0.000567663)))
fy <- c(y1,y2)
fline <- (-0.022167198640804+-0.048726847459291*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-6.4393,col=4,lty=2)
abline(v=0.6032,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.6229,2.001)   #  <-- change to alter plot dims
leg <- c(-2,1.6701)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8504,1.9251)
y2 <- c(1.9257,1.8055)
y3 <- c(2.001,1.6859)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.030050340622283+-0.048726847459291*z)+(1.9602*sqrt(0.000517225+(2*z*0.000146068975076)+((z^2)*0.000567663)))
y2 <- (-0.030050340622283+-0.048726847459291*z)-(1.9602*sqrt(0.000517225+(2*z*0.000146068975076)+((z^2)*0.000567663)))
fy <- c(y1,y2)
fline <- (-0.030050340622283+-0.048726847459291*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-9.9165,col=4,lty=2)
abline(v=0.5661,col=4,lty=2)

All covariates

### Asian ### 
#Instrumental  x moral 
# 51.806^2
# .025^2
# .026^2
# 0.026455648597971^2

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.049812771300131+-0.071733509660943*z)+(1.9602*sqrt(0.000625+(2*z*0.000213611919503)+((z^2)*0.0006999013)))
y2 <- (-0.049812771300131+-0.071733509660943*z)-(1.9602*sqrt(0.000625+(2*z*0.000213611919503)+((z^2)*0.0006999013)))
fy <- c(y1,y2)
fline <- (-0.049812771300131+-0.071733509660943*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moral',ylab='Simple Slope',main='Confidence Bands for Asian Participants', 
     sub= "Model with all covariates")
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.2264,col=4,lty=2)
abline(v=-0.0146,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-122.1646,-121.5812)   #  <-- change to alter plot dims
leg <- c(-2,-122.0917)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-121.8056,-121.7179)
y2 <- c(-121.6934,-121.8927)
y3 <- c(-121.5812,-122.0674)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.049812771300131+-0.071733509660943*z)+(1.9602*sqrt(0.000617262+(2*z*0.000213611919503)+((z^2)*0.000699901)))
y2 <- (-0.049812771300131+-0.071733509660943*z)-(1.9602*sqrt(0.000617262+(2*z*0.000213611919503)+((z^2)*0.000699901)))
fy <- c(y1,y2)
fline <- (-0.049812771300131+-0.071733509660943*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.2209,col=4,lty=2)
abline(v=-0.0201,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-122.1312,-121.6369)   #  <-- change to alter plot dims
leg <- c(-2,-122.0694)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-121.8242,-121.6622)
y2 <- c(-121.7305,-121.8555)
y3 <- c(-121.6369,-122.0488)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.031244760159169+-0.071733509660943*z)+(1.9602*sqrt(0.000699409+(2*z*0.000092661424857)+((z^2)*0.000699901)))
y2 <- (-0.031244760159169+-0.071733509660943*z)-(1.9602*sqrt(0.000699409+(2*z*0.000092661424857)+((z^2)*0.000699901)))
fy <- c(y1,y2)
fline <- (-0.031244760159169+-0.071733509660943*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.9013,col=4,lty=2)
abline(v=0.3664,col=4,lty=2)

##### White ######
#Instrumental  x moral 
# 16.570172787277148^2
# 0.02727730591327^2
# 0.029156796518783^2
# 0.030083434263382^2

z1=-2  #supply lower bound for w2 here
z2=2  #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.058810453573026+-0.077933591478901*z)+(1.9602*sqrt(0.0007440514+(2*z*0.000210752565508)+((z^2)*0.000905013)))
y2 <- (-0.058810453573026+-0.077933591478901*z)-(1.9602*sqrt(0.0007440514+(2*z*0.000210752565508)+((z^2)*0.000905013)))
fy <- c(y1,y2)
fline <- (-0.058810453573026+-0.077933591478901*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moral',ylab='Simple Slope',main='Confidence Bands for White participants', 
     sub="Model with all covariates")
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.8251,col=4,lty=2)
abline(v=-0.0818,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-6.5666,-5.9102)   #  <-- change to alter plot dims
leg <- c(-2,-6.4846)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-6.1225,-6.046)
y2 <- c(-6.0164,-6.2516)
y3 <- c(-5.9102,-6.4572)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.058810453573026+-0.077933591478901*z)+(1.9602*sqrt(0.000744051+(2*z*0.000210752565508)+((z^2)*0.000905013)))
y2 <- (-0.058810453573026+-0.077933591478901*z)-(1.9602*sqrt(0.000744051+(2*z*0.000210752565508)+((z^2)*0.000905013)))
fy <- c(y1,y2)
fline <- (-0.058810453573026+-0.077933591478901*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.8251,col=4,lty=2)
abline(v=-0.0818,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-6.5503,-5.9374)   #  <-- change to alter plot dims
leg <- c(-2,-6.4737)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-6.1315,-6.0188)
y2 <- c(-6.0345,-6.2335)
y3 <- c(-5.9374,-6.4482)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.049751800835947+-0.077933591478901*z)+(1.9602*sqrt(0.000850119+(2*z*0.000216616531576)+((z^2)*0.000905013)))
y2 <- (-0.049751800835947+-0.077933591478901*z)-(1.9602*sqrt(0.000850119+(2*z*0.000216616531576)+((z^2)*0.000905013)))
fy <- c(y1,y2)
fline <- (-0.049751800835947+-0.077933591478901*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.4691,col=4,lty=2)
abline(v=0.1234,col=4,lty=2)

##### Latinx ######
#Instrumental  x moral 

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(31.614,32.3003)   #  <-- change to alter plot dims
leg <- c(-2,31.6997)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(32.0063,32.0518)
y2 <- c(32.1533,31.8901)
y3 <- c(32.3003,31.7283)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.065804040970812+-0.077187956032932*z)+(1.9602*sqrt(0.001397008+(2*z*0.000361428377438)+((z^2)*0.001723698)))
y2 <- (-0.065804040970812+-0.077187956032932*z)-(1.9602*sqrt(0.001397008+(2*z*0.000361428377438)+((z^2)*0.001723698)))
fy <- c(y1,y2)
fline <- (-0.065804040970812+-0.077187956032932*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=0.1424,col=4,lty=2)
abline(v=10.956,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(31.6987,32.2272)   #  <-- change to alter plot dims
leg <- c(-2,31.7648)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(31.9478,32.2272)
y2 <- c(32.0364,32.007)
y3 <- c(32.125,31.7868)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.007351157453784+-0.077187956032932*z)+(1.9602*sqrt(0.001602595+(2*z*0.000384031628992)+((z^2)*0.001723698)))
y2 <- (-0.007351157453784+-0.077187956032932*z)-(1.9602*sqrt(0.001602595+(2*z*0.000384031628992)+((z^2)*0.001723698)))
fy <- c(y1,y2)
fline <- (-0.007351157453784+-0.077187956032932*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

Imputed covariates

### Asian ### 
#Instrumental  x moral
# 29.5391792683096^2
# 0.017902362570038^2
# 0.020667094356348^2
# 0.021724301732672^2
z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.051428193605124+-0.046507864098218*z)+(1.9602*sqrt(0.0003204946+(2*z*0.000137705663667)+((z^2)*0.0004719453)))
y2 <- (-0.051428193605124+-0.046507864098218*z)-(1.9602*sqrt(0.0003204946+(2*z*0.000137705663667)+((z^2)*0.0004719453)))
fy <- c(y1,y2)
fline <- (-0.051428193605124+-0.046507864098218*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moral',ylab='Simple Slope',main='Confidence Bands for Asian participants', 
     sub= "Model with imputed covariates")
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-10.262,col=4,lty=2)
abline(v=-0.394,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-95.3431,-94.873)   #  <-- change to alter plot dims
leg <- c(-2,-95.2844)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-94.9728,-94.9925)
y2 <- c(-94.9229,-95.1286)
y3 <- c(-94.873,-95.2648)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.051428193605124+-0.046507864098218*z)+(1.9602*sqrt(0.000320495+(2*z*0.000137705663667)+((z^2)*0.000471945)))
y2 <- (-0.051428193605124+-0.046507864098218*z)-(1.9602*sqrt(0.000320495+(2*z*0.000137705663667)+((z^2)*0.000471945)))
fy <- c(y1,y2)
fline <- (-0.051428193605124+-0.046507864098218*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-10.262,col=4,lty=2)
abline(v=-0.394,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-95.3282,-94.8979)   #  <-- change to alter plot dims
leg <- c(-2,-95.2744)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-94.9811,-94.9676)
y2 <- c(-94.9395,-95.1121)
y3 <- c(-94.8979,-95.2565)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.04313124191535+-0.046507864098218*z)+(1.9602*sqrt(0.000427129+(2*z*0.000108787420361)+((z^2)*0.000471945)))
y2 <- (-0.04313124191535+-0.046507864098218*z)-(1.9602*sqrt(0.000427129+(2*z*0.000108787420361)+((z^2)*0.000471945)))
fy <- c(y1,y2)
fline <- (-0.04313124191535+-0.046507864098218*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-9.0146,col=4,lty=2)
abline(v=-0.0695,col=4,lty=2)

##### White ######
#Instrumental  x moral 
# 11.959276395221181^2
# 0.021771809019935^2
# 0.024146676771613^2
# 0.02616276070789^2

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.058972999074934+-0.058041129871917*z)+(1.9602*sqrt(0.0004740117+(2*z*0.000201316088697)+((z^2)*0.00068449)))
y2 <- (-0.058972999074934+-0.058041129871917*z)-(1.9602*sqrt(0.0004740117+(2*z*0.000201316088697)+((z^2)*0.00068449)))
fy <- c(y1,y2)
fline <- (-0.058972999074934+-0.058041129871917*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='Moral',ylab='Simple Slope',main='Confidence Bands for White participants', 
     sub= "Model with imputed covariates")
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-6.8451,col=4,lty=2)
abline(v=-0.3276,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-8.187,-7.6254)   #  <-- change to alter plot dims
leg <- c(-2,-8.1168)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-7.7702,-7.774)
y2 <- c(-7.6978,-7.9337)
y3 <- c(-7.6254,-8.0934)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.058972999074934+-0.058041129871917*z)+(1.9602*sqrt(0.000474012+(2*z*0.000201316088697)+((z^2)*0.00068449)))
y2 <- (-0.058972999074934+-0.058041129871917*z)-(1.9602*sqrt(0.000474012+(2*z*0.000201316088697)+((z^2)*0.00068449)))
fy <- c(y1,y2)
fline <- (-0.058972999074934+-0.058041129871917*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-6.8451,col=4,lty=2)
abline(v=-0.3276,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-8.1595,-7.6713)   #  <-- change to alter plot dims
leg <- c(-2,-8.0985)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-7.7855,-7.728)
y2 <- c(-7.7284,-7.9031)
y3 <- c(-7.6713,-8.0781)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.043659139822721+-0.058041129871917*z)+(1.9602*sqrt(0.000583062+(2*z*0.000190666911056)+((z^2)*0.00068449)))
y2 <- (-0.043659139822721+-0.058041129871917*z)-(1.9602*sqrt(0.000583062+(2*z*0.000190666911056)+((z^2)*0.00068449)))
fy <- c(y1,y2)
fline <- (-0.043659139822721+-0.058041129871917*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-4.9681,col=4,lty=2)
abline(v=0.0911,col=4,lty=2)

Race interactions (All groups, AVG DI)

None

##### Black ######
#race  x instrumental

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.5711,1.7728)   #  <-- change to alter plot dims
leg <- c(-2,1.5963)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.6353,1.7728)
y2 <- c(1.6685,1.6888)
y3 <- c(1.7017,1.6047)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.005055161749362+-0.029315673437488*z)+(1.9602*sqrt(0.000189014+(2*z*-0.00000510755792)+((z^2)*0.000038691)))
y2 <- (0.005055161749362+-0.029315673437488*z)-(1.9602*sqrt(0.000189014+(2*z*-0.00000510755792)+((z^2)*0.000038691)))
fy <- c(y1,y2)
fline <- (0.005055161749362+-0.029315673437488*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-0.8284,col=4,lty=2)
abline(v=1.1901,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.5304,1.7932)   #  <-- change to alter plot dims
leg <- c(-2,1.5633)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.6658,1.6814)
y2 <- c(1.7295,1.6278)
y3 <- c(1.7932,1.5742)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.025428187574734+-0.029315673437488*z)+(1.9602*sqrt(0.000073651+(2*z*-0.000015677005096)+((z^2)*0.000038691)))
y2 <- (-0.025428187574734+-0.029315673437488*z)-(1.9602*sqrt(0.000073651+(2*z*-0.000015677005096)+((z^2)*0.000038691)))
fy <- c(y1,y2)
fline <- (-0.025428187574734+-0.029315673437488*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.013,col=4,lty=2)
abline(v=-0.2541,col=4,lty=2)

##### Asian ######
#race  x instrumental

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.4178,1.8476)   #  <-- change to alter plot dims
leg <- c(-2,1.4715)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.4894,1.5716)
y2 <- c(1.6685,1.6888)
y3 <- c(1.8476,1.8059)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.005055161749362+-0.015478922106229*z)+(1.9602*sqrt(0.000189014+(2*z*-0.000004221028783)+((z^2)*0.000040991)))
y2 <- (0.005055161749362+-0.015478922106229*z)-(1.9602*sqrt(0.000189014+(2*z*-0.000004221028783)+((z^2)*0.000040991)))
fy <- c(y1,y2)
fline <- (0.005055161749362+-0.015478922106229*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.2621,col=4,lty=2)
abline(v=3.7732,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.2155,2.0008)   #  <-- change to alter plot dims
leg <- c(-2,1.3136)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.3464,2.0008)
y2 <- c(1.3824,1.9749)
y3 <- c(1.4184,1.949)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.148134885736081+-0.015478922106229*z)+(1.9602*sqrt(0.000061017+(2*z*-0.000006783206549)+((z^2)*0.000040991)))
y2 <- (0.148134885736081+-0.015478922106229*z)-(1.9602*sqrt(0.000061017+(2*z*-0.000006783206549)+((z^2)*0.000040991)))
fy <- c(y1,y2)
fline <- (0.148134885736081+-0.015478922106229*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=5.2963,col=4,lty=2)
abline(v=49.9297,col=4,lty=2)

Pre-registered

# Instrumental X Moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.729,1.9981)   #  <-- change to alter plot dims
leg <- c(-2,1.7627)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8923,1.9981)
y2 <- c(1.944,1.886)
y3 <- c(1.9957,1.7739)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.014510455591654+-0.040953635508422*z)+(1.9602*sqrt(0.000502119+(2*z*0.000218455971617)+((z^2)*0.000856291)))
y2 <- (-0.014510455591654+-0.040953635508422*z)-(1.9602*sqrt(0.000502119+(2*z*0.000218455971617)+((z^2)*0.000856291)))
fy <- c(y1,y2)
fline <- (-0.014510455591654+-0.040953635508422*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.7013,2.0427)   #  <-- change to alter plot dims
leg <- c(-2,1.744)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.908,1.951)
y2 <- c(1.9754,1.8546)
y3 <- c(2.0427,1.7582)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.030181618660384+-0.040953635508422*z)+(1.9602*sqrt(0.000724321+(2*z*0.000261018377047)+((z^2)*0.000856291)))
y2 <- (-0.030181618660384+-0.040953635508422*z)-(1.9602*sqrt(0.000724321+(2*z*0.000261018377047)+((z^2)*0.000856291)))
fy <- c(y1,y2)
fline <- (-0.030181618660384+-0.040953635508422*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

##### Black ######
#race  x instrumental

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.8227,1.9637)   #  <-- change to alter plot dims
leg <- c(-2,1.8403)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.9243,1.9258)
y2 <- c(1.944,1.886)
y3 <- c(1.9637,1.8462)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.014510455591654+-0.01487020038375*z)+(1.9602*sqrt(0.000502119+(2*z*-0.000008422024334)+((z^2)*0.000079418)))
y2 <- (-0.014510455591654+-0.01487020038375*z)-(1.9602*sqrt(0.000502119+(2*z*-0.000008422024334)+((z^2)*0.000079418)))
fy <- c(y1,y2)
fline <- (-0.014510455591654+-0.01487020038375*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.8307,1.9503)   #  <-- change to alter plot dims
leg <- c(-2,1.8456)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.9199,1.9391)
y2 <- c(1.9351,1.8949)
y3 <- c(1.9503,1.8506)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.010064736697483+-0.01487020038375*z)+(1.9602*sqrt(0.000092455+(2*z*-0.000004590685773)+((z^2)*0.000079418)))
y2 <- (-0.010064736697483+-0.01487020038375*z)-(1.9602*sqrt(0.000092455+(2*z*-0.000004590685773)+((z^2)*0.000079418)))
fy <- c(y1,y2)
fline <- (-0.010064736697483+-0.01487020038375*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=1.0204,col=4,lty=2)
abline(v=2.9617,col=4,lty=2)

##### Latinx ######
#race  x instrumental

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.7722,1.9672)   #  <-- change to alter plot dims
leg <- c(-2,1.7966)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.933,1.8047)
y2 <- c(1.944,1.886)
y3 <- c(1.9551,1.9672)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.014510455591654+0.017545675758108*z)+(1.9602*sqrt(0.000502119+(2*z*-0.000008325556089)+((z^2)*0.000101167)))
y2 <- (-0.014510455591654+0.017545675758108*z)-(1.9602*sqrt(0.000502119+(2*z*-0.000008325556089)+((z^2)*0.000101167)))
fy <- c(y1,y2)
fline <- (-0.014510455591654+0.017545675758108*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.7222,2.0278)   #  <-- change to alter plot dims
leg <- c(-2,1.7604)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8723,1.9867)
y2 <- c(1.8227,2.0073)
y3 <- c(1.7731,2.0278)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.046139683103389+0.017545675758108*z)+(1.9602*sqrt(0.000103249+(2*z*-0.000001821465432)+((z^2)*0.000101167)))
y2 <- (0.046139683103389+0.017545675758108*z)-(1.9602*sqrt(0.000103249+(2*z*-0.000001821465432)+((z^2)*0.000101167)))
fy <- c(y1,y2)
fline <- (0.046139683103389+0.017545675758108*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.0101,col=4,lty=2)
abline(v=21.2049,col=4,lty=2)

All covariates

# Instrumental x Moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-14.3448,-13.7043)   #  <-- change to alter plot dims
leg <- c(-2,-14.2647)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-13.9071,-13.8535)
y2 <- c(-13.8057,-14.0458)
y3 <- c(-13.7043,-14.2381)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.060021234170912+-0.073412621896668*z)+(1.9602*sqrt(0.000689369+(2*z*0.000195446848919)+((z^2)*0.000838244)))
y2 <- (-0.060021234170912+-0.073412621896668*z)-(1.9602*sqrt(0.000689369+(2*z*0.000195446848919)+((z^2)*0.000838244)))
fy <- c(y1,y2)
fline <- (-0.060021234170912+-0.073412621896668*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-3.2352,col=4,lty=2)
abline(v=-0.1359,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-14.3186,-13.748)   #  <-- change to alter plot dims
leg <- c(-2,-14.2473)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-13.9217,-13.8098)
y2 <- c(-13.8349,-14.0166)
y3 <- c(-13.748,-14.2235)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.045447225078187+-0.073412621896668*z)+(1.9602*sqrt(0.000786797+(2*z*0.000200961564573)+((z^2)*0.000838244)))
y2 <- (-0.045447225078187+-0.073412621896668*z)-(1.9602*sqrt(0.000786797+(2*z*0.000200961564573)+((z^2)*0.000838244)))
fy <- c(y1,y2)
fline <- (-0.045447225078187+-0.073412621896668*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.5388,col=4,lty=2)
abline(v=0.1739,col=4,lty=2)

Imputed

# Instrumental x Moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-27.1728,-26.6377)   #  <-- change to alter plot dims
leg <- c(-2,-27.1059)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-26.767,-26.7839)
y2 <- c(-26.7024,-26.9338)
y3 <- c(-26.6377,-27.0836)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.057856383913575+-0.053619518718697*z)+(1.9602*sqrt(0.000441496+(2*z*0.000187427741087)+((z^2)*0.000637478)))
y2 <- (-0.057856383913575+-0.053619518718697*z)-(1.9602*sqrt(0.000441496+(2*z*0.000187427741087)+((z^2)*0.000637478)))
fy <- c(y1,y2)
fline <- (-0.057856383913575+-0.053619518718697*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-10.8348,col=4,lty=2)
abline(v=-0.358,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-27.1454,-26.6835)   #  <-- change to alter plot dims
leg <- c(-2,-27.0876)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-26.7822,-26.7382)
y2 <- c(-26.7328,-26.9033)
y3 <- c(-26.6835,-27.0684)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.042617507215769+-0.053619518718697*z)+(1.9602*sqrt(0.000542924+(2*z*0.000177663438461)+((z^2)*0.000637478)))
y2 <- (-0.042617507215769+-0.053619518718697*z)-(1.9602*sqrt(0.000542924+(2*z*0.000177663438461)+((z^2)*0.000637478)))
fy <- c(y1,y2)
fline <- (-0.042617507215769+-0.053619518718697*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-7.613,col=4,lty=2)
abline(v=0.0833,col=4,lty=2)

FIRST DI MODEL

No Covariates

PreRegistered

## Asian
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.7074,2.0111)   #  <-- change to alter plot dims
leg <- c(-2,1.7454)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8507,2.0111)
y2 <- c(1.9269,1.8846)
y3 <- c(2.0031,1.7581)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.010578734045118+-0.050681542268119*z)+(1.9602*sqrt(0.000329982+(2*z*0.000137100973491)+((z^2)*0.000580821)))
y2 <- (-0.010578734045118+-0.050681542268119*z)-(1.9602*sqrt(0.000329982+(2*z*0.000137100973491)+((z^2)*0.000580821)))
fy <- c(y1,y2)
fline <- (-0.010578734045118+-0.050681542268119*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.8803,col=4,lty=2)
abline(v=1.8248,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.6828,2.0469)   #  <-- change to alter plot dims
leg <- c(-2,1.7283)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.8653,1.9674)
y2 <- c(1.9561,1.8554)
y3 <- c(2.0469,1.7435)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.025165939795889+-0.050681542268119*z)+(1.9602*sqrt(0.000531089+(2*z*0.000147650033158)+((z^2)*0.000580821)))
y2 <- (-0.025165939795889+-0.050681542268119*z)-(1.9602*sqrt(0.000531089+(2*z*0.000147650033158)+((z^2)*0.000580821)))
fy <- c(y1,y2)
fline <- (-0.025165939795889+-0.050681542268119*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-5.0336,col=4,lty=2)
abline(v=0.8299,col=4,lty=2)

## Latinx
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(2.0499,2.3883)   #  <-- change to alter plot dims
leg <- c(-2,2.0922)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.1325,2.3883)
y2 <- c(2.2288,2.2473)
y3 <- c(2.3251,2.1063)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.004631015218456+-0.059317511639751*z)+(1.9602*sqrt(0.000552654+(2*z*0.000219925000231)+((z^2)*0.000966041)))
y2 <- (0.004631015218456+-0.059317511639751*z)-(1.9602*sqrt(0.000552654+(2*z*0.000219925000231)+((z^2)*0.000966041)))
fy <- c(y1,y2)
fline <- (0.004631015218456+-0.059317511639751*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-10.5556,col=4,lty=2)
abline(v=-1.0303,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(2.014,2.406)   #  <-- change to alter plot dims
leg <- c(-2,2.063)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.1595,2.3073)
y2 <- c(2.2828,2.1933)
y3 <- c(2.406,2.0793)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.022356900910759+-0.059317511639751*z)+(1.9602*sqrt(0.000826396+(2*z*0.000274679331073)+((z^2)*0.000966041)))
y2 <- (-0.022356900910759+-0.059317511639751*z)-(1.9602*sqrt(0.000826396+(2*z*0.000274679331073)+((z^2)*0.000966041)))
fy <- c(y1,y2)
fline <- (-0.022356900910759+-0.059317511639751*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

All Covariates

## Asian
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-151.7002,-151.2743)   #  <-- change to alter plot dims
leg <- c(-2,-151.647)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-151.4534,-151.3428)
y2 <- c(-151.3639,-151.486)
y3 <- c(-151.2743,-151.6292)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.030540659468736+-0.058192269705877*z)+(1.9602*sqrt(0.00029792+(2*z*0.000123849312208)+((z^2)*0.000325076)))
y2 <- (-0.030540659468736+-0.058192269705877*z)-(1.9602*sqrt(0.00029792+(2*z*0.000123849312208)+((z^2)*0.000325076)))
fy <- c(y1,y2)
fline <- (-0.030540659468736+-0.058192269705877*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.2944,col=4,lty=2)
abline(v=0.0766,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-151.6935,-151.2855)   #  <-- change to alter plot dims
leg <- c(-2,-151.6425)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-151.4572,-151.3317)
y2 <- c(-151.3713,-151.4786)
y3 <- c(-151.2855,-151.6255)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.026814981348789+-0.058192269705877*z)+(1.9602*sqrt(0.000338738+(2*z*0.000004713790553)+((z^2)*0.000325076)))
y2 <- (-0.026814981348789+-0.058192269705877*z)-(1.9602*sqrt(0.000338738+(2*z*0.000004713790553)+((z^2)*0.000325076)))
fy <- c(y1,y2)
fline <- (-0.026814981348789+-0.058192269705877*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.6123,col=4,lty=2)
abline(v=0.169,col=4,lty=2)

## Latinx
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-31.2913,-30.671)   #  <-- change to alter plot dims
leg <- c(-2,-31.2138)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-30.9734,-30.8447)
y2 <- c(-30.8222,-31.0163)
y3 <- c(-30.671,-31.1879)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.048527228529286+-0.080703973593617*z)+(1.9602*sqrt(0.000780555+(2*z*0.000197403798463)+((z^2)*0.000963839)))
y2 <- (-0.048527228529286+-0.080703973593617*z)-(1.9602*sqrt(0.000780555+(2*z*0.000197403798463)+((z^2)*0.000963839)))
fy <- c(y1,y2)
fline <- (-0.048527228529286+-0.080703973593617*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-2.3456,col=4,lty=2)
abline(v=0.0978,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-31.2336,-30.7297)   #  <-- change to alter plot dims
leg <- c(-2,-31.1706)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-31.0117,-30.7297)
y2 <- c(-30.8989,-30.9397)
y3 <- c(-30.786,-31.1496)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.010198580296869+-0.080703973593617*z)+(1.9602*sqrt(0.00088754+(2*z*0.00018354223071)+((z^2)*0.000963839)))
y2 <- (-0.010198580296869+-0.080703973593617*z)-(1.9602*sqrt(0.00088754+(2*z*0.00018354223071)+((z^2)*0.000963839)))
fy <- c(y1,y2)
fline <- (-0.010198580296869+-0.080703973593617*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.1275,col=4,lty=2)
abline(v=1.0436,col=4,lty=2)

## White
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-46.6352,-46.1172)   #  <-- change to alter plot dims
leg <- c(-2,-46.5704)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-46.2979,-46.2003)
y2 <- c(-46.2076,-46.3746)
y3 <- c(-46.1172,-46.5488)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.041750710087427+-0.066164670059188*z)+(1.9602*sqrt(0.000334774+(2*z*0.00009239235514)+((z^2)*0.000416907)))
y2 <- (-0.041750710087427+-0.066164670059188*z)-(1.9602*sqrt(0.000334774+(2*z*0.00009239235514)+((z^2)*0.000416907)))
fy <- c(y1,y2)
fline <- (-0.041750710087427+-0.066164670059188*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.6338,col=4,lty=2)
abline(v=-0.1007,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-46.6355,-46.1166)   #  <-- change to alter plot dims
leg <- c(-2,-46.5706)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-46.2978,-46.2008)
y2 <- c(-46.2072,-46.3749)
y3 <- c(-46.1166,-46.549)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.041935666574239+-0.066164670059188*z)+(1.9602*sqrt(0.000390048+(2*z*0.000100266489861)+((z^2)*0.000416907)))
y2 <- (-0.041935666574239+-0.066164670059188*z)-(1.9602*sqrt(0.000390048+(2*z*0.000100266489861)+((z^2)*0.000416907)))
fy <- c(y1,y2)
fline <- (-0.041935666574239+-0.066164670059188*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.6653,col=4,lty=2)
abline(v=-0.0562,col=4,lty=2)

Imputed

## Asian
# instrumental x moral
xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-107.8735,-107.4712)   #  <-- change to alter plot dims
leg <- c(-2,-107.8232)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-107.5761,-107.5489)
y2 <- c(-107.5237,-107.6777)
y3 <- c(-107.4712,-107.8065)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.03850314898013+-0.045325462704016*z)+(1.9602*sqrt(0.000270843+(2*z*0.000115522766568)+((z^2)*0.000397895)))
y2 <- (-0.03850314898013+-0.045325462704016*z)-(1.9602*sqrt(0.000270843+(2*z*0.000115522766568)+((z^2)*0.000397895)))
fy <- c(y1,y2)
fline <- (-0.03850314898013+-0.045325462704016*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-4.7761,col=4,lty=2)
abline(v=-0.176,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-107.8729,-107.4722)   #  <-- change to alter plot dims
leg <- c(-2,-107.8228)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-107.5765,-107.5478)
y2 <- c(-107.5243,-107.677)
y3 <- c(-107.4722,-107.8061)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.03816084419195+-0.045325462704016*z)+(1.9602*sqrt(0.000366483+(2*z*0.00008621544534)+((z^2)*0.000397895)))
y2 <- (-0.03816084419195+-0.045325462704016*z)-(1.9602*sqrt(0.000366483+(2*z*0.00008621544534)+((z^2)*0.000397895)))
fy <- c(y1,y2)
fline <- (-0.03816084419195+-0.045325462704016*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-5.3044,col=4,lty=2)
abline(v=-0.0173,col=4,lty=2)

## Latinx
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-122.0914,-121.6706)   #  <-- change to alter plot dims
leg <- c(-2,-122.0388)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-121.8207,-121.7328)
y2 <- c(-121.7456,-121.877)
y3 <- c(-121.6706,-122.0213)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.032851471386859+-0.054823449386367*z)+(1.9602*sqrt(0.000494294+(2*z*0.00020247550685)+((z^2)*0.000750575)))
y2 <- (-0.032851471386859+-0.054823449386367*z)-(1.9602*sqrt(0.000494294+(2*z*0.00020247550685)+((z^2)*0.000750575)))
fy <- c(y1,y2)
fline <- (-0.032851471386859+-0.054823449386367*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-17.212,col=4,lty=2)
abline(v=0.3917,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-122.0946,-121.6653)   #  <-- change to alter plot dims
leg <- c(-2,-122.0409)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-121.8189,-121.738)
y2 <- c(-121.7421,-121.8805)
y3 <- c(-121.6653,-122.023)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.034608084724527+-0.054823449386367*z)+(1.9602*sqrt(0.000632073+(2*z*0.000191027798077)+((z^2)*0.000750575)))
y2 <- (-0.034608084724527+-0.054823449386367*z)-(1.9602*sqrt(0.000632073+(2*z*0.000191027798077)+((z^2)*0.000750575)))
fy <- c(y1,y2)
fline <- (-0.034608084724527+-0.054823449386367*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-19.642,col=4,lty=2)
abline(v=0.5152,col=4,lty=2)

## White
# Instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-37.2848,-36.7903)   #  <-- change to alter plot dims
leg <- c(-2,-37.223)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-36.9203,-36.9007)
y2 <- c(-36.8553,-37.0516)
y3 <- c(-36.7903,-37.2024)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.049071167461488+-0.053962314811807*z)+(1.9602*sqrt(0.000333883+(2*z*0.000141238135924)+((z^2)*0.00048541)))
y2 <- (-0.049071167461488+-0.053962314811807*z)-(1.9602*sqrt(0.000333883+(2*z*0.000141238135924)+((z^2)*0.00048541)))
fy <- c(y1,y2)
fline <- (-0.049071167461488+-0.053962314811807*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-3.7345,col=4,lty=2)
abline(v=-0.2878,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-37.2737,-36.8088)   #  <-- change to alter plot dims
leg <- c(-2,-37.2156)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-36.9265,-36.8822)
y2 <- c(-36.8676,-37.0392)
y3 <- c(-36.8088,-37.1962)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.042904836455703+-0.053962314811807*z)+(1.9602*sqrt(0.000412528+(2*z*0.000135466197238)+((z^2)*0.00048541)))
y2 <- (-0.042904836455703+-0.053962314811807*z)-(1.9602*sqrt(0.000412528+(2*z*0.000135466197238)+((z^2)*0.00048541)))
fy <- c(y1,y2)
fline <- (-0.042904836455703+-0.053962314811807*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-3.3561,col=4,lty=2)
abline(v=-0.0728,col=4,lty=2)

Race Interactions (All groups, First DI)

None

## race x instrumental
# Black

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.6351,1.8657)   #  <-- change to alter plot dims
leg <- c(-2,1.6639)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.6735,1.8657)
y2 <- c(1.7235,1.7856)
y3 <- c(1.7734,1.7056)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.015545733111374+-0.032484104535714*z)+(1.9602*sqrt(0.000163151+(2*z*-0.000005534437564)+((z^2)*0.000041866)))
y2 <- (0.015545733111374+-0.032484104535714*z)-(1.9602*sqrt(0.000163151+(2*z*-0.000005534437564)+((z^2)*0.000041866)))
fy <- c(y1,y2)
fline <- (0.015545733111374+-0.032484104535714*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-0.3096,col=4,lty=2)
abline(v=1.3913,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.6351,1.8657)   #  <-- change to alter plot dims
leg <- c(-2,1.6639)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.6735,1.8657)
y2 <- c(1.7235,1.7856)
y3 <- c(1.7734,1.7056)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.015545733111374+-0.032484104535714*z)+(1.9602*sqrt(0.000163151+(2*z*-0.000005534437564)+((z^2)*0.000041866)))
y2 <- (0.015545733111374+-0.032484104535714*z)-(1.9602*sqrt(0.000163151+(2*z*-0.000005534437564)+((z^2)*0.000041866)))
fy <- c(y1,y2)
fline <- (0.015545733111374+-0.032484104535714*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-0.3096,col=4,lty=2)
abline(v=1.3913,col=4,lty=2)

## race x instrumental
# Asian

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.4677,1.9119)   #  <-- change to alter plot dims
leg <- c(-2,1.5233)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.5418,1.6594)
y2 <- c(1.7235,1.7856)
y3 <- c(1.9051,1.9119)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.015545733111374+-0.013850130449993*z)+(1.9602*sqrt(0.000163151+(2*z*-0.000004577295524)+((z^2)*0.000044366)))
y2 <- (0.015545733111374+-0.013850130449993*z)-(1.9602*sqrt(0.000163151+(2*z*-0.000004577295524)+((z^2)*0.000044366)))
fy <- c(y1,y2)
fline <- (0.015545733111374+-0.013850130449993*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-0.9276,col=4,lty=2)
abline(v=19.4438,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.2691,2.0747)   #  <-- change to alter plot dims
leg <- c(-2,1.3698)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(1.4033,2.0747)
y2 <- c(1.4466,2.0625)
y3 <- c(1.4898,2.0504)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.153986276181411+-0.013850130449993*z)+(1.9602*sqrt(0.00006604+(2*z*-0.000007338893081)+((z^2)*0.000044366)))
y2 <- (0.153986276181411+-0.013850130449993*z)-(1.9602*sqrt(0.00006604+(2*z*-0.000007338893081)+((z^2)*0.000044366)))
fy <- c(y1,y2)
fline <- (0.153986276181411+-0.013850130449993*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=5.7404,col=4,lty=2)
abline(v=191.3429,col=4,lty=2)

Pre-registered

## instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.8973,2.1488)   #  <-- change to alter plot dims
leg <- c(-2,1.9288)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.0101,2.1488)
y2 <- c(2.0603,2.044)
y3 <- c(2.1106,1.9392)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.004073626634868+-0.038757740572593*z)+(1.9602*sqrt(0.000452594+(2*z*0.000196314976078)+((z^2)*0.000772783)))
y2 <- (-0.004073626634868+-0.038757740572593*z)-(1.9602*sqrt(0.000452594+(2*z*0.000196314976078)+((z^2)*0.000772783)))
fy <- c(y1,y2)
fline <- (-0.004073626634868+-0.038757740572593*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.8632,2.1801)   #  <-- change to alter plot dims
leg <- c(-2,1.9028)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.0333,2.0792)
y2 <- c(2.1067,1.9976)
y3 <- c(2.1801,1.916)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.027264827481771+-0.038757740572593*z)+(1.9602*sqrt(0.000653566+(2*z*0.000235903586241)+((z^2)*0.000772783)))
y2 <- (-0.027264827481771+-0.038757740572593*z)-(1.9602*sqrt(0.000653566+(2*z*0.000235903586241)+((z^2)*0.000772783)))
fy <- c(y1,y2)
fline <- (-0.027264827481771+-0.038757740572593*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

## Race x instrumental
# Black

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(1.9948,2.0947)   #  <-- change to alter plot dims
leg <- c(-2,2.0073)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.026,2.0766)
y2 <- c(2.0603,2.044)
y3 <- c(2.0947,2.0115)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.004073626634868+-0.016726663990555*z)+(1.9602*sqrt(0.000452594+(2*z*-0.000009114874613)+((z^2)*0.000085955)))
y2 <- (-0.004073626634868+-0.016726663990555*z)-(1.9602*sqrt(0.000452594+(2*z*-0.000009114874613)+((z^2)*0.000085955)))
fy <- c(y1,y2)
fline <- (-0.004073626634868+-0.016726663990555*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(2.0014,2.0915)   #  <-- change to alter plot dims
leg <- c(-2,2.0127)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(2.021,2.0915)
y2 <- c(2.0504,2.054)
y3 <- c(2.0797,2.0164)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (0.000897062379384+-0.016726663990555*z)+(1.9602*sqrt(0.000100063+(2*z*-0.000004969481028)+((z^2)*0.000085955)))
y2 <- (0.000897062379384+-0.016726663990555*z)-(1.9602*sqrt(0.000100063+(2*z*-0.000004969481028)+((z^2)*0.000085955)))
fy <- c(y1,y2)
fline <- (0.000897062379384+-0.016726663990555*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=NaN,col=4,lty=2)
abline(v=NaN,col=4,lty=2)

All

## instrumental x moral
xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-54.7559,-54.2582)   #  <-- change to alter plot dims
leg <- c(-2,-54.6937)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-54.4252,-54.3494)
y2 <- c(-54.3417,-54.5112)
y3 <- c(-54.2582,-54.6729)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.042362829367119+-0.061315286492093*z)+(1.9602*sqrt(0.000297296+(2*z*0.000082068306812)+((z^2)*0.000370806)))
y2 <- (-0.042362829367119+-0.061315286492093*z)-(1.9602*sqrt(0.000297296+(2*z*0.000082068306812)+((z^2)*0.000370806)))
fy <- c(y1,y2)
fline <- (-0.042362829367119+-0.061315286492093*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.7997,col=4,lty=2)
abline(v=-0.1552,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-54.7501,-54.2679)   #  <-- change to alter plot dims
leg <- c(-2,-54.6898)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-54.4285,-54.3397)
y2 <- c(-54.3482,-54.5047)
y3 <- c(-54.2679,-54.6697)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.039126635301705+-0.061315286492093*z)+(1.9602*sqrt(0.000346079+(2*z*0.000089521095584)+((z^2)*0.000370806)))
y2 <- (-0.039126635301705+-0.061315286492093*z)-(1.9602*sqrt(0.000346079+(2*z*0.000089521095584)+((z^2)*0.000370806)))
fy <- c(y1,y2)
fline <- (-0.039126635301705+-0.061315286492093*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-1.71,col=4,lty=2)
abline(v=-0.0504,col=4,lty=2)

Imputed

## instrumental x moral

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-55.6511,-55.1887)   #  <-- change to alter plot dims
leg <- c(-2,-55.5933)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-55.2992,-55.2942)
y2 <- c(-55.244,-55.4341)
y3 <- c(-55.1887,-55.574)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.047539139982592+-0.048789250534861*z)+(1.9602*sqrt(0.00030189+(2*z*0.000127559118657)+((z^2)*0.000439141)))
y2 <- (-0.047539139982592+-0.048789250534861*z)-(1.9602*sqrt(0.00030189+(2*z*0.000127559118657)+((z^2)*0.000439141)))
fy <- c(y1,y2)
fline <- (-0.047539139982592+-0.048789250534861*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-4.9587,col=4,lty=2)
abline(v=-0.3201,col=4,lty=2)

xx <- c(-2,2)   #  <-- change to alter plot dims
yy <- c(-55.6416,-55.2044)   #  <-- change to alter plot dims
leg <- c(-2,-55.587)   #  <-- change to alter legend location
x <- c(-2,2)   #  <-- x-coords for lines
y1 <- c(-55.3045,-55.2785)
y2 <- c(-55.2544,-55.4236)
y3 <- c(-55.2044,-55.5688)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='x1',ylab='Y',main='HLM 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('W2(1)','W2(2)','W2(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

z1=-2  #supply lower bound for w2 here
z2=2   #supply upper bound for w2 here
z <- seq(z1,z2,length=1000)
fz <- c(z,z)
y1 <- (-0.042300600511189+-0.048789250534861*z)+(1.9602*sqrt(0.000373031+(2*z*0.000122691457506)+((z^2)*0.000439141)))
y2 <- (-0.042300600511189+-0.048789250534861*z)-(1.9602*sqrt(0.000373031+(2*z*0.000122691457506)+((z^2)*0.000439141)))
fy <- c(y1,y2)
fline <- (-0.042300600511189+-0.048789250534861*z)
plot(fz,fy,type='p',pch='.',font=2,font.lab=2,col=2,xlab='w2',ylab='Simple Slope',main='Confidence Bands')
lines(z,fline)
f0 <- array(0,c(1000))
lines(z,f0,col=8)
abline(v=-4.4806,col=4,lty=2)
abline(v=-0.1147,col=4,lty=2)