# A Slightly Less Simple Function for Forest Plots

This post is a sequel to my wildly popular – ok, not wildly popular, but very useful for me – post documenting my function to create a forest plot.  Recently, I made models using cross-validation via R’s caret package.  The `train()` function does the heavy lifting, but the problem is that its output is a of the `train` class, which I presume was created for `caret`.  This is a problem because my function assumes it is receiving a generic model, so it fails.

The new function, pasted at the end of this post, makes two modifications.  The main modification is to take the coefficients from `model\$finalModel\$coefficients`.  You of course could pass those directly to the `model` argument in `forestPlot`, but that would require more typing, and I am lazy.  The second update, which is not because of `train`, it to add `na.rm=TRUE` to `min()` and `max()`.  Those arguments should have been in the original function, oh well.  The third change is simply to change the x-axis label, which is hard coded.  That should probably be an argument to pass to the function; oh well.

In response to demand from Alex Hughes, I have included an example of the output.  Initial results suggest that violence, small groups, children, and seeing more individuals protest correlates with larger protests in the future, while the presence of police and large groups correlates with smaller protests. ```'''
This formula has been updated to handle results of train class, from caret package.
The purpose of this script is to define an R function that makes pretty forest plots based on regression results.
Function arguments:
1. model = model results from caret
2. coefficients = character vector of coefficients to keep. Easiest to pass as names(model\$coefficients)
3. coef_labels = labels for y-axis on plot
4. se = user can pass custom standard errors
5. zscore = how user defines width of confidence intervals
6. outfile = name, including filepath, where the output will be saved. Must end in .pdf
7. cex = magnification for text

NB: This plotting works best when the coefficients have similar values, i.e. consider standardizing your variables.
'''

forestPlot <- function(model=NULL, coefficients=NULL, coef_labels=NULL, se = NULL, zscore = 1.96, outpath, cex){
coef <- model\$finalModel\$coefficients # Get coefficients from model

if(is.null(coefficients) == TRUE){ # Generate coefficient names
coefficients <- names(model\$finalModel\$coefficients)
}

coef <- coef[coefficients] # Keep user specific coefficients
coef <- as.numeric(coef)

if(is.null(se) == TRUE){ # If no standard errors given, take from model
stdev <- summary(model)\$coefficients[,2] # Get coefficient standard deviations
stdev <- stdev[coefficients]
stdev <- as.numeric(stdev)
}

if(is.null(se) == FALSE){ # If standard errors given, use those
stdev <- sd
stdev <- stdev[coefficients]
stdev <- as.numeric(stdev)
}

lower <- coef - zscore*stdev # Lower value of confidence interval
upper <- coef + zscore*stdev # Upper value of confidence interval

vars <- length(coef) # Will be used to create y values

if(is.null(coef_labels) == TRUE){
coef_labels <- coefficients
}

# Define x minimum to plot y-axis labels
label_pos <- min(lower, na.rm=TRUE)*-1.05

#Plot all results
pdf(outpath)
par(mar=c(7,10,3,1)) # Need more space on left side of plot (6) for var. labels
plot(x=coef, y=vars:1, xlim=c(min(lower, na.rm=TRUE)*1.1,max(upper, na.rm=TRUE)*1.5), pch = 20, xlab='Standardized Coefficient', bty='n', ylab='', yaxt='n', xaxt='n', cex.lab = cex)
for(i in 1:vars){
lines(x=c(lower[i], upper[i]), y = rep(vars+1-i, each=2), lwd=2)
}
axis(1, cex.axis=cex)
axis(2, at=vars:1, labels=coef_labels, las=1, lwd=0, pos=-label_pos, outer=TRUE, cex.axis=cex) # pos is x-axis location of labels
abline(v=0, lty=2)
dev.off()
}```

This site uses Akismet to reduce spam. Learn how your comment data is processed.
{"ft_lang_label":"__label__en","ft_lang_prob":0.79262984,"math_prob":0.90602773,"size":3729,"snap":"2023-14-2023-23","text_gpt3_token_len":943,"char_repetition_ratio":0.12778524,"word_repetition_ratio":0.010714286,"special_character_ratio":0.256369,"punctuation_ratio":0.1556157,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.997397,"pos_list":[0,1,2],"im_url_duplicate_count":[null,5,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2023-06-09T23:46:26Z\",\"WARC-Record-ID\":\"<urn:uuid:a79cf675-0ce1-4d4c-ba1a-cb3553b083a5>\",\"Content-Length\":\"93884\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:c5af13c4-6564-40d0-aa81-019c463bf8bf>\",\"WARC-Concurrent-To\":\"<urn:uuid:603cd48d-9caa-443f-a8e2-84b99559fa41>\",\"WARC-IP-Address\":\"192.0.78.24\",\"WARC-Target-URI\":\"https://zacharyst.com/2018/05/04/a-slightly-less-simple-function-for-forest-plots/\",\"WARC-Payload-Digest\":\"sha1:GARDBJO6A6G2KAOH2KGH5Q7ANHEUZQ4E\",\"WARC-Block-Digest\":\"sha1:QJ4CME53BMMVXDB6OS33WGRSRXYNSRV6\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2023/CC-MAIN-2023-23/CC-MAIN-2023-23_segments_1685224656869.87_warc_CC-MAIN-20230609233952-20230610023952-00720.warc.gz\"}"}