Hwk 3
Marissa Laurel-Wilson
10/5/2020
library (haven)df1<- read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0001/da27021p1.dta")df1w <- read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0025/da27021p25.dta")wv1 <- merge(df1, df1w, by="aid")wv3<-read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0003/da27021p3.dta")wv2 <- merge(wv1, wv3, by="aid")wave4a <- read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0012/da27021p12.dta")wave4wgt <- read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0028/da27021p28.dta")wave4b <- read_dta("/projects/add_health/data/11613344/ICPSR_27021/DS0017/da27021p17.dta")wave4 <- merge(wave4a, wave4b, by="aid")wave4totalwgt <- merge(wave4, wave4wgt, by="aid")addhealth <-merge(wave4totalwgt, wv2, by="aid")summary(wave4wgt)## aid psuscid region GSWGT4
## Length:15701 Length:15701 Min. :1.000 Min. : 26.55
## Class :character Class :character 1st Qu.:2.000 1st Qu.: 409.15
## Mode :character Mode :character Median :3.000 Median : 1641.44
## Mean :2.412 Mean : 1969.89
## 3rd Qu.:3.000 3rd Qu.: 2988.77
## Max. :4.000 Max. :16323.66
## NA's :358 NA's :6280
## GSWGT4_2
## Min. : 20.57
## 1st Qu.: 249.31
## Median : 1204.89
## Mean : 1484.29
## 3rd Qu.: 2238.11
## Max. :18342.36
## NA's :901summary(addhealth$psuscid.x)## Length Class Mode
## 11898 character characterlibrary(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, unionlibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerlibrary(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(questionr)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(foreign)
library(survey)
library(survival)
library(eha)summary(wv1$H1SU2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 7.000 7.000 6.137 7.000 8.000summary(addhealth$H4TO35)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.00 3.00 29.46 97.00 98.00summary(addhealth$alcohol)## Length Class Mode
## 0 NULL NULLComplete Case
addhealth2<-addhealth%>%
select(GSWGT4,region.x,psuscid.x,aid,iyear,BIO_SEX,H4ID5H,iyear,imonth,iday,H1FV6,H4DS18,H4SE1,H4SE2,IMONTH4,IYEAR4,H1SU1,H4OD1M, H4OD1Y, H1FS6, H1GI1M, H1GI1Y, H3OD1Y,H3OD1M, IYEAR3, IMONTH3, H3TO130, H1SU2, H3TO131, H4SE2, H4TO35, H1TO15)
addhealth2<-addhealth2%>%
filter(complete.cases(.))Suicidal Ideation
addhealth2$suicidewave1<-recode(addhealth2$H1SU1, recodes="0=0; 1=1; 6=NA; 8=NA; 9=NA")
is.numeric(addhealth2$suicidewave1) ## [1] TRUEtable(addhealth2$suicidewave1)##
## 0 1
## 6869 1222addhealth2$suicidewave3<-recode(addhealth2$H3TO130, recodes="0=0; 1=1; 6=NA; 8=NA; 9=NA")
is.numeric(addhealth2$suicidewave3) ## [1] TRUEtable(addhealth2$suicidewave3)##
## 0 1 7
## 7565 381 100addhealth2$alcohol<-recode(addhealth2$H4TO35, recodes="0=0; 1:6=1; else=NA")
is.numeric(addhealth2$alcohol) ## [1] TRUEtable(addhealth2$alcohol)##
## 0 1
## 757 5117Age for Wave 1 and Wave 4
addhealth2$birthmonthw1 <- factor(ifelse(addhealth2$H1GI1M=="1",01,
ifelse(addhealth2$H1GI1M=="2", 02,
ifelse(addhealth2$H1GI1M=="3", 03,
ifelse(addhealth2$H1GI1M=="4", 04,
ifelse(addhealth2$H1GI1M=="5", 05,
ifelse(addhealth2$H1GI1M=="6", 06,
ifelse(addhealth2$H1GI1M=="7", 07,
ifelse(addhealth2$H1GI1M=="8", 08,
ifelse(addhealth2$H1GI1M=="9r", 09,
ifelse(addhealth2$H1GI1M=="10", 10,
ifelse(addhealth2$H1GI1M=="11", 11,
ifelse(addhealth2$H1GI1M=="12",12,
ifelse(addhealth2$H1GI1M=="96", 6,NA))))))))))))))addhealth2$birthyearw1 <- factor(ifelse(addhealth2$H1GI1Y=="74",1974,
ifelse(addhealth2$H1GI1Y=="75",1975,
ifelse(addhealth2$H1GI1Y=="76",1976,
ifelse(addhealth2$H1GI1Y=="77",1977,
ifelse(addhealth2$H1GI1Y=="78",1978,
ifelse(addhealth2$H1GI1Y=="79",1979,
ifelse(addhealth2$H1GI1Y=="80",1980,
ifelse(addhealth2$H1GI1Y=="81",1981,
ifelse(addhealth2$H1GI1Y=="82",1982,
ifelse(addhealth2$H1GI1Y=="83",1983,1979)))))))))))addhealth2$birthdatew1 <- as.yearmon(paste(addhealth2$birthyearw1, addhealth2$birthmonthw1), "%Y %m")
addhealth2$iyearfix <- Recode(addhealth2$iyear,recodes="'95'='1995'")
addhealth2$interviewdatew1 <- as.yearmon(paste(addhealth2$iyearfix, addhealth2$imonth), "%Y %m")
addhealth2$age3<-(as.numeric(round(addhealth2$interviewdatew1 - addhealth2$birthdatew1)))
addhealth2$age<-Recode(addhealth2$age3, recode="12=12; 13=13; 14=14; 15=15; 16=16; 17=17; 18=18; 19=19; 20=20; else=NA")
table(addhealth2$age)##
## 12 13 14 15 16 17 18 19 20
## 16 469 944 1135 1770 1578 1218 270 30addhealth2$birthdatew3<-as.yearmon(paste(addhealth2$H3OD1Y,addhealth2$H3OD1M),"%Y %m")
addhealth2$interviewdatew3<-as.yearmon(paste(addhealth2$IYEAR3,addhealth2$IMONTH3 ),"%Y %m")
addhealth2$agew3<-(as.numeric(round(addhealth2$interviewdatew3 - addhealth2$birthdatew3)))table(addhealth2$agew3)##
## 18 19 20 21 22 23 24 25 26 27
## 11 251 867 1176 1755 1721 1685 538 127 14Wave 4
addhealth2$month<- Recode(addhealth2$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")addhealth2$year <- Recode(addhealth2$H4OD1Y,recodes="1974=1974;1975=1975;1976=1976;1977=1977;1978=1978;1979=1979;1980=1980;1981=1981;1982=1982;1983=1983")addhealth2$date <- as.yearmon(paste(addhealth2$year, addhealth2$month), "%Y %m")addhealth2$intmon <- Recode(addhealth2$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")
addhealth2$intyear <- Recode(addhealth2$IYEAR4,recodes="2007=2007;2008=2008;2009=2009")
addhealth2$intdate<-as.yearmon(paste(addhealth2$intyear,addhealth2$intmon),"%Y %m")addhealth2$agewv4<-(as.numeric(round(addhealth2$intdate - addhealth2$date )))
table(addhealth2$agewv4)##
## 24 25 26 27 28 29 30 31 32 33 34
## 4 34 605 934 1563 1655 2027 974 309 33 7table(addhealth2$age)##
## 12 13 14 15 16 17 18 19 20
## 16 469 944 1135 1770 1578 1218 270 30variables for having no suicide thoughts and then having suicide thoughts
addhealth2$yessuicidethoughtW4<-Recode(addhealth2$H4SE1,recodes="0='no thoughts';1='yes thoughts';else=NA")
addhealth2$nosuicidethoughtW1<-Recode(addhealth2$H1SU1,recodes="0='no thoughts';1='yes thoughts';else=NA")1. Event Variable
variables for having a no suicide thoughts to having suicide thoughts
addhealth2$suicideevent<-ifelse(addhealth2$nosuicidethoughtW1=='no thoughts' & addhealth2$yessuicidethoughtW4=='yes thoughts', 1, 0)
table((addhealth2$suicideevent))##
## 0 1
## 7685 419summary(addhealth2$suicideevent)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 0.0000 0.0517 0.0000 1.0000 412. Duration
addhealth2$suicideevent.age<-ifelse(addhealth2$suicideevent==1, (addhealth2$age+addhealth2$agewv4)/2,
addhealth2$agewv4)summary((addhealth2$age+addhealth2$agewv4)/2)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 18.0 21.5 22.5 22.5 23.5 26.5 715addhealth2$depressed<-Recode(addhealth2$H1FS6, recodes="0=0; 1:3=1; 6=NA; 8=NA")
table(addhealth2$depressed)##
## 0 1
## 4394 3739addhealth2$alcoholwv1<-Recode(addhealth2$H1TO15, recodes="7=0; 1:6=1; else=NA")
table(addhealth2$alcoholwv1)##
## 0 1
## 807 4091summary(addhealth2$alcoholwv1)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 1.000 1.000 0.835 1.000 1.000 3247Outcome Variable: Suicide thoughts using ADD health data, Wave 1 and Wave 4
Event - having suicide thoughts and/or attempt Duration - from wave 1 to wave 4 Censoring - if someone is depressed and uses alcohol Predictors - biological sex, age
Using Proportional Hazards model because measuring Hazard not time.
addhealthmini<-subset(addhealth2, is.na(H1SU1)==F&is.na(H3TO130)==F&is.na(age)==F&is.na(agew3)==F&is.na(agewv4)==F&H1SU1!=1)addhealth2$W1.suicide.attempt<-Recode(addhealth2$H1SU2, recodes="0='a.no'; 1:4='b.yes'; else=NA")
addhealth2$W3.suicide.thought<-ifelse(addhealth2$H3TO130==1& addhealth2$H1SU1==0, 1, ifelse(addhealth2$H4SE1==1& addhealth2$H3TO130==0, 1, 0))
addhealth2$transition.w1<-ifelse(addhealth2$W1.suicide.attempt=="a.no"&addhealth2$W3.suicide.thought==1, 1, 0)
table(addhealth2$W3.suicide.thought)##
## 0 1
## 7370 775addhealth2$W1.suicide.thought<-Recode(addhealth2$H1SU1, recodes="0='a.no'; 1:4='b.yes'; else=NA")
addhealth2$W3.suicide.thought<-ifelse(addhealth2$H3TO130==1& addhealth2$H1SU1==0, 1, ifelse(addhealth2$H4SE1==1& addhealth2$H3TO130==0, 1, 0))
addhealth2$transition.w1<-ifelse(addhealth2$W1.suicide.thought=="a.no"&addhealth2$W3.suicide.thought==1, 1, 0)
table(addhealth2$W3.suicide.thought)##
## 0 1
## 7370 775table(addhealth2$W3.suicide.thought)##
## 0 1
## 7370 775library(muhaz)library (eha)addhealthb<-addhealth2%>%
select(age,agew3, agewv4, depressed, BIO_SEX, W3.suicide.thought, H1SU1, H3TO130, H4SE1, alcoholwv1, GSWGT4, region.x, psuscid.x )
#filter complete cases
addhealthb<-addhealthb%>%
filter(complete.cases(.))adlong<-reshape(data.frame(addhealthb), idvar = 'aid', varying=list(c('age', 'agew3'), c('agew3','agewv4')),
v.names = c('age_enter', 'age_exit'),
times = 1:2, direction='long')
adlong<-adlong[order(adlong$aid, adlong$time),]
head(adlong)## depressed BIO_SEX W3.suicide.thought H1SU1 H3TO130 H4SE1 alcoholwv1
## 1.1 0 2 0 0 0 0 1
## 1.2 0 2 0 0 0 0 1
## 2.1 0 2 0 0 0 0 1
## 2.2 0 2 0 0 0 0 1
## 3.1 0 2 0 0 0 0 1
## 3.2 0 2 0 0 0 0 1
## GSWGT4 region.x psuscid.x time age_enter age_exit aid
## 1.1 443.6169 3 044 1 14 20 1
## 1.2 443.6169 3 044 2 20 26 1
## 2.1 443.6169 3 044 1 14 20 2
## 2.2 443.6169 3 044 2 20 26 2
## 3.1 443.6169 3 044 1 14 20 3
## 3.2 443.6169 3 044 2 20 26 3adlong$povtran <-NA
adlong$povtran[adlong$H1SU1==0&adlong$H3TO130==1&adlong$time==1]<-1
adlong$povtran[adlong$H3TO130==0&adlong$H4SE1==1&adlong$time==2]<-1adlong$povtran[adlong$H1SU1==0&adlong$H3TO130==0&adlong$time==1]<-0
adlong$povtran[adlong$H3TO130==0&adlong$H4SE1==0&adlong$time==2]<-0failed1<-which(is.na(adlong$povtran)==T)
adlong<-adlong[-failed1,]adlong$age1r<-round(adlong$age_enter, 0)
adlong$age2r<-round(adlong$age_exit,0)
head(adlong, n=10)## depressed BIO_SEX W3.suicide.thought H1SU1 H3TO130 H4SE1 alcoholwv1
## 1.1 0 2 0 0 0 0 1
## 1.2 0 2 0 0 0 0 1
## 2.1 0 2 0 0 0 0 1
## 2.2 0 2 0 0 0 0 1
## 3.1 0 2 0 0 0 0 1
## 3.2 0 2 0 0 0 0 1
## 4.1 1 1 0 0 0 0 1
## 4.2 1 1 0 0 0 0 1
## 5.1 1 1 0 0 0 0 1
## 5.2 1 1 0 0 0 0 1
## GSWGT4 region.x psuscid.x time age_enter age_exit aid povtran age1r
## 1.1 443.6169 3 044 1 14 20 1 0 14
## 1.2 443.6169 3 044 2 20 26 1 0 20
## 2.1 443.6169 3 044 1 14 20 2 0 14
## 2.2 443.6169 3 044 2 20 26 2 0 20
## 3.1 443.6169 3 044 1 14 20 3 0 14
## 3.2 443.6169 3 044 2 20 26 3 0 20
## 4.1 1845.0308 2 064 1 15 21 4 0 15
## 4.2 1845.0308 2 064 2 21 27 4 0 21
## 5.1 1845.0308 2 064 1 15 21 5 0 15
## 5.2 1845.0308 2 064 2 21 27 5 0 21
## age2r
## 1.1 20
## 1.2 26
## 2.1 20
## 2.2 26
## 3.1 20
## 3.2 26
## 4.1 21
## 4.2 27
## 5.1 21
## 5.2 27summary(addhealth$alcoholwv1)## Length Class Mode
## 0 NULL NULLlibrary(survminer)## Loading required package: ggplot2## Loading required package: ggpubradlong <-adlong%>% filter(complete.cases(age_enter, age_exit, povtran, depressed, alcoholwv1, age, agew3, agewv4, BIO_SEX, W3.suicide.thought, H1SU1, H3TO130, H4SE1,GSWGT4,region.x,psuscid.x))
library(survey)GSWGT4,region.x,psuscid.x
des2<-svydesign(ids=~psuscid.x,
strata = ~region.x,
weights=~GSWGT4,
data=adlong,
nest=T)fitcox<-coxreg(Surv(time = age_enter, time2 = age_exit, event = povtran)~1+depressed+alcoholwv1, data=adlong)summary(fitcox)## Call:
## coxreg(formula = Surv(time = age_enter, time2 = age_exit, event = povtran) ~
## 1 + depressed + alcoholwv1, data = adlong)
##
## Covariate Mean Coef Rel.Risk S.E. Wald p
## depressed 0.481 0.489 1.631 0.097 0.000
## alcoholwv1 0.834 0.020 1.020 0.132 0.882
##
## Events 457
## Total time at risk 49490
## Max. log. likelihood -3523.6
## LR test statistic 26.16
## Degrees of freedom 2
## Overall p-value 2.08926e-06fitcox2<-coxreg(Surv(time = time,event = povtran)~age_enter-1,data = adlong)
summary(fitcox2)## Call:
## coxreg(formula = Surv(time = time, event = povtran) ~ age_enter -
## 1, data = adlong)
##
## Covariate Mean Coef Rel.Risk S.E. Wald p
## age_enter 20.755 -0.235 0.790 0.022 0.000
##
## Events 457
## Total time at risk 11793
## Max. log. likelihood -3813.1
## LR test statistic 124.61
## Degrees of freedom 1
## Overall p-value 0plot(survfit(fitcox2), ylim = c(.6,1), xlim = c(6,9), main="Survivorship function for Cox Regression-Suicide Thoughts")
plot(survfit(fitcox), ylim = c(.6,1), xlim = c(6,9), main="Survivorship function for Cox Regression-Suicide Thoughts")
Include all main effects in model
library(survey)des2<-svydesign(ids=~psuscid.x,
strata = ~region.x,
weights=~GSWGT4,
data=adlong,
nest=T)fit1b<-svycoxph(Surv(time=time,event=povtran)~BIO_SEX+depressed+alcoholwv1+W3.suicide.thought,design=des2)## Warning in fitter(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 4 ; coefficient may be infinite.summary(fit1b)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (130) clusters.
## svydesign(ids = ~psuscid.x, strata = ~region.x, weights = ~GSWGT4,
## data = adlong, nest = T)## Call:
## svycoxph(formula = Surv(time = time, event = povtran) ~ BIO_SEX +
## depressed + alcoholwv1 + W3.suicide.thought, design = des2)
##
## n= 7668, number of events= 457
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## BIO_SEX 6.668e-02 1.069e+00 9.427e-02 1.095e-01 0.609 0.543
## depressed 1.064e-01 1.112e+00 9.574e-02 9.481e-02 1.122 0.262
## alcoholwv1 4.720e-02 1.048e+00 1.375e-01 1.019e-01 0.463 0.643
## W3.suicide.thought 2.420e+01 3.226e+10 1.132e+03 7.749e-02 312.245 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## BIO_SEX 1.069e+00 9.355e-01 8.625e-01 1.325e+00
## depressed 1.112e+00 8.991e-01 9.236e-01 1.339e+00
## alcoholwv1 1.048e+00 9.539e-01 8.585e-01 1.280e+00
## W3.suicide.thought 3.226e+10 3.100e-11 2.771e+10 3.755e+10
##
## Concordance= 0.985 (se = 0.002 )
## Likelihood ratio test= NA on 4 df, p=NA
## Wald test = 111624 on 4 df, p=<2e-16
## Score (logrank) test = NA on 4 df, p=NA
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).fit.interact1<-svycoxph(Surv(time=time,event=povtran)~BIO_SEX+depressed+alcoholwv1+W3.suicide.thought,
design=des2)## Warning in fitter(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 4 ; coefficient may be infinite.summary(fit.interact1)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (130) clusters.
## svydesign(ids = ~psuscid.x, strata = ~region.x, weights = ~GSWGT4,
## data = adlong, nest = T)## Call:
## svycoxph(formula = Surv(time = time, event = povtran) ~ BIO_SEX +
## depressed + alcoholwv1 + W3.suicide.thought, design = des2)
##
## n= 7668, number of events= 457
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## BIO_SEX 6.668e-02 1.069e+00 9.427e-02 1.095e-01 0.609 0.543
## depressed 1.064e-01 1.112e+00 9.574e-02 9.481e-02 1.122 0.262
## alcoholwv1 4.720e-02 1.048e+00 1.375e-01 1.019e-01 0.463 0.643
## W3.suicide.thought 2.420e+01 3.226e+10 1.132e+03 7.749e-02 312.245 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## BIO_SEX 1.069e+00 9.355e-01 8.625e-01 1.325e+00
## depressed 1.112e+00 8.991e-01 9.236e-01 1.339e+00
## alcoholwv1 1.048e+00 9.539e-01 8.585e-01 1.280e+00
## W3.suicide.thought 3.226e+10 3.100e-11 2.771e+10 3.755e+10
##
## Concordance= 0.985 (se = 0.002 )
## Likelihood ratio test= NA on 4 df, p=NA
## Wald test = 111624 on 4 df, p=<2e-16
## Score (logrank) test = NA on 4 df, p=NA
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).#hazard for alcohol
fit.kaplan<-survfit(Surv(time = time , event=povtran)~alcoholwv1,data=adlong)
summary(fit.kaplan)## Call: survfit(formula = Surv(time = time, event = povtran) ~ alcoholwv1,
## data = adlong)
##
## alcoholwv1=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1269 22 0.983 0.00366 0.976 0.990
## 2 660 46 0.914 0.01032 0.894 0.935
##
## alcoholwv1=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 6399 115 0.982 0.00166 0.979 0.985
## 2 3465 274 0.904 0.00475 0.895 0.914ggsurvplot(fit.kaplan,risk.table = T,fun="cumhaz",title="Cumulative Hazard Function for suicide thoughts by alcohol usage")## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
#hazard for biological sex
fit.kaplan.biosex<-survfit(Surv(time = time , event=povtran)~BIO_SEX,data=adlong)
summary(fit.kaplan.biosex)## Call: survfit(formula = Surv(time = time, event = povtran) ~ BIO_SEX,
## data = adlong)
##
## BIO_SEX=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 2871 49 0.983 0.00242 0.978 0.988
## 2 1496 116 0.907 0.00715 0.893 0.921
##
## BIO_SEX=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 4797 88 0.982 0.00194 0.978 0.985
## 2 2629 204 0.905 0.00542 0.895 0.916ggsurvplot(fit.kaplan.biosex,risk.table = T,fun="cumhaz",title="Cumulative Hazard Function for suicide thoughts by biological sex")## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
#hazard for depression
fit.kaplan.depressed<-survfit(Surv(time = time , event=povtran)~depressed,data=adlong)
summary(fit.kaplan.depressed)## Call: survfit(formula = Surv(time = time, event = povtran) ~ depressed,
## data = adlong)
##
## depressed=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 3984 51 0.987 0.00178 0.984 0.991
## 2 2053 123 0.928 0.00544 0.917 0.939
##
## depressed=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 3684 86 0.977 0.00249 0.972 0.982
## 2 2072 197 0.884 0.00668 0.871 0.897ggsurvplot(fit.kaplan.depressed,risk.table = T,fun="cumhaz",title="Cumulative Hazard Function for suicide thoughts by depression")## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.