HW 2 B. Solomon
Bryan Solomon
2023-09-11
library(purrr)
library(haven)
library(janitor)##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(plyr)##
## Attaching package: 'plyr'## The following object is masked from 'package:purrr':
##
## compactlibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discardlibrary(sur)
library(summarytools)
library(Rmisc)## Loading required package: latticelibrary(car)## Loading required package: carData##
## Attaching package: 'carData'## The following objects are masked from 'package:sur':
##
## Anscombe, States##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somelibrary(usmap)
library(ggthemes)
library(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ readr 2.1.4 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::arrange() masks plyr::arrange()
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ plyr::compact() masks purrr::compact()
## ✖ dplyr::count() masks plyr::count()
## ✖ dplyr::desc() masks plyr::desc()
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::failwith() masks plyr::failwith()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::id() masks plyr::id()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::mutate() masks plyr::mutate()
## ✖ car::recode() masks dplyr::recode()
## ✖ dplyr::rename() masks plyr::rename()
## ✖ car::some() masks purrr::some()
## ✖ dplyr::summarise() masks plyr::summarise()
## ✖ dplyr::summarize() masks plyr::summarize()
## ✖ tibble::view() masks summarytools::view()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(tidycensus)
library("ggpubr")##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:plyr':
##
## mutatelibrary(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary("writexl")
library(readxl)
library("tinytex")
library("lubridate")
library(splitstackshape)
library(tinytex)
library(psych)##
## Attaching package: 'psych'
##
## The following object is masked from 'package:car':
##
## logit
##
## The following object is masked from 'package:sur':
##
## skew
##
## The following objects are masked from 'package:scales':
##
## alpha, rescale
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(glmmTMB)
library(readxl)
library(openxlsx)library(readr)
X2017_2019_Male_Data_2 <-read_csv("C:/Users/bryan/OneDrive/Desktop/2017_2019_Male_Data_2.csv")## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)## Rows: 5206 Columns: 3009
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (1939): caseid, rscrninf, rscrage, rscrhisp, rscrrace, age_a, age_r, age...
## lgl (1070): sedbclc3, sedbclc4, sedwhlc3, sedwhlc4, sedconlc4, p2currwife, p...
##
## ℹ 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.#Subset data
shortdata <- X2017_2019_Male_Data_2[, c("age_r", "hisp", "fmarit", "hisagem", "agemarrn")]#The shortdata set was creased to cut down from the entire NSFG male respondent data set. The variables that were pulled were age_r, hisp, fmarit, hisagem, and agemarrn.
#censored in life table, interviwed but not married #create variables for marriage status
dat <- shortdata
# Define ageevent for evmarried based on your criteria (replace with the appropriate variable)
evmarried <- filter(dat, fmarit %in% c(1, 2, 3, 4))
evmarried$ageevent <- evmarried$agemarrn # Replace with the correct variable
# Define ageevent for nevmarried based on your criteria (replace with the appropriate variable)
nevmarried <- filter(dat, fmarit == 5)
nevmarried$ageevent <- nevmarried$age_r # Replace with the correct variable
# Combine data
atriskW4 <- rbind(evmarried, nevmarried)
# Now you can create your tables
tabyl(evmarried$ageevent)## evmarried$ageevent n percent valid_percent
## 16 3 0.0015298317 0.004862237
## 17 10 0.0050994391 0.016207455
## 18 28 0.0142784294 0.045380875
## 19 42 0.0214176441 0.068071313
## 20 49 0.0249872514 0.079416532
## 21 80 0.0407955125 0.129659643
## 22 46 0.0234574197 0.074554295
## 23 48 0.0244773075 0.077795786
## 24 58 0.0295767466 0.094003241
## 25 51 0.0260071392 0.082658023
## 26 25 0.0127485977 0.040518639
## 27 31 0.0158082611 0.050243112
## 28 25 0.0127485977 0.040518639
## 29 30 0.0152983172 0.048622366
## 30 13 0.0066292708 0.021069692
## 31 15 0.0076491586 0.024311183
## 32 9 0.0045894952 0.014586710
## 33 11 0.0056093830 0.017828201
## 34 9 0.0045894952 0.014586710
## 35 6 0.0030596634 0.009724473
## 36 6 0.0030596634 0.009724473
## 37 6 0.0030596634 0.009724473
## 38 5 0.0025497195 0.008103728
## 40 3 0.0015298317 0.004862237
## 41 2 0.0010198878 0.003241491
## 42 1 0.0005099439 0.001620746
## 43 1 0.0005099439 0.001620746
## 45 1 0.0005099439 0.001620746
## 98 3 0.0015298317 0.004862237
## NA 1344 0.6853646099 NAtabyl(nevmarried$ageevent)## nevmarried$ageevent n percent
## 15 186 0.057425131
## 16 186 0.057425131
## 17 223 0.068848410
## 18 229 0.070700834
## 19 206 0.063599877
## 20 115 0.035504785
## 21 128 0.039518370
## 22 127 0.039209633
## 23 120 0.037048472
## 24 128 0.039518370
## 25 121 0.037357209
## 26 104 0.032108676
## 27 117 0.036122260
## 28 129 0.039827107
## 29 89 0.027477617
## 30 109 0.033652362
## 31 94 0.029021303
## 32 99 0.030564989
## 33 89 0.027477617
## 34 61 0.018832973
## 35 62 0.019141710
## 36 51 0.015745600
## 37 46 0.014201914
## 38 49 0.015128126
## 39 50 0.015436863
## 40 41 0.012658228
## 41 35 0.010805804
## 42 42 0.012966965
## 43 27 0.008335906
## 44 27 0.008335906
## 45 26 0.008027169
## 46 32 0.009879592
## 47 33 0.010188330
## 48 28 0.008644643
## 49 30 0.009262118atriskW4 <- shortdataCreate a new variable ‘evmarW4’ based on ‘fmarit’
atriskW4$evmarW4 <- ifelse(atriskW4$fmarit %in% c(1, 2, 3, 4), 1, 0)#The table shows a summary of the distribution of ‘ageevement’ values within the ‘evmarried’ subset of the data. There are 1,344 missing observations represented with NA in the ‘evmarried’ variable. The ages with the highest values are 18 to 25.
#evmarW4
atriskW4<-rbind(evmarried,nevmarried)
tabyl(evmarried$ageevent)## evmarried$ageevent n percent valid_percent
## 16 3 0.0015298317 0.004862237
## 17 10 0.0050994391 0.016207455
## 18 28 0.0142784294 0.045380875
## 19 42 0.0214176441 0.068071313
## 20 49 0.0249872514 0.079416532
## 21 80 0.0407955125 0.129659643
## 22 46 0.0234574197 0.074554295
## 23 48 0.0244773075 0.077795786
## 24 58 0.0295767466 0.094003241
## 25 51 0.0260071392 0.082658023
## 26 25 0.0127485977 0.040518639
## 27 31 0.0158082611 0.050243112
## 28 25 0.0127485977 0.040518639
## 29 30 0.0152983172 0.048622366
## 30 13 0.0066292708 0.021069692
## 31 15 0.0076491586 0.024311183
## 32 9 0.0045894952 0.014586710
## 33 11 0.0056093830 0.017828201
## 34 9 0.0045894952 0.014586710
## 35 6 0.0030596634 0.009724473
## 36 6 0.0030596634 0.009724473
## 37 6 0.0030596634 0.009724473
## 38 5 0.0025497195 0.008103728
## 40 3 0.0015298317 0.004862237
## 41 2 0.0010198878 0.003241491
## 42 1 0.0005099439 0.001620746
## 43 1 0.0005099439 0.001620746
## 45 1 0.0005099439 0.001620746
## 98 3 0.0015298317 0.004862237
## NA 1344 0.6853646099 NA#this is number of censored for each age, match and use formulas for life tble
tabyl(nevmarried$ageevent)## nevmarried$ageevent n percent
## 15 186 0.057425131
## 16 186 0.057425131
## 17 223 0.068848410
## 18 229 0.070700834
## 19 206 0.063599877
## 20 115 0.035504785
## 21 128 0.039518370
## 22 127 0.039209633
## 23 120 0.037048472
## 24 128 0.039518370
## 25 121 0.037357209
## 26 104 0.032108676
## 27 117 0.036122260
## 28 129 0.039827107
## 29 89 0.027477617
## 30 109 0.033652362
## 31 94 0.029021303
## 32 99 0.030564989
## 33 89 0.027477617
## 34 61 0.018832973
## 35 62 0.019141710
## 36 51 0.015745600
## 37 46 0.014201914
## 38 49 0.015128126
## 39 50 0.015436863
## 40 41 0.012658228
## 41 35 0.010805804
## 42 42 0.012966965
## 43 27 0.008335906
## 44 27 0.008335906
## 45 26 0.008027169
## 46 32 0.009879592
## 47 33 0.010188330
## 48 28 0.008644643
## 49 30 0.009262118atriskW4$censor<- recode(atriskW4$fmarit, "5=1; else=0")
atriskW4$married<- recode(atriskW4$fmarit, "1:4=1; else = 0")
table(atriskW4$ageevent,atriskW4$fmarit)##
## 1 2 3 4 5
## 15 0 0 0 0 186
## 16 2 0 1 0 186
## 17 5 0 4 1 223
## 18 13 3 11 1 229
## 19 18 0 17 7 206
## 20 24 0 20 5 115
## 21 37 2 37 4 128
## 22 15 0 27 4 127
## 23 19 1 23 5 120
## 24 28 2 20 8 128
## 25 9 1 30 11 121
## 26 6 2 17 0 104
## 27 10 0 14 7 117
## 28 2 0 17 6 129
## 29 10 1 13 6 89
## 30 2 0 8 3 109
## 31 3 0 9 3 94
## 32 2 0 6 1 99
## 33 2 1 5 3 89
## 34 2 1 4 2 61
## 35 0 0 6 0 62
## 36 1 0 5 0 51
## 37 1 0 4 1 46
## 38 0 0 3 2 49
## 39 0 0 0 0 50
## 40 1 0 1 1 41
## 41 0 1 1 0 35
## 42 0 0 1 0 42
## 43 0 0 1 0 27
## 44 0 0 0 0 27
## 45 0 0 1 0 26
## 46 0 0 0 0 32
## 47 0 0 0 0 33
## 48 0 0 0 0 28
## 49 0 0 0 0 30
## 98 1 0 1 1 0#The table above is a contingency table showing the distribution of ‘ageevent’ values based on the ‘fmarit’ or marital status variable. The columns in this table are split into five categories that correspond with the different marital status categories in the original code book from NSFG 2017 to 2019 for male respondents.
atriskW4$censor<- recode(atriskW4$fmarit, "5=1; else=0")
table(atriskW4$ageevent,atriskW4$fmarit)##
## 1 2 3 4 5
## 15 0 0 0 0 186
## 16 2 0 1 0 186
## 17 5 0 4 1 223
## 18 13 3 11 1 229
## 19 18 0 17 7 206
## 20 24 0 20 5 115
## 21 37 2 37 4 128
## 22 15 0 27 4 127
## 23 19 1 23 5 120
## 24 28 2 20 8 128
## 25 9 1 30 11 121
## 26 6 2 17 0 104
## 27 10 0 14 7 117
## 28 2 0 17 6 129
## 29 10 1 13 6 89
## 30 2 0 8 3 109
## 31 3 0 9 3 94
## 32 2 0 6 1 99
## 33 2 1 5 3 89
## 34 2 1 4 2 61
## 35 0 0 6 0 62
## 36 1 0 5 0 51
## 37 1 0 4 1 46
## 38 0 0 3 2 49
## 39 0 0 0 0 50
## 40 1 0 1 1 41
## 41 0 1 1 0 35
## 42 0 0 1 0 42
## 43 0 0 1 0 27
## 44 0 0 0 0 27
## 45 0 0 1 0 26
## 46 0 0 0 0 32
## 47 0 0 0 0 33
## 48 0 0 0 0 28
## 49 0 0 0 0 30
## 98 1 0 1 1 0# Recode variables
atriskW4$censor <- recode(atriskW4$fmarit, "5=1; else=0")
atriskW4$married <- recode(atriskW4$fmarit, "1:4=1; else=0")
# Create tables
table(atriskW4$ageevent, atriskW4$fmarit)##
## 1 2 3 4 5
## 15 0 0 0 0 186
## 16 2 0 1 0 186
## 17 5 0 4 1 223
## 18 13 3 11 1 229
## 19 18 0 17 7 206
## 20 24 0 20 5 115
## 21 37 2 37 4 128
## 22 15 0 27 4 127
## 23 19 1 23 5 120
## 24 28 2 20 8 128
## 25 9 1 30 11 121
## 26 6 2 17 0 104
## 27 10 0 14 7 117
## 28 2 0 17 6 129
## 29 10 1 13 6 89
## 30 2 0 8 3 109
## 31 3 0 9 3 94
## 32 2 0 6 1 99
## 33 2 1 5 3 89
## 34 2 1 4 2 61
## 35 0 0 6 0 62
## 36 1 0 5 0 51
## 37 1 0 4 1 46
## 38 0 0 3 2 49
## 39 0 0 0 0 50
## 40 1 0 1 1 41
## 41 0 1 1 0 35
## 42 0 0 1 0 42
## 43 0 0 1 0 27
## 44 0 0 0 0 27
## 45 0 0 1 0 26
## 46 0 0 0 0 32
## 47 0 0 0 0 33
## 48 0 0 0 0 28
## 49 0 0 0 0 30
## 98 1 0 1 1 0# Load survival libraries
library(tidycmprsk)
library(survival)
# Fit survival model
s1 <- survfit(Surv(fmarit) ~ 1, data = atriskW4)
str(s1)## List of 16
## $ n : int 5200
## $ time : num [1:5] 1 2 3 4 5
## $ n.risk : num [1:5] 5200 3735 3716 3368 3239
## $ n.event : num [1:5] 1465 19 348 129 3239
## $ n.censor : num [1:5] 0 0 0 0 0
## $ surv : num [1:5] 0.718 0.715 0.648 0.623 0
## $ std.err : num [1:5] 0.00869 0.00876 0.01023 0.01079 Inf
## $ cumhaz : num [1:5] 0.282 0.287 0.38 0.419 1.419
## $ std.chaz : num [1:5] 0.00736 0.00745 0.00899 0.0096 0.02002
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type: chr "log"
## $ lower : num [1:5] 0.706 0.702 0.635 0.61 NA
## $ upper : num [1:5] 0.731 0.727 0.661 0.636 NA
## $ call : language survfit(formula = Surv(fmarit) ~ 1, data = atriskW4)
## - attr(*, "class")= chr "survfit"#This output provides a summary ofthe results of the survival analysis, including the probabilities over time, stand error, cumulative hazard estimates and cofidence intervals. It helps understand and visualize the survival data.
library(dplyr)
# Filter individuals who had been married (fmarit == 1, 2, 3, or 4)
married_data <- shortdata %>%
filter(fmarit %in% c(1, 2, 3, 4))
# Filter individuals who had not been married (fmarit == 5)
not_married_data <- shortdata %>%
filter(fmarit == 5)# Create a variable for age at first marriage.
dat <- shortdata %>%
mutate(agemarry1 = pmin(hisagem, agemarrn, na.rm=TRUE))
head(dat)## # A tibble: 6 × 6
## age_r hisp fmarit hisagem agemarrn agemarry1
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 31 5 5 NA NA NA
## 2 17 5 5 NA NA NA
## 3 16 1 5 NA NA NA
## 4 49 5 3 NA 18 18
## 5 39 5 5 NA NA NA
## 6 37 5 1 27 NA 27nrow(dat)## [1] 5206dat<-as.data.frame(dat)
tabyl(dat$agemarry1)## dat$agemarry1 n percent valid_percent
## 16 6 0.0011525163 0.003210273
## 17 11 0.0021129466 0.005885500
## 18 42 0.0080676143 0.022471910
## 19 70 0.0134460238 0.037453184
## 20 84 0.0161352286 0.044943820
## 21 136 0.0261237034 0.072766185
## 22 121 0.0232424126 0.064740503
## 23 131 0.0251632731 0.070090958
## 24 139 0.0266999616 0.074371322
## 25 163 0.0313100269 0.087212413
## 26 136 0.0261237034 0.072766185
## 27 117 0.0224740684 0.062600321
## 28 109 0.0209373799 0.058319957
## 29 117 0.0224740684 0.062600321
## 30 88 0.0169035728 0.047084002
## 31 70 0.0134460238 0.037453184
## 32 63 0.0121014214 0.033707865
## 33 57 0.0109489051 0.030497592
## 34 47 0.0090280446 0.025147138
## 35 29 0.0055704956 0.015516319
## 36 31 0.0059546677 0.016586410
## 37 27 0.0051863235 0.014446228
## 38 15 0.0028812908 0.008025682
## 39 7 0.0013446024 0.003745318
## 40 16 0.0030733769 0.008560728
## 41 9 0.0017287745 0.004815409
## 42 8 0.0015366884 0.004280364
## 43 5 0.0009604303 0.002675227
## 44 4 0.0007683442 0.002140182
## 45 4 0.0007683442 0.002140182
## 46 3 0.0005762582 0.001605136
## 98 4 0.0007683442 0.002140182
## NA 3337 0.6409911640 NA#number of each of marriage, the kinds of fmarit variable
table(dat$agemarry1, dat$fmarit)##
## 0 1 2 3 4 5
## 16 0 5 0 1 0 0
## 17 0 6 0 4 1 0
## 18 0 27 3 11 1 0
## 19 0 46 0 17 7 0
## 20 0 59 0 20 5 0
## 21 0 93 2 37 4 0
## 22 0 90 0 27 4 0
## 23 0 102 1 23 5 0
## 24 0 109 2 20 8 0
## 25 0 121 1 30 11 0
## 26 0 117 2 17 0 0
## 27 0 96 0 14 7 0
## 28 0 86 0 17 6 0
## 29 0 97 1 13 6 0
## 30 0 77 0 8 3 0
## 31 0 58 0 9 3 0
## 32 0 56 0 6 1 0
## 33 0 48 1 5 3 0
## 34 0 40 1 4 2 0
## 35 0 23 0 6 0 0
## 36 0 26 0 5 0 0
## 37 0 22 0 4 1 0
## 38 0 10 0 3 2 0
## 39 0 7 0 0 0 0
## 40 0 14 0 1 1 0
## 41 0 7 1 1 0 0
## 42 0 7 0 1 0 0
## 43 0 4 0 1 0 0
## 44 0 4 0 0 0 0
## 45 0 3 0 1 0 0
## 46 0 3 0 0 0 0
## 98 0 2 0 1 1 0length(atriskW4)## [1] 8library(tidycmprsk)
library(survival)
s1 <- survfit(Surv(fmarit) ~ 1, data = atriskW4)
str(s1)## List of 16
## $ n : int 5200
## $ time : num [1:5] 1 2 3 4 5
## $ n.risk : num [1:5] 5200 3735 3716 3368 3239
## $ n.event : num [1:5] 1465 19 348 129 3239
## $ n.censor : num [1:5] 0 0 0 0 0
## $ surv : num [1:5] 0.718 0.715 0.648 0.623 0
## $ std.err : num [1:5] 0.00869 0.00876 0.01023 0.01079 Inf
## $ cumhaz : num [1:5] 0.282 0.287 0.38 0.419 1.419
## $ std.chaz : num [1:5] 0.00736 0.00745 0.00899 0.0096 0.02002
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type: chr "log"
## $ lower : num [1:5] 0.706 0.702 0.635 0.61 NA
## $ upper : num [1:5] 0.731 0.727 0.661 0.636 NA
## $ call : language survfit(formula = Surv(fmarit) ~ 1, data = atriskW4)
## - attr(*, "class")= chr "survfit"#This is a survival analysis of the ‘fmarit’ variable. It is estimating survival probabilities over time. The total number of observations in the dataset are 5,200. The number at risk ‘n. risk’ represent the subjects at risk, which decreases over time as people get married.
#Survival Curve for “Age” a d Cumulative Survival Probability”
library(ggsurvfit)##
## Attaching package: 'ggsurvfit'## The following object is masked from 'package:psych':
##
## %+%survfit2(Surv(ageevent) ~ 1, data = atriskW4) %>%
ggsurvfit() +
labs(
x = "Age",
y = "Cumulative survival probability"
) +
add_confidence_interval() +
add_risktable()
summary(survfit(Surv(ageevent) ~ 1, data = atriskW4), times = 25)## Call: survfit(formula = Surv(ageevent) ~ 1, data = atriskW4)
##
## 1344 observations deleted due to missingness
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 25 1844 2184 0.434 0.00798 0.418 0.45#The results of the survival analysis indicate that 1344 observations were deleted due to missingness since survival survival analysis require complete data for valid results. There are 1,844 individuals who are at risk and 2,184 events that occurred up to time of the survival analysis.
phtest <- survfit(Surv(hisagem, fmarit) ~ hisp, data = atriskW4)## Warning in Surv(hisagem, fmarit): Invalid status value, converted to NA#https://rdrr.io/cran/survival/man/plot.survfit.html
# now running some other plots
plot(phtest, fun = "S", xlab = "age",
ylab = "cumulative survival", main = "cumulative survival curves by hispanic")
#A plot of cumulative survial curves by Hispanic ethnicity was done to
examine the probability of individuals reamining unmarried varried
across different age groups.
plot(phtest, fun = "cumhaz", xlab = "age",
ylab = "cumulative hazard", main = "cumulative hazard curves by hispanic")
#A cumulative hazard curve was done to assess how the cumulative hazard
of individuals experiencing marriage varies across different ages.
plot(phtest, fun = "cloglog", xlab = "age",
ylab = "log-log survival", main = "log-log curves by hispanic")