Exercise 4
Mengni Chen
Sep-09-2026
Warming up
Use pairfam wave1 dataset, please check the following variables in the questionnaire and codebook:
- age,
- sex_gen,
- sat6,
- sat1i1,
- sat1i4
- relstat,
- nkidsbio,
- cohort,
- val1i7
Question
- Which variable measures satisfaction with family?
- Which variable measures attitude towards marriage?
Answer for a.
sat1i4
Answer for b.
val1i7
No. 1
Question
1.1 Use pairfam wave1 dataset, compile a dataset consisting of variables stated in Warming up: 1) age, 2) sex_gen, 3) sat6, 4) sat1i1, 5) sat1i4 6) relstat, 7) nkidsbio, 8) cohort, 9) val1i7
1.2 Check the variables to see any of them having invalid answers, and the number of missingness for each variable.
1.3 Drop observations with missing values.
1.4 What is the original sample size and what is the sample size after cleaning?
Answer 1.1
library(tidyverse) # Add the tidyverse package to my current library.
library(haven) # Import data.
library(janitor) #cleaning data
library(estimatr) # Allows us to estimate (cluster-)robust standard errors.
library(texreg) # Allows us to make nicely-formatted Html & Latex regression tables.pairfam <- read_dta("anchor1_50percent_Eng.dta")
dataset1 <- select(pairfam,
age,
sex_gen,
sat6,
sat1i1,
sat1i4,
relstat,
nkidsbio,
cohort,
val1i7)Answer 1.2
tabyl(dataset1,age)## age n percent
## 14 41 0.006611837
## 15 708 0.114175133
## 16 722 0.116432833
## 17 667 0.107563296
## 18 35 0.005644251
## 24 24 0.003870343
## 25 577 0.093049508
## 26 678 0.109337204
## 27 647 0.104338010
## 28 87 0.014029995
## 34 22 0.003547815
## 35 502 0.080954685
## 36 618 0.099661345
## 37 772 0.124496049
## 38 101 0.016287696tabyl(dataset1,sex_gen)## sex_gen n percent
## 1 3029 0.4884696
## 2 3172 0.5115304tabyl(dataset1,sat6) # 5 cases with invalid answer## sat6 n percent
## -2 2 0.0003225286
## -1 3 0.0004837929
## 0 26 0.0041928721
## 1 18 0.0029027576
## 2 45 0.0072568940
## 3 110 0.0177390743
## 4 133 0.0214481535
## 5 395 0.0636994033
## 6 508 0.0819222706
## 7 1178 0.1899693598
## 8 1877 0.3026931140
## 9 1157 0.1865828092
## 10 749 0.1207869698tabyl(dataset1,sat1i1) #23 cases with invalid answer## sat1i1 n percent
## -2 13 0.002096436
## -1 10 0.001612643
## 0 103 0.016610224
## 1 50 0.008063216
## 2 89 0.014352524
## 3 208 0.033542977
## 4 219 0.035316884
## 5 593 0.095629737
## 6 494 0.079664570
## 7 1154 0.186099016
## 8 1490 0.240283825
## 9 773 0.124657313
## 10 1005 0.162070634tabyl(dataset1,sat1i4) #9 cases with invalid answer## sat1i4 n percent
## -2 3 0.0004837929
## -1 6 0.0009675859
## 0 28 0.0045154007
## 1 21 0.0033865506
## 2 55 0.0088695372
## 3 70 0.0112885019
## 4 80 0.0129011450
## 5 237 0.0382196420
## 6 271 0.0437026286
## 7 486 0.0783744557
## 8 901 0.1452991453
## 9 1181 0.1904531527
## 10 2862 0.4615384615tabyl(dataset1,relstat) #34 cases with invalid answer## relstat n percent
## -7 34 0.0054829866
## 1 2448 0.3947750363
## 2 1012 0.1631994840
## 3 660 0.1064344461
## 4 1735 0.2797935817
## 5 23 0.0037090792
## 6 146 0.0235445896
## 7 63 0.0101596517
## 8 76 0.0122560877
## 9 3 0.0004837929
## 10 1 0.0001612643#or
tabyl(as_factor(dataset1$relstat))## as_factor(dataset1$relstat) n percent
## -7 Incomplete data 34 0.0054829866
## 1 Never married single 2448 0.3947750363
## 2 Never married LAT 1012 0.1631994840
## 3 Never married COHAB 660 0.1064344461
## 4 Married COHAB 1735 0.2797935817
## 5 Married noncohabiting 23 0.0037090792
## 6 Divorced/separated single 146 0.0235445896
## 7 Divorced/separated LAT 63 0.0101596517
## 8 Divorced/separated COHAB 76 0.0122560877
## 9 Widowed single 3 0.0004837929
## 10 Widowed LAT 1 0.0001612643
## 11 Widowed COHAB 0 0.0000000000tabyl(dataset1,nkidsbio) #4 cases with invalid answer## nkidsbio n percent
## -7 4 0.0006450572
## 0 4136 0.6669891953
## 1 839 0.1353007579
## 2 868 0.1399774230
## 3 277 0.0446702145
## 4 54 0.0087082729
## 5 17 0.0027414933
## 6 4 0.0006450572
## 10 2 0.0003225286tabyl(dataset1,cohort) ## cohort n percent
## 1 2173 0.3504274
## 2 2013 0.3246251
## 3 2015 0.3249476tabyl(as_factor(dataset1$val1i7)) #63cases with invalid answer## as_factor(dataset1$val1i7) n percent
## -5 Inconsistent value 0 0.000000000
## -4 Filter error / Incorrect entry 0 0.000000000
## -3 Does not apply 0 0.000000000
## -2 No answer 10 0.001612643
## -1 Don't know 53 0.008547009
## 1 Disagree completely 913 0.147234317
## 2 800 0.129011450
## 3 1328 0.214159007
## 4 1152 0.185776488
## 5 Agree completely 1945 0.313659087Answer 1.3
dataset2 <- dataset1 %>%
transmute(
age,
sex=as_factor(sex_gen) %>% fct_drop(),
#treat sex as a categorical, and drop unused levels
sat6=case_when(sat6<0 ~ as.numeric(NA),
TRUE ~ as.numeric(sat6)),
#specify when sat should be considered missing
work_sat=case_when(sat1i1<0 ~ as.numeric(NA),
TRUE ~ as.numeric(sat1i1)),
#specify when work_sat should be considered missing
fam_sat=case_when(sat1i4<0 ~ as.numeric(NA),
TRUE ~ as.numeric(sat1i4)),
#specify when work_sat should be considered missing
relstat=as_factor(relstat), #treat relationship status as categorical
relstat1=case_when(relstat=="-7 Incomplete data" ~ as.character(NA),
TRUE ~ as.character(relstat)) %>%
#specify when relstat1 should be missing
as_factor() %>% fct_drop(),#make relstat1 as a factor again & drop unused levels
nkidsbio=case_when(nkidsbio<0 ~ as.numeric(NA),
TRUE ~ as.numeric(nkidsbio)),
#specify when nkidsbio should be missing
cohort=as_factor(cohort) %>% fct_drop(),
#treat sex as a categorical, and drop unused levels
mar_att=as_factor(val1i7),
mar_att1=case_when(mar_att %in% c("-2 No answer","-1 Don't know") ~ as.character(NA),
TRUE ~ as.character(mar_att)%>%
#specify when val1i5 should be missing
as_factor() %>% fct_drop()#make mar_att1 as a factor again & drop unused levels
)
) %>%
drop_na() #remove all observations that are missingAnswer 1.4
count(dataset1) ## # A tibble: 1 × 1
## n
## <int>
## 1 6201count(dataset2)## # A tibble: 1 × 1
## n
## <int>
## 1 6065original sample size is 6201; the clean sampled size is 6065.
No. 2
Question
Generate a new variable “marital” for relationship status with categories of “Never married”, “Married”, “Divorced”(including separated), “Widowed”, regardless of whether individuals are cohabiting, noncohabiting, LAT.
Generate a new variable “parenthood” based on nkidsbio with categories of “Have kids”, “No kids”.
Answer
tabyl(dataset2,relstat1)#check the distribution first.## relstat1 n percent
## 1 Never married single 2403 0.3962077494
## 4 Married COHAB 1703 0.2807914262
## 2 Never married LAT 1000 0.1648804617
## 3 Never married COHAB 654 0.1078318219
## 8 Divorced/separated COHAB 75 0.0123660346
## 9 Widowed single 3 0.0004946414
## 6 Divorced/separated single 141 0.0232481451
## 7 Divorced/separated LAT 62 0.0102225886
## 5 Married noncohabiting 23 0.0037922506
## 10 Widowed LAT 1 0.0001648805dataset3 <- dataset2 %>%
mutate(
marital=case_when(
relstat1 %in% c("1 Never married single",
"2 Never married LAT",
"3 Never married COHAB") ~ "Nevermarried",
# when relstat1 has any of the three situations, I assign "Nevermarried" to new variable "marital"
relstat1 %in% c("4 Married COHAB",
"5 Married noncohabiting") ~ 'Married',
# when relstat1 has any of the two situations, I assign "Married" to new variable "marital
relstat1 %in% c("6 Divorced/separated single",
"7 Divorced/separated LAT",
"8 Divorced/separated COHAB") ~ 'Divorced',
# when relstat1 has any of the three situations, I assign "Divorced" to new variable "marital"
relstat1 %in% c("9 Widowed single","10 Widowed LAT") ~ 'Widowed'
# when relstat1 has any of the two situations, I assign "Widow" to new variable "marital"
) %>% as_factor(),#make marital a categorical variable
parenthood=case_when(
nkidsbio>0 ~ "Have kids",#when nkidsbio should be "Have kids"
nkidsbio==0 ~ "No kids" #when nkidsbio should be "No kids"
) %>%as_factor() #make parenthood a categorical variable
)No. 3
Question
Estimate the average life satisfaction and average satifaction with family for the married and the divorced individuals.
Answer
dataset4 <- dataset3%>%
filter(marital!= "Widowed" & marital!= "Nevermarried")
#drop the widowed and never-married
mar_sat <- dataset4 %>%
group_by(marital) %>%
summarise(
mean_sat6=mean(sat6), #calculate the mean of sat6 by marital1
mean_famsat=mean(fam_sat) )
mar_sat## # A tibble: 2 × 3
## marital mean_sat6 mean_famsat
## <fct> <dbl> <dbl>
## 1 Married 7.80 8.98
## 2 Divorced 6.68 7.63No. 4
Question
Generate a table to show the distribution of attitudes toward marriage for the married and the divorced individuals.
Then use a bar chart to show the distributions of attitudes the married and the divorced individuals.
Answer a.
dataset4$marital <- fct_drop(dataset4$marital) #drop unused levels
distribute_table<- tabyl(dataset4, mar_att1, marital) %>%
adorn_totals("row") %>%#add row %
adorn_percentages("col") %>% #add col
adorn_pct_formatting() #format the percentage
distribute_table## mar_att1 Married Divorced
## 1 Disagree completely 12.1% 39.9%
## 2 12.0% 15.8%
## 3 20.9% 22.7%
## 4 16.9% 9.4%
## 5 Agree completely 38.1% 12.2%
## Total 100.0% 100.0%Answer a.
plot1<- ggplot(data = dataset4,
mapping=aes(x = mar_att1,
y = ..prop..,# ask R to show % rather than absolute count
group = marital,#calculate the % of sat6 within each group of relstat_new
fill = marital)#color the bars different by relstat_new
)+
geom_bar()+
facet_wrap(~ marital)+
coord_flip()+ #swap the coordinate system to make a horizontal bar, so that the labels are clear
labs(title="Attitude: Marriage is a lifelong union that should not be broken")
plot1
No. 5
Question
I want to use OLS regression to examine whether the divorced people have lower life satisfaction than the married ones. Please use standard error robust OLS to do the modelling
5.1 Model1: regress life satisfaction on age and marital status.
5.2 What is your interpretation the coefficient of age and marital?
5.3 output the result to a html file.
Answer 5.1
regression1 <- lm_robust(data = dataset4,
formula = sat6 ~ age + marital )
summary(regression1)##
## Call:
## lm_robust(formula = sat6 ~ age + marital, data = dataset4)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 7.951356 0.299248 26.5711 1.705e-133 7.36449 8.53823 2001
## age -0.004529 0.008814 -0.5139 6.074e-01 -0.02181 0.01276 2001
## maritalDivorced -1.119185 0.134549 -8.3180 1.632e-16 -1.38306 -0.85531 2001
##
## Multiple R-squared: 0.04723 , Adjusted R-squared: 0.04628
## F-statistic: 35.11 on 2 and 2001 DF, p-value: 1.029e-15Answer 5.2
Interpreting the age coefficient: When the marital status is controlled, with one year increase in age, the life satisfaction decreased by 0.0045. But this age effect is not statistically significant at 5%.
Interpreting the marital coefficient: When the age is controlled, a divorced individuals has lower life satisfaction than a married individual by amount of 1.1192 score. This effect of marital status is statistically significant.
Answer 5.3
# I export it into a html file. You can do whatever you want, say, word, or excel.
texreg::htmlreg(regression1,
include.ci = FALSE,
digits = 4,
file = "output1.html") #html## The table was written to the file 'output1.html'.
No. 6
Question
6.1 Please change the reference category of marital, using “divorced” as the reference in the regression conducted in Question 5.
6.2 Add a new variable, “parenthood” to your model in Question 6.1. Is there any change in the coefficents of age and marital?
Answer 6.1
dataset4$marital <- fct_relevel(dataset4$marital, "Divorced")
regression1a <- lm_robust(data = dataset4,
formula = sat6 ~ age + marital )
summary(regression1a)##
## Call:
## lm_robust(formula = sat6 ~ age + marital, data = dataset4)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 6.832171 0.330459 20.6748 3.190e-86 6.18409 7.48025 2001
## age -0.004529 0.008814 -0.5139 6.074e-01 -0.02181 0.01276 2001
## maritalMarried 1.119185 0.134549 8.3180 1.632e-16 0.85531 1.38306 2001
##
## Multiple R-squared: 0.04723 , Adjusted R-squared: 0.04628
## F-statistic: 35.11 on 2 and 2001 DF, p-value: 1.029e-15Answer 6.2
regression2 <- lm_robust(data = dataset4,
formula = sat6 ~ age + marital + parenthood )
summary(regression2)##
## Call:
## lm_robust(formula = sat6 ~ age + marital + parenthood, data = dataset4)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 6.831537 0.330525 20.6688 3.568e-86 6.18333 7.47975
## age -0.005542 0.009384 -0.5906 5.549e-01 -0.02394 0.01286
## maritalMarried 1.117001 0.134529 8.3031 1.844e-16 0.85317 1.38083
## parenthoodHave kids 0.044252 0.112518 0.3933 6.942e-01 -0.17641 0.26492
## DF
## (Intercept) 2000
## age 2000
## maritalMarried 2000
## parenthoodHave kids 2000
##
## Multiple R-squared: 0.04731 , Adjusted R-squared: 0.04588
## F-statistic: 23.39 on 3 and 2000 DF, p-value: 7.086e-15The coefficient of “married” gets slightly smaller. The coefficient of age gets slightly more negative.
No. 7
Question
Using standard error robust OLS to do the following models
7.1 Model1: regress life satisfaction with family on age, marital status, parenthood status, and sex. Using “Married” as the reference group.
7.2 Export results of Model 1 to a html file and rename the names of coefficients in a clearer way to understand: specifically customized the names to be “Intercept”, “Age”, “Divorced (Ref.=Married)”, “Have kids (Ref.=no kid)”,“Female (Ref.=Male)”
Answer 7.1
dataset4$marital <- fct_relevel(dataset4$marital, "Married")
model1 <- lm_robust(data = dataset4,
formula = fam_sat ~ age + marital + parenthood + sex )
summary(model1 )##
## Call:
## lm_robust(formula = fam_sat ~ age + marital + parenthood + sex,
## data = dataset4)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 9.36688 0.314726 29.762 1.937e-161 8.74966 9.984107
## age -0.02684 0.009239 -2.905 3.714e-03 -0.04496 -0.008719
## maritalDivorced -1.32305 0.161993 -8.167 5.524e-16 -1.64074 -1.005355
## parenthoodHave kids 0.53755 0.123867 4.340 1.498e-05 0.29463 0.780474
## sex2 Female 0.11903 0.080757 1.474 1.406e-01 -0.03934 0.277411
## DF
## (Intercept) 1999
## age 1999
## maritalDivorced 1999
## parenthoodHave kids 1999
## sex2 Female 1999
##
## Multiple R-squared: 0.08208 , Adjusted R-squared: 0.08024
## F-statistic: 23.07 on 4 and 1999 DF, p-value: < 2.2e-16Answer 7.2
texreg::htmlreg(
list(model1),
custom.coef.names = c("Intercept",
"Age",
"Divorced (Ref.=Married)",
"Have kids (Ref.=no kid)",
"Female (Ref.=Male)"
),
include.ci = FALSE, digits = 4,
file = "models_rename.html") #html## The table was written to the file 'models_rename.html'.