Recovering Marginal Effects and Standard Errors from Interaction Terms in R

[This article was first published on Journal - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

When I fit models with interactions, I often want to recover not only the interaction effect but also the marginal effect (the main effect + the interaction) and of course the standard errors. There are a couple of ways to do this in R but I ended writing my own function (essentially a wrapper around the deltaMethod() function) to fit my needs.

In this case, I have a model where a continuous variable has been interacted with a discrete variable. I want to create a data.frame that stores the main effect, the marginal effects of the interactions, and their standard errors.

I also want the standard errors for the marginal effects on the interaction to match the standard errors on the main effect which in this case were HAC (Heteroskedastic and Autocorrelation Consistent ie. Newey-West).

The function below is a little sloppy as I built it iteratively while I was working through a specific problem, but it is generalized to any lm() type object. The basic use case is that you have 1 or mor models, each with 1 or more interactions and the discrete terms in the interactions have more than two levels (it will work with two levels, but in that case the function might be overkill). In my case I fit a series of regions specific panel models where two of the main effects were interacted with a year term. I wanted have a generalizable function where I could just feed it the model object, variable name, and covariance estimator, and return a data.frame with estimates and standard errors. In the next post I’ll demonstrate the function with simulated data and show how to visualize the marginal in ggplot.

 Note: The code below assumes you have already loaded car (for deltaMethod()), sandwich (for vcovHAC), and lmtest (for coeftest()).

 

##Calls this function which is a wrapper for coeftest from the sandwhich package
sehac<-function(fit,vcov=vcovHAC){ #Convenience function for HAC standard erros
	coeftest(fit,vcov)
}
 
funinteff<-function(mod,var,vcov=vcovHAC){ #mod is an lm() object, var is the name of the main effect that was interacted, vcov is the type of variance covariance method you want to use 
	#Extract Coefficient names create 'beta names' to feed to deltaMethod()
	cnames<-coef(mod)
	pnams<-data.frame('b'=paste('b',0:(length(cnames)-1),sep=""),'est'=cnames) #assign parameter names so that deltaMethod does not throw an error
 
	#Extract the specific parameters of interest
	vars<-grep(var,names(cnames),value=T)
	var1<-vars[1]
	intvars<-vars[2:length(vars)]
	bi<-pnams[var1,'b']
 
	#--Create Data Frame to store Main Effect
	int<-sehac(mod,vcov=vcov)[var1,c('Estimate','Std. Error')]
	int<-as.data.frame(t(int))
	names(int)<-c('Estimate','SE')
	row.names(int)<-var1
 
	#Loop through and store the results in a data.frame
	for(i in 1:length(intvars)){
		bint<-pnams[intvars[i],'b']
		eq<-paste(bi,bint,sep="+")
		interac<-deltaMethod(mod,eq,parameterNames=pnams[,1],vcov=vcov)
		row.names(interac)<-intvars[i]
 
		int<-rbind(int,interac)
	}
	return(int)
}

 

To leave a comment for the author, please follow the link and comment on their blog: Journal - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)