DEM 7223 Homework 1 Introduction to Event History Analysis and Functions of Survival Time Add Health Suicide Ideation Ricardo Lowe
PART I
Using Public Use Add Health data files from 1995 to 2008, I aim to investigate age-specific exposure to suicide ideation (event variable).
First, I load my libraries:
library(haven)
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, unionlibrary(survminer)## Loading required package: ggplot2## Loading required package: ggpubrlibrary(dplyr)
library(muhaz)
library(haven)
library(survival)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeThen, I pull in the data files:
#ADD HEALTH 1995 (WI), 96 (WII), 01 (WIII) and 08 (WIV)
add95 <- read_sav("~/Desktop/ICPSR_21600/DS0001/21600-0001-Data.sav")
add96 <- read_sav("~/Desktop/ICPSR_21600-2/DS0005/21600-0005-Data.sav")
add01 <- read_sav("~/Desktop/ICPSR_21600 - 3/DS0008/21600-0008-Data.sav")
add08 <- read_sav("~/Desktop/ICPSR_21600-040/DS0022/21600-0022-Data.sav")I join the files
#Join 96 and 95 since the former does not have Race variable
add96_2 <- left_join(add95, add96, by = c('AID'))The clean it up a bit so that select variables can be manipulated.
#Year in 95 is only two-digit, so making a 4 digit variable
add95$DIG <- ifelse(add95$AID>0,19," ")
add95$AGE_ <- as.character(add95$H1GI1Y)
add95$AGE1 <- as.numeric(paste(add95$DIG,add95$AGE_,sep = ""))## Warning: NAs introduced by coercion#TEST
a9N5 <- add95 %>% select(AID,DIG,AGE1)I rename variables
#Rename Variables in Datasets to prep for Stack
add95 <- rename(add95, YEAR = IYEAR, BIRTHYR = AGE1, RACE = H1GI9, SEX = BIO_SEX, SUI = H1SU1, SUA = H1SU2)
add96 <- rename(add96_2, YEAR = IYEAR2, BIRTHYR = H2GI1Y, RACE = H1GI9, SEX = BIO_SEX2, SUI = H2SU1, SUA = H2SU2)
add01 <- rename(add01, YEAR = IYEAR3, BIRTHYR = H3OD1Y, RACE = H3IR4, SEX = BIO_SEX3, SUI = H3TO130, SUA = H3TO131)
add08 <- rename(add08, YEAR = IYEAR4, BIRTHYR = H4OD1Y, RACE = H4IR4, SEX = BIO_SEX4, SUI = H4SE1, SUA = H4SE2)Year for 1995 and 1996 were intially two digits, so I corrected them
add95$YEAR <- ifelse(add95$YEAR == 95,1995,"")
add96$YEAR <- ifelse(add96$YEAR == 96,1996,"")Now i select the variables I need fro each dataset
#-----------------id---year---race-sex-ideation-------attempt
a95 <- add95 %>% select(AID,YEAR,BIRTHYR,RACE,SEX,SUA,SUI)
a96 <- add96 %>% select(AID,YEAR,BIRTHYR,RACE,SEX,SUA,SUI)
a01 <- add01 %>% select(AID,YEAR,BIRTHYR,RACE,SEX,SUA,SUI)
a08 <- add08 %>% select(AID,YEAR,BIRTHYR,RACE,SEX,SUA,SUI)And stack all the datasets into one. Then rearrage columns.
sui <- rbind(a95,a01,a08)
X <- sui %>%
arrange(desc(AID))
X$AGE <- as.numeric(X$YEAR) - as.numeric(X$BIRTHYR)
X$SUA <- NULL
EV <- X[, c(1, 2, 3, 7, 4,5,6)]
EV$PERIOD <- 2008 - as.numeric(EV$YEAR)
EV## # A tibble: 16,500 x 8
## AID YEAR BIRTHYR AGE RACE SEX SUI PERIOD
## <chr> <chr> <dbl> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl>
## 1 997199… 1995 1981 14 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 13
## 2 997199… 2008 1981 27 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 0
## 3 997199… 1995 1980 15 1 [(1) White] 2 [(2) F… 0 [(0) No (sk… 13
## 4 997199… 2001 1980 21 1 [(1) White] 2 [(2) F… 1 [(1) Yes] 7
## 5 997199… 2008 1980 28 1 [(1) White] 2 [(2) F… 0 [(0) No (sk… 0
## 6 997199… 1995 1981 14 1 [(1) White] 1 [(1) M… 1 [(1) Yes] 13
## 7 997199… 2008 1981 27 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 0
## 8 997199… 1995 1981 14 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 13
## 9 997199… 2001 1981 20 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 7
## 10 997199… 2008 1981 27 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 0
## # … with 16,490 more rowsNow I order by AID (ID number)
EV2 <- EV %>%
arrange(desc(AID))
EV2## # A tibble: 16,500 x 8
## AID YEAR BIRTHYR AGE RACE SEX SUI PERIOD
## <chr> <chr> <dbl> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl>
## 1 997199… 1995 1981 14 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 13
## 2 997199… 2008 1981 27 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 0
## 3 997199… 1995 1980 15 1 [(1) White] 2 [(2) F… 0 [(0) No (sk… 13
## 4 997199… 2001 1980 21 1 [(1) White] 2 [(2) F… 1 [(1) Yes] 7
## 5 997199… 2008 1980 28 1 [(1) White] 2 [(2) F… 0 [(0) No (sk… 0
## 6 997199… 1995 1981 14 1 [(1) White] 1 [(1) M… 1 [(1) Yes] 13
## 7 997199… 2008 1981 27 1 [(1) White] 1 [(1) M… 0 [(0) No (sk… 0
## 8 997199… 1995 1981 14 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 13
## 9 997199… 2001 1981 20 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 7
## 10 997199… 2008 1981 27 2 [(2) Black/Afr… 1 [(1) M… 0 [(0) No (sk… 0
## # … with 16,490 more rowsIn the Add Health data, suicide ideation is recorded as ‘yes’ or ‘no’ (I excluded other options like DK or Refusal). I created a variable that provides the age of the respondent at time of ideation. If the respondent did not contemplate suicide by age 33, then the variable is SUI == 0 and the age at ideation is censored.
If the age at death is censored, then the age at ideation (censored age at ideation) is the year of the interview - year of birth. If the respondent contemplated suicide by age of 33, then the age at ideation in years is the variable sui.age2. Here I code this:
EV2$sui.age2<-ifelse(EV2$SUI==1,EV2$AGE,"")
EV2$sui.age<-ifelse(EV2$SUI==0,
EV2$AGE,
EV2$sui.age2)
EV2$sui.age2<-as.numeric(EV2$sui.age2,na.rm=TRUE)
EV2$sui.age<-as.numeric(EV2$sui.age,na.rm=TRUE)
#censoring indicator for death by age 1, in months (12 months)
EV2$d.event<-ifelse(is.na(EV2$sui.age2)==T|EV2$sui.age2>33,0,1)
EV2$d.eventfac<-factor(EV2$d.event); levels(EV2$d.eventfac)<-c("No Suicide Ideation", "Suicide Ideation")
table(EV2$d.eventfac)##
## No Suicide Ideation Suicide Ideation
## 15056 1444We see 1444 suicide ideators among the 15056 respondents across all 4 waves. RACE
summary(EV2$AGE)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 12.00 17.00 21.00 21.89 27.00 34.00 4library(survival)
head(EV2[,c("sui.age","d.event")], n=20)## # A tibble: 20 x 2
## sui.age d.event
## <dbl> <dbl>
## 1 14 0
## 2 27 0
## 3 15 0
## 4 21 1
## 5 28 0
## 6 14 1
## 7 27 0
## 8 14 0
## 9 20 0
## 10 27 0
## 11 13 0
## 12 19 0
## 13 26 0
## 14 14 0
## 15 20 0
## 16 27 0
## 17 14 0
## 18 20 0
## 19 27 0
## 20 15 0head(Surv(EV2$sui.age, EV2$d.event), n=20)## [1] 14+ 27+ 15+ 21 28+ 14 27+ 14+ 20+ 27+ 13+ 19+ 26+ 14+ 20+ 27+ 14+ 20+ 27+
## [20] 15+mort<-survfit(Surv(sui.age, d.event)~1, data=EV2,conf.type="none")
plot(mort, ylim=c(.7,1), xlim=c(13,33), main="Survival Function for Suicide Ideation Among Respondents")
1000*(1-summary(mort)$surv[16])## [1] 123.6626summary(mort)## Call: survfit(formula = Surv(sui.age, d.event) ~ 1, data = EV2, conf.type = "none")
##
## 241 observations deleted due to missingness
## time n.risk n.event survival std.err
## 13 16251 49 0.997 0.000430
## 14 15677 95 0.991 0.000751
## 15 14783 142 0.981 0.001089
## 16 13712 150 0.971 0.001386
## 17 12565 185 0.956 0.001718
## 18 11402 144 0.944 0.001969
## 19 10259 85 0.936 0.002128
## 20 9513 62 0.930 0.002251
## 21 8840 51 0.925 0.002360
## 22 8046 67 0.917 0.002521
## 23 7232 34 0.913 0.002615
## 24 6401 30 0.909 0.002717
## 25 5636 22 0.905 0.002810
## 26 5211 40 0.898 0.002995
## 27 4633 47 0.889 0.003247
## 28 3895 56 0.876 0.003622
## 29 3007 56 0.860 0.004159
## 30 2112 62 0.835 0.005126
## 31 1227 48 0.802 0.006753
## 32 356 16 0.766 0.010917
## 33 53 3 0.723 0.026407This Kaplan-Meier estimate of the survival function shows/ the number of respondents at risk of ideation and the number actually contemplating suivide. We see the highest number of suicide ideations occurred between 14 and 18 years of age.
EV2 <- EV2 %>% filter(!is.na(sui.age))
library(muhaz)
hazard<-kphaz.fit(time=EV2$sui.age, status=EV2$d.event, method = "product-limit")
kphaz.plot(hazard, main="Hazard function plot")
data.frame(hazard)## time haz var
## 1 13.5 0.006250236 4.113146e-07
## 2 14.5 0.009981525 7.019656e-07
## 3 15.5 0.011469560 8.775582e-07
## 4 16.5 0.015430919 1.288091e-06
## 5 17.5 0.013346521 1.238277e-06
## 6 18.5 0.008581702 8.667703e-07
## 7 19.5 0.006750514 7.352982e-07
## 8 20.5 0.006078859 7.251578e-07
## 9 21.5 0.008724834 1.137240e-06
## 10 22.5 0.005046990 7.504069e-07
## 11 23.5 0.004968199 8.238959e-07
## 12 24.5 0.004074120 7.547632e-07
## 13 25.5 0.008163162 1.667444e-06
## 14 26.5 0.011098821 2.627641e-06
## 15 27.5 0.016484103 4.879146e-06
## 16 28.5 0.022318668 8.990932e-06
## 17 29.5 0.039678076 2.601705e-05
## 18 30.5 0.073649752 1.297388e-04
## 19 31.5 0.117409476 1.180331e-03
## 20 32.5 0.062944351 1.325387e-03The hazard is highest among respondents aged between 15 and 18. There is also a large jump in risk between ages 31 and 32. Now we have our S(t) and h(t) functions. We can derive the other functions of survival time from these but integrating (summing) and differentiating these functions.
#cumulative hazard
plot(cumsum(hazard$haz)~hazard$time,
main = "Cumulative Hazard function",
ylab="H(t)",xlab="Time in Years",
type="l",xlim=c(13,33), lwd=2,col=3)
#Survival function, I just store this in an object so I can use it
surv<-mort
#here is a cheap version of the pdf
ft<- -diff(mort$surv)
plot(ft, xlim=c(13.0,33.0),
type="s",
ylab="f(t)",xlab="Time in Years",
main="Probability Density Function")
#here is the cumulative distribution function
Ft<-cumsum(ft)
plot(Ft, xlim=c(13.0,33.0), type="s", ylab="F(t)",xlab="Time in Years", main="Cumulative Distribution Function")
PART II
We will now test for differences in survival by characteristics of the student. First we will examine whether the survival chances are the same for Black students, compared to white students.
Here, 0=White, and 1=Black.
EV2$race<-Recode(EV2$RACE, recodes ="1 = 0; 2=1; else=NA")
fit1<-survfit(Surv(sui.age2, d.event)~race, data=EV2)
fit1## Call: survfit(formula = Surv(sui.age2, d.event) ~ race, data = EV2)
##
## 14930 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## race=0 1039 1039 18 18 19
## race=1 290 290 18 18 20ggsurvplot(fit1, xlim=c(13,33), conf.int=T,
title="Survival Function for Suicide Ideation - White vs. Black Students", ylim=c(0, 1))
summary(fit1)## Call: survfit(formula = Surv(sui.age2, d.event) ~ race, data = EV2)
##
## 14930 observations deleted due to missingness
## race=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 13 1039 30 0.971126 0.005195 0.960997 0.98136
## 14 1009 67 0.906641 0.009026 0.889122 0.92451
## 15 942 100 0.810395 0.012161 0.786907 0.83458
## 16 842 111 0.703561 0.014168 0.676333 0.73189
## 17 731 131 0.577478 0.015324 0.548211 0.60831
## 18 600 97 0.484119 0.015504 0.454666 0.51548
## 19 503 63 0.423484 0.015329 0.394481 0.45462
## 20 440 45 0.380173 0.015060 0.351773 0.41087
## 21 395 38 0.343600 0.014733 0.315903 0.37372
## 22 357 48 0.297401 0.014181 0.270866 0.32654
## 23 309 25 0.273340 0.013826 0.247541 0.30183
## 24 284 24 0.250241 0.013438 0.225241 0.27801
## 25 260 14 0.236766 0.013188 0.212279 0.26408
## 26 246 31 0.206930 0.012568 0.183707 0.23309
## 27 215 35 0.173244 0.011741 0.151694 0.19785
## 28 180 45 0.129933 0.010431 0.111015 0.15207
## 29 135 43 0.088547 0.008813 0.072853 0.10762
## 30 92 40 0.050048 0.006765 0.038401 0.06523
## 31 52 39 0.012512 0.003448 0.007290 0.02147
## 32 13 12 0.000962 0.000962 0.000136 0.00683
## 33 1 1 0.000000 NaN NA NA
##
## race=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 13 290 14 0.9517 0.01259 0.92737 0.9767
## 14 276 20 0.8828 0.01889 0.84650 0.9206
## 15 256 25 0.7966 0.02364 0.75154 0.8443
## 16 231 22 0.7207 0.02635 0.67086 0.7742
## 17 209 35 0.6000 0.02877 0.54618 0.6591
## 18 174 31 0.4931 0.02936 0.43879 0.5541
## 19 143 13 0.4483 0.02920 0.39454 0.5093
## 20 130 14 0.4000 0.02877 0.34741 0.4606
## 21 116 11 0.3621 0.02822 0.31077 0.4218
## 22 105 14 0.3138 0.02725 0.26468 0.3720
## 23 91 7 0.2897 0.02664 0.24188 0.3469
## 24 84 6 0.2690 0.02604 0.22248 0.3252
## 25 78 7 0.2448 0.02525 0.20002 0.2997
## 26 71 9 0.2138 0.02407 0.17145 0.2666
## 27 62 11 0.1759 0.02236 0.13708 0.2256
## 28 51 11 0.1379 0.02025 0.10344 0.1839
## 29 40 9 0.1069 0.01814 0.07665 0.1491
## 30 31 16 0.0517 0.01301 0.03160 0.0847
## 31 15 9 0.0207 0.00836 0.00937 0.0457
## 32 6 4 0.0069 0.00486 0.00173 0.0274
## 33 2 2 0.0000 NaN NA NAsurvdiff(Surv(sui.age2, d.event)~race, data=EV2)## Call:
## survdiff(formula = Surv(sui.age2, d.event) ~ race, data = EV2)
##
## n=1329, 14930 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## race=0 1039 1039 1030 0.0822 0.448
## race=1 290 290 299 0.2829 0.448
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5In this case, we see no difference in survival status based on race. How about sex?
0=Male,1=Female
EV2$sex<-Recode(EV2$SEX, recodes ="1 = 0; 2=1; else=NA")
fit2<-survfit(Surv(sui.age2, d.event)~sex, data=EV2)
fit2## Call: survfit(formula = Surv(sui.age2, d.event) ~ sex, data = EV2)
##
## 14815 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## sex=0 570 570 19 18 19
## sex=1 874 874 18 18 18ggsurvplot(fit2, xlim=c(13,33), conf.int=T,
title="Survival Function for Suicide Ideation - Male vs. Female Students", ylim=c(0, 1))
summary(fit2)## Call: survfit(formula = Surv(sui.age2, d.event) ~ sex, data = EV2)
##
## 14815 observations deleted due to missingness
## sex=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 13 570 20 0.9649 0.00771 0.94992 0.9801
## 14 550 27 0.9175 0.01152 0.89524 0.9404
## 15 523 41 0.8456 0.01513 0.81647 0.8758
## 16 482 55 0.7491 0.01816 0.71437 0.7856
## 17 427 70 0.6263 0.02026 0.58783 0.6673
## 18 357 61 0.5193 0.02093 0.47986 0.5620
## 19 296 40 0.4491 0.02083 0.41009 0.4919
## 20 256 22 0.4105 0.02060 0.37206 0.4530
## 21 234 23 0.3702 0.02022 0.33258 0.4120
## 22 211 30 0.3175 0.01950 0.28154 0.3582
## 23 181 18 0.2860 0.01893 0.25117 0.3256
## 24 163 16 0.2579 0.01832 0.22437 0.2964
## 25 147 12 0.2368 0.01781 0.20439 0.2744
## 26 135 13 0.2140 0.01718 0.18288 0.2505
## 27 122 20 0.1789 0.01606 0.15009 0.2134
## 28 102 20 0.1439 0.01470 0.11775 0.1758
## 29 82 31 0.0895 0.01196 0.06886 0.1163
## 30 51 23 0.0491 0.00905 0.03423 0.0705
## 31 28 21 0.0123 0.00461 0.00588 0.0256
## 32 7 7 0.0000 NaN NA NA
##
## sex=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 13 874 29 0.96682 0.00606 0.95502 0.9788
## 14 845 68 0.88902 0.01063 0.86843 0.9101
## 15 777 101 0.77346 0.01416 0.74620 0.8017
## 16 676 95 0.66476 0.01597 0.63419 0.6968
## 17 581 115 0.53318 0.01688 0.50111 0.5673
## 18 466 83 0.43822 0.01678 0.40653 0.4724
## 19 383 45 0.38673 0.01647 0.35575 0.4204
## 20 338 40 0.34096 0.01603 0.31094 0.3739
## 21 298 28 0.30892 0.01563 0.27976 0.3411
## 22 270 37 0.26659 0.01496 0.23883 0.2976
## 23 233 16 0.24828 0.01461 0.22123 0.2786
## 24 217 14 0.23227 0.01428 0.20589 0.2620
## 25 203 10 0.22082 0.01403 0.19497 0.2501
## 26 193 27 0.18993 0.01327 0.16563 0.2178
## 27 166 27 0.15904 0.01237 0.13655 0.1852
## 28 139 36 0.11785 0.01091 0.09830 0.1413
## 29 103 25 0.08924 0.00964 0.07221 0.1103
## 30 78 39 0.04462 0.00698 0.03283 0.0606
## 31 39 27 0.01373 0.00394 0.00783 0.0241
## 32 12 9 0.00343 0.00198 0.00111 0.0106
## 33 3 3 0.00000 NaN NA NAsurvdiff(Surv(sui.age2, d.event)~sex, data=EV2)## Call:
## survdiff(formula = Surv(sui.age2, d.event) ~ sex, data = EV2)
##
## n=1444, 14815 observations deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=0 570 570 605 2.04 4.28
## sex=1 874 874 839 1.47 4.28
##
## Chisq= 4.3 on 1 degrees of freedom, p= 0.04Significant!
prop.table(table(EV2$d.event, EV2$sex), margin = 2)##
## 0 1
## 0 0.92500987 0.89905290
## 1 0.07499013 0.10094710chisq.test(table(EV2$d.event, EV2$sex))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(EV2$d.event, EV2$sex)
## X-squared = 33.379, df = 1, p-value = 7.583e-09quantile(fit2, probs=.05)## $quantile
## 5
## sex=0 14
## sex=1 14
##
## $lower
## 5
## sex=0 13
## sex=1 14
##
## $upper
## 5
## sex=0 14
## sex=1 14EV2 <- EV2 %>% filter(!is.na(sui.age2))
haz2<-kphaz.fit(EV2$sui.age2, EV2$d.event, EV2$sex)
haz2## $time
## [1] 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5
## [16] 28.5 29.5 30.5 31.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5
## [31] 24.5 25.5 26.5 27.5 28.5 29.5 30.5 31.5 32.5
##
## $haz
## [1] 0.05028991 0.08155608 0.12102658 0.17881883 0.18708800 0.14491839
## [7] 0.08967295 0.10323042 0.15296901 0.10444237 0.10298446 0.08485618
## [13] 0.10086010 0.17824704 0.21706235 0.47120690 0.59164214 1.33431390
## [19] 2.59285714 0.08384451 0.13915117 0.15132144 0.22035288 0.19591831
## [25] 0.12481545 0.12575406 0.09849773 0.14708983 0.07098311 0.06653272
## [31] 0.05038839 0.15028181 0.17693007 0.29849124 0.27647012 0.68677801
## [37] 1.15033236 1.26987734 1.83333333
##
## $var
## [1] 9.368917e-05 1.623190e-04 2.666421e-04 4.580208e-04 5.754771e-04
## [6] 5.259525e-04 3.657553e-04 4.637377e-04 7.815044e-04 6.065609e-04
## [11] 6.634462e-04 6.004053e-04 7.831797e-04 1.592800e-03 2.365044e-03
## [16] 7.295802e-03 1.566726e-02 9.805289e-02 1.511797e+00 1.034415e-04
## [21] 1.920229e-04 2.414937e-04 4.239319e-04 4.639386e-04 3.466471e-04
## [26] 3.958730e-04 3.467728e-04 5.857952e-04 3.150443e-04 3.163020e-04
## [31] 2.539521e-04 8.380408e-04 1.162440e-03 2.493333e-03 3.076921e-03
## [36] 1.257642e-02 5.464232e-02 2.038655e-01 1.361111e+00
##
## $strata
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1plot(y=haz2$haz[1:12], x=haz2$time[1:12], col=1, lty=1, type="s")
lines(y=haz2$haz[13:24], x=haz2$time[13:24], col=2, lty=1, type="s")
plot(y=haz2$haz[1:33], x=haz2$time[1:33], col=1, lty=1, type="s")
lines(y=haz2$haz[13:40], x=haz2$time[13:40], col=2, lty=1, type="s")