Section 11.2 Fitting a non-spatial regression model with R
The standard approach traditionally was to start with a non-spatial linear regression model and then to test whether or not this baseline model needs to be extended with spatial effects [197]. Given that we need to check autocorrelation in residuals, part of the workflow for modelling through regression spatial data will involve fitting a non-spatial regression model. In order to fit the model we use the
stats::lm() function using the formula specification (Y ~ X). Typically you want to store your regression model in a "variable", let’s call it "fit1_90". Here we model homicide rate in the 1990s (HR90) using, to start with, just two variables of interest: RD90 — Resource Deprivation/Affluence, and SOUTH — Counties in the southern region scored 1.
ncovr$SOUTH_f <- as.factor(ncovr$SOUTH) # specify this is a categorical variable
fit1_90 <- lm(HR90 ~ RD90 + SOUTH_f, data=ncovr)
You will see in your R Studio global environment space that there is a new object called
fit1_90 with 13 elements on it. We can get a sense for what this object is and includes using the functions we introduced in previous weeks:
class(fit1_90)
attributes(fit1_90)
R is telling us that this is an object of class
lm and that it includes a number of attributes. One of the beauties of R is that you are producing all the results from running the model, putting them in an object, and then giving you the opportunity for using them later on. If you want to simply see the basic results from running the model, you can use the standard summary() function.
summary(fit1_90)
## ## Call: ## lm(formula = HR90 ~ RD90 + SOUTH_f, data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -16.48 -3.00 -0.58 2.22 68.15 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.727 0.139 33.9 <2e-16 *** ## RD90 2.965 0.111 26.8 <2e-16 *** ## SOUTH_f1 3.181 0.222 14.3 <2e-16 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 5.3 on 3082 degrees of freedom ## Multiple R-squared: 0.365, Adjusted R-squared: 0.364 ## F-statistic: 884 on 2 and 3082 DF, p-value: <2e-16
The coefficients suggest that, according to our fitted model, when comparing two counties with the same level of resource deprivation, the homicide rate on average will be about 3 points higher in the Southern counties than in the Northern counties. And when comparing counties "adjusting" for their location, the homicide rate will be also around 3 points higher for every one unit increase in the measure of resource deprivation. Although it is tempting to talk of the "effect" of these variables in the outcome (homicide rate), this can mislead you and your reader about the capacity of your analysis to establish causality, so it is safer to interpret regression coefficients as comparisons [190].
With more than one input, you need to ask yourself whether all of the regression coefficients are zero. This hypothesis is tested with an F test. Again we are assuming the residuals are normally distributed, though with large samples the F statistics approximates the F distribution. You see the F test printed at the bottom of the summary output and the associated p value, which in this case is way below the conventional .05 that we use to declare statistical significance and reject the null hypothesis. At least one of our inputs must be related to our response variable. Notice that the table printed also reports a t test for each of the predictors. These are testing whether each of these predictors is associated with the response variable when adjusting for the other variables in the model. They report the "partial effect of adding that variable to the model" [198]. In this case we can see that all inputs seem to be significantly associated with the output.
The interpretation of regression coefficients is sensitive to the scale of measurement of the predictors. This means one cannot compare the magnitude of the coefficients to compare the relevance of variables to predict the response variable. One way of dealing with this is by rescaling the input variables. A common method involves subtracting the mean and dividing by the standard deviation of each numerical input. The coefficient in these models is the expected difference in the response variable, comparing units that differ by one standard deviation in the predictor while adjusting for other predictors in the model. Instead, [189] has proposed dividing each numeric variables by two times its standard deviation, so that the generic comparison is with inputs equal to plus/minus one standard deviation. As Gelman explains, the resulting coefficients are then comparable to untransformed binary predictors. The implementation of this approach in the
arm package subtracts the mean of each binary input, while it subtracts the mean and divides by two standard deviations for every numeric input. The way we would obtain these rescaled inputs uses the standardize() function of the arm package, which takes as an argument the name of the stored fit model.
standardize(fit1_90)
## ## Call: ## lm(formula = HR90 ~ z.RD90 + c.SOUTH_f, data = ncovr) ## ## Coefficients: ## (Intercept) z.RD90 c.SOUTH_f ## 6.18 5.93 3.18
Notice the main change affects the numerical predictors. The unstandardised coefficients are influenced by the degree of variability in your predictors, which means that typically they will be larger for your binary inputs. With unstandardised coefficients, you are comparing complete change in one variable (whether one is a Southern county or not) with one-unit changes in your numerical variable, which may not amount to much change. So, by putting in a comparable scale, you avoid this problem. Standardising in the way described here will help you to make fairer comparisons. These standardised coefficients are comparable in a way that the unstandardised coefficients are not. We can now see what inputs have a comparatively stronger effect. It is very important to realise, though, that one should not compare standardised coefficients across different models.
In the social sciences there is a great interest in what are called conditional hypotheses or interactions. Many of our theories do not assume simply additive effects but multiplicative effects. For example, you may think that a particular crime prevention programme may work in some environments but not in others. The interest in this kind of conditional hypothesis is growing. One of the assumptions of the regression model is that the relationship between the response variable and your predictors is additive. That is, if you have two predictors \(x_1\) and \(x_2\text{,}\) regression assumes that the effect of \(x_1\) on \(y\) is the same at all levels of \(x_2\text{.}\) If that is not the case, you are then violating one of the assumptions of regression. This is in fact one of the most important assumptions of regression, even if researchers often overlook it.
One way of extending our model to accommodate for interaction effects is to add more terms to our model, a third predictor \(x_3\text{,}\) where \(x_3\) is simply the product of multiplying \(x_1\) by \(x_2\text{.}\) Notice we keep a term for each of the main effects (the original predictors) as well as a new term for the interaction effect. "Analysts should include all constitutive terms when specifying multiplicative interaction models except in very rare circumstances" [199].
How do we do this in R? One way is to use the following notation in the formula argument. Notice how we have added a third term
RD90:SOUTH_f, which is asking R to test the conditional hypothesis that resource deprivation may have a different impact on homicide for southern and northern counties.
fit2_90 <- lm(HR90 ~ RD90 + SOUTH_f + RD90:SOUTH_f , data=ncovr)
# which is equivalent to:
# fit_2 <- lm(HR90 ~ RD90 * SOUTH_f , data=ncovr)
summary(fit2_90)
## ## Call: ## lm(formula = HR90 ~ RD90 + SOUTH_f + RD90:SOUTH_f, data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.05 -3.00 -0.57 2.23 68.14 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.548 0.159 28.68 <2e-16 *** ## RD90 2.581 0.196 13.15 <2e-16 *** ## SOUTH_f1 3.261 0.225 14.51 <2e-16 *** ## RD90:SOUTH_f1 0.562 0.238 2.37 0.018 * ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 5.29 on 3081 degrees of freedom ## Multiple R-squared: 0.366, Adjusted R-squared: 0.365 ## F-statistic: 592 on 3 and 3081 DF, p-value: <2e-16
You see here that essentially you have only two inputs (resource deprivation and south) but several regression coefficients. [188] suggest reserving the term input for the variables encoding the information and to use the term predictor to refer to each of the terms in the model. So here we have two inputs and three predictors (one for SOUTH, another for resource deprivation, and a final one for the interaction effect). In this case the test for the interaction effect is significant, which suggests there may be such an interaction. Let’s visualise the results with the
effects package:
plot(allEffects(fit2_90), ask=FALSE)
Notice that essentially what we are doing is running two regression lines and testing whether the slope is different for the two groups. The intercept is different; we know that Southern counties are more violent, but what we are testing here is whether the level of homicide goes up in a steeper fashion (and in the same direction) for one or the other group as the level of resource deprivation goes up. We see that’s the case here. The estimated lines are almost parallel, but the slope is a bit steeper in the Southern counties. In Southern counties, resource deprivation seems to have more of an impact on homicide than in Northern counties. This is related to one of the issues we raised earlier regarding spatial regression, that of heterogeneity. Here what we see is that there is some evidence of non-stationarity in the relationship between our predictors and our outcome.
A word of warning about interactions: the moment you introduce an interaction effect, the meaning of the coefficients for the other predictors changes (what is often referred to as the "main effects" as opposed to the interaction effect). You cannot retain the interpretation we introduced earlier. Now, for example, the coefficient for the
SOUTH variable relates to the marginal effect of this variable when RD90 equals zero. The typical table of results helps you to understand whether the effects are significant but offers little of interest that will help you to meaningfully interpret what the effects are. For this, is better to use some of the graphical displays we have covered.
Essentially what happens is that the regression coefficients that get printed are interpretable only for certain groups. So now:
-
The intercept still represents the predicted score of homicide for Southern counties and have a score of 0 in resource deprivation (as before).
-
The coefficient of SOUTH_f1 now can be thought of as the difference between the predicted score of homicide rate for Northern counties that have a score of 0 in resource deprivation and Southern counties that have a score of 0 in resource deprivation.
-
The coefficient of RD90 now becomes the comparison of mean homicide rate for Southern counties that differ by one point in resource deprivation.
-
The coefficient for the interaction term represents the difference in the slope for RD90 comparing Southern and Northern counties, the difference in the slope of the two lines that we visualised above.
Models with interaction terms are too often misinterpreted. We strongly recommend you read [199] to understand some of the issues involved.
