Prenatal-to-3 Coding Excercise

Author

R. Luttinen

Data exercise using ACS PUMS data from 2019 and 2021

-4 data files per sample (2 for household, 2 for individual)

-stcode= state ID

-serialno= household ID

1.

#load data

library(readr)

#2019

#household files
ptone2019hh <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/hh19ne_1.csv")
Rows: 59919 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): serialno, stcode
dbl (5): hwt, n_persons, hhtype, hhinc, n_kids

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pttwo2019hh <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/hh19ne_2.csv")
Rows: 216048 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): serialno
dbl (6): stcode, hwt, n_persons, hhtype, hhinc, n_kids

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#person files
ptone2019p <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/pp19ne_1.csv")
Rows: 119650 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): SERIALNO, STCODE
dbl (8): SPORDER, PWT, AGEP, SEX, PPINC, IS_CHILD, POVPIP, RACEETH3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pttwo2019p <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/pp19ne_2.csv")
Rows: 447327 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): SERIALNO
dbl (9): SPORDER, STCODE, PWT, AGEP, SEX, PPINC, IS_CHILD, POVPIP, RACEETH3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#2021

#household files
ptone2021hh <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/hh21ne_1.csv")
Rows: 61368 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): serialno
dbl (6): stcode, hwt, n_persons, hhtype, hhinc, n_kids

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pttwo2021hh <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/hh21ne_2.csv")
Rows: 220070 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): serialno
dbl (6): stcode, hwt, n_persons, hhtype, hhinc, n_kids

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#person files

ptone2021p <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/pp21ne_1.csv")
Rows: 123656 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): serialno, stcode
dbl (8): sporder, pwt, agep, sex, ppinc, is_child, povpip, raceeth3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pttwo2021p <- read_csv("C:/Users/Rebecca/Downloads/DataExercisePacket/rawdata/pp21ne_2.csv")
Rows: 454249 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): serialno
dbl (9): sporder, stcode, pwt, agep, sex, ppinc, is_child, povpip, raceeth3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

2.

#combine household datasets into one

householdfile2019<-rbind(ptone2019hh, pttwo2019hh)

householdfile2021<-rbind(ptone2021hh, pttwo2021hh)

householdfileall<-rbind(householdfile2019,householdfile2021)

#show  layout

head(householdfileall)
# A tibble: 6 × 7
  serialno      stcode   hwt n_persons hhtype hhinc n_kids
  <chr>         <chr>  <dbl>     <dbl>  <dbl> <dbl>  <dbl>
1 2019GQ0000032 09         0         1     NA    NA     NA
2 2019GQ0000222 09         0         1     NA    NA     NA
3 2019GQ0000384 09         0         1     NA    NA     NA
4 2019GQ0000494 09         0         1     NA    NA     NA
5 2019GQ0000680 09         0         1     NA    NA     NA
6 2019GQ0000699 09         0         1     NA    NA     NA
#make a year identifier using the year listed in the serial code

library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
householdfileall$year<- substr(householdfileall$serialno, 1, 4)

class(householdfileall$year)
[1] "character"
#convert to numeric

householdfileall$year<-as.numeric(householdfileall$year)
#combine person datasets into one

personfile2019<-rbind(ptone2019p, pttwo2019p)

#make all the variable names for the 2019 person file lowercase to make appending possible


personfile2019 <- personfile2019 %>% 
  rename_with(tolower)


personfile2021<-rbind(ptone2021p, pttwo2021p)

personfileall<-rbind(personfile2019,personfile2021)

#show  layout

head(personfileall)
# A tibble: 6 × 10
  serialno      sporder stcode   pwt  agep   sex ppinc is_child povpip raceeth3
  <chr>           <dbl> <chr>  <dbl> <dbl> <dbl> <dbl>    <dbl>  <dbl>    <dbl>
1 2019GQ0000032       1 09         9    20     1     0       NA     NA        1
2 2019GQ0000222       1 09        53    20     1  1200       NA     NA        1
3 2019GQ0000384       1 09        44    21     2  2000       NA     NA        1
4 2019GQ0000494       1 09        11    87     2     0       NA     NA        1
5 2019GQ0000680       1 09        11    95     2     0       NA     NA        1
6 2019GQ0000699       1 09        79    21     2  1000       NA     NA        1
#make a year identifier using the year listed in the serial code

library(dplyr)

personfileall$year<- substr(personfileall$serialno, 1, 4)

class(householdfileall$year)
[1] "numeric"
#convert to numeric

personfileall$year<-as.numeric(personfileall$year)

3.

#remove any observations that are not families

class(householdfileall$hhtype)
[1] "numeric"
householdfileallfilter<-householdfileall%>%
  filter(hhtype<4)

#double check

unique(householdfileallfilter$hhtype)
[1] 1 3 2

4.

#standardize income: $1.00 in 2019 is equivalent to $1.06 in 2021


#perform conversion for 2021; use case_when

householdfileallfilter<-householdfileallfilter%>%
  mutate(adjustedincome=case_when(year==2021~ hhinc/1.06,
                                  TRUE ~ hhinc))

5.

Note: number six is missing, instead number 5 has two parts.

#cap the adjusted income variable at the 99th percentile, opting to ignore missing vluaes

householdfileallfilter_capped<-householdfileallfilter%>%
  filter(adjustedincome<= quantile(adjustedincome, 0.99, na.rm=TRUE))

#did family income change bewtween 2019 an 2021?

householdfileallfilter_capped2019<-householdfileallfilter_capped%>%
  filter(year==2019)


householdfileallfilter_capped2021<-householdfileallfilter_capped%>%
  filter(year==2021)

mean(householdfileallfilter_capped2019$adjustedincome)
[1] 124075.8
mean(householdfileallfilter_capped2021$adjustedincome)
[1] 121037.6
#check for statistical significance using a t-test

ttest<-t.test(adjustedincome~year, data=householdfileallfilter_capped)

ttest

    Welch Two Sample t-test

data:  adjustedincome by year
t = 7.757, df = 289738, p-value = 8.723e-15
alternative hypothesis: true difference in means between group 2019 and group 2021 is not equal to 0
95 percent confidence interval:
 2270.543 3805.876
sample estimates:
mean in group 2019 mean in group 2021 
          124075.8           121037.6 
124075.8- 121037.6 
[1] 3038.2

The t-test confirms that the change in mean adjusted income from 2019 to 2021 is significant. The mean income was $124,075.80 in 2019 and $121,037.60 in 2021, a drop of 3,038 dollars.

Note: I’m going to continue working with the dataframe that dropped any income values above the 99th percentile.

7.

#merge adjusted income variable with the person-level data so each person has their adjusted household income

#create dataframe with only household identifier and adjusted income

hhandid<-householdfileallfilter_capped%>%
  select(serialno, adjustedincome)


personfileallincome<-merge(personfileall, hhandid, by= ("serialno"))

#note: this dataframe is not going to ONLY have information for individuals living in a household considered a family household that have an income not above the 99th percentile of the original family income data. 

#show  layout

head(personfileallincome)
       serialno sporder stcode pwt agep sex ppinc is_child povpip raceeth3 year
1 2019HU0000007       3     09 145    6   2    NA        1    501        1 2019
2 2019HU0000007       1     09 139   35   2 85000        0    501        1 2019
3 2019HU0000007       2     09 140   38   1 45000        0    501        1 2019
4 2019HU0000007       4     09 212    4   2    NA        1    501        1 2019
5 2019HU0000008       1     36 121   36   2 58000        0    501       NA 2019
6 2019HU0000008       3     36 156    4   1    NA        1    501        1 2019
  adjustedincome
1         130000
2         130000
3         130000
4         130000
5         273000
6         273000

8.

#make a variable that indicates whether a person lives in a household the earns less than 150% of the federal poverty level

#check distribution

summary(personfileallincome$povpip)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0   252.0   449.0   371.6   501.0   501.0    1870 
personfileallincome<-personfileallincome%>%
  mutate(below150=(povpip<150))

#this makes TRUE/FALSE variable in which TRUE indicates that individuals live in a household that earns less than 150% of FPL

#see distribution

library(janitor)

Attaching package: 'janitor'
The following objects are masked from 'package:stats':

    chisq.test, fisher.test
tabyl(personfileallincome, below150)
 below150      n     percent valid_percent
    FALSE 764437 0.872437823     0.8743038
     TRUE 109901 0.125427981     0.1256962
       NA   1870 0.002134196            NA
#note, I haven't dropped any NA's from this dataframe as the instructions have not indicated whether I should. Eventually I would do this during my data cleaning process. This table shows that 1870 values for my newly created varaible are NA's.

9.

#among children aged 3 or younger, how did the percent of children living below 150% poverty change between 2019 and 2021?
#filter for only observations of children aged 3 or younger


personfileallincomebelow3<-personfileallincome%>%
  filter(agep<=3)

#the is_child variable is equal to 0 then the specific is not a child, but if it is equal to 1 it is a child.


tabyl(personfileallincomebelow3, is_child)
 is_child     n   percent
        0  4790 0.1221253
        1 34432 0.8778747
#since the question refers to children only, I'm going to filter for is_child==1 because some of these cases may be adults with children

childrenonly<-personfileallincomebelow3%>%
  filter(is_child==1)
#create visualization of how the percent of children living below poverty changed from 2019 and 2021



visualinfo<-childrenonly%>%
  group_by(year) %>%
  summarize(
    total_children = n(),
    below150 = sum(below150, na.rm = TRUE),
    percentlivingbelow150 = mean(below150/ total_children, na.rm = TRUE) * 100
  )

head(visualinfo)
# A tibble: 2 × 4
   year total_children below150 percentlivingbelow150
  <dbl>          <int>    <int>                 <dbl>
1  2019          17177     3609                  21.0
2  2021          17255     3527                  20.4
#make a bar graph

library(ggplot2)


ggplot(visualinfo, aes(x = factor(year), y = percentlivingbelow150, fill = factor(year))) +
  geom_col(width = 0.5, color = "black", alpha = 0.8) +
  geom_text(aes(label = paste0(round(percentlivingbelow150, 1), "%")), 
            vjust = -0.5, fontface = "bold", size = 5) +
  scale_fill_manual(values = c("2019" = "orange", "2021" = "blue")) +
  # Sets the axis limit slightly higher than the max percentage for spacing
  scale_y_continuous(labels = function(x) paste0(x, "%"), limits = c(0, 30)) +
  labs(
    title = "% of children living in a household with an income below 150% of the federal poverty line",
    x = "Year",
    y = "% of children"
  ) +
  theme_minimal(base_size = 8) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )

10.

#for each state estiamte the proportion of white, black and hispanic children aged 3 or younger who lived below 150% of poverty in 2021 using survey-weights

library(srvyr)

Attaching package: 'srvyr'
The following object is masked from 'package:stats':

    filter
#apply survey design

survey_design <- childrenonly %>%
  as_survey_design(weights = pwt)



state_level_pov <- survey_design %>%
  group_by(stcode, raceeth3) %>%
  summarise(
    prop_below_150 = survey_mean(below150, na.rm = TRUE)
  ) %>% 
  mutate(
    poverty_rate_pct = prop_below_150 * 100
  )


print(state_level_pov)
# A tibble: 36 × 5
# Groups:   stcode [9]
   stcode raceeth3 prop_below_150 prop_below_150_se poverty_rate_pct
   <chr>     <dbl>          <dbl>             <dbl>            <dbl>
 1 09            1          0.104           0.0126              10.4
 2 09            2          0.422           0.0524              42.2
 3 09            3          0.429           0.0319              42.9
 4 09           NA          0.185           0.0308              18.5
 5 23            1          0.250           0.0234              25.0
 6 23            2          0.972           0.0290              97.2
 7 23            3          0.302           0.168               30.2
 8 23           NA          0.231           0.0664              23.1
 9 25            1          0.100           0.00829             10.0
10 25            2          0.335           0.0464              33.5
# ℹ 26 more rows
#note, I still haven't dropped NA's so they are included in this list.

11.

#make a variable that shows whether a household contains a person aged 3 or younger


#use the orignal all person with adjusted income here

personfileallincomegroupedbyhh<-personfileallincome%>%
  mutate(child3oryounger= agep<=3)