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
library(magrittr)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(foreign)
#download datasets - AddHealth
wave1<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0001/da27021p1.dta")
wave3<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0003/da27021p3.dta")
wave4<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0012/da27021p12.dta")
wave4a<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0013/da27021p13.dta")
wave5<-read.xport("/projects/add_health/data/wave_5/wave5.xpt")
wave5age<-read.xport("/projects/add_health/data/wave_5/Constructed Files/Wave V Constructed Age/agew5.xpt")
#wave3section19<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0006/da27021p6.dta")
#download biomarkers Wave 4
w4glu <- read.xport("/projects/add_health/data/wave_4/biomarkers/glu_a1c.xpt")
w4crp <- read.xport("/projects/add_health/data/wave_4/biomarkers/crp_ebv.xpt")
w4lipids <- read.xport("/projects/add_health/data/wave_4/biomarkers/lipids.xpt")
#weightW4 <- read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0028/da27021p28.dta")
#weightW4<-read.xport("/projects/add_health/data/wave_4/Wave IV Grand Sample Weights/weights4/weights4.xpt")
#download biomarkers wave 5
w5brenal<-read.xport("/projects/add_health/data/wave_5/brenal5.xpt")
w5blipids<-read.xport("/projects/add_health/data/wave_5/blipids5.xpt")
w5bcardio<-read.xport("/projects/add_health/data/wave_5/bcardio5.xpt")
w5bglua1c<-read.xport("/projects/add_health/data/wave_5/bglua1c5.xpt")
w5banthro<-read.xport("/projects/add_health/data/wave_5/banthro5.xpt")
w5bcrp<-read.xport("/projects/add_health/data/wave_5/bcrp5.xpt")
w5bdemo<-read.xport("/projects/add_health/data/wave_5/bdemo5.xpt")
w5bweight<-read.xport("/projects/add_health/data/wave_5/bweight5.xpt")
#w5bweight=read.dta("/projects/add_health/data/wave_5/wave5weight.dta")
friendship <- read.dta("/projects/add_health/data/11613344/ICPSR_27033/DS0001/schoolnetworkvariables.dta")
w1school<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0019/da27021p19.dta")
w1school$s11<-w1school$S11
w1school$s17<-w1school$S17
w1school$s25<-w1school$S25
w1school$s26<-w1school$S26
w1school$s44<-w1school$S44
w1school$s46a<-w1school$S46A
w1school$s46d<-w1school$S46D

w1school2<-w1school%>%
  dplyr::select(aid,s11,s17,s25,s26,s44,s46a,s46d)
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
names(wave5)[1] <- "aid"
names(wave5age)[1] <- "aid"
names(w4glu)[1] <- "aid"
names(w4crp)[1] <- "aid"
names(w4lipids)[1] <- "aid"
#names(weightW4)[1] <- "aid"

names(w5brenal)[1] <- "aid"
names(w5blipids)[1] <- "aid"
names(w5bcardio)[1] <- "aid"
names(w5bglua1c)[1] <- "aid"
names(w5banthro)[1] <- "aid"
names(w5bcrp)[1] <- "aid"
names(w5bdemo)[1] <- "aid"
names(w5bweight)[1] <- "aid"
#change uppercase to lower case
#wave5$aid<-wave5$AID
#wave5age$aid<-wave5age$AID
#w4glu$aid<-w4glu$AID
#w4crp$aid<-w4crp$AID
#w4lipids$aid<-w4lipids$AID
#remove variables
#wave5 = subset(wave5, select = -c(AID) )
#w4glu = subset(w4glu, select = -c(AID) )
#w4crp = subset(w4crp, select = -c(AID) )
#w4lipids = subset(w4lipids, select = -c(AID) )
#merge datasets
a<-merge(wave1,wave3,by="aid")
bb<-merge(a,wave4,by="aid")
b<-merge(bb,wave4a,by="aid")
c<-merge(b,wave5,by="aid")
d<-merge(c,w4glu,by="aid")
e<-merge(d,w4crp,by="aid")
f<-merge(e,w4lipids,by="aid")
addhealtha<-merge(f,w5bweight,by="aid")
addhealthb<-merge(addhealtha,wave5age,by="aid")
addhealthc<-merge(addhealthb,w5brenal,by="aid")
addhealthd<-merge(addhealthc,w5blipids,by="aid")
addhealthe<-merge(addhealthd,w5bcardio,by="aid")
addhealthf<-merge(addhealthe,w5bglua1c,by="aid")
addhealthg<-merge(addhealthf,w5banthro,by="aid")
addhealthh<-merge(addhealthg,w5bcrp,by="aid")
addhealthi<-merge(addhealthh,w5bdemo,by="aid")
addhealthj<-merge(addhealthi,friendship,by="aid")
addhealth<-merge(addhealthj,w1school2,by="aid")
addhealth$livesMother<-Recode(addhealth$s11,recodes="0=0;1=1;else=NA")
addhealth$livesFather<-Recode(addhealth$s17,recodes="0=0;1=1;else=NA")

addhealth$livesParent<-ifelse(addhealth$livesMother==1&addhealth$livesFather==1,"lives_2_parents",
                             ifelse(addhealth$livesMother==1&addhealth$livesFather==0,"lives_mom",
                                    ifelse(addhealth$livesMother==0&addhealth$livesFather==1,"lives_dad",
                                           ifelse(addhealth$livesMother==0&addhealth$livesFather==0,"lives_other",NA))))

addhealth$livesParent2<-ifelse(addhealth$livesMother==1&addhealth$livesFather==1,"lives_2_parents",
                             ifelse(addhealth$livesMother==1&addhealth$livesFather==0,"lives_1_parent",
                                    ifelse(addhealth$livesMother==0&addhealth$livesFather==1,"lives_1_parent",
                                           ifelse(addhealth$livesMother==0&addhealth$livesFather==0,"lives_other",NA))))

addhealth$adopt<-ifelse(addhealth$s25==0,"no",
                       ifelse(addhealth$s25==1&addhealth$s26==0,"adoptNONbioparent",
                              ifelse(addhealth$s25==1&addhealth$s26==1,"adoptbioparent",NA)))
addhealth$schoolClubJoin<-addhealth$s44
addhealth$teacherTrouble<-Recode(addhealth$s46a,recodes="0:1=0;2:4=1;else=NA")
addhealth$studentTrouble<-Recode(addhealth$s46d,recodes="0:1=0;2:4=1;else=NA")
table(addhealth$schoolClubJoin)
## 
##    0    1 
## 4980  907
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)
#age wave 4
#birth month wave 4
#addhealth$birthmonthw4<-Recode(addhealth$H4OD1M,recodes="1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12")
#birth year wave 4
#addhealth$birthyearw4<-addhealth$H4OD1Y
#addhealth$birthyearw4 <- Recode(addhealth$H4OD1Y,recodes="1974=1974;1975=1975;1976=1976;1977=1977;1978=1978;1979=1979;1980=1980;1981=1981;1982=1982;1983=1983")
#combine year and month for birth date wave 4
#addhealth$birthdatew4 <- as.yearmon(paste(addhealth$birthyearw4, addhealth$birthmonthw4), "%Y %m")
#interview date w4
#addhealth$interviewmonthw4 <- Recode(addhealth$IMONTH4,recodes="1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12")
#addhealth$interviewyearw4 <- Recode(addhealth$IYEAR4,recodes="2007=2007;2008=2008;2009=2009")
#addhealth$interviewdatew4<-as.yearmon(paste(addhealth$interviewyearw4,addhealth$interviewmonthw4),"%Y %m")
#age - derived from date of birth subtracted from wave interview date 
#addhealth$ageW4<-(as.numeric(round(addhealth$interviewdatew4 - addhealth$birthdatew4)))
#table(addhealth$ageW4)
#Date of birth variables; Wave 1: H1GI1M (month), H1GI1Y (year); Wave 3: H3OD1M (month), H3OD1Y (year); Wave 4: H4OD1M (month), H4OD1Y (year)
addhealth$birthmonthw4 <- Recode(addhealth$H4OD1M,recodes="1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12")
addhealth$birthyearw1 <- factor(ifelse(addhealth$H1GI1Y=="74",1974,
                                             ifelse(addhealth$H1GI1Y=="75",1975,
                                                           ifelse(addhealth$H1GI1Y=="76",1976,
                                                                         ifelse(addhealth$H1GI1Y=="77",1977,
                                                                                       ifelse(addhealth$H1GI1Y=="78",1978,
                                                                                                     ifelse(addhealth$H1GI1Y=="79",1979,
                                                                                                                   ifelse(addhealth$H1GI1Y=="80",1980,
                                                                                                                                ifelse(addhealth$H1GI1Y=="81",1981,
                                                                                                                                               ifelse(addhealth$H1GI1Y=="82",1982,
                                                                                                                                                             ifelse(addhealth$H1GI1Y=="83",1983,1979)))))))))))
addhealth$birthmonthw1 <- Recode(addhealth$H1GI1M,recodes="1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12;else=NA")
addhealth$birthyearw4<-addhealth$H4OD1Y
addhealth$birthyearw4 <- Recode(addhealth$H4OD1Y,recodes="1974=1974;1975=1975;1976=1976;1977=1977;1978=1978;1979=1979;1980=1980;1981=1981;1982=1982;1983=1983")
addhealth$interviewmonthw4 <- Recode(addhealth$IMONTH4,recodes="1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12")
table(addhealth$interviewmonthw4)
## 
##    1    2    3    4    5    6    7    8    9   10   11   12 
##  277 1056 1141 1147  870  450  332  241  147   98   81   47
addhealth$interviewyearw4 <- Recode(addhealth$IYEAR4,recodes="2007=2007;2008=2008;2009=2009")
addhealth$iyearfix <- Recode(addhealth$iyear,recodes="'94'='1994';'95'='1995'")
addhealth$birthdatew4 <- as.yearmon(paste(addhealth$birthyearw4, addhealth$birthmonthw4), "%Y %m")
addhealth$interviewdatew4<-as.yearmon(paste(addhealth$interviewyearw4,addhealth$interviewmonthw4),"%Y %m")
addhealth$birthdatew3<-as.yearmon(paste(addhealth$H3OD1Y,addhealth$H3OD1M),"%Y %m")
addhealth$interviewdatew3<-as.yearmon(paste(addhealth$IYEAR3,addhealth$IMONTH3 ),"%Y %m")
addhealth$birthdatew1 <- as.yearmon(paste(addhealth$birthyearw1, addhealth$birthmonthw1), "%Y %m")
addhealth$interviewdatew1 <- as.yearmon(paste(addhealth$iyearfix, addhealth$imonth), "%Y %m")

addhealth$ageW4<-(as.numeric(round(addhealth$interviewdatew4 - addhealth$birthdatew4)))
addhealth$ageW3<-(as.numeric(round(addhealth$interviewdatew3 - addhealth$birthdatew3)))
addhealth$ageW1<-(as.numeric(round(addhealth$interviewdatew1 - addhealth$birthdatew1)))
addhealth$ageW5<-addhealth$AGE
#sex assigned at birth
addhealth$sex<-Recode(addhealth$H5OD2A,recodes="1='Male';2='Female';else=NA")
addhealth$sexW5<-factor(addhealth$sex,levels=c("Male","Female"))

#Sexual Orientation
# Wave 3
addhealth$sxorW3<-Recode(addhealth$H3SE13, recodes="1:2='Straight'; 3='Bisexual';4:5='Gay';else=NA")
addhealth$sexualOrientationW3 <- factor(addhealth$sxorW3, levels = c("Straight", "Gay", "Bisexual"))
# Wave 4
addhealth$sxorW4<-Recode(addhealth$H4SE31, recodes="1:2='Straight'; 3='Bisexual';4:5='Gay';else=NA")
addhealth$sexualOrientationW4 <- factor(addhealth$sxorW4, levels = c("Straight", "Gay", "Bisexual"))
# Wave 5
addhealth$sxorW5<-Recode(addhealth$H5SE15,recodes="1:2='Straight';3='Bisexual';4:5='Gay';else=NA")
addhealth$sexualOrientationW5 <- factor(addhealth$sxorW5, levels = c("Straight", "Gay", "Bisexual"))
#education wave 5
addhealth$educW5<-Recode(addhealth$H5OD11,recodes="1:2=0;3:9=1;10:16=2;else=NA")

addhealth$educW5<-as.factor(addhealth$educW5)
#employment status wave 1
#addhealth$H1EE3[addhealth$H1EE3 == 6] <- NA
#addhealth$H1EE3[addhealth$H1EE3 == 8] <- NA
#addhealth$employmentW1<-addhealth$H1EE3
#employment status wave 3
#addhealth$H3LM1[addhealth$H3LM1 == 6] <- NA
#addhealth$H3LM1[addhealth$H3LM1 == 8] <- NA
#addhealth$H3LM2[addhealth$H3LM2 == 6] <- NA
#addhealth$H3LM2[addhealth$H3LM2 == 7] <- NA
#addhealth$H3LM2[addhealth$H3LM2 == 8] <- NA
#addhealth$H3LM2[addhealth$H3LM2 == 9] <- NA
#addhealth$H3LM7[addhealth$H3LM7 == 7] <- NA
#addhealth$employmentW3<-ifelse(addhealth$H3LM1<2&addhealth$H3LM7>0,1,
#                               ifelse(addhealth$H3LM1<2&addhealth$H3LM7<1,0,
#                                      ifelse(addhealth$H3LM2<2&addhealth$H3LM7>0,1,
#                                             ifelse(addhealth$H3LM2<2&addhealth$H3LM7<1,0,NA))))
#employment status wave 4
#addhealth$H4LM1[addhealth$H4LM1 == 6] <- NA
#addhealth$H4LM2[addhealth$H4LM2 == 6] <- NA
#addhealth$H4LM2[addhealth$H4LM2 == 7] <- NA

#addhealth$employmentW4<-ifelse(addhealth$H4LM1<2&addhealth$H4LM14<11,0,
#                                ifelse(addhealth$H4LM2<2&addhealth$H4LM14<11,0,1))

#employment status Wave 5
#addhealth$employmentW5<-ifelse(addhealth$H5LM5==1,1,0)

#addhealth$employmentW1<-as.factor(addhealth$employmentW1)
#addhealth$employmentW3<-as.factor(addhealth$employmentW3)
#addhealth$employmentW4<-as.factor(addhealth$employmentW4)
#addhealth$employmentW5<-as.factor(addhealth$employmentW5)
#marital status Wave 1
addhealth$H1GI15[addhealth$H1GI15 == 6] <- NA
addhealth$H1GI15[addhealth$H1GI15 == 7] <- NA
addhealth$H1GI15[addhealth$H1GI15 == 8] <- NA

addhealth$maritalW1<-ifelse(addhealth$ageW1<15,0,
                             ifelse(addhealth$ageW1>14&addhealth$H1GI15<1,0,
                                    ifelse(addhealth$ageW1>14&addhealth$H1GI15>0,1,NA)))

#marital status Wave 3 NOTE 0=never married, 1=married, 2=widow, 3=divorce
addhealth$H3MR1[addhealth$H3MR1 == 6] <- NA
addhealth$H3MR1[addhealth$H3MR1 == 8] <- NA
addhealth$H3MR1[addhealth$H3MR1 == 9] <- NA
addhealth$H3MR3_A[addhealth$H3MR3_A == 7] <- NA
addhealth$H3MR3_A[addhealth$H3MR3_A == 8] <- NA
addhealth$maritalW3<-ifelse(addhealth$H3MR1==0,0,
                            ifelse(addhealth$H3MR1>0&addhealth$H3MR3_A==1,1,2))
#marital status wave 4
addhealth$H4TR1[addhealth$H4TR1 == 6] <- NA
addhealth$H4TR1[addhealth$H4TR1 == 8] <- NA
addhealth$maritalW4<-ifelse(addhealth$H4TR1==0,0,
                            ifelse(addhealth$H4TR1>0&addhealth$H4TR13==1,1,
                                   ifelse(addhealth$H4TR1>0&addhealth$H4TR13!=1,2,NA)))
#marital status wave 5 
addhealth$maritalW5<-Recode(addhealth$H5HR1,recodes="1=1;2:4=2;5=0;else=NA")


#addhealth$maritalW1<-as.factor(addhealth$maritalW1)
#addhealth$maritalW3<-as.factor(addhealth$maritalW3)
#addhealth$maritalW4<-as.factor(addhealth$maritalW4)
#addhealth$maritalW5<-as.factor(addhealth$maritalW5)
#table(addhealth$maritalW5)
#combine race to one variable.  This is for white, black, hispanic, asian, native american, pacific islander, and other.
addhealth$raceW5all <- factor(ifelse(addhealth$H5OD4A==1,"a_white", 
                                            ifelse(addhealth$H5OD4B==1, "b_black", 
                                            ifelse(addhealth$H5OD4C==1, "c_hispanic",
                                            ifelse(addhealth$H5OD4D==1,"d_asian",
                                            ifelse(addhealth$H5OD4F==1,"e_nativeamer",
                                            ifelse(addhealth$H5OD4E==1,"f_pacificisl",
                                            ifelse(addhealth$H5OD4G==1,"g_other",NA))))))))
#combine race to one variable.  This is for white, black, hispanic, asian, native american, pacific islander, and other.
addhealth$raceW5min <- factor(ifelse(addhealth$H5OD4A==1,"a_white", 
                                            ifelse(addhealth$H5OD4B==1, "b_black", 
                                            ifelse(addhealth$H5OD4C==1, "c_hispanic",NA))))
library(StatMeasures)
#biomarkers wave 5
#1 systolic bp
addhealth$H5SBP<-Recode(addhealth$H5SBP,recodes="9999=NA;9997=NA;9996=NA")
addhealth$systolicbpW5<-ifelse(addhealth$H5SBP>=140,1,0)
#2 diastolic bp
addhealth$H5DBP<-Recode(addhealth$H5DBP,recodes="9999=NA;9997=NA;9996=NA")
addhealth$diastolicbpW5<-ifelse(addhealth$H5DBP>=90,1,0)
#3 pulse rate
addhealth$H5PR<-Recode(addhealth$H5PR,recodes="9999=NA;9997=NA;9996=NA")
addhealth$pulserateW5<-ifelse(addhealth$H5PR>=85,1,0)
#4 pulse pressure
addhealth$H5PP<-Recode(addhealth$H5PP,recodes="9999=NA;9997=NA;9996=NA")
addhealth$pulsepressureW5<-ifelse(addhealth$H5PP>=60,1,0)
#5 TOTAL CHOLESTEROL W5
addhealth$tcW5<-decile(addhealth$H5TC)
addhealth$totalcholW5<-ifelse(addhealth$tcW5>=8,1,0)
table(addhealth$totalcholW5)
## 
##    0    1 
## 3741 1519
#6 HDL, HIGH-DENSITY LIPOPROTEIN CHOLESTEROL
addhealth$H5HDL<-Recode(addhealth$H5HDL,recodes="9999=NA")
addhealth$hdlW5<-decile(addhealth$H5HDL)
addhealth$hidensity_lipW5<-ifelse(addhealth$hdlW5<=2,1,0)
#7 TRIGLYCERIDES
addhealth$TG<-decile(addhealth$H5TG)
addhealth$triglycerideW5<-ifelse(addhealth$TG>=8,1,0)
#8 HEMOGLOBIN
addhealth$hemoglobinW5<-ifelse(addhealth$H5HBA1C>=6.4,1,0)
table(addhealth$hemoglobinW5)
## 
##    0    1 
## 4901  246
#9 glucose
addhealth$GLUCOSE<-Recode(addhealth$H5GLUCOS,recodes="9999=NA")
addhealth$glucoseW5<-ifelse(addhealth$GLUCOSE>=200,1,0)
table(addhealth$glucoseW5)
## 
##    0    1 
## 5123  132
#10 BMI
addhealth$H5BMI<-Recode(addhealth$H5BMI,recodes="9999=NA;9996=NA;9997=NA;9994=NA")
addhealth$bmiW5<-ifelse(addhealth$H5BMI>=30,1,0)
table(addhealth$bmiW5)
## 
##    0    1 
## 2982 2780
#11 waist circumference 
addhealth$waistW5<-Recode(addhealth$H5WSTCLS,recodes="1=0;2=1")
#12 high-sensitivity C-reactive protein test. 
addhealth$hsCRPW5<-ifelse(addhealth$H5CRP>=3,1,0)
table(addhealth$hsCRPW5)
## 
##    0    1 
## 3124 1938
# COMBINE BIOMARKERS
addhealth$scoreW5<-addhealth$systolicbpW5+addhealth$diastolicbpW5+addhealth$pulserateW5+addhealth$pulsepressureW5+addhealth$totalcholW5+addhealth$hidensity_lipW5+addhealth$triglycerideW5+addhealth$hemoglobinW5+addhealth$glucoseW5+addhealth$bmiW5+addhealth$waistW5+addhealth$hsCRPW5
table(addhealth$scoreW5)
## 
##   0   1   2   3   4   5   6   7   8   9  10 
## 784 687 761 788 608 418 337 169  76  32  20
#biomarkers Wave 4
#1 systolic bp
addhealth$H4SBP<-Recode(addhealth$H4SBP,recodes="999=NA;997=NA;996=NA")
addhealth$systolicbp<-ifelse(addhealth$H4SBP>=140,1,0)
#2 diastolic bp
addhealth$H4DBP<-Recode(addhealth$H4DBP,recodes="999=NA;997=NA;996=NA")
addhealth$diastolicbp<-ifelse(addhealth$H4DBP>=90,1,0)
#3 pulse rate
addhealth$H4PR<-Recode(addhealth$H4PR,recodes="999=NA;997=NA;996=NA")
addhealth$pulserate<-ifelse(addhealth$H4PR>=85,1,0)
#4 pulse pressure
addhealth$H4PP<-Recode(addhealth$H4PP,recodes="999=NA;997=NA;996=NA")
addhealth$pulsepressure<-ifelse(addhealth$H4PP>=60,1,0)
#5 TOTAL CHOLESTEROL
addhealth$TC<-Recode(addhealth$TC,recodes="99=NA")
addhealth$totalchol<-ifelse(addhealth$TC>=8,1,0)
#6 HDL, HIGH-DENSITY LIPOPROTEIN CHOLESTEROL
addhealth$HDL<-Recode(addhealth$HDL,recodes="99=NA")
addhealth$hidensity_lip<-ifelse(addhealth$HDL<=2,1,0)

#7 TRIGLYCERIDES
addhealth$TG<-Recode(addhealth$TG,recodes="99=NA")
addhealth$triglyceride<-ifelse(addhealth$TG>=8,1,0)

#8 HEMOGLOBIN
addhealth$HBA1C<-Recode(addhealth$HBA1C,recodes="99=NA")
addhealth$hemoglobin<-ifelse(addhealth$HBA1C>=6.4,1,0)

#9 glucose
addhealth$GLUCOSE<-Recode(addhealth$GLUCOSE,recodes="999=NA")
addhealth$glucose<-ifelse(addhealth$GLUCOSE>=200,1,0)

#10 BMI
addhealth$H4BMI<-Recode(addhealth$H4BMI,recodes="999=NA;888=NA;889=NA;996=NA;997=NA")
addhealth$BMI<-ifelse(addhealth$H4BMI>=30,1,0)

# 11 high-sensitivity C-reactive protein test. 
addhealth$CRPa<-Recode(addhealth$CRP,recodes="999=NA;998=NA")
addhealth$hsCRP<-ifelse(addhealth$CRPa>=3,1,0)

#EPSTEIN BARR
#quantile(addhealth$EBVa, c(.80), na.rm = T)
#addhealth$EBVa<-Recode(addhealth$EBV,recodes="9999=NA;9998=NA")
#addhealth$EpsteinBar<-ifelse(addhealth$EBVa>=243,1,0)

#12 waist circumference 
addhealth$H4WAIST[addhealth$H4WAIST == 996] <- NA
addhealth$H4WAIST[addhealth$H4WAIST == 997] <- NA
addhealth$H4WAIST[addhealth$H4WAIST == 998] <- NA
addhealth$waist<-ifelse(addhealth$H5OD2A==2&addhealth$H4WAIST>88,1,
                        ifelse(addhealth$H5OD2A==2&addhealth$H4WAIST<89,0,
                        ifelse(addhealth$H5OD2A==1&addhealth$H4WAIST>102,1,
                               ifelse(addhealth$H5OD2A==1&addhealth$H4WAIST<103,0,NA))))

#combine biomarkers
addhealth$scoreW4<-addhealth$hsCRP+addhealth$BMI+addhealth$triglyceride+addhealth$hemoglobin+addhealth$hidensity_lip+addhealth$totalchol+addhealth$pulserate+addhealth$systolicbp+addhealth$pulsepressure+addhealth$diastolicbp+addhealth$glucose+addhealth$waist
summary(addhealth$HHincomeW3)
## Length  Class   Mode 
##      0   NULL   NULL
#household income wave 3
#addhealth$HHincomeW3<-Recode(addhealth$H3EC6)
addhealth$HHincomeW3<-ifelse(addhealth$H3EC6>0&addhealth$H3EC6<30000,'a_<30K',
                             ifelse(addhealth$H3EC6>29999&addhealth$H3EC6<50000,'b_30k<50k',
                                    ifelse(addhealth$H3EC6>49999&addhealth$H3EC6<100000,'c_50k>100k',
                                           ifelse(addhealth$H3EC6>99999,'d_>100k',NA))))
addhealth$HHincomeW3<-as.factor(addhealth$HHincomeW3)
#household income wave 5
addhealth$HHincomeW5<- Recode(addhealth$H5EC2,recodes="1:6='a_<30K';7:8='b_30k<50k';9:10='c_50k>100k';11:13='d_>100k';else=NA")

#household income wave 4
addhealth$HHincomeW4<- Recode(addhealth$H4EC1,recodes="1:6='a_<30K';7:8='b_30k<50k';9:10='c_50k>100k';11:12='d_>100k';else=NA")
#percieved closeness mom WAVE 1
addhealth$close.momW1<-Recode(addhealth$H1WP9,recodes="1=1;2=2;3=3;4=4;5=5;else=NA")
#percieved closeness mom WAVE 3
addhealth$close.momW3<-Recode(addhealth$H3WP20,recodes="5=1;4=2;3=3;2=4;1=5;else=NA")
#percieved closeness mom WAVE 4
addhealth$close.momW4<-Recode(addhealth$H4WP24,recodes="1=1;2=2;3=3;4=4;5=5;else=NA")
#percieved closeness mom WAVE 5
addhealth$close.momW5<-Recode(addhealth$H5WP14,recodes="1=1;2=2;3=3;4=4;5=5;else=NA")

#parent SUPPORT/SUPERVISION involvement 
addhealth$curfew<-Recode(addhealth$H1WP1,recodes="0=0;1=1;else=NA")
addhealth$hang<-Recode(addhealth$H1WP2,recodes="0=0;1=1;else=NA")
addhealth$wear<-Recode(addhealth$H1WP3,recodes="0=0;1=1;else=NA")
addhealth$tv<-Recode(addhealth$H1WP4,recodes="0=0;1=1;else=NA")
addhealth$prog<-Recode(addhealth$H1WP5,recodes="0=0;1=1;else=NA")
addhealth$bed<-Recode(addhealth$H1WP6,recodes="0=0;1=1;else=NA")
addhealth$eat<-Recode(addhealth$H1WP7,recodes="0=0;1=1;else=NA")

addhealth$parentSupervise<-addhealth$curfew+addhealth$hang+addhealth$wear+addhealth$tv+addhealth$prog+addhealth$bed+addhealth$eat

#mothers education
addhealth$mom.educ<-Recode(addhealth$H1RM1,recodes="1:3=0;4:7=1;8:9=2;else=NA")
addhealth$mom.educ2<-Recode(addhealth$H1RM1,recodes="1:3=0;4:9=1;else=NA")
#mothers nativity
addhealth$mom.nativity<-Recode(addhealth$H1RM2,recodes="0=0;1=1;else=NA")
#mother employed
addhealth$mom.employed<-Recode(addhealth$H1RM4,recodes="16=0;1:15=1;else=NA")
#mother accecpt public assistance
addhealth$mom.pubass<-Recode(addhealth$H1RM9,recodes="0=0;1=1;else=NA")
addhealth$motherInequity<-addhealth$mom.educ2+addhealth$mom.nativity+addhealth$mom.employed+addhealth$mom.pubass
addhealth$relw1<-Recode(addhealth$H1RE1,recodes="0=0;1:28=1;else=NA")

addhealth$religionW1<-ifelse(addhealth$relw1==0,0,
                              ifelse(addhealth$relw1<2&addhealth$H1RE3==4,0,
                                  ifelse(addhealth$relw1<2&addhealth$H1RE3==3,0,
                                      ifelse(addhealth$relw1<2&addhealth$H1RE3==2,0,
                                          ifelse(addhealth$relw1<2&addhealth$H1RE3==1,1,NA)))))
                                      
                                  

addhealth$religionW3<-Recode(addhealth$H3RE24,recodes="0:4=0;5:6=1;else=NA")                                        

addhealth$religionW4<-Recode(addhealth$H4RE7,recodes="0:3=0;4:5=1;else=NA")                                                
                                                            
addhealth$religionW5<-Recode(addhealth$H5RE2,recodes="1:4=0;5:6=1")

#school prestige friends
addhealth$pre<-addhealth$prxprest
addhealth$schoolSeg<-addhealth$SEG1RCE5
#education level of 0=less high school, 1=highschool, 
#education wave 3
addhealth$educW3<-ifelse(addhealth$H3ED2==1,1,
                         ifelse(addhealth$H3ED3==1,1,
                                ifelse(addhealth$H3ED5==1,2,
                                       ifelse(addhealth$H3ED6==1&addhealth$H3ED5!=1,2,
                                              ifelse(addhealth$H3ED7==1&addhealth$H3ED5!=1,2,
                                                     ifelse(addhealth$H3ED8==1&addhealth$H3ED5!=1,2,
                                                            ifelse(addhealth$H3ED9==1,0,NA)))))))
#education wave 4
addhealth$educW4<-Recode(addhealth$H4ED2,recodes="1:2=0;3:6=1;7:13=2;else=NA")
#education wave 5
addhealth$educW5<-Recode(addhealth$H5OD11,recodes="1:2=0;3:9=1;10:16=2;else=NA")
#education wave 1
addhealth$educW1<-addhealth$H1GI18

addhealth$educW1<-as.factor(addhealth$educW1)
addhealth$educW3<-as.factor(addhealth$educW3)
addhealth$educW4<-as.factor(addhealth$educW4)
addhealth$educW5<-as.factor(addhealth$educW5)
#employment
#employment binary scaled
addhealth$currently.employed<-Recode(addhealth$H5LM5,recodes="1=1;2:3=0;else=NA")
#employer provides insurance, benefits, sick leave
addhealth$H5LM13A<-Recode(addhealth$H5LM13A,recodes="0=0;1=1;else=NA")
addhealth$H5LM13B<-Recode(addhealth$H5LM13B,recodes="0=0;1=1;else=NA")
addhealth$H5LM13C<-Recode(addhealth$H5LM13C,recodes="0=0;1=1;else=NA")
#scale benefits
addhealth$benefitscore<-addhealth$H5LM13A+addhealth$H5LM13B+addhealth$H5LM13C
addhealth$beneScore<-Recode(addhealth$benefitscore,recodes="0=0;1:2=1;3=2")
table(addhealth$benefitscore)
## 
##    0    1    2    3 
##  604  380  351 3406
table(addhealth$beneScore)
## 
##    0    1    2 
##  604  731 3406
addhealth$employment<-ifelse(addhealth$currently.employed==0,"unemployed",
                             ifelse(addhealth$currently.employed==1&addhealth$beneScore==0,"empNoBene",
                                    ifelse(addhealth$currently.employed==1&addhealth$beneScore==1,"empPartBene",
                                           ifelse(addhealth$currently.employed==1&addhealth$beneScore==2,"empFullBene",NA))))
table(addhealth$currently.employed)
## 
##    0    1 
## 1084 4783
table(addhealth$benefitscore)
## 
##    0    1    2    3 
##  604  380  351 3406
table(addhealth$beneScore)
## 
##    0    1    2 
##  604  731 3406
addhealth$employment<-as.factor(addhealth$employment)
summary(addhealth$employment)
## empFullBene   empNoBene empPartBene  unemployed        NA's 
##        3406         604         731        1084          62
#install.packages("magrittr")
#library(data.table)
library(magrittr)
addhealthSub<-addhealth%>%
  dplyr::select(aid,W5BIOWGT,raceW5min,sexW5,scoreW4,scoreW5,ageW1,ageW3,ageW4,ageW5,sexualOrientationW3,sexualOrientationW4,sexualOrientationW5,educW1,educW3,educW4,educW5,employment,maritalW1,maritalW3,maritalW4,maritalW5,religionW1,religionW3,religionW4,religionW5,HHincomeW3,HHincomeW4,HHincomeW5,livesParent2,close.momW1,close.momW3,close.momW4,close.momW5,parentSupervise,mom.educ,mom.educ2,mom.nativity,mom.employed,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,scoreW4,scoreW5)
out<-melt(setDT(addhealthSub), id = "aid",
          measure.vars = list(age=c("ageW1","ageW3","ageW4","ageW5"),
                             sexualorientation=c("sexualOrientationW3","sexualOrientationW3", "sexualOrientationW4", "sexualOrientationW5"), 
                              education=c("educW1","educW3", "educW4", "educW5"),
                              hhincome=c("HHincomeW3","HHincomeW4", "HHincomeW4", "HHincomeW5"),
                              marriage=c("maritalW1","maritalW3","maritalW4","maritalW5"),
                              religion=c("religionW1","religionW3","religionW4","religionW5"),
                              closeMom=c("close.momW1","close.momW3","close.momW4","close.momW5"),
                             score=c("scoreW4","scoreW4","scoreW4","scoreW5")))%>%
  setorder(aid)

#library(dplyr)

addhealthLong<-addhealthSub%>%
  dplyr::select(aid, W5BIOWGT,raceW5min,sexW5,livesParent2,parentSupervise,mom.educ,mom.educ2,mom.nativity,mom.employed,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,employment)%>%
  left_join(., out, "aid")


addhealthLong$wave<-Recode(addhealthLong$variable,recodes="1=1;2=3;3=4;4=5")

addhealthS1long <- addhealthLong[,c(27,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,19,20,21,22,23,24,25,26)]
addhealthS1long$sexualorientation <- factor(addhealthS1long$sexualorientation, levels = c("Straight", "Gay", "Bisexual"))
addhealthS1long$raceXsexo<-interaction(addhealthS1long$raceW5min,addhealthS1long$sexualorientation,drop=TRUE)

addhealthS1long$education<-as.factor(addhealthS1long$education)
addhealthS1long$hhincome<-as.factor(addhealthS1long$hhincome)
#addhealthS1long$employment<-as.factor(addhealthS1long$employment)
addhealthS1long$marriage<-as.factor(addhealthS1long$marriage)
addhealthS1long$mom.educ<-as.factor(addhealthS1long$mom.educ)
addhealthS1long$mom.educ2<-as.factor(addhealthS1long$mom.educ2)
addhealthS1long$mom.employed<-as.factor(addhealthS1long$mom.employed)
addhealthS1long$mom.pubass<-as.factor(addhealthS1long$mom.pubass)
addhealthS1long$mom.nativity<-as.factor(addhealthS1long$mom.nativity)
addhealthS1long$religion<-as.factor(addhealthS1long$religion)

addhealthS1long$livesParent2<-as.factor(addhealthS1long$livesParent2)
addhealthS1long$schoolClubJoin<-as.factor(addhealthS1long$schoolClubJoin)
addhealthS1long$teacherTrouble<-as.factor(addhealthS1long$teacherTrouble)
addhealthS1long$studentTrouble<-as.factor(addhealthS1long$studentTrouble)
addhealthS1long$parentSupervise<-as.factor(addhealthS1long$parentSupervise)
addhealthS1long$closeMom<-as.factor(addhealthS1long$closeMom)
addhealthexp<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo)

#library(tidyr)

addhealthBOMB<-addhealthexp%>%
  filter(complete.cases(.))
addhealthBOMB$livesParent2<-factor(addhealthBOMB$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthBOMB$closeMom<-as.numeric(addhealthBOMB$closeMom)
addhealthBOMB$closeMom2<-Recode(addhealthBOMB$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthBOMB$parentSupervise<-as.numeric(addhealthBOMB$parentSupervise)
#table(addhealthBOMB$parentSupervise)
addhealthBOMB$parentSupervise2<-Recode(addhealthBOMB$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthBOMB$teacherTrouble2<-as.numeric(addhealthBOMB$teacherTrouble)
addhealthBOMB$studentTrouble2<-as.numeric(addhealthBOMB$studentTrouble)
addhealthBOMB$schoolconflict<-addhealthBOMB$studentTrouble2+addhealthBOMB$teacherTrouble2
addhealthBOMB$schoolconflict<-as.factor(addhealthBOMB$schoolconflict)
#install.packages("geepack")
library(geepack)
#wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,,,,,,,,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,

fit<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBOMB, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fit)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage, data = addhealthBOMB, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                                Estimate    Std.err   Wald Pr(>|W|)    
## (Intercept)                   3.4649189  0.6168953 31.547 1.95e-08 ***
## age                          -0.0875666  0.0379377  5.328  0.02099 *  
## I(age^2)                      0.0020074  0.0007525  7.116  0.00764 ** 
## sexW5Female                   0.1581419  0.1211654  1.703  0.19183    
## raceXsexob_black.Straight     0.3806230  0.1539909  6.109  0.01345 *  
## raceXsexoc_hispanic.Straight  0.2365964  0.2456997  0.927  0.33557    
## raceXsexoa_white.Gay         -0.3981976  0.1595432  6.229  0.01257 *  
## raceXsexob_black.Gay          0.1747051  0.2207857  0.626  0.42878    
## raceXsexoc_hispanic.Gay       0.8971712  0.6982601  1.651  0.19884    
## raceXsexoa_white.Bisexual     0.3043039  0.1867727  2.655  0.10326    
## raceXsexob_black.Bisexual    -0.1947599  0.3981531  0.239  0.62473    
## raceXsexoc_hispanic.Bisexual  0.2228296  0.3274918  0.463  0.49624    
## education1                   -0.1437687  0.2037183  0.498  0.48036    
## education2                   -0.4305850  0.2211951  3.789  0.05158 .  
## employmentempNoBene           0.0930432  0.1954456  0.227  0.63403    
## employmentempPartBene         0.1428496  0.1765695  0.655  0.41850    
## employmentunemployed          0.4651464  0.1819475  6.536  0.01057 *  
## hhincomeb_30k<50k            -0.3015014  0.1581869  3.633  0.05665 .  
## hhincomec_50k>100k           -0.0190452  0.1552813  0.015  0.90238    
## hhincomed_>100k              -0.0955309  0.1868809  0.261  0.60922    
## marriage1                     0.0170306  0.1052989  0.026  0.87151    
## marriage2                     0.2332339  0.2060918  1.281  0.25776    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)    3.602  0.2629
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha   0.7911  0.1541
## Number of clusters:   2063  Maximum cluster size: 450
#education+employment+hhincome+marriage
#fit2<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBOMB, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fit2)

fit3<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBOMB, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fit3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage + livesParent2 + 
##     mom.educ + mom.pubass + closeMom2 + parentSupervise2 + religion + 
##     schoolconflict + schoolClubJoin + scale(schoolSeg) + scale(pre), 
##     data = addhealthBOMB, weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), 
##     id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                               Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                   3.761026  0.724491 26.95  2.1e-07 ***
## age                          -0.094735  0.041195  5.29   0.0215 *  
## I(age^2)                      0.002105  0.000795  7.01   0.0081 ** 
## sexW5Female                   0.120919  0.118919  1.03   0.3092    
## raceXsexob_black.Straight     0.218361  0.162938  1.80   0.1802    
## raceXsexoc_hispanic.Straight -0.034799  0.268646  0.02   0.8969    
## raceXsexoa_white.Gay         -0.397398  0.162860  5.95   0.0147 *  
## raceXsexob_black.Gay         -0.038733  0.228033  0.03   0.8651    
## raceXsexoc_hispanic.Gay       0.499350  0.629164  0.63   0.4274    
## raceXsexoa_white.Bisexual     0.315399  0.178738  3.11   0.0776 .  
## raceXsexob_black.Bisexual    -0.425689  0.408923  1.08   0.2979    
## raceXsexoc_hispanic.Bisexual -0.133019  0.329553  0.16   0.6865    
## education1                   -0.166293  0.202525  0.67   0.4116    
## education2                   -0.450156  0.224251  4.03   0.0447 *  
## employmentempNoBene           0.084546  0.192836  0.19   0.6611    
## employmentempPartBene         0.070546  0.171022  0.17   0.6800    
## employmentunemployed          0.415848  0.177423  5.49   0.0191 *  
## hhincomeb_30k<50k            -0.279200  0.149203  3.50   0.0613 .  
## hhincomec_50k>100k           -0.020737  0.152293  0.02   0.8917    
## hhincomed_>100k              -0.091313  0.181855  0.25   0.6156    
## marriage1                     0.027869  0.108661  0.07   0.7976    
## marriage2                     0.248431  0.210721  1.39   0.2384    
## livesParent2lives_1_parent    0.238479  0.160819  2.20   0.1381    
## livesParent2lives_other      -0.109856  0.332476  0.11   0.7411    
## mom.educ1                    -0.360076  0.168951  4.54   0.0331 *  
## mom.educ2                    -0.488375  0.187645  6.77   0.0093 ** 
## mom.pubass1                   0.357081  0.232316  2.36   0.1243    
## closeMom2b_close              0.005083  0.094281  0.00   0.9570    
## parentSupervise2b_moderate    0.095166  0.238347  0.16   0.6897    
## parentSupervise2c_strict      0.348284  0.235723  2.18   0.1395    
## religion1                    -0.146896  0.122950  1.43   0.2322    
## schoolconflict3               0.050046  0.132491  0.14   0.7056    
## schoolconflict4              -0.033908  0.153692  0.05   0.8254    
## schoolClubJoin1               0.293135  0.185467  2.50   0.1140    
## scale(schoolSeg)              0.095455  0.059914  2.54   0.1111    
## scale(pre)                   -0.086799  0.043651  3.95   0.0468 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)      3.6   0.255
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.834   0.132
## Number of clusters:   2063  Maximum cluster size: 450
#fit4<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBOMB, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fit4)
library(huxtable)
## 
## Attaching package: 'huxtable'
## The following object is masked from 'package:StatMeasures':
## 
##     contents
## The following object is masked from 'package:dplyr':
## 
##     add_rownames
allmodels<-huxreg(fit,fit3,statistics = character(0)) 
allmodels
(1)(2)
(Intercept)3.465 ***3.761 ***
(0.617)   (0.724)   
age-0.088 *  -0.095 *  
(0.038)   (0.041)   
I(age^2)0.002 ** 0.002 ** 
(0.001)   (0.001)   
sexW5Female0.158    0.121    
(0.121)   (0.119)   
raceXsexob_black.Straight0.381 *  0.218    
(0.154)   (0.163)   
raceXsexoc_hispanic.Straight0.237    -0.035    
(0.246)   (0.269)   
raceXsexoa_white.Gay-0.398 *  -0.397 *  
(0.160)   (0.163)   
raceXsexob_black.Gay0.175    -0.039    
(0.221)   (0.228)   
raceXsexoc_hispanic.Gay0.897    0.499    
(0.698)   (0.629)   
raceXsexoa_white.Bisexual0.304    0.315    
(0.187)   (0.179)   
raceXsexob_black.Bisexual-0.195    -0.426    
(0.398)   (0.409)   
raceXsexoc_hispanic.Bisexual0.223    -0.133    
(0.327)   (0.330)   
education1-0.144    -0.166    
(0.204)   (0.203)   
education2-0.431    -0.450 *  
(0.221)   (0.224)   
employmentempNoBene0.093    0.085    
(0.195)   (0.193)   
employmentempPartBene0.143    0.071    
(0.177)   (0.171)   
employmentunemployed0.465 *  0.416 *  
(0.182)   (0.177)   
hhincomeb_30k<50k-0.302    -0.279    
(0.158)   (0.149)   
hhincomec_50k>100k-0.019    -0.021    
(0.155)   (0.152)   
hhincomed_>100k-0.096    -0.091    
(0.187)   (0.182)   
marriage10.017    0.028    
(0.105)   (0.109)   
marriage20.233    0.248    
(0.206)   (0.211)   
livesParent2lives_1_parent        0.238    
        (0.161)   
livesParent2lives_other        -0.110    
        (0.332)   
mom.educ1        -0.360 *  
        (0.169)   
mom.educ2        -0.488 ** 
        (0.188)   
mom.pubass1        0.357    
        (0.232)   
closeMom2b_close        0.005    
        (0.094)   
parentSupervise2b_moderate        0.095    
        (0.238)   
parentSupervise2c_strict        0.348    
        (0.236)   
religion1        -0.147    
        (0.123)   
schoolconflict3        0.050    
        (0.132)   
schoolconflict4        -0.034    
        (0.154)   
schoolClubJoin1        0.293    
        (0.185)   
scale(schoolSeg)        0.095    
        (0.060)   
scale(pre)        -0.087 *  
        (0.044)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
nobs(fit3)
## [1] 24657
addhealthwh<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

addhealthwhi<-addhealthwh%>%
  subset(raceW5min=='a_white')

addhealthWhite<-addhealthwhi%>%
  filter(complete.cases(.))
addhealthWhite$livesParent2<-factor(addhealthWhite$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthWhite$closeMom<-as.numeric(addhealthWhite$closeMom)
addhealthWhite$closeMom2<-Recode(addhealthWhite$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthWhite$parentSupervise<-as.numeric(addhealthWhite$parentSupervise)
#table(addhealthWhite$parentSupervise)
addhealthWhite$parentSupervise2<-Recode(addhealthWhite$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthWhite$teacherTrouble2<-as.numeric(addhealthWhite$teacherTrouble)
addhealthWhite$studentTrouble2<-as.numeric(addhealthWhite$studentTrouble)
addhealthWhite$schoolconflict<-addhealthWhite$studentTrouble2+addhealthWhite$teacherTrouble2
addhealthWhite$schoolconflict<-as.factor(addhealthWhite$schoolconflict)
addhealthWhite$raceXsexo<-as.factor(addhealthWhite$raceXsexo)
addhealthWhite$raceXsexo <- droplevels(addhealthWhite$raceXsexo)
table(addhealthWhite$raceXsexo)
## 
## a_white.Straight      a_white.Gay a_white.Bisexual 
##            15312              490              609
fitWh<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthWhite, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitWh)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage, data = addhealthWhite, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                            Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                2.810140  0.567034 24.56  7.2e-07 ***
## age                       -0.031307  0.030707  1.04    0.308    
## I(age^2)                   0.000898  0.000649  1.91    0.166    
## sexW5Female                0.100711  0.135129  0.56    0.456    
## raceXsexoa_white.Gay      -0.350094  0.174113  4.04    0.044 *  
## raceXsexoa_white.Bisexual  0.223013  0.187780  1.41    0.235    
## education1                -0.177601  0.269628  0.43    0.510    
## education2                -0.532091  0.285953  3.46    0.063 .  
## employmentempNoBene        0.083545  0.224884  0.14    0.710    
## employmentempPartBene      0.135111  0.194388  0.48    0.487    
## employmentunemployed       0.444388  0.201968  4.84    0.028 *  
## hhincomeb_30k<50k         -0.150562  0.180686  0.69    0.405    
## hhincomec_50k>100k         0.160655  0.145681  1.22    0.270    
## hhincomed_>100k            0.071093  0.180332  0.16    0.693    
## marriage1                 -0.021352  0.127952  0.03    0.867    
## marriage2                  0.128264  0.259592  0.24    0.621    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.63   0.321
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.557  0.0845
## Number of clusters:   1532  Maximum cluster size: 392
#education+employment+hhincome+marriage
#fit2<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthWhite, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fit2)

fitWh3<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthWhite, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitWh3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage + livesParent2 + 
##     mom.educ + mom.pubass + closeMom2 + parentSupervise2 + religion + 
##     schoolconflict + schoolClubJoin + scale(schoolSeg) + scale(pre), 
##     data = addhealthWhite, weights = W5BIOWGT/mean(W5BIOWGT, 
##         na.rm = T), id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                             Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                 2.925793  0.637485 21.06  4.4e-06 ***
## age                        -0.031069  0.030807  1.02   0.3132    
## I(age^2)                    0.000887  0.000655  1.84   0.1755    
## sexW5Female                 0.079847  0.132689  0.36   0.5473    
## raceXsexoa_white.Gay       -0.366043  0.177494  4.25   0.0392 *  
## raceXsexoa_white.Bisexual   0.232861  0.189703  1.51   0.2196    
## education1                 -0.159683  0.268463  0.35   0.5520    
## education2                 -0.488616  0.285778  2.92   0.0873 .  
## employmentempNoBene         0.114681  0.220798  0.27   0.6035    
## employmentempPartBene       0.056384  0.189687  0.09   0.7663    
## employmentunemployed        0.394450  0.198933  3.93   0.0474 *  
## hhincomeb_30k<50k          -0.149004  0.181718  0.67   0.4122    
## hhincomec_50k>100k          0.166551  0.145649  1.31   0.2528    
## hhincomed_>100k             0.077205  0.181304  0.18   0.6702    
## marriage1                  -0.018830  0.128482  0.02   0.8835    
## marriage2                   0.129446  0.264252  0.24   0.6242    
## livesParent2lives_1_parent  0.140914  0.196087  0.52   0.4724    
## livesParent2lives_other     0.124124  0.394499  0.10   0.7530    
## mom.educ1                  -0.432673  0.198002  4.78   0.0289 *  
## mom.educ2                  -0.569843  0.212356  7.20   0.0073 ** 
## mom.pubass1                 0.536387  0.321916  2.78   0.0957 .  
## closeMom2b_close            0.012924  0.116646  0.01   0.9118    
## parentSupervise2b_moderate  0.097124  0.274416  0.13   0.7234    
## parentSupervise2c_strict    0.370353  0.268900  1.90   0.1684    
## religion1                   0.006613  0.069894  0.01   0.9246    
## schoolconflict3            -0.026086  0.151280  0.03   0.8631    
## schoolconflict4             0.020142  0.174150  0.01   0.9079    
## schoolClubJoin1             0.304592  0.211260  2.08   0.1494    
## scale(schoolSeg)            0.115825  0.060799  3.63   0.0568 .  
## scale(pre)                 -0.108676  0.039955  7.40   0.0065 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.52   0.301
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.559   0.111
## Number of clusters:   1532  Maximum cluster size: 392
#fit4<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthWhite, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fit4)
#library(huxtable)
allmodelsWh<-huxreg(fitWh,fitWh3,statistics = character(0)) 
allmodelsWh
(1)(2)
(Intercept)2.810 ***2.926 ***
(0.567)   (0.637)   
age-0.031    -0.031    
(0.031)   (0.031)   
I(age^2)0.001    0.001    
(0.001)   (0.001)   
sexW5Female0.101    0.080    
(0.135)   (0.133)   
raceXsexoa_white.Gay-0.350 *  -0.366 *  
(0.174)   (0.177)   
raceXsexoa_white.Bisexual0.223    0.233    
(0.188)   (0.190)   
education1-0.178    -0.160    
(0.270)   (0.268)   
education2-0.532    -0.489    
(0.286)   (0.286)   
employmentempNoBene0.084    0.115    
(0.225)   (0.221)   
employmentempPartBene0.135    0.056    
(0.194)   (0.190)   
employmentunemployed0.444 *  0.394 *  
(0.202)   (0.199)   
hhincomeb_30k<50k-0.151    -0.149    
(0.181)   (0.182)   
hhincomec_50k>100k0.161    0.167    
(0.146)   (0.146)   
hhincomed_>100k0.071    0.077    
(0.180)   (0.181)   
marriage1-0.021    -0.019    
(0.128)   (0.128)   
marriage20.128    0.129    
(0.260)   (0.264)   
livesParent2lives_1_parent        0.141    
        (0.196)   
livesParent2lives_other        0.124    
        (0.394)   
mom.educ1        -0.433 *  
        (0.198)   
mom.educ2        -0.570 ** 
        (0.212)   
mom.pubass1        0.536    
        (0.322)   
closeMom2b_close        0.013    
        (0.117)   
parentSupervise2b_moderate        0.097    
        (0.274)   
parentSupervise2c_strict        0.370    
        (0.269)   
religion1        0.007    
        (0.070)   
schoolconflict3        -0.026    
        (0.151)   
schoolconflict4        0.020    
        (0.174)   
schoolClubJoin1        0.305    
        (0.211)   
scale(schoolSeg)        0.116    
        (0.061)   
scale(pre)        -0.109 ** 
        (0.040)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
addhealthM<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

addhealthMa<-addhealthM%>%
  subset(sexW5=='Male')


addhealthMale<-addhealthMa%>%
  filter(complete.cases(.))
addhealthMale$livesParent2<-factor(addhealthMale$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthMale$closeMom<-as.numeric(addhealthMale$closeMom)
addhealthMale$closeMom2<-Recode(addhealthMale$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthMale$parentSupervise<-as.numeric(addhealthMale$parentSupervise)
#table(addhealthMale$parentSupervise)
addhealthMale$parentSupervise2<-Recode(addhealthMale$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthMale$teacherTrouble2<-as.numeric(addhealthMale$teacherTrouble)
addhealthMale$studentTrouble2<-as.numeric(addhealthMale$studentTrouble)
addhealthMale$schoolconflict<-addhealthMale$studentTrouble2+addhealthMale$teacherTrouble2
addhealthMale$schoolconflict<-as.factor(addhealthMale$schoolconflict)
addhealthMale$raceXsexo<-factor(addhealthMale$raceXsexo)
table(addhealthMale$raceXsexo)
## 
##    a_white.Straight    b_black.Straight c_hispanic.Straight         a_white.Gay 
##                5846                2887                 643                 188 
##         b_black.Gay      c_hispanic.Gay    a_white.Bisexual 
##                  92                  13                  48
fitM<-geeglm(score~age+I(age^2)+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthMale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitM)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + raceXsexo + education + 
##     employment + hhincome + marriage, data = addhealthMale, weights = W5BIOWGT/mean(W5BIOWGT, 
##     na.rm = T), id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                              Estimate  Std.err  Wald Pr(>|W|)    
## (Intercept)                   4.22685  0.87293 23.45  1.3e-06 ***
## age                          -0.15694  0.06021  6.79   0.0091 ** 
## I(age^2)                      0.00344  0.00117  8.63   0.0033 ** 
## raceXsexob_black.Straight    -0.01165  0.30750  0.00   0.9698    
## raceXsexoc_hispanic.Straight  0.15088  0.41525  0.13   0.7163    
## raceXsexoa_white.Gay         -0.50843  0.34139  2.22   0.1364    
## raceXsexob_black.Gay         -0.09098  0.71681  0.02   0.8990    
## raceXsexoc_hispanic.Gay       1.01029  0.93637  1.16   0.2806    
## raceXsexoa_white.Bisexual     0.31074  0.47721  0.42   0.5150    
## education1                    0.19883  0.13938  2.03   0.1537    
## education2                   -0.14912  0.15676  0.90   0.3415    
## employmentempNoBene          -0.15896  0.30748  0.27   0.6052    
## employmentempPartBene         0.02251  0.29305  0.01   0.9388    
## employmentunemployed          0.14147  0.40020  0.12   0.7237    
## hhincomeb_30k<50k            -0.42724  0.23161  3.40   0.0651 .  
## hhincomec_50k>100k           -0.10488  0.24239  0.19   0.6652    
## hhincomed_>100k              -0.37208  0.25952  2.06   0.1516    
## marriage1                    -0.09106  0.18998  0.23   0.6317    
## marriage2                     0.11396  0.26340  0.19   0.6653    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.21   0.306
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.703  0.0991
## Number of clusters:   737  Maximum cluster size: 450
#education+employment+hhincome+marriage
#fitM2<-geeglm(score~age+I(age^2)+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthMale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitM2)

fitM3<-geeglm(score~age+I(age^2)+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthMale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitM3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + raceXsexo + education + 
##     employment + hhincome + marriage + livesParent2 + mom.educ + 
##     mom.pubass + closeMom2 + parentSupervise2 + religion + schoolconflict + 
##     schoolClubJoin + scale(schoolSeg) + scale(pre), data = addhealthMale, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                              Estimate  Std.err  Wald Pr(>|W|)    
## (Intercept)                   5.30907  1.14633 21.45  3.6e-06 ***
## age                          -0.17267  0.06241  7.65   0.0057 ** 
## I(age^2)                      0.00364  0.00117  9.61   0.0019 ** 
## raceXsexob_black.Straight    -0.11123  0.31145  0.13   0.7210    
## raceXsexoc_hispanic.Straight  0.13942  0.43207  0.10   0.7469    
## raceXsexoa_white.Gay         -0.55034  0.31293  3.09   0.0786 .  
## raceXsexob_black.Gay         -0.15838  0.60509  0.07   0.7935    
## raceXsexoc_hispanic.Gay       0.84909  0.87434  0.94   0.3315    
## raceXsexoa_white.Bisexual     0.21492  0.37331  0.33   0.5648    
## education1                    0.12981  0.16484  0.62   0.4310    
## education2                   -0.25042  0.19903  1.58   0.2083    
## employmentempNoBene          -0.15489  0.30346  0.26   0.6098    
## employmentempPartBene         0.05821  0.28806  0.04   0.8399    
## employmentunemployed          0.18728  0.36023  0.27   0.6031    
## hhincomeb_30k<50k            -0.32866  0.20729  2.51   0.1128    
## hhincomec_50k>100k           -0.09808  0.22357  0.19   0.6609    
## hhincomed_>100k              -0.37283  0.23419  2.53   0.1114    
## marriage1                    -0.04189  0.20299  0.04   0.8365    
## marriage2                     0.17150  0.27428  0.39   0.5318    
## livesParent2lives_1_parent   -0.07800  0.28660  0.07   0.7855    
## livesParent2lives_other      -0.31027  0.71035  0.19   0.6623    
## mom.educ1                    -0.57358  0.33184  2.99   0.0839 .  
## mom.educ2                    -0.57345  0.35003  2.68   0.1014    
## mom.pubass1                  -0.18391  0.33165  0.31   0.5792    
## closeMom2b_close              0.01080  0.13751  0.01   0.9374    
## parentSupervise2b_moderate   -0.42140  0.41298  1.04   0.3075    
## parentSupervise2c_strict      0.11430  0.41283  0.08   0.7819    
## religion1                    -0.29426  0.19562  2.26   0.1325    
## schoolconflict3              -0.05308  0.21739  0.06   0.8071    
## schoolconflict4              -0.22930  0.25373  0.82   0.3662    
## schoolClubJoin1               0.30088  0.32016  0.88   0.3473    
## scale(schoolSeg)              0.19923  0.09243  4.65   0.0311 *  
## scale(pre)                    0.03139  0.09229  0.12   0.7337    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.17   0.352
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.771   0.172
## Number of clusters:   737  Maximum cluster size: 450
#fitM4<-geeglm(score~age+I(age^2)+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthMale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitM4)
#library(huxtable)
allmodelsM<-huxreg(fitM,fitM3,statistics = character(0)) 
allmodelsM
(1)(2)
(Intercept)4.227 ***5.309 ***
(0.873)   (1.146)   
age-0.157 ** -0.173 ** 
(0.060)   (0.062)   
I(age^2)0.003 ** 0.004 ** 
(0.001)   (0.001)   
raceXsexob_black.Straight-0.012    -0.111    
(0.307)   (0.311)   
raceXsexoc_hispanic.Straight0.151    0.139    
(0.415)   (0.432)   
raceXsexoa_white.Gay-0.508    -0.550    
(0.341)   (0.313)   
raceXsexob_black.Gay-0.091    -0.158    
(0.717)   (0.605)   
raceXsexoc_hispanic.Gay1.010    0.849    
(0.936)   (0.874)   
raceXsexoa_white.Bisexual0.311    0.215    
(0.477)   (0.373)   
education10.199    0.130    
(0.139)   (0.165)   
education2-0.149    -0.250    
(0.157)   (0.199)   
employmentempNoBene-0.159    -0.155    
(0.307)   (0.303)   
employmentempPartBene0.023    0.058    
(0.293)   (0.288)   
employmentunemployed0.141    0.187    
(0.400)   (0.360)   
hhincomeb_30k<50k-0.427    -0.329    
(0.232)   (0.207)   
hhincomec_50k>100k-0.105    -0.098    
(0.242)   (0.224)   
hhincomed_>100k-0.372    -0.373    
(0.260)   (0.234)   
marriage1-0.091    -0.042    
(0.190)   (0.203)   
marriage20.114    0.171    
(0.263)   (0.274)   
livesParent2lives_1_parent        -0.078    
        (0.287)   
livesParent2lives_other        -0.310    
        (0.710)   
mom.educ1        -0.574    
        (0.332)   
mom.educ2        -0.573    
        (0.350)   
mom.pubass1        -0.184    
        (0.332)   
closeMom2b_close        0.011    
        (0.138)   
parentSupervise2b_moderate        -0.421    
        (0.413)   
parentSupervise2c_strict        0.114    
        (0.413)   
religion1        -0.294    
        (0.196)   
schoolconflict3        -0.053    
        (0.217)   
schoolconflict4        -0.229    
        (0.254)   
schoolClubJoin1        0.301    
        (0.320)   
scale(schoolSeg)        0.199 *  
        (0.092)   
scale(pre)        0.031    
        (0.092)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
addhealthFem<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

addhealthFema<-addhealthFem%>%
  subset(sexW5=='Female')


addhealthFemale<-addhealthFema%>%
  filter(complete.cases(.))
addhealthFemale$livesParent2<-factor(addhealthFemale$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthFemale$closeMom<-as.numeric(addhealthFemale$closeMom)
addhealthFemale$closeMom2<-Recode(addhealthFemale$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthFemale$parentSupervise<-as.numeric(addhealthFemale$parentSupervise)
#table(addhealthFemale$parentSupervise)
addhealthFemale$parentSupervise2<-Recode(addhealthFemale$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthFemale$teacherTrouble2<-as.numeric(addhealthFemale$teacherTrouble)
addhealthFemale$studentTrouble2<-as.numeric(addhealthFemale$studentTrouble)
addhealthFemale$schoolconflict<-addhealthFemale$studentTrouble2+addhealthFemale$teacherTrouble2
addhealthFemale$schoolconflict<-as.factor(addhealthFemale$schoolconflict)
addhealthFemale$raceXsexo<-factor(addhealthFemale$raceXsexo)
table(addhealthFemale$raceXsexo)
## 
##    a_white.Straight    b_black.Straight c_hispanic.Straight         a_white.Gay 
##                9466                3282                1081                 302 
##         b_black.Gay      c_hispanic.Gay    a_white.Bisexual    b_black.Bisexual 
##                 120                  30                 561                  70 
## c_hispanic.Bisexual 
##                  28
fitF<-geeglm(score~age+I(age^2)+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthFemale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitF)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + raceXsexo + education + 
##     employment + hhincome + marriage, data = addhealthFemale, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                               Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                   3.078619  0.639292 23.19  1.5e-06 ***
## age                          -0.026879  0.032320  0.69  0.40561    
## I(age^2)                      0.000722  0.000679  1.13  0.28794    
## raceXsexob_black.Straight     0.643001  0.160909 15.97  6.4e-05 ***
## raceXsexoc_hispanic.Straight  0.259608  0.260898  0.99  0.31971    
## raceXsexoa_white.Gay         -0.389102  0.286384  1.85  0.17425    
## raceXsexob_black.Gay          0.651857  0.309356  4.44  0.03511 *  
## raceXsexoc_hispanic.Gay       2.241957  0.448893 24.94  5.9e-07 ***
## raceXsexoa_white.Bisexual     0.268287  0.206382  1.69  0.19362    
## raceXsexob_black.Bisexual     0.159300  0.386485  0.17  0.68021    
## raceXsexoc_hispanic.Bisexual  0.108363  0.332551  0.11  0.74453    
## education1                   -0.636224  0.370124  2.95  0.08562 .  
## education2                   -0.884396  0.408509  4.69  0.03039 *  
## employmentempNoBene           0.319713  0.246741  1.68  0.19506    
## employmentempPartBene         0.287783  0.226474  1.61  0.20383    
## employmentunemployed          0.679409  0.190497 12.72  0.00036 ***
## hhincomeb_30k<50k            -0.076757  0.173049  0.20  0.65736    
## hhincomec_50k>100k            0.201772  0.141433  2.04  0.15369    
## hhincomed_>100k               0.258528  0.192910  1.80  0.18020    
## marriage1                     0.102510  0.117089  0.77  0.38131    
## marriage2                     0.313539  0.302708  1.07  0.30030    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.79   0.361
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha     0.65  0.0897
## Number of clusters:   1326  Maximum cluster size: 324
#education+employment+hhincome+marriage
#fitF2<-geeglm(score~age+I(age^2)+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthFemale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitF2)

fitF3<-geeglm(score~age+I(age^2)+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthFemale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitF3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + raceXsexo + education + 
##     employment + hhincome + marriage + livesParent2 + mom.educ + 
##     mom.pubass + closeMom2 + parentSupervise2 + religion + schoolconflict + 
##     schoolClubJoin + scale(schoolSeg) + scale(pre), data = addhealthFemale, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                               Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                   2.498749  0.687935 13.19  0.00028 ***
## age                          -0.014099  0.033144  0.18  0.67055    
## I(age^2)                      0.000477  0.000701  0.46  0.49609    
## raceXsexob_black.Straight     0.433277  0.180378  5.77  0.01630 *  
## raceXsexoc_hispanic.Straight -0.038537  0.308913  0.02  0.90072    
## raceXsexoa_white.Gay         -0.499256  0.291505  2.93  0.08677 .  
## raceXsexob_black.Gay          0.737682  0.253153  8.49  0.00357 ** 
## raceXsexoc_hispanic.Gay       1.817272  0.847717  4.60  0.03206 *  
## raceXsexoa_white.Bisexual     0.300520  0.216586  1.93  0.16528    
## raceXsexob_black.Bisexual    -0.539488  0.459133  1.38  0.23999    
## raceXsexoc_hispanic.Bisexual -0.078915  0.406979  0.04  0.84625    
## education1                   -0.672594  0.383559  3.07  0.07951 .  
## education2                   -0.842235  0.427835  3.88  0.04900 *  
## employmentempNoBene           0.344285  0.237596  2.10  0.14733    
## employmentempPartBene         0.087788  0.205074  0.18  0.66859    
## employmentunemployed          0.534909  0.197520  7.33  0.00677 ** 
## hhincomeb_30k<50k            -0.067419  0.179373  0.14  0.70702    
## hhincomec_50k>100k            0.246666  0.144715  2.91  0.08829 .  
## hhincomed_>100k               0.341232  0.208063  2.69  0.10100    
## marriage1                     0.127253  0.121200  1.10  0.29374    
## marriage2                     0.338644  0.309575  1.20  0.27400    
## livesParent2lives_1_parent    0.435899  0.165264  6.96  0.00835 ** 
## livesParent2lives_other      -0.272784  0.375222  0.53  0.46723    
## mom.educ1                    -0.172042  0.194194  0.78  0.37566    
## mom.educ2                    -0.462178  0.220485  4.39  0.03607 *  
## mom.pubass1                   0.742987  0.263850  7.93  0.00486 ** 
## closeMom2b_close             -0.092632  0.125681  0.54  0.46110    
## parentSupervise2b_moderate    0.565146  0.237971  5.64  0.01756 *  
## parentSupervise2c_strict      0.631585  0.234093  7.28  0.00698 ** 
## religion1                     0.051074  0.082478  0.38  0.53575    
## schoolconflict3               0.209267  0.165940  1.59  0.20727    
## schoolconflict4               0.084288  0.180958  0.22  0.64137    
## schoolClubJoin1               0.272321  0.207750  1.72  0.18992    
## scale(schoolSeg)             -0.009204  0.073734  0.02  0.90066    
## scale(pre)                   -0.148057  0.043487 11.59  0.00066 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.96   0.386
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha     1.11   0.257
## Number of clusters:   1326  Maximum cluster size: 324
#fitF4<-geeglm(score~age+I(age^2)+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthFemale, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitF4)
#library(huxtable)
allmodelsF<-huxreg(fitF,fitF3,statistics = character(0)) 
allmodelsF
(1)(2)
(Intercept)3.079 ***2.499 ***
(0.639)   (0.688)   
age-0.027    -0.014    
(0.032)   (0.033)   
I(age^2)0.001    0.000    
(0.001)   (0.001)   
raceXsexob_black.Straight0.643 ***0.433 *  
(0.161)   (0.180)   
raceXsexoc_hispanic.Straight0.260    -0.039    
(0.261)   (0.309)   
raceXsexoa_white.Gay-0.389    -0.499    
(0.286)   (0.292)   
raceXsexob_black.Gay0.652 *  0.738 ** 
(0.309)   (0.253)   
raceXsexoc_hispanic.Gay2.242 ***1.817 *  
(0.449)   (0.848)   
raceXsexoa_white.Bisexual0.268    0.301    
(0.206)   (0.217)   
raceXsexob_black.Bisexual0.159    -0.539    
(0.386)   (0.459)   
raceXsexoc_hispanic.Bisexual0.108    -0.079    
(0.333)   (0.407)   
education1-0.636    -0.673    
(0.370)   (0.384)   
education2-0.884 *  -0.842 *  
(0.409)   (0.428)   
employmentempNoBene0.320    0.344    
(0.247)   (0.238)   
employmentempPartBene0.288    0.088    
(0.226)   (0.205)   
employmentunemployed0.679 ***0.535 ** 
(0.190)   (0.198)   
hhincomeb_30k<50k-0.077    -0.067    
(0.173)   (0.179)   
hhincomec_50k>100k0.202    0.247    
(0.141)   (0.145)   
hhincomed_>100k0.259    0.341    
(0.193)   (0.208)   
marriage10.103    0.127    
(0.117)   (0.121)   
marriage20.314    0.339    
(0.303)   (0.310)   
livesParent2lives_1_parent        0.436 ** 
        (0.165)   
livesParent2lives_other        -0.273    
        (0.375)   
mom.educ1        -0.172    
        (0.194)   
mom.educ2        -0.462 *  
        (0.220)   
mom.pubass1        0.743 ** 
        (0.264)   
closeMom2b_close        -0.093    
        (0.126)   
parentSupervise2b_moderate        0.565 *  
        (0.238)   
parentSupervise2c_strict        0.632 ** 
        (0.234)   
religion1        0.051    
        (0.082)   
schoolconflict3        0.209    
        (0.166)   
schoolconflict4        0.084    
        (0.181)   
schoolClubJoin1        0.272    
        (0.208)   
scale(schoolSeg)        -0.009    
        (0.074)   
scale(pre)        -0.148 ***
        (0.043)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
addhealthBl<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

addhealthBla<-addhealthBl%>%
  subset(raceW5min=='b_black')

addhealthBlack<-addhealthBla%>%
  filter(complete.cases(.))
addhealthBlack$livesParent2<-factor(addhealthBlack$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthBlack$closeMom<-as.numeric(addhealthBlack$closeMom)
addhealthBlack$closeMom2<-Recode(addhealthBlack$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthBlack$parentSupervise<-as.numeric(addhealthBlack$parentSupervise)
#table(addhealthBlack$parentSupervise)
addhealthBlack$parentSupervise2<-Recode(addhealthBlack$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthBlack$teacherTrouble2<-as.numeric(addhealthBlack$teacherTrouble)
addhealthBlack$studentTrouble2<-as.numeric(addhealthBlack$studentTrouble)
addhealthBlack$schoolconflict<-addhealthBlack$studentTrouble2+addhealthBlack$teacherTrouble2
addhealthBlack$schoolconflict<-as.factor(addhealthBlack$schoolconflict)
addhealthBlack$raceXsexo<-factor(addhealthBlack$raceXsexo)
table(addhealthBlack$raceXsexo)
## 
## b_black.Straight      b_black.Gay b_black.Bisexual 
##             6169              212               70
fitB<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBlack, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitB)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage, data = addhealthBlack, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                           Estimate  Std.err  Wald Pr(>|W|)    
## (Intercept)                5.66537  1.37253 17.04  3.7e-05 ***
## age                       -0.25985  0.08813  8.69   0.0032 ** 
## I(age^2)                   0.00531  0.00164 10.50   0.0012 ** 
## sexW5Female                0.64609  0.31170  4.30   0.0382 *  
## raceXsexob_black.Gay      -0.33687  0.19282  3.05   0.0806 .  
## raceXsexob_black.Bisexual -0.56838  0.47390  1.44   0.2304    
## education1                -0.17947  0.10226  3.08   0.0793 .  
## education2                -0.15006  0.16925  0.79   0.3753    
## employmentempNoBene        0.10722  0.40546  0.07   0.7914    
## employmentempPartBene      0.49690  0.52581  0.89   0.3446    
## employmentunemployed       0.83273  0.46445  3.21   0.0730 .  
## hhincomeb_30k<50k         -0.45002  0.26110  2.97   0.0848 .  
## hhincomec_50k>100k        -0.44935  0.28290  2.52   0.1122    
## hhincomed_>100k           -0.44841  0.37636  1.42   0.2335    
## marriage1                  0.20755  0.23772  0.76   0.3826    
## marriage2                  0.40471  0.38797  1.09   0.2969    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.52   0.631
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.991   0.262
## Number of clusters:   367  Maximum cluster size: 450
#education+employment+hhincome+marriage
#fitB2<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBlack, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitB2)

fitB3<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBlack, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitB3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage + livesParent2 + 
##     mom.educ + mom.pubass + closeMom2 + parentSupervise2 + religion + 
##     schoolconflict + schoolClubJoin + scale(schoolSeg) + scale(pre), 
##     data = addhealthBlack, weights = W5BIOWGT/mean(W5BIOWGT, 
##         na.rm = T), id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                            Estimate  Std.err  Wald Pr(>|W|)    
## (Intercept)                 5.77135  1.44674 15.91  6.6e-05 ***
## age                        -0.25706  0.08754  8.62   0.0033 ** 
## I(age^2)                    0.00515  0.00157 10.80   0.0010 ** 
## sexW5Female                 0.53514  0.30411  3.10   0.0785 .  
## raceXsexob_black.Gay       -0.45596  0.20778  4.82   0.0282 *  
## raceXsexob_black.Bisexual  -0.56432  0.46337  1.48   0.2233    
## education1                 -0.33999  0.17437  3.80   0.0512 .  
## education2                 -0.32926  0.28623  1.32   0.2500    
## employmentempNoBene        -0.19516  0.41289  0.22   0.6364    
## employmentempPartBene       0.34731  0.53148  0.43   0.5134    
## employmentunemployed        0.64803  0.49787  1.69   0.1931    
## hhincomeb_30k<50k          -0.41942  0.23688  3.13   0.0766 .  
## hhincomec_50k>100k         -0.44730  0.26349  2.88   0.0896 .  
## hhincomed_>100k            -0.42287  0.35270  1.44   0.2306    
## marriage1                   0.22065  0.24514  0.81   0.3681    
## marriage2                   0.43447  0.37411  1.35   0.2455    
## livesParent2lives_1_parent  0.39896  0.28253  1.99   0.1579    
## livesParent2lives_other    -0.79943  0.53078  2.27   0.1320    
## mom.educ1                  -0.46354  0.44020  1.11   0.2923    
## mom.educ2                  -0.25148  0.48506  0.27   0.6041    
## mom.pubass1                 0.25789  0.38737  0.44   0.5056    
## closeMom2b_close            0.01210  0.20447  0.00   0.9528    
## parentSupervise2b_moderate  0.22281  0.48763  0.21   0.6477    
## parentSupervise2c_strict    0.32550  0.47791  0.46   0.4958    
## religion1                  -0.27117  0.21848  1.54   0.2145    
## schoolconflict3             0.37831  0.33492  1.28   0.2587    
## schoolconflict4             0.13873  0.34572  0.16   0.6882    
## schoolClubJoin1             0.95543  0.46439  4.23   0.0396 *  
## scale(schoolSeg)           -0.08413  0.14201  0.35   0.5536    
## scale(pre)                 -0.02634  0.14798  0.03   0.8587    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.64   0.693
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.964   0.407
## Number of clusters:   367  Maximum cluster size: 450
#fitB4<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthBlack, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitB4)
#library(huxtable)
allmodelsB<-huxreg(fitB,fitB3,statistics = character(0)) 
allmodelsB
(1)(2)
(Intercept)5.665 ***5.771 ***
(1.373)   (1.447)   
age-0.260 ** -0.257 ** 
(0.088)   (0.088)   
I(age^2)0.005 ** 0.005 ** 
(0.002)   (0.002)   
sexW5Female0.646 *  0.535    
(0.312)   (0.304)   
raceXsexob_black.Gay-0.337    -0.456 *  
(0.193)   (0.208)   
raceXsexob_black.Bisexual-0.568    -0.564    
(0.474)   (0.463)   
education1-0.179    -0.340    
(0.102)   (0.174)   
education2-0.150    -0.329    
(0.169)   (0.286)   
employmentempNoBene0.107    -0.195    
(0.405)   (0.413)   
employmentempPartBene0.497    0.347    
(0.526)   (0.531)   
employmentunemployed0.833    0.648    
(0.464)   (0.498)   
hhincomeb_30k<50k-0.450    -0.419    
(0.261)   (0.237)   
hhincomec_50k>100k-0.449    -0.447    
(0.283)   (0.263)   
hhincomed_>100k-0.448    -0.423    
(0.376)   (0.353)   
marriage10.208    0.221    
(0.238)   (0.245)   
marriage20.405    0.434    
(0.388)   (0.374)   
livesParent2lives_1_parent        0.399    
        (0.283)   
livesParent2lives_other        -0.799    
        (0.531)   
mom.educ1        -0.464    
        (0.440)   
mom.educ2        -0.251    
        (0.485)   
mom.pubass1        0.258    
        (0.387)   
closeMom2b_close        0.012    
        (0.204)   
parentSupervise2b_moderate        0.223    
        (0.488)   
parentSupervise2c_strict        0.326    
        (0.478)   
religion1        -0.271    
        (0.218)   
schoolconflict3        0.378    
        (0.335)   
schoolconflict4        0.139    
        (0.346)   
schoolClubJoin1        0.955 *  
        (0.464)   
scale(schoolSeg)        -0.084    
        (0.142)   
scale(pre)        -0.026    
        (0.148)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
addhealthHi<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

addhealthHis<-addhealthHi%>%
  subset(raceW5min=='c_hispanic')

addhealthHispanic<-addhealthHis%>%
  filter(complete.cases(.))
addhealthHispanic$livesParent2<-factor(addhealthHispanic$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthHispanic$closeMom<-as.numeric(addhealthHispanic$closeMom)
addhealthHispanic$closeMom2<-Recode(addhealthHispanic$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthHispanic$parentSupervise<-as.numeric(addhealthHispanic$parentSupervise)
#table(addhealthHispanic$parentSupervise)
addhealthHispanic$parentSupervise2<-Recode(addhealthHispanic$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthHispanic$teacherTrouble2<-as.numeric(addhealthHispanic$teacherTrouble)
addhealthHispanic$studentTrouble2<-as.numeric(addhealthHispanic$studentTrouble)
addhealthHispanic$schoolconflict<-addhealthHispanic$studentTrouble2+addhealthHispanic$teacherTrouble2
addhealthHispanic$schoolconflict<-as.factor(addhealthHispanic$schoolconflict)
addhealthHispanic$raceXsexo<-as.factor(addhealthHispanic$raceXsexo)
addhealthHispanic$raceXsexo <- droplevels(addhealthHispanic$raceXsexo)
table(addhealthHispanic$raceXsexo)
## 
## c_hispanic.Straight      c_hispanic.Gay c_hispanic.Bisexual 
##                1724                  43                  28
fitH<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthHispanic, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitH)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage, data = addhealthHispanic, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                               Estimate   Std.err Wald Pr(>|W|)   
## (Intercept)                   2.454657  0.876583 7.84   0.0051 **
## age                           0.023512  0.044766 0.28   0.5994   
## I(age^2)                     -0.000475  0.000913 0.27   0.6031   
## sexW5Female                   0.045089  0.516764 0.01   0.9305   
## raceXsexoc_hispanic.Gay       0.565429  0.491550 1.32   0.2500   
## raceXsexoc_hispanic.Bisexual  0.035445  0.209377 0.03   0.8656   
## education1                    0.272502  0.179368 2.31   0.1287   
## education2                   -0.017844  0.292044 0.00   0.9513   
## employmentempNoBene          -0.144671  0.777441 0.03   0.8524   
## employmentempPartBene        -0.865062  0.536248 2.60   0.1067   
## employmentunemployed          0.670840  0.852456 0.62   0.4313   
## hhincomeb_30k<50k            -0.192146  0.224673 0.73   0.3924   
## hhincomec_50k>100k           -0.103171  0.136918 0.57   0.4511   
## hhincomed_>100k              -0.314455  0.176674 3.17   0.0751 . 
## marriage1                     0.161644  0.385867 0.18   0.6753   
## marriage2                     0.203663  0.263185 0.60   0.4390   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.18   0.447
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.855    0.25
## Number of clusters:   164  Maximum cluster size: 144
#education+employment+hhincome+marriage
#fitH2<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthHispanic, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitH2)

fitH3<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthHispanic, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitH3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage + livesParent2 + 
##     mom.educ + mom.pubass + closeMom2 + parentSupervise2 + religion + 
##     schoolconflict + schoolClubJoin + scale(schoolSeg) + scale(pre), 
##     data = addhealthHispanic, weights = W5BIOWGT/mean(W5BIOWGT, 
##         na.rm = T), id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                               Estimate   Std.err  Wald Pr(>|W|)    
## (Intercept)                   3.873844  1.078085 12.91  0.00033 ***
## age                           0.028908  0.044596  0.42  0.51684    
## I(age^2)                     -0.000671  0.000923  0.53  0.46770    
## sexW5Female                  -0.575478  0.441900  1.70  0.19282    
## raceXsexoc_hispanic.Gay       0.726368  0.582213  1.56  0.21218    
## raceXsexoc_hispanic.Bisexual  0.350310  0.466421  0.56  0.45262    
## education1                    0.144758  0.217765  0.44  0.50622    
## education2                   -0.203820  0.279874  0.53  0.46646    
## employmentempNoBene           0.372174  0.580288  0.41  0.52129    
## employmentempPartBene        -1.108990  0.837685  1.75  0.18554    
## employmentunemployed          1.362607  0.726149  3.52  0.06059 .  
## hhincomeb_30k<50k            -0.125897  0.221052  0.32  0.56899    
## hhincomec_50k>100k           -0.059879  0.154492  0.15  0.69832    
## hhincomed_>100k              -0.271077  0.182578  2.20  0.13762    
## marriage1                     0.213604  0.379001  0.32  0.57303    
## marriage2                     0.246290  0.258663  0.91  0.34101    
## livesParent2lives_1_parent   -0.523378  0.531810  0.97  0.32504    
## livesParent2lives_other      -1.140102  0.899234  1.61  0.20485    
## mom.educ1                     0.916162  0.505433  3.29  0.06989 .  
## mom.educ2                    -0.755974  0.603713  1.57  0.21049    
## mom.pubass1                  -0.579308  0.612153  0.90  0.34397    
## closeMom2b_close             -0.251500  0.248975  1.02  0.31243    
## parentSupervise2b_moderate   -1.394671  0.718773  3.76  0.05234 .  
## parentSupervise2c_strict     -0.771564  0.754622  1.05  0.30657    
## religion1                    -0.096078  0.119911  0.64  0.42299    
## schoolconflict3               1.178578  0.539219  4.78  0.02884 *  
## schoolconflict4              -0.785504  0.600584  1.71  0.19091    
## schoolClubJoin1              -0.137420  0.479221  0.08  0.77430    
## scale(schoolSeg)              0.155736  0.265302  0.34  0.55720    
## scale(pre)                    0.212052  0.149876  2.00  0.15711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     2.99    0.63
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.723   0.467
## Number of clusters:   164  Maximum cluster size: 144
#fitH4<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoi#n+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthHispanic, #weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitH4)
#library(huxtable)
allmodelsH<-huxreg(fitH,fitH3,statistics = character(0)) 
allmodelsH
(1)(2)
(Intercept)2.455 **3.874 ***
(0.877)  (1.078)   
age0.024   0.029    
(0.045)  (0.045)   
I(age^2)-0.000   -0.001    
(0.001)  (0.001)   
sexW5Female0.045   -0.575    
(0.517)  (0.442)   
raceXsexoc_hispanic.Gay0.565   0.726    
(0.492)  (0.582)   
raceXsexoc_hispanic.Bisexual0.035   0.350    
(0.209)  (0.466)   
education10.273   0.145    
(0.179)  (0.218)   
education2-0.018   -0.204    
(0.292)  (0.280)   
employmentempNoBene-0.145   0.372    
(0.777)  (0.580)   
employmentempPartBene-0.865   -1.109    
(0.536)  (0.838)   
employmentunemployed0.671   1.363    
(0.852)  (0.726)   
hhincomeb_30k<50k-0.192   -0.126    
(0.225)  (0.221)   
hhincomec_50k>100k-0.103   -0.060    
(0.137)  (0.154)   
hhincomed_>100k-0.314   -0.271    
(0.177)  (0.183)   
marriage10.162   0.214    
(0.386)  (0.379)   
marriage20.204   0.246    
(0.263)  (0.259)   
livesParent2lives_1_parent       -0.523    
       (0.532)   
livesParent2lives_other       -1.140    
       (0.899)   
mom.educ1       0.916    
       (0.505)   
mom.educ2       -0.756    
       (0.604)   
mom.pubass1       -0.579    
       (0.612)   
closeMom2b_close       -0.251    
       (0.249)   
parentSupervise2b_moderate       -1.395    
       (0.719)   
parentSupervise2c_strict       -0.772    
       (0.755)   
religion1       -0.096    
       (0.120)   
schoolconflict3       1.179 *  
       (0.539)   
schoolconflict4       -0.786    
       (0.601)   
schoolClubJoin1       -0.137    
       (0.479)   
scale(schoolSeg)       0.156    
       (0.265)   
scale(pre)       0.212    
       (0.150)   
*** p < 0.001; ** p < 0.01; * p < 0.05.
addhealthL<-addhealthS1long%>%
  dplyr::select(wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,livesParent2,parentSupervise,mom.educ2,mom.educ,mom.employed,mom.nativity,mom.pubass,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,sexualorientation,religion,closeMom,score,raceXsexo,raceW5min)

#library(tidyr)
addhealthLG<-addhealthL%>%
  subset(sexualorientation=='Gay'|sexualorientation=="Bisexual")
addhealthLG$sexualorientation<-factor(addhealthLG$sexualorientation)

addhealthLG$raceXsexo<-interaction(addhealthLG$raceW5min,addhealthLG$sexualorientation)
addhealthLG$raceXsexo<-factor(addhealthLG$raceXsexo)

addhealthLGBTQ<-addhealthLG%>%
  filter(complete.cases(.))
addhealthLGBTQ$livesParent2<-factor(addhealthLGBTQ$livesParent2,levels=c("lives_2_parents","lives_1_parent","lives_other"))
addhealthLGBTQ$closeMom<-as.numeric(addhealthLGBTQ$closeMom)
addhealthLGBTQ$closeMom2<-Recode(addhealthLGBTQ$closeMom,recodes="1:3='a_notClose';4:5='b_close'")
addhealthLGBTQ$parentSupervise<-as.numeric(addhealthLGBTQ$parentSupervise)
#table(addhealthLGBTQ$parentSupervise)
addhealthLGBTQ$parentSupervise2<-Recode(addhealthLGBTQ$parentSupervise,recodes="1:3='a_easy';4:6='b_moderate';7:8='c_strict'")
addhealthLGBTQ$teacherTrouble2<-as.numeric(addhealthLGBTQ$teacherTrouble)
addhealthLGBTQ$studentTrouble2<-as.numeric(addhealthLGBTQ$studentTrouble)
addhealthLGBTQ$schoolconflict<-addhealthLGBTQ$studentTrouble2+addhealthLGBTQ$teacherTrouble2
addhealthLGBTQ$schoolconflict<-as.factor(addhealthLGBTQ$schoolconflict)
#install.packages("geepack")
#library(geepack)
#wave,aid,W5BIOWGT,age,sexW5,education,hhincome,employment,marriage,,,,,,,,pre,teacherTrouble,studentTrouble,schoolClubJoin,schoolSeg,

fitL<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthLGBTQ, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitL)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage, data = addhealthLGBTQ, 
##     weights = W5BIOWGT/mean(W5BIOWGT, na.rm = T), id = aid, waves = wave, 
##     corstr = "exchangeable")
## 
##  Coefficients:
##                              Estimate  Std.err Wald Pr(>|W|)   
## (Intercept)                   1.57212  2.28545 0.47   0.4915   
## age                           0.07130  0.14311 0.25   0.6183   
## I(age^2)                     -0.00127  0.00253 0.25   0.6147   
## sexW5Female                   0.00792  0.61050 0.00   0.9896   
## raceXsexob_black.Gay          1.70967  0.87584 3.81   0.0509 . 
## raceXsexoc_hispanic.Gay       1.43045  1.66016 0.74   0.3889   
## raceXsexoa_white.Bisexual     0.52741  0.29871 3.12   0.0775 . 
## raceXsexob_black.Bisexual     3.35179  1.18462 8.01   0.0047 **
## raceXsexoc_hispanic.Bisexual  1.20103  1.48668 0.65   0.4192   
## education1                   -0.55057  1.12586 0.24   0.6248   
## education2                   -1.08734  1.20407 0.82   0.3665   
## employmentempNoBene          -0.66919  0.86869 0.59   0.4411   
## employmentempPartBene         0.68322  0.73119 0.87   0.3501   
## employmentunemployed         -1.60702  0.71602 5.04   0.0248 * 
## hhincomeb_30k<50k            -0.23942  0.34515 0.48   0.4879   
## hhincomec_50k>100k            0.66625  0.25613 6.77   0.0093 **
## hhincomed_>100k               0.06250  0.36483 0.03   0.8640   
## marriage1                     0.48378  0.54472 0.79   0.3745   
## marriage2                    -0.42727  0.42582 1.01   0.3157   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     3.58   0.849
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.464   0.251
## Number of clusters:   107  Maximum cluster size: 108
#education+employment+hhincome+marriage
#fitL2<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthLGBTQ, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitL2)

fitL3<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+education+employment+hhincome+marriage+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+scale(schoolSeg)+scale(pre),id=aid , wave = wave, corstr ="exchangeable",   data=addhealthLGBTQ, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
summary(fitL3)
## 
## Call:
## geeglm(formula = score ~ age + I(age^2) + sexW5 + raceXsexo + 
##     education + employment + hhincome + marriage + livesParent2 + 
##     mom.educ + mom.pubass + closeMom2 + parentSupervise2 + religion + 
##     schoolconflict + schoolClubJoin + scale(schoolSeg) + scale(pre), 
##     data = addhealthLGBTQ, weights = W5BIOWGT/mean(W5BIOWGT, 
##         na.rm = T), id = aid, waves = wave, corstr = "exchangeable")
## 
##  Coefficients:
##                              Estimate  Std.err  Wald Pr(>|W|)    
## (Intercept)                   1.42176  2.32674  0.37   0.5412    
## age                           0.11172  0.13767  0.66   0.4171    
## I(age^2)                     -0.00157  0.00247  0.40   0.5249    
## sexW5Female                  -0.14146  0.51334  0.08   0.7829    
## raceXsexob_black.Gay          1.06082  0.90949  1.36   0.2435    
## raceXsexoc_hispanic.Gay       1.48145  1.47239  1.01   0.3143    
## raceXsexoa_white.Bisexual     0.25413  0.35223  0.52   0.4706    
## raceXsexob_black.Bisexual     3.91513  1.24795  9.84   0.0017 ** 
## raceXsexoc_hispanic.Bisexual -0.87056  1.30667  0.44   0.5053    
## education1                   -1.36058  1.52254  0.80   0.3715    
## education2                   -1.93949  1.59006  1.49   0.2226    
## employmentempNoBene          -1.75401  0.95039  3.41   0.0650 .  
## employmentempPartBene         0.76040  0.51717  2.16   0.1415    
## employmentunemployed         -0.31129  0.73838  0.18   0.6733    
## hhincomeb_30k<50k            -0.12108  0.37570  0.10   0.7472    
## hhincomec_50k>100k            0.90375  0.20384 19.66  9.3e-06 ***
## hhincomed_>100k               0.21928  0.32643  0.45   0.5017    
## marriage1                     0.28394  0.35736  0.63   0.4269    
## marriage2                    -0.31794  0.35697  0.79   0.3731    
## livesParent2lives_1_parent   -0.35521  0.56644  0.39   0.5306    
## livesParent2lives_other       3.33862  1.03256 10.45   0.0012 ** 
## mom.educ1                     1.10980  0.67836  2.68   0.1018    
## mom.educ2                     0.86709  0.69770  1.54   0.2139    
## mom.pubass1                  -0.41933  1.04078  0.16   0.6870    
## closeMom2b_close              0.31025  0.24392  1.62   0.2034    
## parentSupervise2b_moderate   -2.08377  0.67144  9.63   0.0019 ** 
## parentSupervise2c_strict     -0.05697  0.65492  0.01   0.9307    
## religion1                     1.03278  0.17360 35.39  2.7e-09 ***
## schoolconflict3              -0.11584  0.57974  0.04   0.8416    
## schoolconflict4               0.27792  0.66680  0.17   0.6768    
## schoolClubJoin1              -0.04422  0.92939  0.00   0.9620    
## scale(schoolSeg)             -0.01927  0.27462  0.00   0.9441    
## scale(pre)                   -0.56551  0.24892  5.16   0.0231 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     1.75   0.523
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha    0.251   0.164
## Number of clusters:   107  Maximum cluster size: 108
#fitL4<-geeglm(score~age+I(age^2)+sexW5+raceXsexo+livesParent2+mom.educ+mom.pubass+closeMom2+parentSupervise2+religion+schoolconflict+schoolClubJoin+schoolSeg+pre+education+employment+hhincome+marriage,id=aid , wave = wave, corstr ="exchangeable",   data=addhealthLGBTQ, weights=W5BIOWGT/mean(W5BIOWGT,na.rm=T))
#summary(fitL4)
#library(huxtable)
allmodelsL<-huxreg(fitL,fitL3,statistics = character(0)) 
allmodelsL
(1)(2)
(Intercept)1.572   1.422    
(2.285)  (2.327)   
age0.071   0.112    
(0.143)  (0.138)   
I(age^2)-0.001   -0.002    
(0.003)  (0.002)   
sexW5Female0.008   -0.141    
(0.611)  (0.513)   
raceXsexob_black.Gay1.710   1.061    
(0.876)  (0.909)   
raceXsexoc_hispanic.Gay1.430   1.481    
(1.660)  (1.472)   
raceXsexoa_white.Bisexual0.527   0.254    
(0.299)  (0.352)   
raceXsexob_black.Bisexual3.352 **3.915 ** 
(1.185)  (1.248)   
raceXsexoc_hispanic.Bisexual1.201   -0.871    
(1.487)  (1.307)   
education1-0.551   -1.361    
(1.126)  (1.523)   
education2-1.087   -1.939    
(1.204)  (1.590)   
employmentempNoBene-0.669   -1.754    
(0.869)  (0.950)   
employmentempPartBene0.683   0.760    
(0.731)  (0.517)   
employmentunemployed-1.607 * -0.311    
(0.716)  (0.738)   
hhincomeb_30k<50k-0.239   -0.121    
(0.345)  (0.376)   
hhincomec_50k>100k0.666 **0.904 ***
(0.256)  (0.204)   
hhincomed_>100k0.062   0.219    
(0.365)  (0.326)   
marriage10.484   0.284    
(0.545)  (0.357)   
marriage2-0.427   -0.318    
(0.426)  (0.357)   
livesParent2lives_1_parent       -0.355    
       (0.566)   
livesParent2lives_other       3.339 ** 
       (1.033)   
mom.educ1       1.110    
       (0.678)   
mom.educ2       0.867    
       (0.698)   
mom.pubass1       -0.419    
       (1.041)   
closeMom2b_close       0.310    
       (0.244)   
parentSupervise2b_moderate       -2.084 ** 
       (0.671)   
parentSupervise2c_strict       -0.057    
       (0.655)   
religion1       1.033 ***
       (0.174)   
schoolconflict3       -0.116    
       (0.580)   
schoolconflict4       0.278    
       (0.667)   
schoolClubJoin1       -0.044    
       (0.929)   
scale(schoolSeg)       -0.019    
       (0.275)   
scale(pre)       -0.566 *  
       (0.249)   
*** p < 0.001; ** p < 0.01; * p < 0.05.