#install.packages("SASxport")
#install.packages("foreign")
#install.packages('Hmisc')
#library(haven)
#library(Hmisc)
#library(SASxport)
library(foreign)
#download datasets - AddHealth
wave1<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0001/da27021p1.dta")
wave4<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0012/da27021p12.dta")
#wave4birth<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0016/da27021p16.dta")
#addhealth.wghtW1<-read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0025/da27021p25.dta")
#wave3<-read.dta("/projects/add_health/data/11613344/ICPSR_27021/DS0003/da27021p3.dta")
w4meds <- read.xport("/projects/add_health/data/wave_4/biomarkers/w4meds.xpt")
w4glu <- read.xport("/projects/add_health/data/wave_4/biomarkers/glu_a1c.xpt")
#w4barorecepsensitv<-read.xport("/projects/add_health/data/data12_18/W4SupplementalFiles-Biomarker/Baroreceptor Sensitivity/brs.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")
#w1context <- read.dta("/projects/add_health/data/11613344/ICPSR_27024/DS0001/waveIcontext.dta")
#addhealth<-merge(wave1,wave3,by="aid")
#addhealth1<-merge(wave1,wave3,by="aid")
#addhealth1<-merge(addhealth1,wave4,by="aid")
#merge wave1, wave 4, biomarkers, medication.
addhealth1<-merge(wave4,wave1,by="aid")
#addhealth1<-merge(addhealth1a,wave3,by="aid")
addhealth2<-merge(w4meds,w4glu,by="AID")
addhealth3<-merge(w4crp,w4lipids,by="AID")
#Change uppercase AID to lowercase aid in order to merge by AID/aid.
addhealth2$aid<-addhealth2$AID
addhealth3$aid<-addhealth3$AID
addhealthBiomerge<-merge(addhealth2,addhealth3,by="aid")
#merge into one dataset
addhealthmerge1<-merge(addhealth1,addhealthBiomerge,by="aid")
addhealth<-merge(addhealthmerge1,weightW4,by="aid")
#addhealth<-merge(addhealth,wave4birth,by="aid")
library(car)
## Loading required package: carData
#Date of birth variables; Wave 1: H1GI1M (month), H1GI1Y (year); Wave 3: H3OD1M (month), H3OD1Y (year); Wave 4: H4OD1M (month), H4OD1Y (year)
#addhealth$birthmonthw1 <- factor(ifelse(addhealth$H1GI1M=="1",01,
# ifelse(addhealth$H1GI1M=="2", 02,
# ifelse(addhealth$H1GI1M=="3", 03,
# ifelse(addhealth$H1GI1M=="4", 04,
# ifelse(addhealth$H1GI1M=="5", 05,
# ifelse(addhealth$H1GI1M=="6", 06,
# ifelse(addhealth$H1GI1M=="7", 07,
# ifelse(addhealth$H1GI1M=="8", 08,
# ifelse(addhealth$H1GI1M=="9r", 09,
# ifelse(addhealth$H1GI1M=="10", 10,
# ifelse(addhealth$H1GI1M=="11", 11,
# ifelse(addhealth$H1GI1M=="12",12,
# ifelse(addhealth$H1GI1M=="96", 6,NA))))))))))))))
#table(addhealth$H1GI1M)
#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)))))))))))
#
#table(addhealth$birthyearw1)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
#addhealth$birthdatew1 <- as.yearmon(paste(addhealth$birthyearw1, addhealth$birthmonthw1), "%Y %m")
#table(addhealth$birthdatew1)
#interview date w1
#addhealth$iyearfix <- Recode(addhealth$iyear,recodes="'95'='1995'")
#addhealth$interviewdatew1 <- as.yearmon(paste(addhealth$iyearfix, addhealth$imonth), "%Y %m")
#table(addhealth$iyearfix)
#age - derived from date of birth subtracted from wave interview date
#addhealth$agew1<-(as.numeric(round(addhealth$interviewdatew1 - addhealth$birthdatew1)))
#table(addhealth$agew1)
#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)
##
## 24 25 26 27 28 29 30 31 32 33 34
## 1 84 972 1217 1903 1940 2170 1740 601 59 24
#sexual orientation from wave 3. Reference group is straight people
#addhealth$sexorient <- factor(ifelse(addhealth$H3SE13=="(1) 100% heterosexual (straight)",1,
# ifelse(addhealth$H3SE13=="(2) Mostly #heterosexual.somewhat attracted to people of own",2,
# ifelse(addhealth$H3SE13=="(3) Bisexual-attracted #to men and women equally",3,
# ifelse(addhealth$H3SE13=="(4) Mostly #homosexual.somewhat attracted to opposite sex",4,
# ifelse(addhealth$H3SE13=="(5) 100% #homosexual (gay)",5,NA))))))
#addhealth$sexorient<-Recode(addhealth$H3SE13, recodes="1:2='a_straight'; 3='b_bisexual';4:5='c_LGB';else=NA",as.factor=T)
addhealth$sexorient<-Recode(addhealth$H4SE31,recodes="1:2='a_straight';3='b_bisexual';4:5='c_LGB';else=NA",as.factor=T)
#addhealth$sexorient<-Recode(addhealth$H4SE31,recodes="1:2='a_straight';3='b_bisexual';4:5='c_LGB';else=NA")
#table(addhealth$sexorienta)
table(addhealth$sexorient)
##
## a_straight b_bisexual c_LGB
## 10060 221 285
#sex
addhealth$sex <- ifelse(addhealth$BIO_SEX4==2,'b_female','a_male')
#addhealth$sex_bio<-Recode(addhealth$sex, recodes="1='male'; 0='female'",as.factor=T)
#table(addhealth$sex_bio)
table(addhealth$sex)
##
## a_male b_female
## 3201 7510
#Education. Less then high school is reference group
addhealth$W4_educ<-Recode(addhealth$H4ED2,recodes="1:2='a_less_highschool';3:6='b_HS_grad';7:13='d_college_bach.plus';else=NA")
table(addhealth$W4_educ)
##
## a_less_highschool b_HS_grad d_college_bach.plus
## 711 6271 3727
#biomarkers
#systolic bp
addhealth$H4SBP<-Recode(addhealth$H4SBP,recodes="999=NA;997=NA;996=NA")
addhealth$systolicbp<-ifelse(addhealth$H4SBP>=140,1,0)
table(addhealth$systolicbp)
##
## 0 1
## 9190 1210
#diastolic bp
addhealth$H4DBP<-Recode(addhealth$H4DBP,recodes="999=NA;997=NA;996=NA")
addhealth$diastolicbp<-ifelse(addhealth$H4DBP>=90,1,0)
table(addhealth$diastolicbp)
##
## 0 1
## 8736 1664
#pulse rate
addhealth$H4PR<-Recode(addhealth$H4PR,recodes="999=NA;997=NA;996=NA")
addhealth$pulserate<-ifelse(addhealth$H4PR>=85,1,0)
#pulse pressure
addhealth$pulsepressure<-ifelse(addhealth$H4PP>=60,1,0)
table(addhealth$pulsepressure)
##
## 0 1
## 9953 682
#TOTAL CHOLESTEROL
addhealth$TC<-Recode(addhealth$TC,recodes="99=NA")
addhealth$totalchol<-ifelse(addhealth$TC>=8,1,0)
table(addhealth$totalchol)
##
## 0 1
## 6665 3131
#HDL, HIGH-DENSITY LIPOPROTEIN CHOLESTEROL
addhealth$HDL<-Recode(addhealth$HDL,recodes="99=NA")
addhealth$hidensity_lip<-ifelse(addhealth$HDL<=2,1,0)
table(addhealth$hidensity_lip)
##
## 0 1
## 7885 1779
#TRIGLYCERIDES
addhealth$TG<-Recode(addhealth$TG,recodes="99=NA")
addhealth$triglyceride<-ifelse(addhealth$TG>=8,1,0)
table(addhealth$triglyceride)
##
## 0 1
## 6351 3240
#HEMOGLOBIN
addhealth$HBA1C<-Recode(addhealth$HBA1C,recodes="99=NA")
addhealth$hemoglobin<-ifelse(addhealth$HBA1C>=6.4,1,0)
table(addhealth$hemoglobin)
##
## 0 1
## 9171 873
#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)
table(addhealth$BMI)
##
## 0 1
## 6194 4432
# high-sensitivity C-reactive protein test.
addhealth$CRPa<-Recode(addhealth$CRP,recodes="999=NA;998=NA")
addhealth$hsCRP<-ifelse(addhealth$CRPa>=3,1,0)
table(addhealth$hsCRP)
##
## 0 1
## 4615 5212
#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)
table(addhealth$EpsteinBar)
##
## 0 1
## 8021 1852
addhealth$score<-addhealth$EpsteinBar+addhealth$hsCRP+addhealth$BMI+addhealth$triglyceride+addhealth$hemoglobin+addhealth$hidensity_lip+addhealth$totalchol+addhealth$pulserate+addhealth$systolicbp+addhealth$pulsepressure
table(addhealth$score)
##
## 0 1 2 3 4 5 6 7 8 9
## 1039 1797 2019 1869 1351 704 261 100 13 8
#addhealth.pc<-prcomp(~EpsteinBar+hsCRP+BMI+waistcirm+hemoglobin+hidensity_lip+totalchol+pulserate+systolicbp+pulsepressure,data=addhealth, retx=T)
#screeplot(addhealth.pc, type = "l", main = "Scree Plot")
#abline(h=1)
#summary(addhealth.pc)
#addhealth.pc$rotation
#hist(addhealth.pc$x[,2])
#race variables
#nhwhite
addhealth$nhwhite <- ifelse(addhealth$H1GI6A==1,1,0)
addhealth$nhwhite<-Recode(addhealth$nhwhite, recodes="1='nhwhite'; 0='nonnhwhite'; 6:8=NA",as.factor=T)
#nhblack
addhealth$nhblack <- ifelse(addhealth$H1GI6B==1,1,0)
addhealth$nhblack<-Recode(addhealth$nhblack, recodes="1='nhblack'; 0='nonnhblack'; 6:8=NA",as.factor=T)
#Hispanic
addhealth$hispanic <- ifelse(addhealth$H1GI4==1,1,0)
addhealth$hispanic<-Recode(addhealth$hispanic, recodes="1='hispanic'; 0='nonhispanic'; 6:8=NA",as.factor=T)
#Asian
addhealth$asian <- ifelse(addhealth$H1GI6D==1,1,0)
addhealth$asian<-Recode(addhealth$asian, recodes="1='asian'; 0='nonasian'; 6:8=NA",as.factor=T)
#Native American
#addhealth$native_american <- ifelse(addhealth$H1GI6C==1,1,0)
#addhealth$native_american<-Recode(addhealth$native_american, recodes="1='native_american'; 0='nonnative_american'; 6:8=NA",as.factor=T)
#other
#addhealth$other <- ifelse(addhealth$H1GI6E==1&addhealth$H1GI6C==1,1,0)
#addhealth$other<-Recode(addhealth$other, recodes="1='other'; 0='nonother'; 6:8=NA",as.factor=T)
#combine race to one variable
addhealth$racethnic <- factor(ifelse(addhealth$nhwhite=="nhwhite","a-nhwhite",
ifelse(addhealth$nhblack=="nhblack", "b-nhblack",
ifelse(addhealth$hispanic=="hispanic", "c-hispanic",NA))))
# ifelse(addhealth$asian=="asian","d-asian",NA)))))
# ifelse(addhealth$native_american=="native_american","e-native_american",NA))))))
#ifelse(addhealth$other=="other","f-other",NA)))))))
table(addhealth$racethnic)
##
## a-nhwhite b-nhblack c-hispanic
## 7736 1789 635
#income variable
addhealth$incomeW4<-Recode(addhealth$H4EC1,recodes="1:5='a_<$25,000';6:9='b_$25,000<$75,000';10:11='c_$75,000<$150,000';12='d_>$150,000';else=NA")
#addhealth$treated_lessrespectW4<-Recode(addhealth$H4MH28,recodes="0:1='never/rarely';2:3='sometimes/often';else=NA")
#addhealth$suicideW4think<-Recode(addhealth$H4SE1, recodes="0='no';1='yes';else=NA",as.factor=T)
#addhealth$suicideW4attempt<-Recode(addhealth$H4SE2, recodes="0='no';1='once';2:3='2-4_times';4='5+times';else=NA",as.factor=T)
#Citizenship in Wave 4
addhealth$W4_UScitizenship<-ifelse(addhealth$H4OD4==1,"a_UScitizen_birth",
ifelse(addhealth$H4OD4==0&addhealth$H4OD5==1,"b_UScitizen_naturalized",
ifelse(addhealth$H4OD4==0&addhealth$H4OD5==0,"c_foreignBorn_nonUScitizen",NA)))
#addhealth$US_born<-ifelse(addhealth$H4OD4==1,"a_UScitizen_birth","b_foreign_born_non_UScitizen")
#born in the USA
addhealth$born_US_foreign<-ifelse(addhealth$H4OD4==1,"a_born.USA",
ifelse(addhealth$H4OD4==0,"b_born.foreign",NA))
#number of times married
addhealth$W4.married<-Recode(addhealth$H4TR1,recodes="0='a_never married';1='b_married_once';2:4='c_married_twice+';else=NA")
#Self reported general health
addhealth$health<-Recode(addhealth$H4GH1,recodes="1:3='a_Good.health';4:5='b_Bad.health';else=NA")
table(addhealth$health)
##
## a_Good.health b_Bad.health
## 8741 1970
#allostatic score by group
addhealth$score.grp<-Recode(addhealth$score,recodes="0:1='0';2:3='1';4:10='2';else=NA")
addhealth$score.grp<-as.factor(addhealth$score.grp)
table(addhealth$score.grp)
##
## 0 1 2
## 2836 3888 2437
quantile(addhealth$score,na.rm=T)
## 0% 25% 50% 75% 100%
## 0 1 2 4 9
str(addhealth$score.grp)
## Factor w/ 3 levels "0","1","2": 1 1 1 NA NA 1 1 2 1 3 ...
addhealth$raceXsexorient<-interaction(addhealth$racethnic,addhealth$sexorient)
addhealth$sexXsexorient<-interaction(addhealth$sex,addhealth$sexorient)
table(addhealth$sexXsexorient)
##
## a_male.a_straight b_female.a_straight a_male.b_bisexual b_female.b_bisexual
## 2979 7081 31 190
## a_male.c_LGB b_female.c_LGB
## 142 143
addhealth$seXorientXrace<-interaction(addhealth$sex,addhealth$sexorient,addhealth$racethnic)
table(addhealth$seXorientXrace)
##
## a_male.a_straight.a-nhwhite b_female.a_straight.a-nhwhite
## 2224 5089
## a_male.b_bisexual.a-nhwhite b_female.b_bisexual.a-nhwhite
## 22 149
## a_male.c_LGB.a-nhwhite b_female.c_LGB.a-nhwhite
## 75 90
## a_male.a_straight.b-nhblack b_female.a_straight.b-nhblack
## 446 1203
## a_male.b_bisexual.b-nhblack b_female.b_bisexual.b-nhblack
## 3 27
## a_male.c_LGB.b-nhblack b_female.c_LGB.b-nhblack
## 40 26
## a_male.a_straight.c-hispanic b_female.a_straight.c-hispanic
## 164 424
## a_male.b_bisexual.c-hispanic b_female.b_bisexual.c-hispanic
## 2 8
## a_male.c_LGB.c-hispanic b_female.c_LGB.c-hispanic
## 14 17
#sexorient
#find the most common value
#mcv.lgb<-factor(names(which.max(table(addhealth$sexorient))), levels=levels(addhealth$sexorient))
#mcv.lgb
#impute the cases
#addhealth$lgb.imp<-as.factor(ifelse(is.na(addhealth$sexorient)==T,mcv.lgb, addhealth$sexorient))
#levels(addhealth$lgb.imp)<-levels(addhealth$sexorient)
#prop.table(table(addhealth$sexorient))
#prop.table(table(addhealth$lgb.imp))
#table(addhealth$sexorient)
#mom born in USA
addhealth$mom.born.usa<-Recode(addhealth$PA3,recodes="0='b_mom.born.foreign';1='a_mom.born.usa';else=NA")
#mom race
addhealth$mom.latino<-Recode(addhealth$PA4,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.white<-Recode(addhealth$PA6_1,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.black<-Recode(addhealth$PA6_2,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.asian<-Recode(addhealth$PA6_4,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.nativeamer<-Recode(addhealth$PA6_3,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.other<-Recode(addhealth$PA6_5,recodes="0=0;1=1;6=6;else=NA")
addhealth$mom.race<-ifelse(addhealth$mom.latino==1,"c_latino",
#ifelse(addhealth$mom.asian==1,"d_asian",
ifelse(addhealth$mom.black==1,"b_black",
ifelse(addhealth$mom.white==1,"a_white",NA)))
table(addhealth$mom.race)
##
## a_white b_black c_latino
## 6439 1458 926
#mothers education
addhealth$mom_educ<-Recode(addhealth$PA12,recodes="1:2='a_less.HS';3:7='b_HSgrad/voca';8:9='c_col.grad';10='a_less.HS';else=NA")
#mothers employment
addhealth$mom_emp<-ifelse(addhealth$PA13==1,"b_yes.emp",
ifelse(addhealth$PA13==0,"a_no.emp",NA))
#mom disabled
addhealth$mom_disabiled<-Recode(addhealth$PA18,recodes="0='a_no';1='b_yes';else=NA")
#mom public assistance
addhealth$mom_welfare<-Recode(addhealth$PA21,recodes="0='a_no';1='b_yes';else=NA")
#neighborhood according to parent
addhealth$hood.litter<-Recode(addhealth$PA33,recodes="1=0;2:3=1;else=NA")
addhealth$hood.drugs<-Recode(addhealth$PA34,recodes="1=0;2:3=1;else=NA")
addhealth$hood.move<-Recode(addhealth$PA35,recodes="1=0;2:3=1;else=NA")
addhealth$good.neigh<-addhealth$hood.litter+addhealth$hood.drugs+addhealth$hood.move
addhealth$good.neigh<-as.factor(addhealth$good.neigh)
addhealth$good.neighood<-Recode(addhealth$good.neigh,recodes="0='0';1:2='1';3=3")
table(addhealth$good.neigh)
##
## 0 1 2 3
## 2678 2719 2086 1590
table(addhealth$good.neighood)
##
## 0 1 3
## 2678 4805 1590
#drinking
addhealth$mom.drinking<-Recode(addhealth$PA61,recodes="1='a_non.drinker';2:6='b_drinker';else=NA")
addhealth$mom.smoking<-Recode(addhealth$PA64,recodes="0='a_no';1='b_yes';else=NA")
#respondant smoking drinking
#addhealth$smoke<-Recode(addhealth$H4TO1,recodes="0=0;1=1;else=NA")
#addhealth$have.ever.smoke.regular<-Recode(addhealth$H4TO3,recodes="0=0;1=1;else=NA")
addhealth$smoke<-ifelse(addhealth$H4TO1==0,"a_non.smoker",
ifelse(addhealth$H4TO1==1&addhealth$H4TO3==1,"c_smoker",
ifelse(addhealth$H4TO1==1&addhealth$H4TO3==0,"b_tried.smoking",NA)))
#addhealth$H4TO5[addhealth$H4TO5==96] <- NA
#addhealth$H4TO5[addhealth$H4TO5==98] <- NA
#addhealth$H4TO6[addhealth$H4TO6==996] <- NA
#addhealth$H4TO6[addhealth$H4TO6==997] <- NA
#addhealth$H4TO6[addhealth$H4TO6==998] <- NA
#addhealth$number.smokes.daily<-addhealth$H4TO6
#addhealth$smoke.past.30days<-Recode(addhealth$H4TO5,recodes="0=0;1:30=1;else=NA")
#addhealth$smoke<-ifelse(addhealth$have.smoke==0&addhealth$smoke.daily.currently==0,"a_non.smoker",
# ifelse(addhealth$have.smoke==1&addhealth$smoke.daily.currently==0,"a_non.smoker",
# ifelse(addhealth$have.smoke==1&addhealth$smoke.daily.currently==1,"b_smoker",
# ifelse(addhealth$have.smoke==0&addhealth$smoke.daily.currently==1,"b_smoker",NA))))
#addhealth$smoke<-Recode(addhealth$H4TO5,recodes="0='a_non.smoker';1:30='b_smoker';else=NA")
#addhealth$smoke2<-Recode(addhealth$H4TO5,recodes="0=0;1:30=1;else=NA")
#addhealth$smoke3<-Recode(addhealth$H4TO6,recodes="1:24=0;25:100=1;else=NA")
#addhealth$smoke<-ifelse(addhealth$smoke1==1&addhealth$smoke2==1&addhealth$smoke3==1,"d_heavy.smoker",
# ifelse(addhealth$smoke1==1&addhealth$smoke2==1&addhealth$smoke3==0,"c_moderate.smoker",
# ifelse(addhealth$smoke1==0&addhealth$smoke2==0,"a_non.smoker",
# ifelse(addhealth$smoke1==1&addhealth$smoke2==0,"b_previous.smoker",NA))))
#addhealth$smoke<-addhealth$have.ever.smoke+addhealth$have.ever.smoke.regular+addhealth$smoke.past.30days
#drinking respondant
addhealth$drink.yes.no<-ifelse(addhealth$H4TO33==0,0,
ifelse(addhealth$H4TO33==1,1,NA))
addhealth$drink.yes.12months<-Recode(addhealth$H4TO35,recodes="0:3=1;4:6=2;else=NA")
#addhealth$drink.current<-ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12month==1,1,0)
#addhealth$drink.binge<-Recode(addhealth$H4TO36,recodes="1:4=0;5:18=1;else=NA")
#addhealth$drink<-ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12months==1&addhealth$drink.binge==1,"d_binge.drinker",
# ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12months==1&addhealth$drink.binge==0,"c_current.drinker",
# ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12months==0,"b_previous.drinker",
# ifelse(addhealth$drink.yes.no==0,"a_non.drinker",NA))))
addhealth$drink<-ifelse(addhealth$drink.yes.no==0,"a_non.drinker",
ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12months==1,"b_lite.drinker",
ifelse(addhealth$drink.yes.no==1&addhealth$drink.yes.12months==2,"c_drinker",NA)))
table(addhealth$smoke)
##
## a_non.smoker b_tried.smoking c_smoker
## 3712 2109 4864
table(addhealth$drink)
##
## a_non.drinker b_lite.drinker c_drinker
## 2196 5699 2761
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following object is masked from 'package:car':
##
## logit
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#select variables to use
addhealth.pro<-addhealth%>%
select(psuscid,GSWGT4_2,score.grp,sex,racethnic,mom_educ,mom_emp,mom_disabiled,mom_welfare,good.neigh,mom.drinking,mom.smoking,born_US_foreign,sexorient,W4_educ,raceXsexorient,agew4,W4.married,incomeW4,W4_UScitizenship,health,mom.race,good.neigh,smoke,drink,sexXsexorient,seXorientXrace)
#filter complete cases
#library(tidyverse)
addhealth.pro<-addhealth.pro%>%
filter(complete.cases(.))
#Survey design
library(survey)
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
##
## Attaching package: 'survey'
## The following object is masked from 'package:VGAM':
##
## calibrate
## The following object is masked from 'package:graphics':
##
## dotchart
#addhealth %>% drop_na()
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~psuscid,
weights=~GSWGT4_2,
nest=T,
data=addhealth.pro[!is.na(addhealth.pro$GSWGT4_2),])
fit<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient, family = cumulative(parallel = T, reverse = T),
## data = addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -3.865 -0.4924 0.3024 0.6574 2.307
## logitlink(P[Y>=3]) -2.250 -0.6598 -0.2998 0.5747 3.718
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.73416 0.56310 -1.304 0.19231
## (Intercept):2 -2.55267 0.56516 -4.517 6.28e-06 ***
## agew4 0.06128 0.01923 3.187 0.00144 **
## sexb_female -0.23151 0.07420 -3.120 0.00181 **
## racethnicb-nhblack 0.30565 0.11407 2.679 0.00737 **
## racethnicc-hispanic 0.93497 0.23103 4.047 5.19e-05 ***
## sexorientb_bisexual 0.12734 0.23219 0.548 0.58340
## sexorientc_LGB 0.26262 0.23495 1.118 0.26367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 6305.159 on 5864 degrees of freedom
##
## Log-likelihood: -3152.579 on 5864 degrees of freedom
##
## Number of Fisher scoring iterations: 4
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexb_female racethnicb-nhblack racethnicc-hispanic
## 1.0631977 0.7933364 1.3575049 2.5471281
## sexorientb_bisexual sexorientc_LGB
## 1.1358003 1.3003321
fit2<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit2)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient + W4_educ + born_US_foreign, family = cumulative(parallel = T,
## reverse = T), data = addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2,
## na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -4.026 -0.4194 0.2862 0.6501 2.152
## logitlink(P[Y>=3]) -2.365 -0.6560 -0.2908 0.5589 3.497
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.62114 0.57177 -1.086 0.27732
## (Intercept):2 -2.47326 0.57371 -4.311 1.63e-05 ***
## agew4 0.06238 0.01931 3.229 0.00124 **
## sexb_female -0.18042 0.07480 -2.412 0.01586 *
## racethnicb-nhblack 0.35490 0.11470 3.094 0.00197 **
## racethnicc-hispanic 0.90951 0.23344 3.896 9.78e-05 ***
## sexorientb_bisexual 0.18356 0.23384 0.785 0.43245
## sexorientc_LGB 0.22625 0.23835 0.949 0.34250
## W4_educb_HS_grad -0.03718 0.11478 -0.324 0.74601
## W4_educd_college_bach.plus -0.72738 0.13175 -5.521 3.37e-08 ***
## born_US_foreignb_born.foreign 0.75460 0.38405 1.965 0.04943 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 6237.493 on 5861 degrees of freedom
##
## Log-likelihood: -3118.746 on 5861 degrees of freedom
##
## Number of Fisher scoring iterations: 4
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexb_female
## 1.0643621 0.8349179
## racethnicb-nhblack racethnicc-hispanic
## 1.4260372 2.4831078
## sexorientb_bisexual sexorientc_LGB
## 1.2014906 1.2538944
## W4_educb_HS_grad W4_educd_college_bach.plus
## 0.9635062 0.4831754
## born_US_foreignb_born.foreign
## 2.1267586
fit3<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign+W4.married+incomeW4,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit3)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient + W4_educ + born_US_foreign + W4.married + incomeW4,
## family = cumulative(parallel = T, reverse = T), data = addhealth.pro,
## weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -4.132 -0.4058 0.2885 0.6389 2.028
## logitlink(P[Y>=3]) -2.362 -0.6501 -0.2821 0.5809 3.520
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.460633 0.583533 -0.789 0.42989
## (Intercept):2 -2.333871 0.585245 -3.988 6.67e-05 ***
## agew4 0.061873 0.019957 3.100 0.00193 **
## sexb_female -0.215906 0.075562 -2.857 0.00427 **
## racethnicb-nhblack 0.279307 0.116286 2.402 0.01631 *
## racethnicc-hispanic 0.965688 0.234309 4.121 3.77e-05 ***
## sexorientb_bisexual 0.002472 0.236343 0.010 0.99165
## sexorientc_LGB 0.276237 0.240745 1.147 0.25121
## W4_educb_HS_grad 0.127289 0.120538 1.056 0.29096
## W4_educd_college_bach.plus -0.430115 0.144611 -2.974 0.00294 **
## born_US_foreignb_born.foreign 0.762305 0.383063 1.990 0.04659 *
## W4.marriedb_married_once -0.039862 0.076466 -0.521 0.60215
## W4.marriedc_married_twice+ 0.119815 0.142726 0.839 0.40120
## incomeW4b_$25,000<$75,000 -0.227703 0.091723 -2.483 0.01305 *
## incomeW4c_$75,000<$150,000 -0.585706 0.118772 -4.931 8.17e-07 ***
## incomeW4d_>$150,000 -1.067039 0.218225 -4.890 1.01e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 6193.869 on 5856 degrees of freedom
##
## Log-likelihood: -3096.935 on 5856 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexb_female
## 1.0638275 0.8058111
## racethnicb-nhblack racethnicc-hispanic
## 1.3222133 2.6265952
## sexorientb_bisexual sexorientc_LGB
## 1.0024752 1.3181598
## W4_educb_HS_grad W4_educd_college_bach.plus
## 1.1357452 0.6504345
## born_US_foreignb_born.foreign W4.marriedb_married_once
## 2.1432117 0.9609216
## W4.marriedc_married_twice+ incomeW4b_$25,000<$75,000
## 1.1272888 0.7963606
## incomeW4c_$75,000<$150,000 incomeW4d_>$150,000
## 0.5567124 0.3440257
fit4<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit4)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient + W4_educ + born_US_foreign + W4.married + incomeW4 +
## drink + smoke + health, family = cumulative(parallel = T,
## reverse = T), data = addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2,
## na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -4.805 -0.3989 0.2638 0.5876 2.362
## logitlink(P[Y>=3]) -2.595 -0.6350 -0.2684 0.5999 4.065
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.50330 0.59545 -0.845 0.397973
## (Intercept):2 -2.45910 0.59735 -4.117 3.84e-05 ***
## agew4 0.06899 0.02020 3.415 0.000639 ***
## sexb_female -0.27441 0.07763 -3.535 0.000408 ***
## racethnicb-nhblack -0.02439 0.12061 -0.202 0.839727
## racethnicc-hispanic 0.84955 0.23914 3.553 0.000382 ***
## sexorientb_bisexual -0.13363 0.24099 -0.555 0.579226
## sexorientc_LGB 0.32213 0.24317 1.325 0.185271
## W4_educb_HS_grad 0.20118 0.12402 1.622 0.104762
## W4_educd_college_bach.plus -0.37471 0.14999 -2.498 0.012481 *
## born_US_foreignb_born.foreign 0.70860 0.39069 1.814 0.069721 .
## W4.marriedb_married_once -0.09871 0.07812 -1.263 0.206412
## W4.marriedc_married_twice+ -0.01485 0.14539 -0.102 0.918639
## incomeW4b_$25,000<$75,000 -0.09147 0.09419 -0.971 0.331514
## incomeW4c_$75,000<$150,000 -0.34927 0.12275 -2.845 0.004436 **
## incomeW4d_>$150,000 -0.94538 0.22620 -4.179 2.92e-05 ***
## drinkb_lite.drinker 0.07930 0.09655 0.821 0.411436
## drinkc_drinker -0.24646 0.11748 -2.098 0.035910 *
## smokeb_tried.smoking -0.18696 0.11640 -1.606 0.108231
## smokec_smoker -0.58947 0.08882 -6.637 3.20e-11 ***
## healthb_Bad.health 0.87808 0.09103 9.646 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 6030.878 on 5851 degrees of freedom
##
## Log-likelihood: -3015.439 on 5851 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexb_female
## 1.0714219 0.7600213
## racethnicb-nhblack racethnicc-hispanic
## 0.9759025 2.3385911
## sexorientb_bisexual sexorientc_LGB
## 0.8749114 1.3800634
## W4_educb_HS_grad W4_educd_college_bach.plus
## 1.2228444 0.6874897
## born_US_foreignb_born.foreign W4.marriedb_married_once
## 2.0311411 0.9060069
## W4.marriedc_married_twice+ incomeW4b_$25,000<$75,000
## 0.9852587 0.9125900
## incomeW4c_$75,000<$150,000 incomeW4d_>$150,000
## 0.7052028 0.3885331
## drinkb_lite.drinker drinkc_drinker
## 1.0825314 0.7815590
## smokeb_tried.smoking smokec_smoker
## 0.8294732 0.5546192
## healthb_Bad.health
## 2.4062748
fit5<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit5)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient + W4_educ + born_US_foreign + W4.married + incomeW4 +
## drink + smoke + health + mom.race + mom_educ + mom_emp +
## good.neigh, family = cumulative(parallel = T, reverse = T),
## data = addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -6.073 -0.3967 0.2646 0.5824 2.310
## logitlink(P[Y>=3]) -2.305 -0.6420 -0.2660 0.5634 3.989
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.666111 0.609316 -1.093 0.274302
## (Intercept):2 -2.633491 0.611363 -4.308 1.65e-05 ***
## agew4 0.073687 0.020383 3.615 0.000300 ***
## sexb_female -0.265287 0.078038 -3.399 0.000675 ***
## racethnicb-nhblack -1.490818 0.521751 -2.857 0.004272 **
## racethnicc-hispanic 0.574128 0.280026 2.050 0.040338 *
## sexorientb_bisexual -0.106475 0.242943 -0.438 0.661188
## sexorientc_LGB 0.376658 0.247101 1.524 0.127432
## W4_educb_HS_grad 0.201662 0.126265 1.597 0.110237
## W4_educd_college_bach.plus -0.354121 0.154419 -2.293 0.021834 *
## born_US_foreignb_born.foreign 0.417074 0.418595 0.996 0.319073
## W4.marriedb_married_once -0.108242 0.079233 -1.366 0.171899
## W4.marriedc_married_twice+ -0.005761 0.148033 -0.039 0.968955
## incomeW4b_$25,000<$75,000 -0.113494 0.095724 -1.186 0.235767
## incomeW4c_$75,000<$150,000 -0.346356 0.124588 -2.780 0.005436 **
## incomeW4d_>$150,000 -0.974377 0.229762 -4.241 2.23e-05 ***
## drinkb_lite.drinker 0.124044 0.097626 1.271 0.203868
## drinkc_drinker -0.175298 0.119931 -1.462 0.143834
## smokeb_tried.smoking -0.186026 0.117362 -1.585 0.112954
## smokec_smoker -0.593026 0.089512 -6.625 3.47e-11 ***
## healthb_Bad.health 0.855052 0.091521 9.343 < 2e-16 ***
## mom.raceb_black 1.524265 0.526399 2.896 0.003784 **
## mom.racec_latino 0.335513 0.183811 1.825 0.067954 .
## mom_educb_HSgrad/voca 0.038289 0.104668 0.366 0.714504
## mom_educc_col.grad -0.136362 0.141019 -0.967 0.333556
## mom_empb_yes.emp -0.106211 0.084155 -1.262 0.206914
## good.neigh1 -0.057987 0.099381 -0.583 0.559571
## good.neigh2 0.107635 0.104120 1.034 0.301251
## good.neigh3 0.187959 0.108973 1.725 0.084559 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 6008.573 on 5843 degrees of freedom
##
## Log-likelihood: -3004.286 on 5843 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexb_female
## 1.0764703 0.7669857
## racethnicb-nhblack racethnicc-hispanic
## 0.2251884 1.7755818
## sexorientb_bisexual sexorientc_LGB
## 0.8989972 1.4574054
## W4_educb_HS_grad W4_educd_college_bach.plus
## 1.2234339 0.7017904
## born_US_foreignb_born.foreign W4.marriedb_married_once
## 1.5175142 0.8974103
## W4.marriedc_married_twice+ incomeW4b_$25,000<$75,000
## 0.9942553 0.8927094
## incomeW4c_$75,000<$150,000 incomeW4d_>$150,000
## 0.7072604 0.3774273
## drinkb_lite.drinker drinkc_drinker
## 1.1320660 0.8392071
## smokeb_tried.smoking smokec_smoker
## 0.8302524 0.5526524
## healthb_Bad.health mom.raceb_black
## 2.3514969 4.5917675
## mom.racec_latino mom_educb_HSgrad/voca
## 1.3986573 1.0390314
## mom_educc_col.grad mom_empb_yes.emp
## 0.8725271 0.8992345
## good.neigh1 good.neigh2
## 0.9436626 1.1136408
## good.neigh3
## 1.2067843
fit6<-vglm(as.ordered(score.grp)~agew4+sex+raceXsexorient+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit6)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + raceXsexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh,
## family = cumulative(parallel = T, reverse = T), data = addhealth.pro,
## weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -4.696 -0.3838 0.2642 0.5838 2.289
## logitlink(P[Y>=3]) -2.336 -0.6406 -0.2638 0.5662 3.986
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.581440 0.610093 -0.953 0.340572
## (Intercept):2 -2.553098 0.612041 -4.171 3.03e-05 ***
## agew4 0.070399 0.020414 3.449 0.000564 ***
## sexb_female -0.262079 0.078099 -3.356 0.000792 ***
## raceXsexorientb-nhblack.a_straight -0.943216 0.563818 -1.673 0.094345 .
## raceXsexorientc-hispanic.a_straight 0.595857 0.280294 2.126 0.033518 *
## raceXsexorienta-nhwhite.b_bisexual -0.134224 0.250922 -0.535 0.592703
## raceXsexorientb-nhblack.b_bisexual -0.876550 1.035103 -0.847 0.397093
## raceXsexorientc-hispanic.b_bisexual 11.494274 448.093873 NA NA
## raceXsexorienta-nhwhite.c_LGB 0.644408 0.263756 2.443 0.014558 *
## raceXsexorientb-nhblack.c_LGB -3.072100 1.006043 -3.054 0.002261 **
## raceXsexorientc-hispanic.c_LGB -0.087400 8.174216 -0.011 0.991469
## W4_educb_HS_grad 0.208622 0.126278 1.652 0.098517 .
## W4_educd_college_bach.plus -0.346128 0.154518 -2.240 0.025087 *
## born_US_foreignb_born.foreign 0.386893 0.418199 0.925 0.354893
## W4.marriedb_married_once -0.096964 0.079334 -1.222 0.221623
## W4.marriedc_married_twice+ 0.004092 0.148061 0.028 0.977949
## incomeW4b_$25,000<$75,000 -0.131616 0.095933 -1.372 0.170074
## incomeW4c_$75,000<$150,000 -0.364511 0.124745 -2.922 0.003477 **
## incomeW4d_>$150,000 -0.996977 0.230045 -4.334 1.47e-05 ***
## drinkb_lite.drinker 0.130624 0.097716 1.337 0.181298
## drinkc_drinker -0.170452 0.120035 -1.420 0.155602
## smokeb_tried.smoking -0.190656 0.117412 -1.624 0.104414
## smokec_smoker -0.584809 0.089588 -6.528 6.68e-11 ***
## healthb_Bad.health 0.851158 0.091802 9.272 < 2e-16 ***
## mom.raceb_black 1.018581 0.562151 1.812 0.069996 .
## mom.racec_latino 0.325271 0.183877 1.769 0.076901 .
## mom_educb_HSgrad/voca 0.041120 0.104772 0.392 0.694708
## mom_educc_col.grad -0.131094 0.141107 -0.929 0.352869
## mom_empb_yes.emp -0.115542 0.084284 -1.371 0.170415
## good.neigh1 -0.055951 0.099438 -0.563 0.573657
## good.neigh2 0.117100 0.104220 1.124 0.261190
## good.neigh3 0.191686 0.108999 1.759 0.078644 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 5997.802 on 5839 degrees of freedom
##
## Log-likelihood: -2998.901 on 5839 degrees of freedom
##
## Number of Fisher scoring iterations: 10
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## 'raceXsexorientc-hispanic.b_bisexual'
##
##
## Exponentiated coefficients:
## agew4 sexb_female
## 1.072936e+00 7.694503e-01
## raceXsexorientb-nhblack.a_straight raceXsexorientc-hispanic.a_straight
## 3.893735e-01 1.814586e+00
## raceXsexorienta-nhwhite.b_bisexual raceXsexorientb-nhblack.b_bisexual
## 8.743942e-01 4.162164e-01
## raceXsexorientc-hispanic.b_bisexual raceXsexorienta-nhwhite.c_LGB
## 9.815214e+04 1.904859e+00
## raceXsexorientb-nhblack.c_LGB raceXsexorientc-hispanic.c_LGB
## 4.632379e-02 9.163109e-01
## W4_educb_HS_grad W4_educd_college_bach.plus
## 1.231980e+00 7.074219e-01
## born_US_foreignb_born.foreign W4.marriedb_married_once
## 1.472399e+00 9.075888e-01
## W4.marriedc_married_twice+ incomeW4b_$25,000<$75,000
## 1.004101e+00 8.766772e-01
## incomeW4c_$75,000<$150,000 incomeW4d_>$150,000
## 6.945362e-01 3.689932e-01
## drinkb_lite.drinker drinkc_drinker
## 1.139539e+00 8.432834e-01
## smokeb_tried.smoking smokec_smoker
## 8.264166e-01 5.572125e-01
## healthb_Bad.health mom.raceb_black
## 2.342357e+00 2.769261e+00
## mom.racec_latino mom_educb_HSgrad/voca
## 1.384406e+00 1.041977e+00
## mom_educc_col.grad mom_empb_yes.emp
## 8.771356e-01 8.908829e-01
## good.neigh1 good.neigh2
## 9.455854e-01 1.124232e+00
## good.neigh3
## 1.211290e+00
fit7<-vglm(as.ordered(score.grp)~agew4+sexXsexorient+racethnic+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit7)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sexXsexorient +
## racethnic + W4_educ + born_US_foreign + W4.married + incomeW4 +
## drink + smoke + health + mom.race + mom_educ + mom_emp +
## good.neigh, family = cumulative(parallel = T, reverse = T),
## data = addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -5.634 -0.3941 0.2605 0.5776 3.416
## logitlink(P[Y>=3]) -2.262 -0.6322 -0.2643 0.5668 4.060
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.85425 0.61181 -1.396 0.162635
## (Intercept):2 -2.83740 0.61410 -4.620 3.83e-06 ***
## agew4 0.07817 0.02045 3.822 0.000132 ***
## sexXsexorientb_female.a_straight -0.26929 0.07948 -3.388 0.000703 ***
## sexXsexorienta_male.b_bisexual -2.98048 0.74121 -4.021 5.79e-05 ***
## sexXsexorientb_female.b_bisexual 0.03270 0.26940 0.121 0.903396
## sexXsexorienta_male.c_LGB 1.29225 0.38851 3.326 0.000881 ***
## sexXsexorientb_female.c_LGB -0.51412 0.33500 -1.535 0.124858
## racethnicb-nhblack -1.33172 0.52401 -2.541 0.011041 *
## racethnicc-hispanic 0.60921 0.28042 2.172 0.029820 *
## W4_educb_HS_grad 0.23618 0.12664 1.865 0.062177 .
## W4_educd_college_bach.plus -0.30908 0.15485 -1.996 0.045935 *
## born_US_foreignb_born.foreign 0.25632 0.42080 0.609 0.542435
## W4.marriedb_married_once -0.10783 0.07945 -1.357 0.174724
## W4.marriedc_married_twice+ -0.02830 0.14833 -0.191 0.848679
## incomeW4b_$25,000<$75,000 -0.11250 0.09598 -1.172 0.241156
## incomeW4c_$75,000<$150,000 -0.36306 0.12486 -2.908 0.003640 **
## incomeW4d_>$150,000 -1.02635 0.23130 -4.437 9.11e-06 ***
## drinkb_lite.drinker 0.11827 0.09783 1.209 0.226678
## drinkc_drinker -0.19327 0.12028 -1.607 0.108079
## smokeb_tried.smoking -0.19969 0.11764 -1.697 0.089612 .
## smokec_smoker -0.58815 0.08968 -6.558 5.44e-11 ***
## healthb_Bad.health 0.89257 0.09204 9.698 < 2e-16 ***
## mom.raceb_black 1.35718 0.52831 2.569 0.010202 *
## mom.racec_latino 0.31714 0.18405 1.723 0.084860 .
## mom_educb_HSgrad/voca 0.05094 0.10505 0.485 0.627722
## mom_educc_col.grad -0.13855 0.14140 -0.980 0.327179
## mom_empb_yes.emp -0.07859 0.08449 -0.930 0.352283
## good.neigh1 -0.01767 0.09978 -0.177 0.859472
## good.neigh2 0.11979 0.10431 1.148 0.250792
## good.neigh3 0.18517 0.10919 1.696 0.089907 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 5977.923 on 5841 degrees of freedom
##
## Log-likelihood: -2988.962 on 5841 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## agew4 sexXsexorientb_female.a_straight
## 1.08130712 0.76392411
## sexXsexorienta_male.b_bisexual sexXsexorientb_female.b_bisexual
## 0.05076827 1.03323844
## sexXsexorienta_male.c_LGB sexXsexorientb_female.c_LGB
## 3.64097071 0.59802424
## racethnicb-nhblack racethnicc-hispanic
## 0.26402229 1.83897989
## W4_educb_HS_grad W4_educd_college_bach.plus
## 1.26640855 0.73412344
## born_US_foreignb_born.foreign W4.marriedb_married_once
## 1.29216958 0.89778105
## W4.marriedc_married_twice+ incomeW4b_$25,000<$75,000
## 0.97209541 0.89360003
## incomeW4c_$75,000<$150,000 incomeW4d_>$150,000
## 0.69554682 0.35831320
## drinkb_lite.drinker drinkc_drinker
## 1.12554964 0.82425683
## smokeb_tried.smoking smokec_smoker
## 0.81898152 0.55535190
## healthb_Bad.health mom.raceb_black
## 2.44138871 3.88523670
## mom.racec_latino mom_educb_HSgrad/voca
## 1.37319423 1.05226269
## mom_educc_col.grad mom_empb_yes.emp
## 0.87062331 0.92442280
## good.neigh1 good.neigh2
## 0.98248912 1.12726483
## good.neigh3
## 1.20342360
fit8<-vglm(as.ordered(score.grp)~agew4+seXorientXrace+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T))
summary(fit8)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + seXorientXrace +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh,
## family = cumulative(parallel = T, reverse = T), data = addhealth.pro,
## weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -4.557 -0.3736 0.2587 0.5721 3.544
## logitlink(P[Y>=3]) -2.186 -0.6198 -0.2638 0.5624 4.024
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 -0.88640 0.61481 -1.442
## (Intercept):2 -2.87565 0.61714 -4.660
## agew4 0.07753 0.02051 3.779
## seXorientXraceb_female.a_straight.a-nhwhite -0.21801 0.08435 -2.585
## seXorientXracea_male.b_bisexual.a-nhwhite -2.99925 0.75425 -3.976
## seXorientXraceb_female.b_bisexual.a-nhwhite 0.08675 0.28032 0.309
## seXorientXracea_male.c_LGB.a-nhwhite 1.46023 0.40189 3.633
## seXorientXraceb_female.c_LGB.a-nhwhite -0.17988 0.36430 -0.494
## seXorientXracea_male.a_straight.b-nhblack -0.40523 0.60607 -0.669
## seXorientXraceb_female.a_straight.b-nhblack -1.24349 0.57672 -2.156
## seXorientXraceb_female.b_bisexual.b-nhblack -1.06735 1.03768 -1.029
## seXorientXracea_male.c_LGB.b-nhblack -2.21408 1.92294 -1.151
## seXorientXraceb_female.c_LGB.b-nhblack -3.53934 1.18116 NA
## seXorientXracea_male.a_straight.c-hispanic 0.63069 0.37016 1.704
## seXorientXraceb_female.a_straight.c-hispanic 0.41153 0.37114 1.109
## seXorientXracea_male.b_bisexual.c-hispanic 11.54097 447.96372 NA
## seXorientXraceb_female.c_LGB.c-hispanic -0.15299 8.18135 -0.019
## W4_educb_HS_grad 0.25760 0.12716 2.026
## W4_educd_college_bach.plus -0.29033 0.15571 -1.865
## born_US_foreignb_born.foreign 0.24893 0.42230 0.589
## W4.marriedb_married_once -0.09352 0.07954 -1.176
## W4.marriedc_married_twice+ -0.02109 0.14845 -0.142
## incomeW4b_$25,000<$75,000 -0.14129 0.09639 -1.466
## incomeW4c_$75,000<$150,000 -0.38758 0.12509 -3.098
## incomeW4d_>$150,000 -1.05021 0.23216 -4.524
## drinkb_lite.drinker 0.14384 0.09839 1.462
## drinkc_drinker -0.17134 0.12057 -1.421
## smokeb_tried.smoking -0.19966 0.11776 -1.695
## smokec_smoker -0.58544 0.08982 -6.518
## healthb_Bad.health 0.89497 0.09250 9.675
## mom.raceb_black 0.94755 0.56679 1.672
## mom.racec_latino 0.31279 0.18418 1.698
## mom_educb_HSgrad/voca 0.04609 0.10549 0.437
## mom_educc_col.grad -0.14832 0.14180 -1.046
## mom_empb_yes.emp -0.09564 0.08476 -1.128
## good.neigh1 -0.02059 0.09997 -0.206
## good.neigh2 0.13456 0.10446 1.288
## good.neigh3 0.19791 0.10936 1.810
## Pr(>|z|)
## (Intercept):1 0.149373
## (Intercept):2 3.17e-06 ***
## agew4 0.000157 ***
## seXorientXraceb_female.a_straight.a-nhwhite 0.009752 **
## seXorientXracea_male.b_bisexual.a-nhwhite 6.99e-05 ***
## seXorientXraceb_female.b_bisexual.a-nhwhite 0.756961
## seXorientXracea_male.c_LGB.a-nhwhite 0.000280 ***
## seXorientXraceb_female.c_LGB.a-nhwhite 0.621481
## seXorientXracea_male.a_straight.b-nhblack 0.503730
## seXorientXraceb_female.a_straight.b-nhblack 0.031072 *
## seXorientXraceb_female.b_bisexual.b-nhblack 0.303670
## seXorientXracea_male.c_LGB.b-nhblack 0.249568
## seXorientXraceb_female.c_LGB.b-nhblack NA
## seXorientXracea_male.a_straight.c-hispanic 0.088411 .
## seXorientXraceb_female.a_straight.c-hispanic 0.267504
## seXorientXracea_male.b_bisexual.c-hispanic NA
## seXorientXraceb_female.c_LGB.c-hispanic 0.985081
## W4_educb_HS_grad 0.042781 *
## W4_educd_college_bach.plus 0.062244 .
## born_US_foreignb_born.foreign 0.555552
## W4.marriedb_married_once 0.239705
## W4.marriedc_married_twice+ 0.887024
## incomeW4b_$25,000<$75,000 0.142695
## incomeW4c_$75,000<$150,000 0.001945 **
## incomeW4d_>$150,000 6.08e-06 ***
## drinkb_lite.drinker 0.143755
## drinkc_drinker 0.155277
## smokeb_tried.smoking 0.090000 .
## smokec_smoker 7.12e-11 ***
## healthb_Bad.health < 2e-16 ***
## mom.raceb_black 0.094570 .
## mom.racec_latino 0.089466 .
## mom_educb_HSgrad/voca 0.662134
## mom_educc_col.grad 0.295582
## mom_empb_yes.emp 0.259213
## good.neigh1 0.836851
## good.neigh2 0.197687
## good.neigh3 0.070350 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 5964.185 on 5834 degrees of freedom
##
## Log-likelihood: -2982.093 on 5834 degrees of freedom
##
## Number of Fisher scoring iterations: 10
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## 'seXorientXraceb_female.c_LGB.b-nhblack', 'seXorientXracea_male.b_bisexual.c-hispanic'
##
##
## Exponentiated coefficients:
## agew4
## 1.080613e+00
## seXorientXraceb_female.a_straight.a-nhwhite
## 8.041210e-01
## seXorientXracea_male.b_bisexual.a-nhwhite
## 4.982441e-02
## seXorientXraceb_female.b_bisexual.a-nhwhite
## 1.090626e+00
## seXorientXracea_male.c_LGB.a-nhwhite
## 4.306958e+00
## seXorientXraceb_female.c_LGB.a-nhwhite
## 8.353743e-01
## seXorientXracea_male.a_straight.b-nhblack
## 6.668203e-01
## seXorientXraceb_female.a_straight.b-nhblack
## 2.883771e-01
## seXorientXraceb_female.b_bisexual.b-nhblack
## 3.439170e-01
## seXorientXracea_male.c_LGB.b-nhblack
## 1.092544e-01
## seXorientXraceb_female.c_LGB.b-nhblack
## 2.903243e-02
## seXorientXracea_male.a_straight.c-hispanic
## 1.878909e+00
## seXorientXraceb_female.a_straight.c-hispanic
## 1.509119e+00
## seXorientXracea_male.b_bisexual.c-hispanic
## 1.028446e+05
## seXorientXraceb_female.c_LGB.c-hispanic
## 8.581402e-01
## W4_educb_HS_grad
## 1.293818e+00
## W4_educd_college_bach.plus
## 7.480175e-01
## born_US_foreignb_born.foreign
## 1.282650e+00
## W4.marriedb_married_once
## 9.107176e-01
## W4.marriedc_married_twice+
## 9.791300e-01
## incomeW4b_$25,000<$75,000
## 8.682354e-01
## incomeW4c_$75,000<$150,000
## 6.786954e-01
## incomeW4d_>$150,000
## 3.498650e-01
## drinkb_lite.drinker
## 1.154695e+00
## drinkc_drinker
## 8.425311e-01
## smokeb_tried.smoking
## 8.190126e-01
## smokec_smoker
## 5.568609e-01
## healthb_Bad.health
## 2.447265e+00
## mom.raceb_black
## 2.579370e+00
## mom.racec_latino
## 1.367229e+00
## mom_educb_HSgrad/voca
## 1.047173e+00
## mom_educc_col.grad
## 8.621581e-01
## mom_empb_yes.emp
## 9.087951e-01
## good.neigh1
## 9.796244e-01
## good.neigh2
## 1.144031e+00
## good.neigh3
## 1.218849e+00
AIC(fit)
## [1] 6321.159
AIC(fit2)
## [1] 6259.493
AIC(fit3)
## [1] 6225.869
AIC(fit4)
## [1] 6072.878
AIC(fit5)
## [1] 6066.573
#race x sexorient
AIC(fit6)
## [1] 6063.802
#sex X sexorient
AIC(fit7)
## [1] 6039.923
#sex x sexorient x race
AIC(fit8)
## [1] 6040.185
#fit<-svyolr(score.grp~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,des)
#summary(fit)
#library(VGAM)
#Proportional odds
#fit.vgam<-vglm(as.ordered(score.grp)~sex+racethnic+sexorient+born_US_foreign+W4_educ+incomeW4,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = T, reverse = T)) #<-parallel = T == proportional odds
#summary(fit.vgam)
#continuation Ratio
#fit.crp<-vglm(as.ordered(score.grp)~sex+racethnic,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cratio( parallel = F, reverse = T))
#Non-proportional odds
#fit.vgam2<-vglm(as.ordered(score.grp)~sex+racethnic+sexorient+W4_UScitizenship+W4_educ+incomeW4,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cumulative(parallel = F, reverse = T)) #<-parallel = F == Nonproportional odds
is.parallel(fit5)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
## W4_educ born_US_foreign W4.married incomeW4 drink
## TRUE TRUE TRUE TRUE TRUE
## smoke health mom.race mom_educ mom_emp
## TRUE TRUE TRUE TRUE TRUE
## good.neigh
## TRUE
table(addhealth.pro$raceXsexorient)
##
## a-nhwhite.a_straight b-nhblack.a_straight c-hispanic.a_straight
## 2311 379 95
## a-nhwhite.b_bisexual b-nhblack.b_bisexual c-hispanic.b_bisexual
## 61 7 2
## a-nhwhite.c_LGB b-nhblack.c_LGB c-hispanic.c_LGB
## 59 21 1
table(addhealth.pro$raceXsexorient,addhealth.pro$W4_educ)
##
## a_less_highschool b_HS_grad d_college_bach.plus
## a-nhwhite.a_straight 213 1583 515
## b-nhblack.a_straight 24 246 109
## c-hispanic.a_straight 10 66 19
## a-nhwhite.b_bisexual 17 26 18
## b-nhblack.b_bisexual 1 6 0
## c-hispanic.b_bisexual 0 2 0
## a-nhwhite.c_LGB 2 46 11
## b-nhblack.c_LGB 0 21 0
## c-hispanic.c_LGB 0 0 1
table(addhealth.pro$raceXsexorient,addhealth.pro$sex)
##
## a_male b_female
## a-nhwhite.a_straight 673 1638
## b-nhblack.a_straight 79 300
## c-hispanic.a_straight 35 60
## a-nhwhite.b_bisexual 8 53
## b-nhblack.b_bisexual 0 7
## c-hispanic.b_bisexual 2 0
## a-nhwhite.c_LGB 28 31
## b-nhblack.c_LGB 7 14
## c-hispanic.c_LGB 0 1
#continuation Ratio
fit.crp<-vglm(as.ordered(score.grp)~agew4+sex+racethnic+sexorient+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cratio( parallel = F, reverse = T))
summary(fit.crp)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + sex + racethnic +
## sexorient + W4_educ + born_US_foreign + W4.married + incomeW4 +
## drink + smoke + health + mom.race + mom_educ + mom_emp +
## good.neigh, family = cratio(parallel = F, reverse = T), data = addhealth.pro,
## weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y<2|Y<=2]) -2.566 -0.6293 -5.567e-06 0.3110 11.403
## logitlink(P[Y<3|Y<=3]) -4.259 -0.5655 3.334e-01 0.6332 2.731
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 2.36142 0.84729 2.787 0.005320 **
## (Intercept):2 1.52146 0.73033 2.083 0.037230 *
## agew4:1 -0.13488 0.02846 -4.740 2.14e-06 ***
## agew4:2 -0.02026 0.02434 -0.832 0.405132
## sexb_female:1 0.55492 0.11009 5.040 4.64e-07 ***
## sexb_female:2 0.06067 0.09362 0.648 0.516918
## racethnicb-nhblack:1 1.09624 0.73403 1.493 0.135319
## racethnicb-nhblack:2 1.13192 0.64504 1.755 0.079294 .
## racethnicc-hispanic:1 1.15273 0.46868 2.460 0.013912 *
## racethnicc-hispanic:2 -1.31530 0.32548 -4.041 5.32e-05 ***
## sexorientb_bisexual:1 0.03306 0.32880 0.101 0.919919
## sexorientb_bisexual:2 0.21339 0.28931 0.738 0.460764
## sexorientc_LGB:1 0.04686 0.37567 0.125 0.900742
## sexorientc_LGB:2 -0.60620 0.28022 -2.163 0.030516 *
## W4_educb_HS_grad:1 0.51319 0.18510 2.772 0.005563 **
## W4_educb_HS_grad:2 -0.58486 0.15447 -3.786 0.000153 ***
## W4_educd_college_bach.plus:1 1.16432 0.21897 5.317 1.05e-07 ***
## W4_educd_college_bach.plus:2 -0.21791 0.19242 -1.132 0.257453
## born_US_foreignb_born.foreign:1 0.14213 0.83498 0.170 0.864836
## born_US_foreignb_born.foreign:2 -0.98829 0.47340 -2.088 0.036830 *
## W4.marriedb_married_once:1 0.41723 0.10886 3.833 0.000127 ***
## W4.marriedb_married_once:2 -0.06232 0.09662 -0.645 0.518936
## W4.marriedc_married_twice+:1 -0.10375 0.21189 -0.490 0.624394
## W4.marriedc_married_twice+:2 0.05575 0.17575 0.317 0.751108
## incomeW4b_$25,000<$75,000:1 -0.38004 0.13775 -2.759 0.005798 **
## incomeW4b_$25,000<$75,000:2 0.32399 0.11135 2.910 0.003619 **
## incomeW4c_$75,000<$150,000:1 -0.13233 0.17038 -0.777 0.437333
## incomeW4c_$75,000<$150,000:2 0.45650 0.15218 3.000 0.002702 **
## incomeW4d_>$150,000:1 0.52623 0.29573 1.779 0.075169 .
## incomeW4d_>$150,000:2 1.00804 0.31858 3.164 0.001555 **
## drinkb_lite.drinker:1 -0.17441 0.13961 -1.249 0.211588
## drinkb_lite.drinker:2 -0.08208 0.11432 -0.718 0.472764
## drinkc_drinker:1 0.17250 0.16401 1.052 0.292925
## drinkc_drinker:2 0.19493 0.14786 1.318 0.187392
## smokeb_tried.smoking:1 -0.40362 0.16848 -2.396 0.016592 *
## smokeb_tried.smoking:2 0.43186 0.14284 3.023 0.002499 **
## smokec_smoker:1 0.43725 0.12691 3.445 0.000570 ***
## smokec_smoker:2 0.56634 0.10489 5.399 6.69e-08 ***
## healthb_Bad.health:1 -0.48904 0.14229 -3.437 0.000588 ***
## healthb_Bad.health:2 -0.90054 0.10199 -8.830 < 2e-16 ***
## mom.raceb_black:1 -1.42401 0.74638 -1.908 0.056405 .
## mom.raceb_black:2 -0.97088 0.64561 -1.504 0.132629
## mom.racec_latino:1 -1.57165 0.30109 -5.220 1.79e-07 ***
## mom.racec_latino:2 0.42864 0.24766 1.731 0.083492 .
## mom_educb_HSgrad/voca:1 0.34824 0.15616 2.230 0.025746 *
## mom_educb_HSgrad/voca:2 -0.18103 0.12331 -1.468 0.142072
## mom_educc_col.grad:1 0.44092 0.19784 2.229 0.025837 *
## mom_educc_col.grad:2 0.05195 0.17512 0.297 0.766745
## mom_empb_yes.emp:1 -0.11089 0.12001 -0.924 0.355481
## mom_empb_yes.emp:2 0.18246 0.09869 1.849 0.064476 .
## good.neigh1:1 0.12816 0.13326 0.962 0.336215
## good.neigh1:2 -0.01936 0.12310 -0.157 0.875028
## good.neigh2:1 -0.10266 0.14015 -0.733 0.463857
## good.neigh2:2 -0.14300 0.12777 -1.119 0.263030
## good.neigh3:1 -0.42496 0.15100 -2.814 0.004888 **
## good.neigh3:2 -0.11330 0.13257 -0.855 0.392764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3])
##
## Residual deviance: 5861.449 on 5816 degrees of freedom
##
## Log-likelihood: -2930.725 on 5816 degrees of freedom
##
## Number of Fisher scoring iterations: 8
##
## No Hauck-Donner effect found in any of the estimates
#intersection
fit.crp2<-vglm(as.ordered(score.grp)~agew4+seXorientXrace+W4_educ+born_US_foreign+W4.married+incomeW4+drink+smoke+health+mom.race+mom_educ+mom_emp+good.neigh,addhealth.pro, weights = GSWGT4_2/mean(GSWGT4_2, na.rm=T), family=cratio( parallel = F, reverse = T))
summary(fit.crp2)
##
## Call:
## vglm(formula = as.ordered(score.grp) ~ agew4 + seXorientXrace +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh,
## family = cratio(parallel = F, reverse = T), data = addhealth.pro,
## weights = GSWGT4_2/mean(GSWGT4_2, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y<2|Y<=2]) -3.020 -0.6400 -2.399e-07 0.3051 10.451
## logitlink(P[Y<3|Y<=3]) -4.332 -0.5613 3.328e-01 0.6285 2.467
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept):1 2.318e+00 8.567e-01 2.706
## (Intercept):2 1.925e+00 7.413e-01 2.597
## agew4:1 -1.343e-01 2.871e-02 -4.676
## agew4:2 -2.796e-02 2.460e-02 -1.136
## seXorientXraceb_female.a_straight.a-nhwhite:1 6.121e-01 1.172e-01 5.225
## seXorientXraceb_female.a_straight.a-nhwhite:2 -4.776e-02 1.032e-01 -0.463
## seXorientXracea_male.b_bisexual.a-nhwhite:1 2.115e+00 7.667e-01 2.759
## seXorientXracea_male.b_bisexual.a-nhwhite:2 1.621e+01 6.503e+02 NA
## seXorientXraceb_female.b_bisexual.a-nhwhite:1 2.071e-01 4.157e-01 0.498
## seXorientXraceb_female.b_bisexual.a-nhwhite:2 -2.904e-01 3.192e-01 -0.910
## seXorientXracea_male.c_LGB.a-nhwhite:1 6.284e-01 7.029e-01 0.894
## seXorientXracea_male.c_LGB.a-nhwhite:2 -1.967e+00 4.259e-01 -4.619
## seXorientXraceb_female.c_LGB.a-nhwhite:1 2.276e-02 5.022e-01 0.045
## seXorientXraceb_female.c_LGB.a-nhwhite:2 2.434e-01 4.601e-01 0.529
## seXorientXracea_male.a_straight.b-nhblack:1 9.420e-01 9.276e-01 1.016
## seXorientXracea_male.a_straight.b-nhblack:2 -8.736e-02 7.572e-01 -0.115
## seXorientXraceb_female.a_straight.b-nhblack:1 1.551e+00 8.686e-01 1.785
## seXorientXraceb_female.a_straight.b-nhblack:2 7.925e-01 7.378e-01 1.074
## seXorientXraceb_female.b_bisexual.b-nhblack:1 -1.512e+01 1.672e+03 NA
## seXorientXraceb_female.b_bisexual.b-nhblack:2 1.613e+01 9.965e+02 NA
## seXorientXracea_male.c_LGB.b-nhblack:1 1.692e+00 2.144e+00 0.789
## seXorientXracea_male.c_LGB.b-nhblack:2 1.612e+01 2.019e+03 NA
## seXorientXraceb_female.c_LGB.b-nhblack:1 3.081e+00 1.277e+00 2.412
## seXorientXraceb_female.c_LGB.b-nhblack:2 1.598e+01 8.402e+02 NA
## seXorientXracea_male.a_straight.c-hispanic:1 6.119e-01 6.897e-01 0.887
## seXorientXracea_male.a_straight.c-hispanic:2 -1.426e+00 4.164e-01 -3.423
## seXorientXraceb_female.a_straight.c-hispanic:1 2.140e+00 6.042e-01 3.541
## seXorientXraceb_female.a_straight.c-hispanic:2 -1.417e+00 4.224e-01 -3.355
## seXorientXracea_male.b_bisexual.c-hispanic:1 2.403e+00 1.485e+04 0.000
## seXorientXracea_male.b_bisexual.c-hispanic:2 -1.811e+01 7.425e+03 NA
## seXorientXraceb_female.c_LGB.c-hispanic:1 -1.553e+01 1.479e+04 NA
## seXorientXraceb_female.c_LGB.c-hispanic:2 1.531e+01 9.496e+03 NA
## W4_educb_HS_grad:1 4.808e-01 1.872e-01 2.568
## W4_educb_HS_grad:2 -6.522e-01 1.563e-01 -4.173
## W4_educd_college_bach.plus:1 1.124e+00 2.215e-01 5.076
## W4_educd_college_bach.plus:2 -2.778e-01 1.954e-01 -1.422
## born_US_foreignb_born.foreign:1 1.446e-01 8.351e-01 0.173
## born_US_foreignb_born.foreign:2 -7.415e-01 4.859e-01 -1.526
## W4.marriedb_married_once:1 4.338e-01 1.097e-01 3.954
## W4.marriedb_married_once:2 -9.264e-02 9.753e-02 -0.950
## W4.marriedc_married_twice+:1 -5.247e-02 2.128e-01 -0.247
## W4.marriedc_married_twice+:2 6.790e-02 1.771e-01 0.383
## incomeW4b_$25,000<$75,000:1 -3.483e-01 1.401e-01 -2.486
## incomeW4b_$25,000<$75,000:2 3.556e-01 1.125e-01 3.161
## incomeW4c_$75,000<$150,000:1 -1.095e-01 1.722e-01 -0.636
## incomeW4c_$75,000<$150,000:2 5.067e-01 1.536e-01 3.299
## incomeW4d_>$150,000:1 5.182e-01 2.997e-01 1.729
## incomeW4d_>$150,000:2 1.129e+00 3.271e-01 3.451
## drinkb_lite.drinker:1 -2.047e-01 1.413e-01 -1.449
## drinkb_lite.drinker:2 -1.049e-01 1.156e-01 -0.908
## drinkc_drinker:1 1.666e-01 1.651e-01 1.009
## drinkc_drinker:2 2.098e-01 1.498e-01 1.400
## smokeb_tried.smoking:1 -3.945e-01 1.689e-01 -2.336
## smokeb_tried.smoking:2 4.601e-01 1.445e-01 3.183
## smokec_smoker:1 4.233e-01 1.279e-01 3.310
## smokec_smoker:2 5.644e-01 1.057e-01 5.338
## healthb_Bad.health:1 -5.637e-01 1.460e-01 -3.860
## healthb_Bad.health:2 -9.335e-01 1.033e-01 -9.039
## mom.raceb_black:1 -1.268e+00 8.590e-01 -1.477
## mom.raceb_black:2 -5.094e-01 7.238e-01 -0.704
## mom.racec_latino:1 -1.546e+00 3.010e-01 -5.137
## mom.racec_latino:2 4.694e-01 2.505e-01 1.874
## mom_educb_HSgrad/voca:1 3.654e-01 1.579e-01 2.314
## mom_educb_HSgrad/voca:2 -2.025e-01 1.249e-01 -1.621
## mom_educc_col.grad:1 4.415e-01 1.993e-01 2.215
## mom_educc_col.grad:2 5.848e-02 1.768e-01 0.331
## mom_empb_yes.emp:1 -9.496e-02 1.212e-01 -0.783
## mom_empb_yes.emp:2 1.685e-01 9.975e-02 1.689
## good.neigh1:1 1.126e-01 1.344e-01 0.838
## good.neigh1:2 -4.778e-02 1.247e-01 -0.383
## good.neigh2:1 -1.084e-01 1.407e-01 -0.770
## good.neigh2:2 -1.692e-01 1.285e-01 -1.316
## good.neigh3:1 -4.228e-01 1.517e-01 -2.787
## good.neigh3:2 -1.242e-01 1.335e-01 -0.930
## Pr(>|z|)
## (Intercept):1 0.006809 **
## (Intercept):2 0.009403 **
## agew4:1 2.92e-06 ***
## agew4:2 0.255804
## seXorientXraceb_female.a_straight.a-nhwhite:1 1.74e-07 ***
## seXorientXraceb_female.a_straight.a-nhwhite:2 0.643654
## seXorientXracea_male.b_bisexual.a-nhwhite:1 0.005805 **
## seXorientXracea_male.b_bisexual.a-nhwhite:2 NA
## seXorientXraceb_female.b_bisexual.a-nhwhite:1 0.618373
## seXorientXraceb_female.b_bisexual.a-nhwhite:2 0.362963
## seXorientXracea_male.c_LGB.a-nhwhite:1 0.371344
## seXorientXracea_male.c_LGB.a-nhwhite:2 3.86e-06 ***
## seXorientXraceb_female.c_LGB.a-nhwhite:1 0.963857
## seXorientXraceb_female.c_LGB.a-nhwhite:2 0.596823
## seXorientXracea_male.a_straight.b-nhblack:1 0.309830
## seXorientXracea_male.a_straight.b-nhblack:2 0.908155
## seXorientXraceb_female.a_straight.b-nhblack:1 0.074228 .
## seXorientXraceb_female.a_straight.b-nhblack:2 0.282724
## seXorientXraceb_female.b_bisexual.b-nhblack:1 NA
## seXorientXraceb_female.b_bisexual.b-nhblack:2 NA
## seXorientXracea_male.c_LGB.b-nhblack:1 0.429977
## seXorientXracea_male.c_LGB.b-nhblack:2 NA
## seXorientXraceb_female.c_LGB.b-nhblack:1 0.015873 *
## seXorientXraceb_female.c_LGB.b-nhblack:2 NA
## seXorientXracea_male.a_straight.c-hispanic:1 0.374969
## seXorientXracea_male.a_straight.c-hispanic:2 0.000619 ***
## seXorientXraceb_female.a_straight.c-hispanic:1 0.000398 ***
## seXorientXraceb_female.a_straight.c-hispanic:2 0.000792 ***
## seXorientXracea_male.b_bisexual.c-hispanic:1 0.999871
## seXorientXracea_male.b_bisexual.c-hispanic:2 NA
## seXorientXraceb_female.c_LGB.c-hispanic:1 NA
## seXorientXraceb_female.c_LGB.c-hispanic:2 NA
## W4_educb_HS_grad:1 0.010232 *
## W4_educb_HS_grad:2 3.01e-05 ***
## W4_educd_college_bach.plus:1 3.85e-07 ***
## W4_educd_college_bach.plus:2 0.154975
## born_US_foreignb_born.foreign:1 0.862564
## born_US_foreignb_born.foreign:2 0.126999
## W4.marriedb_married_once:1 7.68e-05 ***
## W4.marriedb_married_once:2 0.342189
## W4.marriedc_married_twice+:1 0.805193
## W4.marriedc_married_twice+:2 0.701455
## incomeW4b_$25,000<$75,000:1 0.012906 *
## incomeW4b_$25,000<$75,000:2 0.001571 **
## incomeW4c_$75,000<$150,000:1 0.524711
## incomeW4c_$75,000<$150,000:2 0.000969 ***
## incomeW4d_>$150,000:1 0.083796 .
## incomeW4d_>$150,000:2 0.000559 ***
## drinkb_lite.drinker:1 0.147466
## drinkb_lite.drinker:2 0.364085
## drinkc_drinker:1 0.312848
## drinkc_drinker:2 0.161465
## smokeb_tried.smoking:1 0.019491 *
## smokeb_tried.smoking:2 0.001457 **
## smokec_smoker:1 0.000933 ***
## smokec_smoker:2 9.38e-08 ***
## healthb_Bad.health:1 0.000113 ***
## healthb_Bad.health:2 < 2e-16 ***
## mom.raceb_black:1 0.139792
## mom.raceb_black:2 0.481566
## mom.racec_latino:1 2.80e-07 ***
## mom.racec_latino:2 0.060992 .
## mom_educb_HSgrad/voca:1 0.020664 *
## mom_educb_HSgrad/voca:2 0.105007
## mom_educc_col.grad:1 0.026754 *
## mom_educc_col.grad:2 0.740877
## mom_empb_yes.emp:1 0.433339
## mom_empb_yes.emp:2 0.091160 .
## good.neigh1:1 0.402209
## good.neigh1:2 0.701523
## good.neigh2:1 0.441172
## good.neigh2:2 0.188101
## good.neigh3:1 0.005324 **
## good.neigh3:2 0.352358
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3])
##
## Residual deviance: 5795.419 on 5798 degrees of freedom
##
## Log-likelihood: -2897.709 on 5798 degrees of freedom
##
## Number of Fisher scoring iterations: 15
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## 'seXorientXracea_male.b_bisexual.a-nhwhite:2', 'seXorientXraceb_female.b_bisexual.b-nhblack:1', 'seXorientXraceb_female.b_bisexual.b-nhblack:2', 'seXorientXracea_male.c_LGB.b-nhblack:2', 'seXorientXraceb_female.c_LGB.b-nhblack:2', 'seXorientXracea_male.b_bisexual.c-hispanic:2', 'seXorientXraceb_female.c_LGB.c-hispanic:1', 'seXorientXraceb_female.c_LGB.c-hispanic:2'
AIC(fit.crp)
## [1] 5973.449
AIC(fit.crp2)
## [1] 5943.419
colSums(depvar(fit.crp))
## 0 1 2
## 799 1232 905
lrtest(fit.crp, fit.crp2)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh
## Model 2: as.ordered(score.grp) ~ agew4 + seXorientXrace + W4_educ + born_US_foreign +
## W4.married + incomeW4 + drink + smoke + health + mom.race +
## mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5816 -2930.7
## 2 5798 -2897.7 -18 66.031 2.096e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(fit, fit5)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient
## Model 2: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5864 -3152.6
## 2 5843 -3004.3 -21 296.59 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(fit,fit6)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient
## Model 2: as.ordered(score.grp) ~ agew4 + sex + raceXsexorient + W4_educ +
## born_US_foreign + W4.married + incomeW4 + drink + smoke +
## health + mom.race + mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5864 -3152.6
## 2 5839 -2998.9 -25 307.36 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(fit,fit7)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient
## Model 2: as.ordered(score.grp) ~ agew4 + sexXsexorient + racethnic + W4_educ +
## born_US_foreign + W4.married + incomeW4 + drink + smoke +
## health + mom.race + mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5864 -3152.6
## 2 5841 -2989.0 -23 327.24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(fit,fit8)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient
## Model 2: as.ordered(score.grp) ~ agew4 + seXorientXrace + W4_educ + born_US_foreign +
## W4.married + incomeW4 + drink + smoke + health + mom.race +
## mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5864 -3152.6
## 2 5834 -2982.1 -30 340.97 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(fit,fit2,fit3,fit4,fit5,fit5,fit7,fit8)
## Likelihood ratio test
##
## Model 1: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient
## Model 2: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign
## Model 3: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4
## Model 4: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health
## Model 5: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh
## Model 6: as.ordered(score.grp) ~ agew4 + sex + racethnic + sexorient +
## W4_educ + born_US_foreign + W4.married + incomeW4 + drink +
## smoke + health + mom.race + mom_educ + mom_emp + good.neigh
## Model 7: as.ordered(score.grp) ~ agew4 + sexXsexorient + racethnic + W4_educ +
## born_US_foreign + W4.married + incomeW4 + drink + smoke +
## health + mom.race + mom_educ + mom_emp + good.neigh
## Model 8: as.ordered(score.grp) ~ agew4 + seXorientXrace + W4_educ + born_US_foreign +
## W4.married + incomeW4 + drink + smoke + health + mom.race +
## mom_educ + mom_emp + good.neigh
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5864 -3152.6
## 2 5861 -3118.8 -3 67.666 1.349e-14 ***
## 3 5856 -3096.9 -5 43.623 2.762e-08 ***
## 4 5851 -3015.4 -5 162.992 < 2.2e-16 ***
## 5 5843 -3004.3 -8 22.305 0.004381 **
## 6 5843 -3004.3 0 0.000 1.000000
## 7 5841 -2989.0 -2 30.649 2.211e-07 ***
## 8 5834 -2982.1 -7 13.738 0.056038 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
is.parallel(fit)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
is.parallel(fit2)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
## W4_educ born_US_foreign
## TRUE TRUE
is.parallel(fit3)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
## W4_educ born_US_foreign W4.married incomeW4
## TRUE TRUE TRUE TRUE
is.parallel(fit4)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
## W4_educ born_US_foreign W4.married incomeW4 drink
## TRUE TRUE TRUE TRUE TRUE
## smoke health
## TRUE TRUE
is.parallel(fit5)
## (Intercept) agew4 sex racethnic sexorient
## FALSE TRUE TRUE TRUE TRUE
## W4_educ born_US_foreign W4.married incomeW4 drink
## TRUE TRUE TRUE TRUE TRUE
## smoke health mom.race mom_educ mom_emp
## TRUE TRUE TRUE TRUE TRUE
## good.neigh
## TRUE
is.parallel(fit6)
## (Intercept) agew4 sex raceXsexorient W4_educ
## FALSE TRUE TRUE TRUE TRUE
## born_US_foreign W4.married incomeW4 drink smoke
## TRUE TRUE TRUE TRUE TRUE
## health mom.race mom_educ mom_emp good.neigh
## TRUE TRUE TRUE TRUE TRUE
is.parallel(fit7)
## (Intercept) agew4 sexXsexorient racethnic W4_educ
## FALSE TRUE TRUE TRUE TRUE
## born_US_foreign W4.married incomeW4 drink smoke
## TRUE TRUE TRUE TRUE TRUE
## health mom.race mom_educ mom_emp good.neigh
## TRUE TRUE TRUE TRUE TRUE
is.parallel(fit8)
## (Intercept) agew4 seXorientXrace W4_educ born_US_foreign
## FALSE TRUE TRUE TRUE TRUE
## W4.married incomeW4 drink smoke health
## TRUE TRUE TRUE TRUE TRUE
## mom.race mom_educ mom_emp good.neigh
## TRUE TRUE TRUE TRUE
sumOrd <- summary(fit8)
(coefOrd <- coef(sumOrd))
## Estimate Std. Error
## (Intercept):1 -0.88639774 0.61480761
## (Intercept):2 -2.87565400 0.61713865
## agew4 0.07752854 0.02051444
## seXorientXraceb_female.a_straight.a-nhwhite -0.21800558 0.08435085
## seXorientXracea_male.b_bisexual.a-nhwhite -2.99925027 0.75424639
## seXorientXraceb_female.b_bisexual.a-nhwhite 0.08675178 0.28032037
## seXorientXracea_male.c_LGB.a-nhwhite 1.46023185 0.40188538
## seXorientXraceb_female.c_LGB.a-nhwhite -0.17987537 0.36430257
## seXorientXracea_male.a_straight.b-nhblack -0.40523469 0.60606546
## seXorientXraceb_female.a_straight.b-nhblack -1.24348638 0.57671597
## seXorientXraceb_female.b_bisexual.b-nhblack -1.06735487 1.03768288
## seXorientXracea_male.c_LGB.b-nhblack -2.21407641 1.92294222
## seXorientXraceb_female.c_LGB.b-nhblack -3.53934182 1.18116131
## seXorientXracea_male.a_straight.c-hispanic 0.63069121 0.37015896
## seXorientXraceb_female.a_straight.c-hispanic 0.41152606 0.37113604
## seXorientXracea_male.b_bisexual.c-hispanic 11.54097473 447.96372431
## seXorientXraceb_female.c_LGB.c-hispanic -0.15298777 8.18135437
## W4_educb_HS_grad 0.25759792 0.12715588
## W4_educd_college_bach.plus -0.29032893 0.15570946
## born_US_foreignb_born.foreign 0.24892858 0.42229877
## W4.marriedb_married_once -0.09352238 0.07954451
## W4.marriedc_married_twice+ -0.02109082 0.14845254
## incomeW4b_$25,000<$75,000 -0.14129245 0.09639082
## incomeW4c_$75,000<$150,000 -0.38758279 0.12508862
## incomeW4d_>$150,000 -1.05020778 0.23216093
## drinkb_lite.drinker 0.14383626 0.09838654
## drinkc_drinker -0.17134473 0.12056875
## smokeb_tried.smoking -0.19965579 0.11776354
## smokec_smoker -0.58543986 0.08981617
## healthb_Bad.health 0.89497107 0.09250001
## mom.raceb_black 0.94754506 0.56679205
## mom.racec_latino 0.31278579 0.18418453
## mom_educb_HSgrad/voca 0.04609403 0.10548589
## mom_educc_col.grad -0.14831665 0.14180015
## mom_empb_yes.emp -0.09563567 0.08476438
## good.neigh1 -0.02058607 0.09996961
## good.neigh2 0.13455783 0.10445666
## good.neigh3 0.19790696 0.10936178
## z value Pr(>|z|)
## (Intercept):1 -1.44174817 1.493734e-01
## (Intercept):2 -4.65965631 3.167378e-06
## agew4 3.77921795 1.573217e-04
## seXorientXraceb_female.a_straight.a-nhwhite -2.58450966 9.751757e-03
## seXorientXracea_male.b_bisexual.a-nhwhite -3.97648607 6.994109e-05
## seXorientXraceb_female.b_bisexual.a-nhwhite 0.30947368 7.569612e-01
## seXorientXracea_male.c_LGB.a-nhwhite 3.63345355 2.796528e-04
## seXorientXraceb_female.c_LGB.a-nhwhite -0.49375268 6.214809e-01
## seXorientXracea_male.a_straight.b-nhblack -0.66863189 5.037303e-01
## seXorientXraceb_female.a_straight.b-nhblack -2.15615040 3.107192e-02
## seXorientXraceb_female.b_bisexual.b-nhblack -1.02859447 3.036703e-01
## seXorientXracea_male.c_LGB.b-nhblack -1.15140039 2.495676e-01
## seXorientXraceb_female.c_LGB.b-nhblack -2.99649319 2.731043e-03
## seXorientXracea_male.a_straight.c-hispanic 1.70383881 8.841121e-02
## seXorientXraceb_female.a_straight.c-hispanic 1.10882806 2.675044e-01
## seXorientXracea_male.b_bisexual.c-hispanic 0.02576319 9.794462e-01
## seXorientXraceb_female.c_LGB.c-hispanic -0.01869957 9.850808e-01
## W4_educb_HS_grad 2.02584356 4.278082e-02
## W4_educd_college_bach.plus -1.86455548 6.224374e-02
## born_US_foreignb_born.foreign 0.58946081 5.555522e-01
## W4.marriedb_married_once -1.17572399 2.397052e-01
## W4.marriedc_married_twice+ -0.14207111 8.870238e-01
## incomeW4b_$25,000<$75,000 -1.46582890 1.426949e-01
## incomeW4c_$75,000<$150,000 -3.09846560 1.945256e-03
## incomeW4d_>$150,000 -4.52361975 6.079090e-06
## drinkb_lite.drinker 1.46195055 1.437548e-01
## drinkc_drinker -1.42113708 1.552769e-01
## smokeb_tried.smoking -1.69539555 9.000041e-02
## smokec_smoker -6.51820106 7.115557e-11
## healthb_Bad.health 9.67536178 3.837320e-22
## mom.raceb_black 1.67176843 9.456999e-02
## mom.racec_latino 1.69821966 8.946631e-02
## mom_educb_HSgrad/voca 0.43696872 6.621340e-01
## mom_educc_col.grad -1.04595552 2.955816e-01
## mom_empb_yes.emp -1.12825304 2.592131e-01
## good.neigh1 -0.20592333 8.368508e-01
## good.neigh2 1.28816894 1.976872e-01
## good.neigh3 1.80965374 7.034950e-02
sumOrd2 <- summary(fit5)
(coefOrd <- coef(sumOrd2))
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.666110615 0.60931619 -1.09321010 2.743016e-01
## (Intercept):2 -2.633490525 0.61136278 -4.30757421 1.650547e-05
## agew4 0.073687424 0.02038317 3.61511045 3.002196e-04
## sexb_female -0.265287146 0.07803793 -3.39946397 6.751807e-04
## racethnicb-nhblack -1.490817689 0.52175073 -2.85733707 4.272119e-03
## racethnicc-hispanic 0.574128125 0.28002594 2.05026762 4.033832e-02
## sexorientb_bisexual -0.106475384 0.24294308 -0.43827297 6.611884e-01
## sexorientc_LGB 0.376657705 0.24710122 1.52430535 1.274325e-01
## W4_educb_HS_grad 0.201661588 0.12626497 1.59713011 1.102367e-01
## W4_educd_college_bach.plus -0.354120549 0.15441923 -2.29324119 2.183412e-02
## born_US_foreignb_born.foreign 0.417073617 0.41859550 0.99636431 3.190732e-01
## W4.marriedb_married_once -0.108242095 0.07923278 -1.36612761 1.718989e-01
## W4.marriedc_married_twice+ -0.005761287 0.14803283 -0.03891898 9.689550e-01
## incomeW4b_$25,000<$75,000 -0.113494133 0.09572437 -1.18563464 2.357666e-01
## incomeW4c_$75,000<$150,000 -0.346356425 0.12458813 -2.78001141 5.435699e-03
## incomeW4d_>$150,000 -0.974377200 0.22976166 -4.24081728 2.227074e-05
## drinkb_lite.drinker 0.124044291 0.09762598 1.27060744 2.038683e-01
## drinkc_drinker -0.175297778 0.11993052 -1.46166114 1.438341e-01
## smokeb_tried.smoking -0.186025509 0.11736213 -1.58505568 1.129537e-01
## smokec_smoker -0.593026060 0.08951204 -6.62509824 3.470173e-11
## healthb_Bad.health 0.855052100 0.09152054 9.34273460 9.387755e-21
## mom.raceb_black 1.524265027 0.52639945 2.89564327 3.783823e-03
## mom.racec_latino 0.335512677 0.18381111 1.82531233 6.795391e-02
## mom_educb_HSgrad/voca 0.038288912 0.10466771 0.36581398 7.145039e-01
## mom_educc_col.grad -0.136361553 0.14101851 -0.96697630 3.335559e-01
## mom_empb_yes.emp -0.106211395 0.08415465 -1.26209772 2.069136e-01
## good.neigh1 -0.057986630 0.09938076 -0.58347946 5.595706e-01
## good.neigh2 0.107634682 0.10412023 1.03375376 3.012513e-01
## good.neigh3 0.187959252 0.10897300 1.72482406 8.455918e-02
sumOrd3 <- summary(fit)
(coefOrd <- coef(sumOrd3))
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.73415764 0.56310038 -1.3037776 1.923094e-01
## (Intercept):2 -2.55267409 0.56516132 -4.5167176 6.280558e-06
## agew4 0.06128111 0.01923115 3.1865547 1.439782e-03
## sexb_female -0.23150793 0.07420239 -3.1199523 1.808803e-03
## racethnicb-nhblack 0.30564837 0.11407145 2.6794468 7.374393e-03
## racethnicc-hispanic 0.93496649 0.23102892 4.0469674 5.188547e-05
## sexorientb_bisexual 0.12733748 0.23218627 0.5484281 5.833980e-01
## sexorientc_LGB 0.26261971 0.23495232 1.1177575 2.636706e-01
table(addhealth$agew4)
##
## 24 25 26 27 28 29 30 31 32 33 34
## 1 84 972 1217 1903 1940 2170 1740 601 59 24