We start by reading the data and recreating setting in groups.
> library(dplyr)
> fpe <- read.table("https://grodri.github.io/datasets/effort.dat") |>
+ mutate(setting_g = cut(setting, breaks=c(min(setting),70,80,max(setting)),
+ right=FALSE, include.lowest=TRUE, labels=c("Low","Medium","High"))) Let us now create a copy of program effort and group it into
categories 0-4, 5-14 and 15+. We’ll call the variable
effort_g for effort in groups.
> fpe <- mutate(fpe, effort_g = cut(effort, breaks=c(min(effort), 5, 15, max(effort)),
+ right=FALSE, include.lowest=TRUE, labels=c("Weak","Moderate","Strong")))Here’s a table showing steeper declines in countries with strong programs, with a smaller difference between weak and moderate:
> group_by(fpe, effort_g) |> summarize(change = mean(change))# A tibble: 3 × 2
effort_g change
<fct> <dbl>
1 Weak 5
2 Moderate 9.33
3 Strong 27.9 Let us fit a model treating both setting and effort as factor variables, with weak programs in low settings as the reference cell
> m2g <- lm(change ~ setting_g + effort_g, data = fpe)
> summary(m2g)
Call:
lm(formula = change ~ setting_g + effort_g, data = fpe)
Residuals:
Min 1Q Median 3Q Max
-7.5347 -4.5883 0.5617 1.7909 11.7845
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.379 3.105 1.732 0.103723
setting_gMedium -1.681 3.855 -0.436 0.669035
setting_gHigh 2.388 4.457 0.536 0.599989
effort_gModerate 3.836 3.575 1.073 0.300138
effort_gStrong 20.672 4.339 4.764 0.000251 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.188 on 15 degrees of freedom
Multiple R-squared: 0.7833, Adjusted R-squared: 0.7255
F-statistic: 13.55 on 4 and 15 DF, p-value: 7.19e-05> anova(m2g)Analysis of Variance Table
Response: change
Df Sum Sq Mean Sq F value Pr(>F)
setting_g 2 1193.79 596.89 15.588 0.0002176 ***
effort_g 2 882.02 441.01 11.517 0.0009320 ***
Residuals 15 574.39 38.29
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Compare these estimates with the results in Table 2.15 and the anova with Table 2.16 in the lecture notes.
Countries with strong family planning programs show steeper declines than countries with weak programs at the same level of social setting, on average 21 percentage points more.
This statement is based on the assumption of additivity, namely that the difference in outcomes across categories of program effort is the same at every level of setting. We will test this assumption below.
As you can see the differences in fertility change by level of effort, among countries with the same level of setting are significant, with an F-ratio of 11.5 on 2 and 15 d.f.
Let us reproduce Table 2.17 in the notes, showing fitted means by
grouped setting and effort. I will use fitted() to generate
fitted values, and then summarize() them by the two
relevant factors
> library(tidyr)
Attaching package: 'tidyr'The following objects are masked from 'package:Matrix':
expand, pack, unpack> fpe <- mutate(fpe, fitted = fitted(m2g))
> group_by(fpe, setting_g, effort_g) |> summarize(fitted = mean(fitted)) |>
+ pivot_wider(names_from=effort_g, values_from=fitted)`summarise()` has grouped output by 'setting_g'. You can override using the `.groups` argument.# A tibble: 3 × 4
# Groups: setting_g [3]
setting_g Weak Moderate Strong
<fct> <dbl> <dbl> <dbl>
1 Low 5.38 9.22 NA
2 Medium 3.70 7.53 24.4
3 High 7.77 11.6 28.4Can you get the missing cell in the upper right corner? How about the margins, can you get those values?
Let us now consider a model with an interaction between social setting and program effort, so differences in fertility decline by effort will vary by setting.
In R we use a colon : to indicate an interaction term,
and an asterisk * to include main effects and interactions,
so A*B is equivalent to A + B + A:B.
> m2gi <- lm(change ~ setting_g*effort_g, data=fpe)
> summary(m2gi)
Call:
lm(formula = change ~ setting_g * effort_g, data = fpe)
Residuals:
Min 1Q Median 3Q Max
-7.3333 -3.0000 0.0000 0.9375 11.6667
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.6667 3.5775 0.745 0.4704
setting_gMedium 3.3333 5.0594 0.659 0.5224
setting_gHigh 6.3333 7.1550 0.885 0.3935
effort_gModerate 8.5833 4.7326 1.814 0.0948 .
effort_gStrong 19.3333 6.6929 2.889 0.0136 *
setting_gMedium:effort_gModerate -14.5833 8.5786 -1.700 0.1149
setting_gHigh:effort_gModerate -6.5833 9.9594 -0.661 0.5211
setting_gMedium:effort_gStrong -0.3333 9.7974 -0.034 0.9734
setting_gHigh:effort_gStrong NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.196 on 12 degrees of freedom
Multiple R-squared: 0.8261, Adjusted R-squared: 0.7247
F-statistic: 8.146 on 7 and 12 DF, p-value: 0.0009208Oops, our software dropped a term. Why? Because there are no countries with strong programs in low settings, so we have only eight groups, but are trying to represent their means using nine parameters, which is obviously one too many. Fortunately this doesn’t affect testing.
> anova(m2gi)Analysis of Variance Table
Response: change
Df Sum Sq Mean Sq F value Pr(>F)
setting_g 2 1193.79 596.89 15.5458 0.0004664 ***
effort_g 2 882.02 441.01 11.4859 0.0016322 **
setting_g:effort_g 3 113.64 37.88 0.9866 0.4317952
Residuals 12 460.75 38.40
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We have no evidence that the differences by effort vary with social setting, with an F just below 1 on 3 and 12 d.f.
This makes the issue of interpreting parameters moot, but it may be worth noting briefly the problems caused by the empty cell. As things stand, the coefficient of moderate effort compares moderate with weak at low setting, but the coefficient of strong effort compares strong with weak at high setting. (Table 2.20 in the notes may help see this point. When the term for high and strong is dropped, the only difference between weak and strong programs at high setting is the coefficient of strong.)
The parametrization I like best for this problem combines the main effects of effort with the interactions, so we obtain differences between strong and weak, and between moderate and weak programs, at each level of setting. This allows us to omit the difference between strong and weak programs at low setting, which is the one we can’t identify. Try the specification below
lm( change ~ setting_g + effort_g:setting_g, data=fpe)You could have obtained the same results in this unit using dummy
variables. For the record, this is how you might build the dummies and
fit the models. For the additive model we need just four indicators, two
for each factor. In R I make sure the dummies are numeric, not boolean,
otherwise the coefficient would have names like
effortStrongTRUE.
d <- as.numeric # for short
mutate(fpe,
settingMed = d(setting_g=="Medium"), settingHigh = d(setting_g=="High"),
effortMod = d(effort_g=="Moderate"), effortStrong = d(effort_g=="Strong"))
lm(change ~ settingMed + settingHigh + effortMod + effortStrong, data=fpe)We need a total of four dummies to represent the interactions, which can be computed simply as the product of the indicators for the main effects.
fpe <- mutate(fpe,
seMedMod = settingMed*effortMod, seMedStrong=settingMed*effortStrong,
seHighMod = settingHigh*effortMod, seHighStrong = settingHigh*effortStrong)
lm(change settingMed + settingHigh + effortMod + effortStrong +
seMedMod + seMedStrong + seHighMod + seHighStrong, data = fpe)Our software will again omit a variable, but you have more control on what to drop. Can you figure out which dummies you would need to show the effects of effort at each level of setting?
Updated fall 2022