This is explained in Ch. 18 of DBDA2E, and an example is given in Figure 18.9 on p. 529 that plots the predicted slopes and their 95% HDIs as a function of the value of the other predictor. The variable "Spend" is \(x_1\) and the variable "PrcntTake" is \(x_2\). The top panel (below) plots the posterior predicted slope on \(x_1\) for values of \(x_2\):

After Figure 18.9, p. 529, of DBDA2E |

**A reader asked me how to compute those predicted slopes in R and make the graph.**I thought I already had the R code in the package of programs, but upon searching discovered that it was not there. (Sorry; with over 700 pages of text and many dozens of programs, some things get overlooked.)

**This blog post supplies the R code! :-)**

First, from the folder of programs, open the file

**Jags-Ymet-XmetMulti-Mrobust-Example.R**Then, at the top of the file, comment out all the data-loading examples and put in the following code for reading in the relevant data and creating an interaction variable. This is explained on p. 528 of DBDA2E:

**# Two predictors with interaction:**

# Read the SAT data file:

# Read the SAT data file:

**myData = read.csv( file="Guber1999data.csv" )**

# Append named column with interaction product:

# Append named column with interaction product:

**myData = cbind( myData , SpendXPrcnt=myData[,"Spend"]*myData[,"PrcntTake"] )**

# Specify the column names:

# Specify the column names:

**yName = "SATT" ; xName = c("Spend","PrcntTake","SpendXPrcnt")**

# Put the saved files into a sub-directory:

# Put the saved files into a sub-directory:

**dirName = "Guber1999data-Jags-Inter"**

if( !dir.exists(dirName) ) { dir.create(dirName) }

fileNameRoot = paste0(dirName,"/Guber1999data-Jags-Inter-")

# Specify the number of saved steps and the thinning (merely to save memory):

if( !dir.exists(dirName) ) { dir.create(dirName) }

fileNameRoot = paste0(dirName,"/Guber1999data-Jags-Inter-")

# Specify the number of saved steps and the thinning (merely to save memory):

**numSavedSteps=15000 ; thinSteps=20 # smaller thinSteps will run faster**

**# Specify the preferred graph file format:**

**graphFileType = "png"**

Then run the MCMC, as already specified in the script. I'm repeating it here for completeness:

**# Load the relevant model into R's working memory:**

**source("Jags-Ymet-XmetMulti-Mrobust.R")**

# Generate the MCMC chain:

# Generate the MCMC chain:

**mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName ,**

numSavedSteps=numSavedSteps , thinSteps=thinSteps ,

saveName=fileNameRoot )

numSavedSteps=numSavedSteps , thinSteps=thinSteps ,

saveName=fileNameRoot )

The script also has commands for displaying MCMC diagnostics, numerical summary, and basic plots of the posterior distribution.

Now the novel R code for computing the slopes and making the graph! At the end of the script, run the following code:

**# Convert mcmc coda object to a matrix:**

**mcmcMat=as.matrix(mcmcCoda,chains=TRUE)**

# Make combs of X1 and X2 values for displaying slopes and their HDIs:

# Make combs of X1 and X2 values for displaying slopes and their HDIs:

**X1Comb = seq(min(myData[,xName[1]]),max(myData[,xName[1]]),length=20)**

X2Comb = seq(min(myData[,xName[2]]),max(myData[,xName[2]]),length=20)

# Create place holders for X1 median slopes and HDIs:

X2Comb = seq(min(myData[,xName[2]]),max(myData[,xName[2]]),length=20)

# Create place holders for X1 median slopes and HDIs:

**X1SlopeMed = rep(0,length(X2Comb))**

X1SlopeHDI = matrix(0,nrow=2,ncol=length(X2Comb))

# Compute the slope on X1 for each value of X2,

# as explained at top of p. 530 of DBDA2E:

X1SlopeHDI = matrix(0,nrow=2,ncol=length(X2Comb))

# Compute the slope on X1 for each value of X2,

# as explained at top of p. 530 of DBDA2E:

**for ( combIdx in 1:length(X2Comb) ) {**

credSlope = mcmcMat[,"beta[1]"] + mcmcMat[,"beta[3]"]*X2Comb[combIdx]

X1SlopeHDI[,combIdx] = HDIofMCMC( credSlope )

X1SlopeMed[combIdx] = median( credSlope )

}

# Same as above but for X2:

credSlope = mcmcMat[,"beta[1]"] + mcmcMat[,"beta[3]"]*X2Comb[combIdx]

X1SlopeHDI[,combIdx] = HDIofMCMC( credSlope )

X1SlopeMed[combIdx] = median( credSlope )

}

# Same as above but for X2:

**X2SlopeHDI = matrix(0,nrow=2,ncol=length(X1Comb))**

X2SlopeMed = rep(0,length(X1Comb))

for ( combIdx in 1:length(X1Comb) ) {

credSlope = mcmcMat[,"beta[2]"] + mcmcMat[,"beta[3]"]*X1Comb[combIdx]

X2SlopeHDI[,combIdx] = HDIofMCMC( credSlope )

X2SlopeMed[combIdx] = median( credSlope )

}

# Make the graph as in lower part of Figure 18.9, p. 529 of DBDA2E:

X2SlopeMed = rep(0,length(X1Comb))

for ( combIdx in 1:length(X1Comb) ) {

credSlope = mcmcMat[,"beta[2]"] + mcmcMat[,"beta[3]"]*X1Comb[combIdx]

X2SlopeHDI[,combIdx] = HDIofMCMC( credSlope )

X2SlopeMed[combIdx] = median( credSlope )

}

# Make the graph as in lower part of Figure 18.9, p. 529 of DBDA2E:

**openGraph()**

layout(matrix(1:2,nrow=2))

par( mar=c(3.5,3.5,2,1) , mgp=c(2,0.7,0) )

# Panel to display slope on X1:

layout(matrix(1:2,nrow=2))

par( mar=c(3.5,3.5,2,1) , mgp=c(2,0.7,0) )

# Panel to display slope on X1:

**plot( X2Comb , X1SlopeMed , ylim=range(X1SlopeHDI) ,**

**pch=19 , cex=1.5 , col="skyblue" ,**

ylab=paste("Slope on",xName[1]) , xlab=paste("Value of",xName[2]) ,

ylab=paste("Slope on",xName[1]) , xlab=paste("Value of",xName[2]) ,

**cex.lab=1.5 )**

segments( x0=X2Comb , y0=X1SlopeHDI[1,] ,x1=X2Comb , y1=X1SlopeHDI[2,] ,

col="skyblue" , lwd=3 )

abline(h=0,lty="dotted")

# Panel to display slope on X2:

segments( x0=X2Comb , y0=X1SlopeHDI[1,] ,x1=X2Comb , y1=X1SlopeHDI[2,] ,

col="skyblue" , lwd=3 )

abline(h=0,lty="dotted")

# Panel to display slope on X2:

**plot( X1Comb , X2SlopeMed , ylim=range(X2SlopeHDI) , pch=19 , cex=1.5 ,**

**col="skyblue" ,**

ylab=paste("Slope on",xName[2]) , xlab=paste("Value of",xName[1]) ,

ylab=paste("Slope on",xName[2]) , xlab=paste("Value of",xName[1]) ,

**cex.lab=1.5 )**

segments( x0=X1Comb , y0=X2SlopeHDI[1,] ,x1=X1Comb , y1=X2SlopeHDI[2,] ,

col="skyblue" , lwd=3 )

abline(h=0,lty="dotted")

# Save it:

segments( x0=X1Comb , y0=X2SlopeHDI[1,] ,x1=X1Comb , y1=X2SlopeHDI[2,] ,

col="skyblue" , lwd=3 )

abline(h=0,lty="dotted")

# Save it:

**saveGraph( file=paste0(fileNameRoot,"Slopes") , type=graphFileType )**

Hopefully the comments in the script (above) are sufficient for you to decipher each step.

## No comments:

## Post a Comment