Assignment 9

library(ipumsr)

setwd("~/Downloads")

ddi <- read_ipums_ddi("nhis_00004.xml")
data <- read_ipums_micro(ddi)

## Use of data from IPUMS NHIS is subject to conditions including that users
## should cite the data appropriately. Use command `ipums_conditions()` for more
## details.

data<- haven::zap_labels(data)

names(data) <- tolower(gsub(pattern = "_",replacement =  "",x =  names(data)))

Data Recoding

data$male <- ifelse(data$sex == 1,1,0)

#Numeric Outcome variable
data$generalhealth<-Recode(data$health,
                               recodes="1=5;2=4;3=3; 4=2; 5=1;  else=NA",
                               as.factor = F)


data$anxiety_disoreder<- car::Recode(data$anxietyev,
                       recodes="1= 0; 2= 1;else=NA")



data$educ<- car::Recode(data$educ,
                     recodes="102:116='1Less than HS'; 201:202='2hsgrad'; 301:503='3More than HS';997:999=NA;000=NA",
                     as.factor=T)
data$educ<-fct_relevel(data$educ,'1Less than HS','2hsgrad','3More than HS') 


# ever had depression
data$hadpression<- car::Recode(data$depressev,
                       recodes="1= 0 ; 2= 1;else=NA")


data$familyincome_num <- data$famtotinc

#had depression or anxiety


data$dep_anx <- paste( data$anxiety_disoreder, data$hadpression, sep ="")

data$dep_anx<- as.numeric(data$dep_anx)

## Warning: NAs introduced by coercion

#race/ethnicity
data$nhwhite<- car::Recode(data$hisprace,
                       recodes="02=1; 99=NA; else=0")

data$other_race_eth<- car::Recode(data$hisprace,
                      recodes="c(1,3,4,5,6,7)=1; 99=NA; else=0")


data$race_eth<-car::Recode(data$hisprace,
recodes="02='NH_White'; c(1,3,4,5,6,7)='other_race_eth';else=NA",
as.factor = T)
data$race_eth<-relevel(data$race_eth,
                          ref = "NH_White")

#Categorial outcome variable

data$dep_anx2<- car::Recode(data$dep_anx,
                       recodes="01:11=1; 00 = 0; else=NA")

data <- data%>%
filter(familyincome_num >=0 & age<=999996)


data <- data%>%
filter(age >=18 & age<=099)

Selecting only needed variables

data2<- data %>% 
  select( strata, sampweight, dep_anx2, race_eth, nhwhite, other_race_eth, educ, generalhealth, age, familyincome_num, male)%>%
  filter(complete.cases(.))

Survey Design

options(survey.lonely.psu = "adjust")

des<-svydesign(ids=~1, strata=~strata, weights=~sampweight, data = data2 )
des

## Stratified Independent Sampling design (with replacement)
## svydesign(ids = ~1, strata = ~strata, weights = ~sampweight, 
##     data = data2)

1) define an outcome, this can be of any form, but use an appropriate distribution for your outcome. Also define at least 1 continuous predictor.

Categorial Outcome: Anxiety/Depression

Continuous Outcome: Family Income

Checking to see the selected continuous variable are linear

Since the above are significant it shows that the age variable is not linear

Plotting to know if the selected continuous variable are linear

plot(datagam)

This gives more evidence that the above is not linear

2) Using the gam() function, estimate a model with only linear terms in your model

datagam_l<-gam(dep_anx2~age+ familyincome_num +generalhealth + educ + male + nhwhite,
              data=data2,
                weights = sampweight/mean(sampweight, na.rm=T),
              family=binomial)

summary(datagam_l)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## dep_anx2 ~ age + familyincome_num + generalhealth + educ + male + 
##     nhwhite
## 
## Parametric coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.700e+00  5.820e-02  29.202  < 2e-16 ***
## age               -1.861e-02  6.114e-04 -30.438  < 2e-16 ***
## familyincome_num  -3.030e-06  2.038e-07 -14.870  < 2e-16 ***
## generalhealth     -6.627e-01  1.094e-02 -60.599  < 2e-16 ***
## educ2hsgrad       -3.501e-02  3.647e-02  -0.960    0.337    
## educ3More than HS  2.054e-01  3.479e-02   5.904 3.55e-09 ***
## male              -7.643e-01  2.163e-02 -35.335  < 2e-16 ***
## nhwhite            8.878e-01  2.418e-02  36.713  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.121   Deviance explained = 11.1%
## UBRE = -0.090304  Scale est. = 1         n = 62931

3) Repeat Step 2, but include a smooth term for at least one continuous variable

Numerical outcome

datagam_s<-gam(dep_anx2~s(age)+ familyincome_num +generalhealth + educ + male + nhwhite,
              data=data2,
                weights = sampweight/mean(sampweight, na.rm=T),
              family=binomial)

summary(datagam_s)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## dep_anx2 ~ s(age) + familyincome_num + generalhealth + educ + 
##     male + nhwhite
## 
## Parametric coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        7.228e-01  4.354e-02  16.601  < 2e-16 ***
## familyincome_num  -3.209e-06  2.079e-07 -15.434  < 2e-16 ***
## generalhealth     -6.574e-01  1.097e-02 -59.901  < 2e-16 ***
## educ2hsgrad       -4.806e-02  3.668e-02  -1.311     0.19    
## educ3More than HS  1.819e-01  3.520e-02   5.167 2.38e-07 ***
## male              -7.672e-01  2.165e-02 -35.437  < 2e-16 ***
## nhwhite            8.979e-01  2.426e-02  37.020  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df Chi.sq p-value    
## s(age) 8.33  8.874  958.5  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.123   Deviance explained = 11.2%
## UBRE = -0.091175  Scale est. = 1         n = 62931

4) Test if the model in step 3 is fitting the data better than the purely linear model.

anova( datagam_l, datagam_s, test="Chisq")

Comparing the two models, the model with a smooth term for the age variable fits the model better. This is evident because of the very small p-value.

5) Produce a plot of the smooth effect from the model in step 3

plot(datagam_s)