install.packages(c("fixest", "plm"))26 Fixed effects
Does the association between parent and respondent schooling persist when we compare respondents within the same region? The Chapter 24 chapter adjusted for age and sex. We now let regions have different intercepts while retaining a common parental-schooling slope. A second application uses repeated observations of the same people to explain what changes when the group is a person rather than a region.
By the end of this chapter, you will be able to:
- Fit and interpret a linear model with group fixed effects.
- Compare specifications on the same observations and check equivalent slope estimates.
- Verify the structure and within-person variation of a panel dataset.
- Fit person and year fixed effects with person-clustered uncertainty.
- Explain the limits of fixed effects, including absorbed variables and causal claims.
Use the table operations from earlier chapters, factors from Chapter 15, and formulas and coefficient tables from Chapter 24. Prepare gss24.rds using the Appendix A and keep it in your working directory. The setup below is complete for a fresh session. Package prefixes such as fixest:: identify a function’s package without attaching it.
The new fitting package is fixest. The later panel application also uses data included in plm. Install both once, rather than every time you run the chapter. This installation chunk is deliberately not evaluated when the book renders.
library(tidyverse)
library(tidystata)
gss24 <- readRDS("gss24.rds")26.1 Compare within regions
The main outcome remains educ, the highest year of school completed. maeduc records mother’s schooling or younger-parent schooling for same-sex parents, following the source definition. Age and sex retain their earlier roles. The new variable, region, identifies the respondent’s current region. It does not necessarily describe where the respondent grew up or attended school.
A group fixed effect gives each group its own intercept. Respondents in the same region share that region’s intercept, while each respondent retains their own predictor values. In our additive model, the parent-schooling, age, and sex coefficients are shared across regions. Different intercepts do not imply different slopes. An interaction would be needed to let the parent-schooling slope vary by region; that is a different model from the one fitted here.
For this GSS comparison, fixed effects are a way to include categorical region adjustment in a linear model. They do not turn the GSS into a panel. Each row is still a different respondent observed once. The later application will use repeated observations and person fixed effects, where the available comparisons change more substantially.
First reconstruct the preceding OLS and logit chapters’ complete-case sample, adding the requirement that region is known. We retain degree even though this chapter does not model it, so the comparison preserves the earlier sample definition. Both models below use the resulting table without further filtering. Male remains the reference sex; converting region to a factor ensures that region codes are treated as categories rather than a numeric trend.
region_data <- gss24 |>
select(id, educ, maeduc, age, sex, degree, region) |>
filter(!is.na(id), !is.na(educ), !is.na(maeduc), !is.na(age),
!is.na(sex), !is.na(degree), !is.na(region)) |>
mutate(
sex = factor(sex, levels = c("male", "female")),
region = factor(region)
)
nrow(region_data)[1] 3329tab(region_data, region)# A tibble: 4 × 6
region n percent cum valid valid_cum
<fct> <int> <dbl> <dbl> <dbl> <dbl>
1 northeast 507 15.2 15.2 15.2 15.2
2 midwest 792 23.8 39.0 23.8 39.0
3 south 1290 38.8 77.8 38.8 77.8
4 west 740 22.2 100 22.2 100
summ(region_data, educ, maeduc, age)# A tibble: 3 × 8
name type n unique min mean sd max
<chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 educ dbl 3329 20 0 14.5 2.82 20
2 maeduc dbl 3329 21 0 12.5 3.44 20
3 age dbl 3329 72 18 49.0 17.5 89
The analysis retains 3329 respondents. The frequency table counts respondents within each of the four recorded regions; it is not a count of repeated observations of people. The numerical summary describes schooling and age in years. These checks establish the sample and units before a model comparison. Excluding unknown values makes the comparison consistent, but does not establish that retained respondents resemble those excluded.
The GSS models in this chapter are unweighted teaching examples with ordinary model-based uncertainty. They do not account for the survey’s sampling design and do not provide survey-adjusted population estimates. The available survey weight is not a substitute for specifying that design. Nor would clustering on only four broad regions provide a reliable general solution to uncertainty here.
26.1.1 Include region indicators
We can fit region fixed effects with the familiar stats::lm(). Add region to the formula alongside the other explanatory terms. R represents the region factor using an omitted reference region and indicators for the other regions. In this sample, northeast is the reference region. Each displayed region coefficient is an intercept difference from the reference, conditional on the other modeled terms.
region_adjusted <- stats::lm(
educ ~ maeduc + age + sex, data = region_data
)
region_lm <- stats::lm(
educ ~ maeduc + age + sex + region, data = region_data
)
broom::tidy(region_lm, conf.int = TRUE)# A tibble: 7 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 10.0 0.277 36.2 1.22e-242 9.50 10.6
2 maeduc 0.299 0.0136 22.0 2.41e-100 0.273 0.326
3 age 0.0184 0.00268 6.88 7.26e- 12 0.0132 0.0237
4 sexfemale -0.0175 0.0921 -0.190 8.49e- 1 -0.198 0.163
5 regionmidwest -0.434 0.150 -2.89 3.85e- 3 -0.727 -0.140
6 regionsouth -0.216 0.138 -1.57 1.17e- 1 -0.487 0.0544
7 regionwest -0.0443 0.152 -0.291 7.71e- 1 -0.342 0.254
stats::nobs(region_adjusted)[1] 3329stats::nobs(region_lm)[1] 3329The output retains the familiar maeduc, age, and sexfemale coefficients, followed by the region contrasts. The omitted region has no separate coefficient because its intercept is already represented by (Intercept). A positive region contrast means higher fitted schooling than the reference region at the same modeled parent schooling, age, and sex. It is not a claim that moving there would increase schooling. Current residence may itself reflect earlier educational and employment choices.
Both models use 3329 observations. To focus on our question, make one coefficient table containing just the parental-schooling rows. bind_rows() stacks tables that have matching column names. Its named inputs become values in the new model column when we supply .id = "model". The familiar tidy() supplies coefficient tables; filter() and select() retain the comparison we want.
region_comparison <- bind_rows(
"Age and sex" = broom::tidy(region_adjusted, conf.int = TRUE),
"Age, sex, and region" = broom::tidy(region_lm, conf.int = TRUE),
.id = "model"
) |>
filter(term == "maeduc") |>
select(model, estimate, std.error, conf.low, conf.high)
region_comparison# A tibble: 2 × 5
model estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 Age and sex 0.299 0.0136 0.272 0.326
2 Age, sex, and region 0.299 0.0136 0.273 0.326
Each row is an estimate from a model, rather than a respondent. The columns contain the parental-schooling slope, its standard error, and the endpoints of its ordinary 95% confidence interval. The slope unit is respondent-schooling years per parent-schooling year.
knitr::kable() formats this small table for reading. digits = 3 rounds displayed numeric values to three decimals without changing the stored estimates. The caption identifies the sample and uncertainty; the automatic table number can be used in the text.
knitr::kable(region_comparison, digits = 3)| model | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|
| Age and sex | 0.299 | 0.014 | 0.272 | 0.326 |
| Age, sex, and region | 0.299 | 0.014 | 0.273 | 0.326 |
The parental-schooling association is essentially unchanged after allowing region intercepts. In Table 26.1, both estimates round to 0.299 additional years of respondent schooling per additional year of parent schooling. More precisely, they are 0.299025 without region and 0.299164 with region. We should report this small change rather than expect fixed effects to produce a dramatic revision.
The region-adjusted slope describes a shared within-region association, conditional on modeled age and sex. It does not hold every family circumstance constant. Broad current regions leave substantial differences within regions, and do not measure childhood institutions directly. A nearly unchanged coefficient neither proves that geography is unimportant nor establishes that the parental-schooling association is causal.
26.1.2 Absorb the intercepts
For a few groups, displaying their indicators is manageable. With hundreds of people or organizations, the group intercepts can crowd out the coefficients of interest. fixest::feols() fits a linear fixed-effects model while absorbing those intercepts: it accounts for them without listing every one among the reported slope coefficients.
In its formula, the single vertical bar separates ordinary explanatory terms from the fixed-effect grouping variable. educ ~ maeduc + age + sex | region keeps schooling as the outcome, puts the three explanatory variables before |, and places region after it. This formula separator is different from the native pipe |>, which passes an object into a function. vcov = "iid" requests ordinary independent, constant-variance model-based uncertainty for this GSS demonstration.
region_fe <- fixest::feols(
educ ~ maeduc + age + sex | region,
data = region_data, vcov = "iid"
)
region_feOLS estimation, Dep. Var.: educ
Observations: 3,329
Fixed-effects: region: 4
Standard-errors: IID
Estimate Std. Error t value Pr(>|t|)
maeduc 0.299164 0.013591 22.011086 < 2.2e-16 ***
age 0.018409 0.002677 6.877281 7.2584e-12 ***
sexfemale -0.017525 0.092127 -0.190229 8.4914e-01
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 2.62519 Adj. R2: 0.130056
Within R2: 0.128879The printed model identifies region as a fixed effect and reports the nonabsorbed coefficients. The parent-schooling estimate should match the lm() model that included the region factor. Absorption changes how we obtain and present that estimate, not the underlying additive comparison.
Use the earlier coefficient-table tools to inspect that equivalence. The maeduc rows should have the same estimate apart from numerical rounding. Check the shared observation count as well. Exact reporting conventions for degrees of freedom or fit measures can differ across software; matching slope estimates is the equivalence demonstrated here.
region_equivalence <- bind_rows(
"lm indicators" = broom::tidy(region_lm),
"feols absorption" = broom::tidy(region_fe),
.id = "method"
) |>
filter(term == "maeduc") |>
select(method, estimate)
region_equivalence# A tibble: 2 × 2
method estimate
<chr> <dbl>
1 lm indicators 0.299
2 feols absorption 0.299
stats::nobs(region_fe)[1] 3329Both implementations use the same respondents and produce the same parental-schooling slope. The fixed-effect model has not estimated a separate parental-schooling slope for every region. It has estimated a common slope while permitting the baseline level of schooling to differ across regions.
NoteOptional: subtract the region means
A bounded numerical demonstration helps explain “within region.” For this demonstration only, omit age and sex and compare a simpler region-adjusted model with a model of deviations from region means. This is not a reconstruction of the full adjusted model above.
Inside mutate(), .by = region calculates each mean within a region and repeats that mean for its members. Subtracting it gives a respondent’s distance above or below the region mean. We create two new columns in a separate table, preserving the original schooling values.
region_centered <- region_data |>
mutate(
educ_within = educ - mean(educ),
maeduc_within = maeduc - mean(maeduc),
.by = region
)
region_centered |>
select(region, educ, educ_within, maeduc, maeduc_within) |>
head()# A tibble: 6 × 5
region educ educ_within maeduc maeduc_within
<fct> <dbl> <dbl> <dbl> <dbl>
1 northeast 16 1.26 14 1.43
2 northeast 14 -0.742 12 -0.568
3 northeast 13 -1.74 11 -1.57
4 northeast 14 -0.742 5 -7.57
5 northeast 13 -1.74 12 -0.568
6 northeast 14 -0.742 0 -12.6
For example, a negative educ_within means fewer years of schooling than the sample mean in that respondent’s region. The same row’s maeduc_within gives the analogous parent-schooling deviation. These are differences in years, not standardized scores or percentages.
In an lm() formula, - 1 removes the intercept. After subtracting both region means, we fit a line through zero to the two deviation columns. The following two slopes match because both describe the simple schooling association after allowing region-specific intercepts.
region_simple <- stats::lm(educ ~ maeduc + region, data = region_data)
region_within <- stats::lm(
educ_within ~ maeduc_within - 1, data = region_centered
)
broom::tidy(region_simple) |>
filter(term == "maeduc")# A tibble: 1 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 maeduc 0.279 0.0133 20.9 1.73e-91
broom::tidy(region_within)# A tibble: 1 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 maeduc_within 0.279 0.0133 20.9 1.35e-91
This demonstrates slope equivalence only. Ordinary uncertainty from the demeaned fit does not automatically count the estimated region intercepts correctly. To reproduce the full adjusted slope by this route, we would also need to transform all the other explanatory terms appropriately. Use the full lm() or feols() fit for reporting, rather than treating this two-column illustration as a replacement.
26.2 Secondary application: a panel
How are wages associated with union membership and marriage when we follow the same people over time? The GSS extract does not follow respondents annually, so this question requires a different dataset. We now use Males, a historical panel supplied by plm. Keeping the application together makes the dataset change explicit: the unit becomes one person in one year, and the outcome becomes the logarithm of hourly wages.
The installed plm help entry for Males describes observations of 545 young men over 1980–1987, with eight annual observations per person. Its source is the wage-panel application associated with Vella and Verbeek (1998), “Whose wages do unions raise? A dynamic model of unionism and wage rate determination for young men,” Journal of Applied Econometrics, 13, 163–183. This restricted historical sample should not be described as all workers or as a contemporary national workforce.
Record the package versions alongside this source reference. utils::packageVersion() takes a quoted package name and returns its installed version. The first result identifies the data package; the second identifies the fixed-effects implementation used for this application. Recording versions helps distinguish this particular reproducible example from later package updates.
utils::packageVersion("plm")[1] '2.6.7'utils::packageVersion("fixest")[1] '0.14.1'utils::data() loads a dataset distributed inside an installed package. Its first argument names the dataset and package = "plm" identifies its source. It creates Males in the session; no external file download is needed for the lesson.
utils::data("Males", package = "plm")as_tibble() converts this data frame to the table form used throughout the book. We assign the result to a new object so the packaged source remains available unchanged.
panel_source <- as_tibble(Males)
panel_source# A tibble: 4,360 × 12
nr year school exper union ethn married health wage industry
<int> <int> <int> <int> <fct> <fct> <fct> <fct> <dbl> <fct>
1 13 1980 14 1 no other no no 1.20 Business_and_Repa…
2 13 1981 14 2 yes other no no 1.85 Personal_Service
3 13 1982 14 3 no other no no 1.34 Business_and_Repa…
4 13 1983 14 4 no other no no 1.43 Business_and_Repa…
5 13 1984 14 5 no other no no 1.57 Personal_Service
6 13 1985 14 6 no other no no 1.70 Business_and_Repa…
7 13 1986 14 7 no other no no -0.720 Business_and_Repa…
8 13 1987 14 8 no other no no 1.67 Business_and_Repa…
9 17 1980 13 4 no other no no 1.68 Trade
10 17 1981 13 5 no other no no 1.52 Trade
# ℹ 4,350 more rows
# ℹ 2 more variables: occupation <fct>, residence <fct>
The source uses nr for person ID and wage for the natural logarithm of hourly wages. In particular, wage is already logged: taking its logarithm again would change the outcome incorrectly. Use familiar rename() to make these meanings visible as person_id and log_wage. year is the observation year; union and married are yes/no factors. school records schooling years, and ethn identifies the source ethnicity categories.
panel_data <- panel_source |>
rename(person_id = nr, log_wage = wage) |>
select(person_id, year, log_wage, union, married, school, ethn)
head(panel_data)# A tibble: 6 × 7
person_id year log_wage union married school ethn
<int> <int> <dbl> <fct> <fct> <int> <fct>
1 13 1980 1.20 no no 14 other
2 13 1981 1.85 yes no 14 other
3 13 1982 1.34 no no 14 other
4 13 1983 1.43 no no 14 other
5 13 1984 1.57 no no 14 other
6 13 1985 1.70 no no 14 other
Each displayed row now has an explicit person identifier and year. Repeated IDs are expected in a panel; duplicate person-year combinations are not. We must check this structure before treating the observations as repeated measures.
26.2.1 Check repeated observations
count(person_id, year) creates one row per observed ID-year combination and counts its records. Filtering for n != 1 shows any combinations that occur more than once. An empty result means every observed combination is unique; it does not yet establish that every expected year exists for every person.
panel_duplicate_keys <- panel_data |>
count(person_id, year) |>
filter(n != 1)
panel_duplicate_keys# A tibble: 0 × 3
# ℹ 3 variables: person_id <int>, year <int>, n <int>
Next count observations per person, naming that count waves, then count how many people have each number of waves. Count years separately to verify the calendar coverage. These are familiar table operations with two different denominators: people in the first result, person-year records in the second.
panel_wave_counts <- panel_data |>
count(person_id, name = "waves") |>
count(waves, name = "persons")
panel_wave_counts# A tibble: 1 × 2
waves persons
<int> <int>
1 8 545
count(panel_data, year)# A tibble: 8 × 2
year n
<int> <int>
1 1980 545
2 1981 545
3 1982 545
4 1983 545
5 1984 545
6 1985 545
7 1986 545
8 1987 545
nrow(panel_data)[1] 4360The results show 545 people, each with eight observations, for 4,360 person-year rows. Every year from 1980 through 1987 has 545 records. Together with the unique keys, these checks establish a balanced annual panel over this period. The 4,360 records are not 4,360 independent people.
Check missing core model values separately. This filter keeps a row if any required value is missing; | inside the filter means logical “or,” a third use distinct from |> and the feols() formula separator. nrow() counts the selected rows. We do not overwrite panel_data or silently drop records.
panel_missing <- panel_data |>
filter(is.na(person_id) | is.na(year) | is.na(log_wage) |
is.na(union) | is.na(married))
nrow(panel_missing)[1] 0Zero rows have missing values in these core variables. Both wage models can therefore use all 4,360 rows. This check concerns recorded values; it does not assess measurement accuracy or how these people entered the source panel.
26.2.2 Find within-person variation
A person fixed effect permits every person to have a different intercept. The explanatory slopes then use changes within people over time. A variable that never changes for a person cannot explain that person’s deviations from their own level.
n_distinct() counts distinct observed values. Inside summarise() with .by = person_id, it tells us how many union states, marital states, or schooling values each person has across their eight records. Because the core status variables are complete, one state means no recorded change and two states means at least one yes/no change.
panel_variation <- panel_data |>
summarise(
union_states = n_distinct(union),
married_states = n_distinct(married),
school_values = n_distinct(school),
.by = person_id
)
count(panel_variation, union_states)# A tibble: 2 × 2
union_states n
<int> <int>
1 1 299
2 2 246
count(panel_variation, married_states)# A tibble: 2 × 2
married_states n
<int> <int>
1 1 235
2 2 310
count(panel_variation, school_values)# A tibble: 1 × 2
school_values n
<int> <int>
1 1 545
The first two results distinguish people who change union or marital status from those whose recorded status is constant. The last shows one schooling value for all 545 people. Person fixed effects therefore absorb schooling in this dataset: there is no within-person schooling change from which to estimate a separate schooling slope. Source ethnicity is also time invariant and is absorbed.
“Absorbed” does not mean that schooling or ethnicity is substantively unimportant. Their stable association with each person’s wage level cannot be separated from that person’s unrestricted intercept in this specification. Similarly, linear experience increases with calendar time within each person. Its linear term is collinear with person and year effects here, so we do not add it to the two-way model. Estimating a coefficient requires variation that remains after all included fixed effects are accounted for.
26.2.3 Fit person and year effects
First fit a pooled wage model with year fixed effects. It compares records across different people as well as within people, while allowing a common intercept shift for every calendar year. The union and married factors use no as their reference categories; unionyes and marriedyes are yes-minus-no contrasts. Check these categories using the familiar frequency tables before fitting.
tab(panel_data, union)# A tibble: 2 × 6
union n percent cum valid valid_cum
<fct> <int> <dbl> <dbl> <dbl> <dbl>
1 no 3296 75.6 75.6 75.6 75.6
2 yes 1064 24.4 100 24.4 100
tab(panel_data, married)# A tibble: 2 × 6
married n percent cum valid valid_cum
<fct> <int> <dbl> <dbl> <dbl> <dbl>
1 no 2446 56.1 56.1 56.1 56.1
2 yes 1914 43.9 100 43.9 100
For panel uncertainty, vcov = ~ person_id requests standard errors clustered by person. The one-sided formula names the cluster identifier; it does not add a predictor. This allows errors from different years for the same person to be dependent, while treating different people as independent clusters. It changes the uncertainty calculation, not the fitted slope coefficients. We have 545 person clusters, unlike the four regions in the GSS example.
panel_year <- fixest::feols(
log_wage ~ union + married | year,
data = panel_data, vcov = ~ person_id
)
panel_yearOLS estimation, Dep. Var.: log_wage
Observations: 4,360
Fixed-effects: year: 8
Standard-errors: Clustered (person_id)
Estimate Std. Error t value Pr(>|t|)
unionyes 0.176175 0.029253 6.02256 3.1625e-09 ***
marriedyes 0.142476 0.027835 5.11862 4.2762e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.501477 Adj. R2: 0.111452
Within R2: 0.0412 The union and marriage coefficients are positive. Because this pooled model does not allow separate person intercepts, these associations can reflect stable differences between people as well as changes in their circumstances. Year effects account for common year shifts, not for every individual wage trajectory.
Now put person_id + year after the formula separator. The plus sign includes both sets of intercepts. Person fixed effects account for stable person-level differences; year effects account for changes shared across the sample in a given year. The union and marriage coefficients use the remaining within-person variation, after those year shifts are taken into account.
panel_person_year <- fixest::feols(
log_wage ~ union + married | person_id + year,
data = panel_data, vcov = ~ person_id
)
panel_person_yearOLS estimation, Dep. Var.: log_wage
Observations: 4,360
Fixed-effects: person_id: 545, year: 8
Standard-errors: Clustered (person_id)
Estimate Std. Error t value Pr(>|t|)
unionyes 0.083370 0.023060 3.61529 0.0003279 ***
marriedyes 0.058337 0.021337 2.73403 0.0064602 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.330217 Adj. R2: 0.559649
Within R2: 0.007635stats::nobs(panel_year)[1] 4360stats::nobs(panel_person_year)[1] 4360Both models use all 4,360 person-year records. The smaller union and marriage coefficients in the second model therefore do not result from changing the sample. Their comparison has changed from a mixture of between-person and within-person information to information within people, net of common year effects.
Combine the two coefficient tables as before. For these models, tidy() obtains the uncertainty specified when each model was fitted, so the standard errors and confidence intervals reflect clustering by person. This output has two rows per model, one for each nonabsorbed explanatory term.
panel_comparison <- bind_rows(
"Year effects" = broom::tidy(panel_year, conf.int = TRUE),
"Person and year effects" = broom::tidy(panel_person_year, conf.int = TRUE),
.id = "model"
) |>
select(model, term, estimate, std.error, conf.low, conf.high)
knitr::kable(panel_comparison, digits = 3)| model | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|
| Year effects | unionyes | 0.176 | 0.029 | 0.119 | 0.234 |
| Year effects | marriedyes | 0.142 | 0.028 | 0.088 | 0.197 |
| Person and year effects | unionyes | 0.083 | 0.023 | 0.038 | 0.129 |
| Person and year effects | marriedyes | 0.058 | 0.021 | 0.016 | 0.100 |
The positive wage associations are smaller when we compare the same people across years. In Table 26.2, the union coefficient changes from about 0.176 to 0.083 log points, and the marriage coefficient from about 0.142 to 0.058 log points. The second model’s union coefficient describes the fitted difference between a person’s union and nonunion records, at the same modeled marital status and after accounting for common year shifts. The marriage coefficient has the corresponding interpretation at the same modeled union status.
The outcome is a natural logarithm, so these coefficients are differences in log hourly wages, called log points. They are not currency amounts or percentage-point differences. We keep that scale throughout this comparison rather than relabel the estimates as exact percentage changes. The interval columns also use log points and quantify model-based uncertainty with person clustering.
These within-person comparisons do not identify causal union or marriage effects on their own. Job changes, health, employment conditions, and other changing circumstances can affect both wages and the explanatory variables. People may enter a union job or marry partly in response to changing prospects. Person effects remove additive differences that stay constant, not changing omitted explanations. Year effects remove common time shifts, not every person’s particular trajectory. Clustering addresses dependence in uncertainty, not these identification problems.
26.3 Report the comparison
In the GSS application, the parental-schooling association is nearly unchanged after allowing current-region intercepts. It remains an unweighted conditional association among the same complete-case respondents, with no survey-design correction. Region fixed effects adjust levels across four groups, but neither reconstruct childhood context nor explain away differences within a region.
In the secondary panel application, union membership and marriage have smaller positive log-wage associations after adding person effects to the year effects. The contrast is informative because the records are unchanged and the comparison is explicit. Repeated observations permit us to separate stable person differences from within-person change. They do not ensure that changes in union membership or marriage occurred independently of other wage determinants.
Report the observation unit, sample, fixed-effect groups, coefficient scale, and uncertainty method together. “Fixed effects” alone does not tell readers whether the groups were regions, people, years, or another unit. It also does not state whether uncertainty was ordinary, clustered, or adjusted for a survey design.
26.4 Review
- Group fixed effects give groups separate intercepts while keeping the specified slopes common.
lm()includes factor indicators;feols()can absorb their intercepts and recover equivalent slopes.- Keep the sample fixed when comparing specifications and check observation counts with
nobs(). count()checks panel keys and waves;n_distinct()identifies available within-person variation.- Person effects absorb time-invariant variables; person and year effects can also absorb common linear trends.
- State coefficient units and the uncertainty method. Fixed effects and clustered standard errors do not establish a causal interpretation.
The optional Programming part develops expressions and objects further. These models can already be reconstructed from the data preparation, formulas, and checks shown here.
26.5 Exercises
- Reconstruct the GSS region sample and fit the age-and-sex model and the model adding region. Check their observation counts, extract the two
maeducslopes, and interpret the change. Explain why region is not necessarily childhood context. - Optional: subtract region means. Reconstruct the GSS sample and compare the simple region-adjusted parental-schooling slope with the slope obtained after subtracting region means from
educandmaeduc. Explain why this does not reproduce the model that also adjusts for age and sex, and why the demeaned fit’s standard error should not be substituted for the full model’s. - Load
Malesindependently. Verify unique person-year keys, eight waves per person, calendar coverage, and complete core model variables. Count how many people have one or two observed union and marriage states. Explain why unchanged schooling cannot have a separate slope with person effects. - Recreate the panel and fit the year-only and person-plus-year fixed-effects wage models, clustering uncertainty by person in both. Compare their observation counts and coefficients. Explain the log-point scale and give one changing circumstance that could prevent a causal interpretation.
Worked answers
NoteCompare region adjustment
Read the saved extract and load the packages in a fresh session. This answer creates its own sample and models.
library(tidyverse)
library(tidystata)
answer_gss24 <- readRDS("gss24.rds")answer_region_data <- answer_gss24 |>
select(id, educ, maeduc, age, sex, degree, region) |>
filter(!is.na(id), !is.na(educ), !is.na(maeduc), !is.na(age),
!is.na(sex), !is.na(degree), !is.na(region)) |>
mutate(
sex = factor(sex, levels = c("male", "female")),
region = factor(region)
)
answer_adjusted <- stats::lm(
educ ~ maeduc + age + sex, data = answer_region_data
)
answer_region <- stats::lm(
educ ~ maeduc + age + sex + region, data = answer_region_data
)
stats::nobs(answer_adjusted)[1] 3329stats::nobs(answer_region)[1] 3329answer_region_comparison <- bind_rows(
"Age and sex" = broom::tidy(answer_adjusted, conf.int = TRUE),
"Age, sex, and region" = broom::tidy(answer_region, conf.int = TRUE),
.id = "model"
) |>
filter(term == "maeduc") |>
select(model, estimate, std.error, conf.low, conf.high)
answer_region_comparison# A tibble: 2 × 5
model estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 Age and sex 0.299 0.0136 0.272 0.326
2 Age, sex, and region 0.299 0.0136 0.273 0.326
Both fits use 3,329 respondents. Parent schooling has essentially the same positive association with respondent schooling in the two specifications, about 0.299 schooling years per parent-schooling year. Adding region permits different regional intercepts; it does not make the parental-schooling slope region specific. Current residence can differ from childhood residence and may reflect later moves connected to education or work. The unchanged slope is a result to report, not evidence that regional institutions never matter. These unweighted models retain ordinary model-based uncertainty and do not identify a causal effect.
NoteOptional: recover a simple within-region slope
Rebuild the source and sample independently, then create deviations within region. The two models below deliberately omit age and sex.
library(tidyverse)
library(tidystata)
answer_gss24 <- readRDS("gss24.rds")answer_region_data <- answer_gss24 |>
select(id, educ, maeduc, age, sex, degree, region) |>
filter(!is.na(id), !is.na(educ), !is.na(maeduc), !is.na(age),
!is.na(sex), !is.na(degree), !is.na(region)) |>
mutate(
sex = factor(sex, levels = c("male", "female")),
region = factor(region)
)
answer_centered <- answer_region_data |>
mutate(
educ_within = educ - mean(educ),
maeduc_within = maeduc - mean(maeduc),
.by = region
)
answer_simple_region <- stats::lm(
educ ~ maeduc + region, data = answer_region_data
)
answer_simple_within <- stats::lm(
educ_within ~ maeduc_within - 1, data = answer_centered
)
broom::tidy(answer_simple_region) |>
filter(term == "maeduc")# A tibble: 1 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 maeduc 0.279 0.0133 20.9 1.73e-91
broom::tidy(answer_simple_within)# A tibble: 1 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 maeduc_within 0.279 0.0133 20.9 1.35e-91
The two slope estimates match apart from numerical rounding. Subtracting the group means removes the region-specific levels from these two numeric variables. This is a check on the simple region-adjusted slope, not on the slope additionally adjusted for age and sex. The other explanatory terms would also need appropriate transformation for that calculation. The demeaned lm() fit does not automatically account for the region intercepts in its residual degrees of freedom. Report uncertainty from the original fixed-effects model rather than borrowing the demonstration’s standard error.
NoteCheck a panel before fitting
This answer starts from the packaged source, so it does not depend on the GSS or any objects from the main application.
library(tidyverse)
library(tidystata)
utils::data("Males", package = "plm")
answer_panel <- as_tibble(Males) |>
rename(person_id = nr, log_wage = wage) |>
select(person_id, year, log_wage, union, married, school, ethn)
answer_panel |>
count(person_id, year) |>
filter(n != 1)# A tibble: 0 × 3
# ℹ 3 variables: person_id <int>, year <int>, n <int>
answer_panel |>
count(person_id, name = "waves") |>
count(waves, name = "persons")# A tibble: 1 × 2
waves persons
<int> <int>
1 8 545
count(answer_panel, year)# A tibble: 8 × 2
year n
<int> <int>
1 1980 545
2 1981 545
3 1982 545
4 1983 545
5 1984 545
6 1985 545
7 1986 545
8 1987 545
answer_panel |>
filter(is.na(person_id) | is.na(year) | is.na(log_wage) |
is.na(union) | is.na(married)) |>
nrow()[1] 0answer_variation <- answer_panel |>
summarise(
union_states = n_distinct(union),
married_states = n_distinct(married),
school_values = n_distinct(school),
.by = person_id
)
count(answer_variation, union_states)# A tibble: 2 × 2
union_states n
<int> <int>
1 1 299
2 2 246
count(answer_variation, married_states)# A tibble: 2 × 2
married_states n
<int> <int>
1 1 235
2 2 310
count(answer_variation, school_values)# A tibble: 1 × 2
school_values n
<int> <int>
1 1 545
The key check returns zero rows. There are 545 people with eight waves each, 545 records in every year from 1980 to 1987, and no missing core model values. The union and marriage tables count people with constant versus changing recorded status; their denominators are 545 people, not 4,360 records. All people have one schooling value. A person intercept can already account for that person’s stable schooling contribution, so its separate coefficient cannot be identified alongside unrestricted person effects. A variable’s absence from the reported slope table does not establish its substantive irrelevance.
NoteCompare within-person wage associations
Load the panel and recreate both models from their source. fixest and broom are called with package prefixes; their functions do not require attached packages.
library(tidyverse)
library(tidystata)
utils::data("Males", package = "plm")
answer_panel <- as_tibble(Males) |>
rename(person_id = nr, log_wage = wage) |>
select(person_id, year, log_wage, union, married, school, ethn)
answer_year <- fixest::feols(
log_wage ~ union + married | year,
data = answer_panel, vcov = ~ person_id
)
answer_person_year <- fixest::feols(
log_wage ~ union + married | person_id + year,
data = answer_panel, vcov = ~ person_id
)
stats::nobs(answer_year)[1] 4360stats::nobs(answer_person_year)[1] 4360answer_panel_comparison <- bind_rows(
"Year effects" = broom::tidy(answer_year, conf.int = TRUE),
"Person and year effects" = broom::tidy(
answer_person_year, conf.int = TRUE
),
.id = "model"
) |>
select(model, term, estimate, std.error, conf.low, conf.high)
answer_panel_comparison# A tibble: 4 × 6
model term estimate std.error conf.low conf.high
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Year effects unionyes 0.176 0.0293 0.119 0.234
2 Year effects marriedyes 0.142 0.0278 0.0878 0.197
3 Person and year effects unionyes 0.0834 0.0231 0.0381 0.129
4 Person and year effects marriedyes 0.0583 0.0213 0.0164 0.100
Both fits retain 4,360 records. The union coefficient falls from about 0.176 to 0.083 log points and the marriage coefficient from about 0.142 to 0.058 log points after adding person effects. The second model compares each person’s changing circumstances, after allowing common year shifts, rather than pooling stable differences between people into those slopes. These are differences in natural log hourly wages, not currency or exact percentage changes. A job change could alter both union membership and wages, for example. Person and year effects do not remove every changing omitted circumstance, and person-clustered standard errors do not resolve that problem. The conclusions concern this historical sample of young men.