Untitled
Steven Weeks
10/31/2019
library(haven)#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")
addhealth.wghtW4<-read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0028/da27021p28.dta")
#addhealth<-merge(wave1,wave3,by="aid")
addhealth1<-merge(wave1,wave3,by="aid")
addhealth1<-merge(addhealth1,wave4,by="aid")
addhealth1<-merge(addhealth1,addhealth.wghtW4,by="aid")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#select variables to use
addhealth<-addhealth1%>%
select(GSWGT4,region,psuscid,BIO_SEX,H3SE13,H1GI6A,H1GI6B,H1GI4,H1GI6D,H1GI6E,H4ED2,iyear,imonth,H1GI1M,H1GI1Y,IYEAR3,IMONTH3,H3OD1Y,H3OD1M,H4OD1M,H4OD1Y,IMONTH4,IYEAR4,H4MH22,H3GH1,H4GH1,H1GI6C,H1FS6,H3SP9,H4DS18,H3DS18F,H4MH22,H4GH1,H4EC1)
#complete cases
addhealth<-addhealth%>%
filter(complete.cases(.))library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode#general health
addhealth$W4.general_health<-Recode(addhealth$H4GH1,recodes="1:3='c_good/vg/excellent';4:5='b_poor/fair';else=NA")
#depression
addhealth$W4.depression<-Recode(addhealth$H4MH22,recodes="0:1='a_never_some';2:3='b_most_always';else=NA")
#Date of birth variables; Wave 1: H1GI1M (month), H1GI1Y (year); Wave 3: H3OD1M (month), H3OD1Y (year); Wave 4: H4OD1M (month), H4OD1Y (year)
#birth year w1
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)))))))))))
#birth month w1
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")
#combine year and month for a birth date wave 1
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericaddhealth$birthdatew1 <- as.yearmon(paste(addhealth$birthyearw1, addhealth$birthmonthw1), "%Y %m")
#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")
#birth date wave 3
addhealth$birthdatew3<-as.yearmon(paste(addhealth$H3OD1Y,addhealth$H3OD1M),"%Y %m")
#interview date w1
addhealth$iyearfix <- Recode(addhealth$iyear,recodes="'94'='1994';'95'='1995'")
addhealth$interviewdatew1 <- as.yearmon(paste(addhealth$iyearfix, addhealth$imonth), "%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")
#interview date w3
addhealth$interviewdatew3<-as.yearmon(paste(addhealth$IYEAR3,addhealth$IMONTH3 ),"%Y %m")
#age - derived from date of birth subtracted from wave interview date
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)))
#sex variable - Male is reference group
addhealth$sex <- ifelse(addhealth$BIO_SEX==2,2,1)
addhealth$sex<-Recode(addhealth$sex, recodes="1='a_male'; 2='b_female'",as.factor=T)
#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)
#race variables - Reference group is nhwhite
#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,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",
ifelse(addhealth$asian=="asian","d-asian",
ifelse(addhealth$native_american=="native_american","e-native_american",
ifelse(addhealth$other=="other","f-other",NA)))))))
#Education. Less then high school is reference group
addhealth$educ<-Recode(addhealth$H4ED2,recodes="1:2='a_less_highschool';3:6='b_highschool_grad';7='c_college_bach';8:13='d_college+';else=NA")
#income variable
addhealth$incomeW4<-Recode(addhealth$H4EC1,recodes="1:5='a_<$25k';6:8='b_$25k>$50k';9:10='c_$50k>$100k';11:12='e_$100k+';else=NA")#transition
addhealth$W3.jumped.beatenup<-Recode(addhealth$H3DS18F, recodes="0='a.no'; 1='b.yes';else=NA")
addhealth$W4.jumped.beatenup<-Recode(addhealth$H4DS18, recodes="0='a.no'; 1='b.yes';else=NA")
addhealth$transition1<-ifelse(addhealth$W3.jumped.beatenup=="a.no",0,1)
addhealth$transition2<-ifelse(addhealth$W4.jumped.beatenup=="b.yes",1,0)
addhealth$transition<-ifelse(addhealth$transition1==0&addhealth$transition2==1,1,0)
table(addhealth$transition) ##
## 0 1
## 8386 1001#age at transition
addhealth$age_transition<-ifelse(addhealth$transition==1,
addhealth$agew3,addhealth$agew4)
table(addhealth$age_transition)##
## 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
## 45 139 138 221 217 178 111 997 1184 1707 1648 1827 740 213 19
## 34
## 3#select variables to use
addhealth<-addhealth%>%
select(psuscid,region,GSWGT4,sex,sexorient,racethnic,educ,incomeW4,transition,age_transition,W4.depression,W4.general_health)
#filter complete cases
addhealth<-addhealth%>%
filter(complete.cases(.))library(survminer)## Loading required package: ggplot2## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang## Loading required package: ggpubr## Loading required package: magrittrlibrary(survey)## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartlibrary(survival)fit.d<-survfit(Surv(time = age_transition,event = transition)~1,data=addhealth)
summary(fit.d)## Call: survfit(formula = Surv(time = age_transition, event = transition) ~
## 1, data = addhealth)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 8676 40 0.995 0.000727 0.994 0.997
## 20 8636 132 0.980 0.001497 0.977 0.983
## 21 8504 121 0.966 0.001939 0.962 0.970
## 22 8383 205 0.943 0.002497 0.938 0.948
## 23 8178 205 0.919 0.002930 0.913 0.925
## 24 7973 159 0.901 0.003212 0.894 0.907
## 25 7812 46 0.895 0.003286 0.889 0.902
## 26 7710 14 0.894 0.003309 0.887 0.900
## 27 6780 1 0.894 0.003311 0.887 0.900ggsurvplot(fit.d,data = addhealth,risk.table = T,title="Survival function for violence transition",ylim=c(.85,1),xlim=c(17,29))
library(muhaz)
haz<-kphaz.fit(time=addhealth$age_transition,status = addhealth$transition,method = "product-limit")
haz## $time
## [1] 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5
##
## $haz
## [1] 0.015402872 0.014330795 0.024758224 0.025386788 0.020145380 0.005925837
## [7] 0.001916821 0.000150184
##
## $var
## [1] 1.797372e-06 1.697316e-06 2.990249e-06 3.144018e-06 2.552518e-06
## [6] 7.633933e-07 2.626961e-07 2.255523e-08kphaz.plot(haz,main="Hazard Function Plot")
ggsurvplot(fit.d,data = addhealth,risk.table = T,fun="cumhaz",title="Cumulative Hazard function for depression transition")
plot(cumsum(haz$haz)~haz$time,
main="Cumulative Hazard Function",
ylab="H(t)",xlab="Time in Years",
type=NULL,lwd=2,col=3)
data.frame(haz)## time haz var
## 1 19.5 0.015402872 1.797372e-06
## 2 20.5 0.014330795 1.697316e-06
## 3 21.5 0.024758224 2.990249e-06
## 4 22.5 0.025386788 3.144018e-06
## 5 23.5 0.020145380 2.552518e-06
## 6 24.5 0.005925837 7.633933e-07
## 7 25.5 0.001916821 2.626961e-07
## 8 26.5 0.000150184 2.255523e-08fit.kaplan<-survfit(Surv(age_transition,transition)~1,data=addhealth)
summary(fit.kaplan)## Call: survfit(formula = Surv(age_transition, transition) ~ 1, data = addhealth)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 8676 40 0.995 0.000727 0.994 0.997
## 20 8636 132 0.980 0.001497 0.977 0.983
## 21 8504 121 0.966 0.001939 0.962 0.970
## 22 8383 205 0.943 0.002497 0.938 0.948
## 23 8178 205 0.919 0.002930 0.913 0.925
## 24 7973 159 0.901 0.003212 0.894 0.907
## 25 7812 46 0.895 0.003286 0.889 0.902
## 26 7710 14 0.894 0.003309 0.887 0.900
## 27 6780 1 0.894 0.003311 0.887 0.900#based on sexual orientation
fit.kaplan.LGB<-survfit(Surv(age_transition,transition)~sexorient,data=addhealth)
summary(fit.kaplan.LGB)## Call: survfit(formula = Surv(age_transition, transition) ~ sexorient,
## data = addhealth)
##
## sexorient=a_straight
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 8401 38 0.995 0.000732 0.994 0.997
## 20 8363 126 0.980 0.001509 0.978 0.983
## 21 8237 113 0.967 0.001948 0.963 0.971
## 22 8124 194 0.944 0.002510 0.939 0.949
## 23 7930 199 0.920 0.002956 0.914 0.926
## 24 7731 149 0.903 0.003236 0.896 0.909
## 25 7580 43 0.897 0.003311 0.891 0.904
## 26 7482 14 0.896 0.003335 0.889 0.902
## 27 6577 1 0.896 0.003337 0.889 0.902
##
## sexorient=b_bisexual
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 143 1 0.993 0.00697 0.979 1.000
## 20 142 3 0.972 0.01379 0.945 0.999
## 21 139 4 0.944 0.01922 0.907 0.982
## 22 135 8 0.888 0.02636 0.838 0.941
## 23 127 2 0.874 0.02774 0.821 0.930
## 24 125 4 0.846 0.03017 0.789 0.907
## 25 121 3 0.825 0.03176 0.765 0.890
##
## sexorient=c_LGB
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 132 1 0.992 0.00755 0.978 1.000
## 20 131 3 0.970 0.01492 0.941 0.999
## 21 128 4 0.939 0.02077 0.900 0.981
## 22 124 3 0.917 0.02406 0.871 0.965
## 23 121 4 0.886 0.02762 0.834 0.942
## 24 117 6 0.841 0.03184 0.781 0.906ggsurvplot(fit.kaplan.LGB,conf.int = T,risk.table = F,title="Survivorship Function for depression Transition",xlab="Wave of Survey")
fit0 <- coxph(Surv(time=age_transition, event = transition) ~ factor(sexorient), data=addhealth,
iter=0, na.action=na.exclude)
o.minus.e <- tapply(resid(fit0), addhealth$sexorient, sum)
obs <- tapply(addhealth$transition, addhealth$sexorient, sum)
cbind(observed=obs, expected= obs- o.minus.e, "o-e"=o.minus.e) ## observed expected o-e
## a_straight 877 894.73079 -17.730788
## b_bisexual 25 14.64107 10.358935
## c_LGB 21 13.62815 7.371854#depression transition based on sexual orientation
fit<-survfit(Surv(time=age_transition, event=transition)~sexorient,addhealth)
summary(fit)## Call: survfit(formula = Surv(time = age_transition, event = transition) ~
## sexorient, data = addhealth)
##
## sexorient=a_straight
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 8401 38 0.995 0.000732 0.994 0.997
## 20 8363 126 0.980 0.001509 0.978 0.983
## 21 8237 113 0.967 0.001948 0.963 0.971
## 22 8124 194 0.944 0.002510 0.939 0.949
## 23 7930 199 0.920 0.002956 0.914 0.926
## 24 7731 149 0.903 0.003236 0.896 0.909
## 25 7580 43 0.897 0.003311 0.891 0.904
## 26 7482 14 0.896 0.003335 0.889 0.902
## 27 6577 1 0.896 0.003337 0.889 0.902
##
## sexorient=b_bisexual
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 143 1 0.993 0.00697 0.979 1.000
## 20 142 3 0.972 0.01379 0.945 0.999
## 21 139 4 0.944 0.01922 0.907 0.982
## 22 135 8 0.888 0.02636 0.838 0.941
## 23 127 2 0.874 0.02774 0.821 0.930
## 24 125 4 0.846 0.03017 0.789 0.907
## 25 121 3 0.825 0.03176 0.765 0.890
##
## sexorient=c_LGB
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 132 1 0.992 0.00755 0.978 1.000
## 20 131 3 0.970 0.01492 0.941 0.999
## 21 128 4 0.939 0.02077 0.900 0.981
## 22 124 3 0.917 0.02406 0.871 0.965
## 23 121 4 0.886 0.02762 0.834 0.942
## 24 117 6 0.841 0.03184 0.781 0.906fit00<-survfit(Surv(time = age_transition,event = transition)~sexorient+sex,addhealth)
summary(fit00)## Call: survfit(formula = Surv(time = age_transition, event = transition) ~
## sexorient + sex, data = addhealth)
##
## sexorient=a_straight, sex=a_male
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 3837 15 0.996 0.00101 0.994 0.998
## 20 3822 53 0.982 0.00213 0.978 0.986
## 21 3769 48 0.970 0.00276 0.964 0.975
## 22 3721 69 0.952 0.00346 0.945 0.959
## 23 3652 83 0.930 0.00411 0.922 0.938
## 24 3569 61 0.914 0.00452 0.905 0.923
## 25 3507 21 0.909 0.00465 0.900 0.918
## 26 3467 9 0.906 0.00470 0.897 0.916
## 27 3093 1 0.906 0.00471 0.897 0.915
##
## sexorient=a_straight, sex=b_female
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 4564 23 0.995 0.00105 0.993 0.997
## 20 4541 73 0.979 0.00212 0.975 0.983
## 21 4468 65 0.965 0.00273 0.959 0.970
## 22 4403 125 0.937 0.00359 0.930 0.944
## 23 4278 116 0.912 0.00420 0.904 0.920
## 24 4162 88 0.893 0.00458 0.884 0.902
## 25 4073 22 0.888 0.00467 0.879 0.897
## 26 4015 5 0.887 0.00469 0.878 0.896
##
## sexorient=b_bisexual, sex=a_male
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 23 1 0.957 0.0425 0.877 1
## 20 22 1 0.913 0.0588 0.805 1
##
## sexorient=b_bisexual, sex=b_female
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 20 120 2 0.983 0.0117 0.961 1.000
## 21 118 4 0.950 0.0199 0.912 0.990
## 22 114 8 0.883 0.0293 0.828 0.943
## 23 106 2 0.867 0.0310 0.808 0.930
## 24 104 4 0.833 0.0340 0.769 0.903
## 25 100 3 0.808 0.0359 0.741 0.882
##
## sexorient=c_LGB, sex=a_male
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 20 84 1 0.988 0.0118 0.965 1.000
## 21 83 3 0.952 0.0232 0.908 0.999
## 22 80 3 0.917 0.0302 0.859 0.978
## 23 77 2 0.893 0.0337 0.829 0.962
## 24 75 5 0.833 0.0407 0.757 0.917
##
## sexorient=c_LGB, sex=b_female
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 19 48 1 0.979 0.0206 0.940 1.000
## 20 47 2 0.938 0.0349 0.871 1.000
## 21 45 1 0.917 0.0399 0.842 0.998
## 23 44 2 0.875 0.0477 0.786 0.974
## 24 42 1 0.854 0.0509 0.760 0.960#Survey design
library(survey)
des2<-svydesign(ids=~psuscid,
strata = ~region,
weights=~GSWGT4,
data=addhealth,
nest=T)fit.cox<-svycoxph(Surv(time=age_transition,event=transition)~sex+sexorient,design=des2)
summary(fit.cox)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (132) clusters.
## svydesign(ids = ~psuscid, strata = ~region, weights = ~GSWGT4,
## data = addhealth, nest = T)## Call:
## svycoxph(formula = Surv(time = age_transition, event = transition) ~
## sex + sexorient, design = des2)
##
## n= 8676, number of events= 923
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sexb_female 0.1587 1.1720 0.0885 1.794 0.0729 .
## sexorientb_bisexual 0.3811 1.4639 0.2878 1.324 0.1855
## sexorientc_LGB 0.4848 1.6239 0.2470 1.963 0.0497 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexb_female 1.172 0.8532 0.9854 1.394
## sexorientb_bisexual 1.464 0.6831 0.8327 2.574
## sexorientc_LGB 1.624 0.6158 1.0007 2.635
##
## Concordance= 0.527 (se = 0.011 )
## Likelihood ratio test= NA on 3 df, p=NA
## Wald test = 8.37 on 3 df, p=0.04
## Score (logrank) test = NA on 3 df, p=NAplot(survfit(fit.cox,conf.int = F),ylab = "S(t)",xlab="Age",ylim=c(.96,1), xlim=c(17,28))
###model with demographic variables: sex, race/ethnicity, sexual orientation
#model with demographic variables: sex, race/ethnicity, sexual orientation
fit1<-svycoxph(Surv(time=age_transition,event=transition)~sex+racethnic+sexorient+educ+incomeW4+W4.general_health+W4.depression,design=des2)
summary(fit1)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (132) clusters.
## svydesign(ids = ~psuscid, strata = ~region, weights = ~GSWGT4,
## data = addhealth, nest = T)## Call:
## svycoxph(formula = Surv(time = age_transition, event = transition) ~
## sex + racethnic + sexorient + educ + incomeW4 + W4.general_health +
## W4.depression, design = des2)
##
## n= 8676, number of events= 923
##
## coef exp(coef) se(coef) z
## sexb_female 0.15561 1.16837 0.08963 1.736
## racethnicb-nhblack 0.30989 1.36328 0.09228 3.358
## racethnicc-hispanic 0.18010 1.19733 0.15772 1.142
## racethnicd-asian 0.12435 1.13241 0.19182 0.648
## racethnice-native_american 0.12003 1.12754 0.57424 0.209
## racethnicf-other -0.11522 0.89117 0.50978 -0.226
## sexorientb_bisexual 0.37548 1.45569 0.29108 1.290
## sexorientc_LGB 0.51331 1.67081 0.25439 2.018
## educb_highschool_grad -0.18095 0.83448 0.18181 -0.995
## educc_college_bach -0.34089 0.71114 0.21561 -1.581
## educd_college+ -0.26620 0.76628 0.22355 -1.191
## incomeW4b_$25k>$50k -0.11890 0.88789 0.12525 -0.949
## incomeW4c_$50k>$100k -0.20583 0.81397 0.11912 -1.728
## incomeW4e_$100k+ -0.07974 0.92335 0.15731 -0.507
## W4.general_healthc_good/vg/excellent 0.07057 1.07312 0.13813 0.511
## W4.depressionb_most_always 0.16695 1.18170 0.16513 1.011
## Pr(>|z|)
## sexb_female 0.082542 .
## racethnicb-nhblack 0.000784 ***
## racethnicc-hispanic 0.253500
## racethnicd-asian 0.516828
## racethnice-native_american 0.834424
## racethnicf-other 0.821184
## sexorientb_bisexual 0.197069
## sexorientc_LGB 0.043616 *
## educb_highschool_grad 0.319620
## educc_college_bach 0.113876
## educd_college+ 0.233736
## incomeW4b_$25k>$50k 0.342448
## incomeW4c_$50k>$100k 0.083995 .
## incomeW4e_$100k+ 0.612227
## W4.general_healthc_good/vg/excellent 0.609443
## W4.depressionb_most_always 0.312003
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## sexb_female 1.1684 0.8559 0.9801
## racethnicb-nhblack 1.3633 0.7335 1.1377
## racethnicc-hispanic 1.1973 0.8352 0.8790
## racethnicd-asian 1.1324 0.8831 0.7775
## racethnice-native_american 1.1275 0.8869 0.3659
## racethnicf-other 0.8912 1.1221 0.3281
## sexorientb_bisexual 1.4557 0.6870 0.8228
## sexorientc_LGB 1.6708 0.5985 1.0148
## educb_highschool_grad 0.8345 1.1984 0.5843
## educc_college_bach 0.7111 1.4062 0.4660
## educd_college+ 0.7663 1.3050 0.4944
## incomeW4b_$25k>$50k 0.8879 1.1263 0.6946
## incomeW4c_$50k>$100k 0.8140 1.2285 0.6445
## incomeW4e_$100k+ 0.9234 1.0830 0.6784
## W4.general_healthc_good/vg/excellent 1.0731 0.9319 0.8186
## W4.depressionb_most_always 1.1817 0.8462 0.8550
## upper .95
## sexb_female 1.393
## racethnicb-nhblack 1.634
## racethnicc-hispanic 1.631
## racethnicd-asian 1.649
## racethnice-native_american 3.475
## racethnicf-other 2.420
## sexorientb_bisexual 2.575
## sexorientc_LGB 2.751
## educb_highschool_grad 1.192
## educc_college_bach 1.085
## educd_college+ 1.188
## incomeW4b_$25k>$50k 1.135
## incomeW4c_$50k>$100k 1.028
## incomeW4e_$100k+ 1.257
## W4.general_healthc_good/vg/excellent 1.407
## W4.depressionb_most_always 1.633
##
## Concordance= 0.559 (se = 0.013 )
## Likelihood ratio test= NA on 16 df, p=NA
## Wald test = 42.09 on 16 df, p=4e-04
## Score (logrank) test = NA on 16 df, p=NAschoenfeld residuals
schoenresid<-resid(fit1,type="schoenfeld")
fit.sr<-lm(schoenresid~des2$variables$age_transition[des2$variables$transition==1])
summary(fit.sr)## Response sexb_female :
##
## Call:
## lm(formula = sexb_female ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5987 -0.5917 0.4049 0.4075 0.4131
##
## Coefficients:
## Estimate
## (Intercept) 0.078754
## des2$variables$age_transition[des2$variables$transition == 1] -0.001099
## Std. Error
## (Intercept) 0.220988
## des2$variables$age_transition[des2$variables$transition == 1] 0.009912
## t value
## (Intercept) 0.356
## des2$variables$age_transition[des2$variables$transition == 1] -0.111
## Pr(>|t|)
## (Intercept) 0.722
## des2$variables$age_transition[des2$variables$transition == 1] 0.912
##
## Residual standard error: 0.4918 on 921 degrees of freedom
## Multiple R-squared: 1.336e-05, Adjusted R-squared: -0.001072
## F-statistic: 0.0123 on 1 and 921 DF, p-value: 0.9117
##
##
## Response racethnicb-nhblack :
##
## Call:
## lm(formula = `racethnicb-nhblack` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2446 -0.2328 -0.2293 -0.2227 0.7815
##
## Coefficients:
## Estimate
## (Intercept) -0.011235
## des2$variables$age_transition[des2$variables$transition == 1] 0.002357
## Std. Error
## (Intercept) 0.189491
## des2$variables$age_transition[des2$variables$transition == 1] 0.008499
## t value
## (Intercept) -0.059
## des2$variables$age_transition[des2$variables$transition == 1] 0.277
## Pr(>|t|)
## (Intercept) 0.953
## des2$variables$age_transition[des2$variables$transition == 1] 0.782
##
## Residual standard error: 0.4217 on 921 degrees of freedom
## Multiple R-squared: 8.35e-05, Adjusted R-squared: -0.001002
## F-statistic: 0.07691 on 1 and 921 DF, p-value: 0.7816
##
##
## Response racethnicc-hispanic :
##
## Call:
## lm(formula = `racethnicc-hispanic` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09946 -0.08883 -0.07855 -0.07316 0.94242
##
## Coefficients:
## Estimate
## (Intercept) 0.131979
## des2$variables$age_transition[des2$variables$transition == 1] -0.005222
## Std. Error
## (Intercept) 0.123583
## des2$variables$age_transition[des2$variables$transition == 1] 0.005543
## t value
## (Intercept) 1.068
## des2$variables$age_transition[des2$variables$transition == 1] -0.942
## Pr(>|t|)
## (Intercept) 0.286
## des2$variables$age_transition[des2$variables$transition == 1] 0.346
##
## Residual standard error: 0.275 on 921 degrees of freedom
## Multiple R-squared: 0.0009627, Adjusted R-squared: -0.000122
## F-statistic: 0.8875 on 1 and 921 DF, p-value: 0.3464
##
##
## Response racethnicd-asian :
##
## Call:
## lm(formula = `racethnicd-asian` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09134 -0.06581 -0.05937 -0.05249 0.96031
##
## Coefficients:
## Estimate
## (Intercept) -0.113619
## des2$variables$age_transition[des2$variables$transition == 1] 0.006429
## Std. Error
## (Intercept) 0.107282
## des2$variables$age_transition[des2$variables$transition == 1] 0.004812
## t value
## (Intercept) -1.059
## des2$variables$age_transition[des2$variables$transition == 1] 1.336
## Pr(>|t|)
## (Intercept) 0.290
## des2$variables$age_transition[des2$variables$transition == 1] 0.182
##
## Residual standard error: 0.2388 on 921 degrees of freedom
## Multiple R-squared: 0.001934, Adjusted R-squared: 0.0008506
## F-statistic: 1.785 on 1 and 921 DF, p-value: 0.1819
##
##
## Response racethnice-native_american :
##
## Call:
## lm(formula = `racethnice-native_american` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.01410 -0.00596 -0.00390 -0.00181 0.99825
##
## Coefficients:
## Estimate
## (Intercept) -0.046700
## des2$variables$age_transition[des2$variables$transition == 1] 0.002058
## Std. Error
## (Intercept) 0.029509
## des2$variables$age_transition[des2$variables$transition == 1] 0.001324
## t value
## (Intercept) -1.583
## des2$variables$age_transition[des2$variables$transition == 1] 1.555
## Pr(>|t|)
## (Intercept) 0.114
## des2$variables$age_transition[des2$variables$transition == 1] 0.120
##
## Residual standard error: 0.06567 on 921 degrees of freedom
## Multiple R-squared: 0.002619, Adjusted R-squared: 0.001536
## F-statistic: 2.419 on 1 and 921 DF, p-value: 0.1203
##
##
## Response racethnicf-other :
##
## Call:
## lm(formula = `racethnicf-other` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.00963 -0.00631 -0.00460 -0.00299 0.99525
##
## Coefficients:
## Estimate
## (Intercept) 0.034149
## des2$variables$age_transition[des2$variables$transition == 1] -0.001604
## Std. Error
## (Intercept) 0.029523
## des2$variables$age_transition[des2$variables$transition == 1] 0.001324
## t value
## (Intercept) 1.157
## des2$variables$age_transition[des2$variables$transition == 1] -1.211
## Pr(>|t|)
## (Intercept) 0.248
## des2$variables$age_transition[des2$variables$transition == 1] 0.226
##
## Residual standard error: 0.06571 on 921 degrees of freedom
## Multiple R-squared: 0.00159, Adjusted R-squared: 0.0005056
## F-statistic: 1.466 on 1 and 921 DF, p-value: 0.2262
##
##
## Response sexorientb_bisexual :
##
## Call:
## lm(formula = sexorientb_bisexual ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07270 -0.03516 -0.02562 -0.01489 0.99297
##
## Coefficients:
## Estimate
## (Intercept) -0.205997
## des2$variables$age_transition[des2$variables$transition == 1] 0.009387
## Std. Error
## (Intercept) 0.072690
## des2$variables$age_transition[des2$variables$transition == 1] 0.003260
## t value
## (Intercept) -2.834
## des2$variables$age_transition[des2$variables$transition == 1] 2.879
## Pr(>|t|)
## (Intercept) 0.00470 **
## des2$variables$age_transition[des2$variables$transition == 1] 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1618 on 921 degrees of freedom
## Multiple R-squared: 0.008919, Adjusted R-squared: 0.007843
## F-statistic: 8.289 on 1 and 921 DF, p-value: 0.004082
##
##
## Response sexorientc_LGB :
##
## Call:
## lm(formula = sexorientc_LGB ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04080 -0.02609 -0.02210 -0.01818 0.98965
##
## Coefficients:
## Estimate
## (Intercept) -0.082234
## des2$variables$age_transition[des2$variables$transition == 1] 0.003677
## Std. Error
## (Intercept) 0.067019
## des2$variables$age_transition[des2$variables$transition == 1] 0.003006
## t value
## (Intercept) -1.227
## des2$variables$age_transition[des2$variables$transition == 1] 1.223
## Pr(>|t|)
## (Intercept) 0.220
## des2$variables$age_transition[des2$variables$transition == 1] 0.222
##
## Residual standard error: 0.1492 on 921 degrees of freedom
## Multiple R-squared: 0.001622, Adjusted R-squared: 0.0005379
## F-statistic: 1.496 on 1 and 921 DF, p-value: 0.2216
##
##
## Response educb_highschool_grad :
##
## Call:
## lm(formula = educb_highschool_grad ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6031 -0.5827 0.4042 0.4161 0.4342
##
## Coefficients:
## Estimate
## (Intercept) 0.082915
## des2$variables$age_transition[des2$variables$transition == 1] -0.004518
## Std. Error
## (Intercept) 0.221521
## des2$variables$age_transition[des2$variables$transition == 1] 0.009936
## t value
## (Intercept) 0.374
## des2$variables$age_transition[des2$variables$transition == 1] -0.455
## Pr(>|t|)
## (Intercept) 0.708
## des2$variables$age_transition[des2$variables$transition == 1] 0.649
##
## Residual standard error: 0.493 on 921 degrees of freedom
## Multiple R-squared: 0.0002245, Adjusted R-squared: -0.000861
## F-statistic: 0.2068 on 1 and 921 DF, p-value: 0.6494
##
##
## Response educc_college_bach :
##
## Call:
## lm(formula = educc_college_bach ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2026 -0.1991 -0.1978 -0.1959 0.8053
##
## Coefficients:
## Estimate
## (Intercept) 0.0376748
## des2$variables$age_transition[des2$variables$transition == 1] -0.0008175
## Std. Error
## (Intercept) 0.1793433
## des2$variables$age_transition[des2$variables$transition == 1] 0.0080441
## t value
## (Intercept) 0.210
## des2$variables$age_transition[des2$variables$transition == 1] -0.102
## Pr(>|t|)
## (Intercept) 0.834
## des2$variables$age_transition[des2$variables$transition == 1] 0.919
##
## Residual standard error: 0.3991 on 921 degrees of freedom
## Multiple R-squared: 1.121e-05, Adjusted R-squared: -0.001075
## F-statistic: 0.01033 on 1 and 921 DF, p-value: 0.9191
##
##
## Response educd_college+ :
##
## Call:
## lm(formula = `educd_college+` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1667 -0.1344 -0.1258 -0.1097 0.8992
##
## Coefficients:
## Estimate
## (Intercept) -0.168477
## des2$variables$age_transition[des2$variables$transition == 1] 0.008238
## Std. Error
## (Intercept) 0.150081
## des2$variables$age_transition[des2$variables$transition == 1] 0.006732
## t value
## (Intercept) -1.123
## des2$variables$age_transition[des2$variables$transition == 1] 1.224
## Pr(>|t|)
## (Intercept) 0.262
## des2$variables$age_transition[des2$variables$transition == 1] 0.221
##
## Residual standard error: 0.334 on 921 degrees of freedom
## Multiple R-squared: 0.001624, Adjusted R-squared: 0.0005395
## F-statistic: 1.498 on 1 and 921 DF, p-value: 0.2213
##
##
## Response incomeW4b_$25k>$50k :
##
## Call:
## lm(formula = `incomeW4b_$25k>$50k` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3029 -0.2830 -0.2741 0.7108 0.7377
##
## Coefficients:
## Estimate
## (Intercept) -0.126101
## des2$variables$age_transition[des2$variables$transition == 1] 0.005041
## Std. Error
## (Intercept) 0.201836
## des2$variables$age_transition[des2$variables$transition == 1] 0.009053
## t value
## (Intercept) -0.625
## des2$variables$age_transition[des2$variables$transition == 1] 0.557
## Pr(>|t|)
## (Intercept) 0.532
## des2$variables$age_transition[des2$variables$transition == 1] 0.578
##
## Residual standard error: 0.4492 on 921 degrees of freedom
## Multiple R-squared: 0.0003365, Adjusted R-squared: -0.0007489
## F-statistic: 0.31 on 1 and 921 DF, p-value: 0.5778
##
##
## Response incomeW4c_$50k>$100k :
##
## Call:
## lm(formula = `incomeW4c_$50k>$100k` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3868 -0.3770 -0.3700 0.6228 0.6367
##
## Coefficients:
## Estimate
## (Intercept) -0.040941
## des2$variables$age_transition[des2$variables$transition == 1] 0.002697
## Std. Error
## (Intercept) 0.217723
## des2$variables$age_transition[des2$variables$transition == 1] 0.009766
## t value
## (Intercept) -0.188
## des2$variables$age_transition[des2$variables$transition == 1] 0.276
## Pr(>|t|)
## (Intercept) 0.851
## des2$variables$age_transition[des2$variables$transition == 1] 0.782
##
## Residual standard error: 0.4846 on 921 degrees of freedom
## Multiple R-squared: 8.281e-05, Adjusted R-squared: -0.001003
## F-statistic: 0.07627 on 1 and 921 DF, p-value: 0.7825
##
##
## Response incomeW4e_$100k+ :
##
## Call:
## lm(formula = `incomeW4e_$100k+` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1734 -0.1636 -0.1564 -0.1520 0.8557
##
## Coefficients:
## Estimate
## (Intercept) 0.102719
## des2$variables$age_transition[des2$variables$transition == 1] -0.003856
## Std. Error
## (Intercept) 0.164578
## des2$variables$age_transition[des2$variables$transition == 1] 0.007382
## t value
## (Intercept) 0.624
## des2$variables$age_transition[des2$variables$transition == 1] -0.522
## Pr(>|t|)
## (Intercept) 0.533
## des2$variables$age_transition[des2$variables$transition == 1] 0.602
##
## Residual standard error: 0.3663 on 921 degrees of freedom
## Multiple R-squared: 0.0002962, Adjusted R-squared: -0.0007893
## F-statistic: 0.2728 on 1 and 921 DF, p-value: 0.6016
##
##
## Response W4.general_healthc_good/vg/excellent :
##
## Call:
## lm(formula = `W4.general_healthc_good/vg/excellent` ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9029 0.0996 0.1016 0.1031 0.1074
##
## Coefficients:
## Estimate
## (Intercept) 0.017207
## des2$variables$age_transition[des2$variables$transition == 1] -0.001161
## Std. Error
## (Intercept) 0.136042
## des2$variables$age_transition[des2$variables$transition == 1] 0.006102
## t value
## (Intercept) 0.126
## des2$variables$age_transition[des2$variables$transition == 1] -0.190
## Pr(>|t|)
## (Intercept) 0.899
## des2$variables$age_transition[des2$variables$transition == 1] 0.849
##
## Residual standard error: 0.3028 on 921 degrees of freedom
## Multiple R-squared: 3.928e-05, Adjusted R-squared: -0.001046
## F-statistic: 0.03618 on 1 and 921 DF, p-value: 0.8492
##
##
## Response W4.depressionb_most_always :
##
## Call:
## lm(formula = W4.depressionb_most_always ~ des2$variables$age_transition[des2$variables$transition ==
## 1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14776 -0.09987 -0.08740 -0.06416 0.93771
##
## Coefficients:
## Estimate
## (Intercept) -0.258588
## des2$variables$age_transition[des2$variables$transition == 1] 0.011943
## Std. Error
## (Intercept) 0.128400
## des2$variables$age_transition[des2$variables$transition == 1] 0.005759
## t value
## (Intercept) -2.014
## des2$variables$age_transition[des2$variables$transition == 1] 2.074
## Pr(>|t|)
## (Intercept) 0.0443 *
## des2$variables$age_transition[des2$variables$transition == 1] 0.0384 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2858 on 921 degrees of freedom
## Multiple R-squared: 0.004647, Adjusted R-squared: 0.003567
## F-statistic: 4.3 on 1 and 921 DF, p-value: 0.03838Grambsch and Therneau test
#Test the fit of model. Grambsch and Therneau formal test using weighted residuals
fit.test<-cox.zph(fit1)
fit.test## rho chisq p
## sexb_female -0.017728 6.22e-01 4.30e-01
## racethnicb-nhblack -0.039982 2.38e+00 1.23e-01
## racethnicc-hispanic 0.030570 1.62e+00 2.04e-01
## racethnicd-asian 0.032832 1.96e+00 1.61e-01
## racethnice-native_american 0.000332 1.82e-04 9.89e-01
## racethnicf-other 0.071279 6.45e+00 1.11e-02
## sexorientb_bisexual -0.022284 1.15e+00 2.84e-01
## sexorientc_LGB -0.002577 1.02e-02 9.20e-01
## educb_highschool_grad -0.086827 1.80e+01 2.24e-05
## educc_college_bach -0.081685 1.56e+01 7.97e-05
## educd_college+ -0.082949 1.53e+01 9.06e-05
## incomeW4b_$25k>$50k 0.054263 5.32e+00 2.10e-02
## incomeW4c_$50k>$100k 0.076500 9.17e+00 2.46e-03
## incomeW4e_$100k+ 0.053468 5.33e+00 2.10e-02
## W4.general_healthc_good/vg/excellent -0.038349 2.63e+00 1.05e-01
## W4.depressionb_most_always 0.055881 6.31e+00 1.20e-02
## GLOBAL NA 4.19e+01 4.12e-04Plot of residuals
#plot of residuals
par(mfrow=c(3,3))
plot(fit.test,df=2)
par(mfrow=c(1,1))
###Martingale residuals
#martingale residual
res.mar<-resid(fit1,type = "martingale")
scatter.smooth(des2$variables$sexorient,res.mar,degree = 2,span = 1)## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 4.2626e-16## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 4.2626e-16## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401
scatter.smooth(des2$variables$sexorient,res.mar,degree=2,span=1,ylab="Martingale Residual",col=1,cex=.5,lpars = list(col="red",lwd=3))## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 4.2626e-16## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : pseudoinverse used at 0.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : reciprocal condition number 4.2626e-16## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## FALSE, : There are other near singularities as well. 4.0401title(main="Martingale residuals for sexual orientation")
stratify by sexual orientation
#stratify by sexual orientation
fit4.strat<-svycoxph(Surv(time=age_transition,event=transition)~strata(sexorient)+sex+racethnic+educ+incomeW4,design=des2)
summary(fit4.strat)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (132) clusters.
## svydesign(ids = ~psuscid, strata = ~region, weights = ~GSWGT4,
## data = addhealth, nest = T)## Call:
## svycoxph(formula = Surv(time = age_transition, event = transition) ~
## strata(sexorient) + sex + racethnic + educ + incomeW4, design = des2)
##
## n= 8676, number of events= 923
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sexb_female 0.15856 1.17182 0.08917 1.778 0.07539 .
## racethnicb-nhblack 0.31117 1.36503 0.09114 3.414 0.00064 ***
## racethnicc-hispanic 0.17723 1.19390 0.15624 1.134 0.25665
## racethnicd-asian 0.12542 1.13362 0.19142 0.655 0.51235
## racethnice-native_american 0.12514 1.13330 0.56571 0.221 0.82493
## racethnicf-other -0.11859 0.88817 0.50886 -0.233 0.81573
## educb_highschool_grad -0.18806 0.82856 0.18206 -1.033 0.30160
## educc_college_bach -0.34916 0.70528 0.21486 -1.625 0.10415
## educd_college+ -0.27496 0.75960 0.22392 -1.228 0.21947
## incomeW4b_$25k>$50k -0.13126 0.87699 0.12223 -1.074 0.28287
## incomeW4c_$50k>$100k -0.22026 0.80231 0.11553 -1.907 0.05658 .
## incomeW4e_$100k+ -0.09170 0.91238 0.15259 -0.601 0.54788
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sexb_female 1.1718 0.8534 0.9839 1.396
## racethnicb-nhblack 1.3650 0.7326 1.1417 1.632
## racethnicc-hispanic 1.1939 0.8376 0.8790 1.622
## racethnicd-asian 1.1336 0.8821 0.7790 1.650
## racethnice-native_american 1.1333 0.8824 0.3740 3.435
## racethnicf-other 0.8882 1.1259 0.3276 2.408
## educb_highschool_grad 0.8286 1.2069 0.5799 1.184
## educc_college_bach 0.7053 1.4179 0.4629 1.075
## educd_college+ 0.7596 1.3165 0.4898 1.178
## incomeW4b_$25k>$50k 0.8770 1.1403 0.6902 1.114
## incomeW4c_$50k>$100k 0.8023 1.2464 0.6397 1.006
## incomeW4e_$100k+ 0.9124 1.0960 0.6765 1.230
##
## Concordance= 0.551 (se = 0.013 )
## Likelihood ratio test= NA on 12 df, p=NA
## Wald test = 33.89 on 12 df, p=7e-04
## Score (logrank) test = NA on 12 df, p=NAstratify by sex
#stratify by sex
fit4.strat<-svycoxph(Surv(time=age_transition,event=transition)~strata(sex)+racethnic+sexorient+educ+incomeW4,design=des2)
summary(fit4.strat)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (132) clusters.
## svydesign(ids = ~psuscid, strata = ~region, weights = ~GSWGT4,
## data = addhealth, nest = T)## Call:
## svycoxph(formula = Surv(time = age_transition, event = transition) ~
## strata(sex) + racethnic + sexorient + educ + incomeW4, design = des2)
##
## n= 8676, number of events= 923
##
## coef exp(coef) se(coef) z Pr(>|z|)
## racethnicb-nhblack 0.31123 1.36510 0.09123 3.411 0.000646 ***
## racethnicc-hispanic 0.17718 1.19384 0.15635 1.133 0.257140
## racethnicd-asian 0.12554 1.13376 0.19144 0.656 0.511968
## racethnice-native_american 0.12461 1.13271 0.56518 0.220 0.825494
## racethnicf-other -0.11794 0.88875 0.50927 -0.232 0.816863
## sexorientb_bisexual 0.37690 1.45775 0.29135 1.294 0.195790
## sexorientc_LGB 0.50852 1.66283 0.25382 2.003 0.045131 *
## educb_highschool_grad -0.18936 0.82749 0.18240 -1.038 0.299186
## educc_college_bach -0.35023 0.70453 0.21508 -1.628 0.103447
## educd_college+ -0.27624 0.75863 0.22421 -1.232 0.217926
## incomeW4b_$25k>$50k -0.13171 0.87660 0.12220 -1.078 0.281120
## incomeW4c_$50k>$100k -0.22047 0.80214 0.11544 -1.910 0.056167 .
## incomeW4e_$100k+ -0.09225 0.91188 0.15250 -0.605 0.545230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## racethnicb-nhblack 1.3651 0.7325 1.1416 1.632
## racethnicc-hispanic 1.1938 0.8376 0.8787 1.622
## racethnicd-asian 1.1338 0.8820 0.7791 1.650
## racethnice-native_american 1.1327 0.8828 0.3741 3.429
## racethnicf-other 0.8888 1.1252 0.3276 2.411
## sexorientb_bisexual 1.4578 0.6860 0.8235 2.580
## sexorientc_LGB 1.6628 0.6014 1.0111 2.735
## educb_highschool_grad 0.8275 1.2085 0.5788 1.183
## educc_college_bach 0.7045 1.4194 0.4622 1.074
## educd_college+ 0.7586 1.3182 0.4889 1.177
## incomeW4b_$25k>$50k 0.8766 1.1408 0.6899 1.114
## incomeW4c_$50k>$100k 0.8021 1.2467 0.6397 1.006
## incomeW4e_$100k+ 0.9119 1.0966 0.6763 1.230
##
## Concordance= 0.55 (se = 0.013 )
## Likelihood ratio test= NA on 13 df, p=NA
## Wald test = 36.33 on 13 df, p=5e-04
## Score (logrank) test = NA on 13 df, p=NA