Homework 7
Deven Murray
2026-11-03
# SIPP data are in person-month format, meaning each record represents one month for a specific person.
# Unique persons are identified using SSUID+PNUM. Unique households are identified using SSUID+ERESIDENCEID. For
# additional guidance on using SIPP data, please see the SIPP Users' Guide at <https://www.census.gov/programs-surveys/sipp/guidance/users-guide.html>
# Note the 'select' statement in the first use of fread(). Most machines do not have enough memory to read
# the entire SIPP file into memory. Use a 'select' statement to read in only the columns you are interested in using.
# If you still encounter an out-of-memory error, you must select less columns or less observations.
# Load the "data.table", "dplyr", and "bit64" libraries. Use install.packages function if you don't already have these installed.
library("data.table")##
## Attaching package: 'data.table'## The following object is masked from 'package:base':
##
## %notin%library("bit64")## ********************************************************
## R-core is collecting use cases for 64-bit integers as they explore native support for these vectors.
##
## See https://stat.ethz.ch/pipermail/r-devel/2026-July/084631.html and reach out to Luke Tierney (luke-tierney@uiowa.edu).
## ********************************************************##
## Attaching package: 'bit64'## The following object is masked from 'package:utils':
##
## hashtab## The following objects are masked from 'package:base':
##
## %in%, :, array, as.factor, as.ordered, colSums, factor, intersect,
## is.double, is.element, match, matrix, order, rank, rowSums,
## setdiff, setequal, table, unionlibrary("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:bit64':
##
## intersect, setdiff, setequal, union## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("car")## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary("lmtest")## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:data.table':
##
## yearmon, yearqtr## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary("sandwich")# Store the text "pu2025.csv" in a name called ds
ds <- c("pu2025.csv")
# The file is pipe-delimited even though it's named .csv
pu <- fread(ds, sep = "|", select = c("SSUID", "ERESIDENCEID", "PNUM", "MONTHCODE", "SPANEL", "WPFINWGT", "ESEEING", "RFOODS", "RFOODR", "TAGE", "ESEX", "THHLDSTATUS"))
# Out of 5,203 variables, 12 variables were selected. The first six in the above line represent ID and structure variables. SSUID is sample unit ID with PNUM identifies a person and ERESIDENCEID identifies a household. PNUM is person number within the sample unit. MONTHCODE is reference month (1–12). SPANEL is panel year (2022-2025). WPFINWGT is the final person weight.The next six represent ESEEING blind or serious difficulty seeing even with glasses/contacts. RFOODS food security status recode: 1 high/marginal, 2 low, 3 very low. RFOODR is the raw food security score a count of affirmative responses (0–6). TAGE is age as of last birthday (capped at 90 for confidentiality). ESEX is sex with 1 Male, 2 Female. THHLDSTATUS represents household status and was added in the midst of coding due to needing to prove why I received 5,423 NAs.
# Make sure all the column names are upper-case
names(pu) <- toupper(names(pu))
# Preview the data
head(pu, 20)## SSUID ERESIDENCEID PNUM MONTHCODE SPANEL WPFINWGT ESEEING RFOODS
## <i64> <int> <int> <int> <int> <num> <int> <int>
## 1: 11428525324 100001 101 1 2024 8290.772 2 1
## 2: 11428525324 100001 101 2 2024 8476.619 2 1
## 3: 11428525324 100001 101 3 2024 9238.426 2 1
## 4: 11428525324 100001 101 4 2024 9371.330 2 1
## 5: 11428525324 100001 101 5 2024 9547.140 2 1
## 6: 11428525324 100001 101 6 2024 9337.827 2 1
## 7: 11428525324 100001 101 7 2024 9096.493 2 1
## 8: 11428525324 100001 101 8 2024 9310.043 2 1
## 9: 11428525324 100001 101 9 2024 9209.630 2 1
## 10: 11428525324 100001 101 10 2024 9150.450 2 1
## 11: 11428525324 100001 101 11 2024 9195.871 2 1
## 12: 11428525324 100001 101 12 2024 9433.126 2 1
## 13: 11428534325 100001 101 1 2025 4473.576 2 3
## 14: 11428534325 100001 101 2 2025 4561.484 2 3
## 15: 11428534325 100001 101 3 2025 4761.070 2 3
## 16: 11428534325 100001 101 4 2025 5111.557 2 3
## 17: 11428534325 100001 101 5 2025 5001.163 2 3
## 18: 11428534325 100001 101 6 2025 4933.283 2 3
## 19: 11428534325 100001 101 7 2025 4894.891 2 3
## 20: 11428534325 100001 101 8 2025 4716.378 2 3
## SSUID ERESIDENCEID PNUM MONTHCODE SPANEL WPFINWGT ESEEING RFOODS
## <i64> <int> <int> <int> <int> <num> <int> <int>
## RFOODR TAGE ESEX THHLDSTATUS
## <int> <int> <int> <int>
## 1: 0 70 1 1
## 2: 0 70 1 1
## 3: 0 70 1 1
## 4: 0 70 1 1
## 5: 0 70 1 1
## 6: 0 70 1 1
## 7: 0 70 1 1
## 8: 0 70 1 1
## 9: 0 70 1 1
## 10: 0 70 1 1
## 11: 0 70 1 1
## 12: 0 70 1 1
## 13: 6 50 1 2
## 14: 6 50 1 2
## 15: 6 50 1 2
## 16: 6 50 1 2
## 17: 6 50 1 2
## 18: 6 50 1 2
## 19: 6 50 1 2
## 20: 6 50 1 2
## RFOODR TAGE ESEX THHLDSTATUS
## <int> <int> <int> <int># Rows are person-months, not people therefore 379,215 rows equals 32,052 people.
nrow(pu)## [1] 379215n_distinct(pu$SSUID, pu$PNUM)## [1] 32052n_distinct(pu$SSUID, pu$ERESIDENCEID)## [1] 15215summary(pu)## SSUID ERESIDENCEID PNUM MONTHCODE
## Min. : 11428525324 Min. :100001 Min. :101.0 Min. : 1.000
## 1st Qu.:33286052753124 1st Qu.:100001 1st Qu.:101.0 1st Qu.: 3.000
## Median :66728586009222 Median :100001 Median :102.0 Median : 6.000
## Mean :56487541761701 Mean :115599 Mean :111.3 Mean : 6.482
## 3rd Qu.:87686017231924 3rd Qu.:100001 3rd Qu.:103.0 3rd Qu.: 9.000
## Max. :88199599599922 Max. :400005 Max. :405.0 Max. :12.000
##
## SPANEL WPFINWGT ESEEING RFOODS
## Min. :2022 Min. : 0 Min. :1.00 Min. :1.000
## 1st Qu.:2022 1st Qu.: 5673 1st Qu.:2.00 1st Qu.:1.000
## Median :2023 Median : 8573 Median :2.00 Median :1.000
## Mean :2023 Mean : 10505 Mean :1.95 Mean :1.157
## 3rd Qu.:2025 3rd Qu.: 13105 3rd Qu.:2.00 3rd Qu.:1.000
## Max. :2025 Max. :133335 Max. :2.00 Max. :3.000
## NAs :5423
## RFOODR TAGE ESEX THHLDSTATUS
## Min. :0.0000 Min. : 0.00 Min. :1.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:25.00 1st Qu.:1.000 1st Qu.:1.000
## Median :0.0000 Median :46.00 Median :2.000 Median :1.000
## Mean :0.4814 Mean :45.31 Mean :1.513 Mean :1.365
## 3rd Qu.:0.0000 3rd Qu.:66.00 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :6.0000 Max. :90.00 Max. :2.000 Max. :6.000
## # Descriptive stats shows 5,423 NAs in ESEEING, these NAs were not asked because they were outside the variable's universe, per the dictionary.No other variables report NA.
na_rows <- pu %>% filter(is.na(ESEEING))
table(na_rows$TAGE)##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 693 31 67 23 32 32 44 30 21 29 42 22 24 8 30 35 58 33 179 282
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 142 140 122 166 138 118 83 74 141 160 71 75 86 103 36 57 56 68 59 90
## 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
## 45 73 27 16 9 51 35 24 5 12 41 21 12 24 20 31 15 30 48 31
## 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## 28 13 20 48 39 11 28 14 20 15 39 24 43 33 66 67 69 33 67 21
## 80 81 82 83 84 85 86 87 88 89 90
## 54 40 5 12 37 58 20 11 86 47 85sum(na_rows$TAGE>= 1)## [1] 4730# 4730 are the NA rows aged 1 or older. Age doesn't explain them.
table(na_rows$THHLDSTATUS, na_rows$TAGE >= 1)##
## FALSE TRUE
## 1 10 0
## 2 677 0
## 3 0 1335
## 4 6 2909
## 5 0 140
## 6 0 346# All 5,423 NAs fall outside ESEEING's universe. 687 failed age only, 4,730 failed status only, 6 failed both.table(pu$SPANEL)##
## 2022 2023 2024 2025
## 120652 69613 75519 113431spanel_2025 <- pu %>% filter(SPANEL == 2025)
nrow(spanel_2025)## [1] 113431# Reference years differ across panels; weights are designed within panel. Kept 2025 for recency, sample from 2025 is 113,431 person-months.n_distinct(spanel_2025$SSUID, spanel_2025$PNUM)## [1] 9490month_1 <- spanel_2025 %>% filter(MONTHCODE == 1)
nrow(month_1)## [1] 9419# 9,490 people in the 2025 panel, one row per person is set, nrow produces 71-person gap.
first_months <- spanel_2025 %>% group_by(SSUID, PNUM) %>% summarise(first_month = min(MONTHCODE), .groups = "drop")
table(first_months$first_month)##
## 1 2 3 5 6 7 8 9 10 11 12
## 9419 5 6 6 13 8 8 7 5 5 8# All 71 missing people have their first record in the data after January, spread from February to December. None of them have a January row, so filtering to January drops them. Data doesn't show why (e.g., birth, moving, etc.)model1 <- lm(RFOODR ~ TAGE + ESEEING + ESEX, data = month_1)
# adding the variable ESEX as a predictor. ESEX is stored as 1 for male, 2 for female.
summary(model1)##
## Call:
## lm(formula = RFOODR ~ TAGE + ESEEING + ESEX, data = month_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2870 -0.5933 -0.4443 -0.3164 5.7352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8050736 0.1506301 11.983 < 2e-16 ***
## TAGE -0.0057293 0.0006019 -9.519 < 2e-16 ***
## ESEEING -0.5504952 0.0702746 -7.833 5.26e-15 ***
## ESEX 0.0305470 0.0283216 1.079 0.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.372 on 9415 degrees of freedom
## Multiple R-squared: 0.01411, Adjusted R-squared: 0.0138
## F-statistic: 44.92 on 3 and 9415 DF, p-value: < 2.2e-16# sex shows no detectable effect on food insecurity score once age and vision are held constant.
plot(model1, which = 1)
# as previously mentioned RFOODR has seven values (0-6). Each stripe is people with the same RFOODR score.
raintest(model1, order.by = ~TAGE, data = month_1)##
## Rainbow test
##
## data: model1
## Rain = 0.83665, df1 = 4710, df2 = 4705, p-value = 1# red line is roughly flat, no bend. Rainbow test = 0.84, p = 1, no evidence against linearity.
sum(month_1$RFOODR == 0) / nrow(month_1)## [1] 0.8220618# The stripes come from RFOODR being whole numbers 0-6, with 82% at zero.dwtest(model1)##
## Durbin-Watson test
##
## data: model1
## DW = 0.90611, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0# dw is low at 0.906, p-value < 2.2e-16, the assumption is violated.
month_1 %>% select(SSUID, ERESIDENCEID, PNUM, RFOODR) %>% head(20)## SSUID ERESIDENCEID PNUM RFOODR
## <i64> <int> <int> <int>
## 1: 11428534325 100001 101 6
## 2: 11428534325 100001 102 6
## 3: 11428534325 100001 103 6
## 4: 11428534325 100001 104 6
## 5: 11428534325 100001 105 6
## 6: 11486013425 100001 101 0
## 7: 11486013425 100001 102 0
## 8: 11486013425 100001 103 0
## 9: 11486013425 100001 104 0
## 10: 11486013425 100001 105 0
## 11: 41803467625 100001 101 6
## 12: 41803467625 100001 102 6
## 13: 41803742025 100001 101 0
## 14: 41803742025 100001 102 0
## 15: 41803755425 100001 101 0
## 16: 41803755425 100001 102 0
## 17: 41803758425 100001 101 0
## 18: 41806825925 100002 101 0
## 19: 41808623325 100001 101 0
## 20: 41808623325 100001 102 0
## SSUID ERESIDENCEID PNUM RFOODR
## <i64> <int> <int> <int>hh <- month_1 %>% group_by(SSUID, ERESIDENCEID) %>% summarise(values = n_distinct(RFOODR), .groups = "drop")
table(hh$values)##
## 1
## 4223# The dataset is sorted by household and RFOODR is a household measure copied to each person. This is shown with 4,223 households all with exactly 1 distinct RFOODR value. Consequently, standard errors are understated.plot(model1, which = 3)
# the red line goes from 0.4 to 1.0. Errors get bigger as the predicted score gets bigger.
bptest(model1)##
## studentized Breusch-Pagan test
##
## data: model1
## BP = 70.418, df = 3, p-value = 3.473e-15# BP = 70.418, df = 3, p-value = 3.473e-15, assumption violated.
plot(model1, which = 2)
# Most people score 0, which puts most residuals in a tight negative band, and then a few hundred people score 4, 5 or 6; normality is violated.shapiro.test(residuals(model1))## Error in `shapiro.test()`:
## ! sample size must be between 3 and 5000# Error in shapiro.test(residuals(model1)) : sample size must be between 3 and 5000vif(model1)## TAGE ESEEING ESEX
## 1.021419 1.019264 1.002735# all about 1.02, the predictors are essentially unrelated to each other
cor(month_1[, c("TAGE", "ESEEING", "ESEX")])## TAGE ESEEING ESEX
## TAGE 1.00000000 -0.13678554 0.05035019
## ESEEING -0.13678554 1.00000000 -0.02062433
## ESEX 0.05035019 -0.02062433 1.00000000# ESEEING is 1 for difficulty and 2 for no difficulty, so a negative correlation with age means older people are more likely to have vision difficulty. The largest correlation is -0.137, no multicollinearityThe model does not meet the assumptions, three of five failed. For linearity, the model passed since there was a flat red line, Rain = 0.837, p = 1. The model failed the independence assumption due to DW = 0.906 and p < 2.2e-16. For homoscedasticity, the model failed due to the red line increasing from 0.4 to 1.0, BP = 70.4, and p = 3.5e-15. For normality, the model failed because the Q-Q tail lifts off. Also, Shapiro can’t run at n = 9,419. Lastly, the model passed multicollinearity with VIF about 1.02 and max r = -0.137. The independence failure means the standard errors in the summary(model1) output are too small, so the t values are too large and the p values too small. The model looks more certain than the data supports.
month_1$hhid <- paste(month_1$SSUID, month_1$ERESIDENCEID)
# paste joins SSUID and ERESIDENCEID into one household label. The clustering function takes a single grouping variable, not a pair.
coeftest(model1, vcov = vcovCL, cluster = month_1$hhid)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.80507360 0.20164109 8.9519 < 2.2e-16 ***
## TAGE -0.00572926 0.00084347 -6.7925 1.169e-11 ***
## ESEEING -0.55049515 0.08999947 -6.1166 9.937e-10 ***
## ESEX 0.03054703 0.02554985 1.1956 0.2319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the estimates don't change, only the standard errors. Same coefficient table with standard errors that treat each household as one group instead of treating each person as independent.
month_1 %>% group_by(hhid) %>% summarise(mixed = n_distinct(ESEX) > 1, .groups = "drop") %>% summarise(mean(mixed))## # A tibble: 1 × 1
## `mean(mixed)`
## <dbl>
## 1 0.593# 59 percent of households contain both sexes, this is likely why ESEX standard error fell while the others roseThe model violates three assumptions: independence, homoscedasticity and normality. The violation of the assumption of independence is the main problem due to the model treating 9,419 rows as independent when only 4,223 outcome values are. This produces large t values and small p values. The first failure comes from 4,223 households having one value each. The other two failures are from the model predicting scores between 0.25 and 1.27. Most people score 0 so their errors are small. The few who score 5 or 6 have huge errors. That lopsidedness breaks both normality and equal spread.
To mitigate the independence violation, I used cluster-robust standard errors grouped by household. The estimates did not change. TAGE’s standard error went from 0.0006019 to 0.0008435, and its t value from -9.52 to -6.79. TAGE and ESEEING are still significant, so the conclusion held. The clustered standard errors are also robust to unequal spread, so they cover homoscedasticity as well.