This example will illustrate how to fit a multistate hazard model using the multinomial logit model. The outcome for the example is “type of non-parental child care” and whether a family changes their particular type of childcare between waves 1 and 5 of the data. The data from the ECLS-K.

#Load required libraries
library(survival)
library(car)
library(survey)
## Loading required package: grid
## Loading required package: Matrix
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
library(lattice)

Using Longitudinal Data

As in the other examples, I illustrate fitting these models to data that are longitudinal, instead of person-duration. In this example, we will examine how to fit the Cox model to a longitudinally collected data set.

First we load our data

load("~/Google Drive/dem7223/class2016/data/eclsk.Rdata")
names(eclsk)<-tolower(names(eclsk))
#get out only the variables I'm going to use for this example
myvars<-c( "childid","gender", "race", "r1_kage","r4age", "r5age", "r6age", "r7age","c1r4mtsc", "c4r4mtsc", "c5r4mtsc", "c6r4mtsc", "c7r4mtsc", "w1povrty","w3povrty", "w5povrty", "w8povrty","wkmomed","w3momed", "w5momed","p1hmemp", "p5hmemp","p1hdemp", "p5hdemp","p1hparnt", "p2hparnt", "p5hparnt", "s2_id", "c1_5fp0", "c15fpstr", "c15fppsu", "p1primpk", "p1agefrs","p1primnw", "p4primnw", "p5primnw")
eclsk<-eclsk[,myvars]


eclsk$age1<-ifelse(eclsk$r1_kage==-9, NA, eclsk$r1_kage/12)
eclsk$age2<-ifelse(eclsk$r4age==-9, NA, eclsk$r4age/12)
#for the later waves, the NCES group the ages into ranges of months, so 1= <105 months, 2=105 to 108 months. So, I fix the age at the midpoint of the interval they give, and make it into years by dividing by 12
eclsk$age3<-recode(eclsk$r5age,recodes="1=105; 2=107; 3=109; 4=112; 5=115; 6=117; -9=NA")/12

eclsk$pov1<-ifelse(eclsk$w1povrty==1,1,0)
eclsk$pov2<-ifelse(eclsk$w3povrty==1,1,0)
eclsk$pov3<-ifelse(eclsk$w5povrty==1,1,0)

eclsk$care1<-recode(eclsk$p1primnw, recodes="0='none'; 1:2='relative'; 3:4='nonrelative'; 5='centerbased'; 6:7='other/multi';9=NA; else=NA", as.factor.result = T, levels = c("none", "relative", "nonrelative", "centerbased", "other/multi") )
eclsk$care2<-recode(eclsk$p4primnw, recodes="0='none'; 1:2='relative'; 3:4='nonrelative'; 5='centerbased'; 6:7='other/multi';9=NA; else=NA", as.factor.result = T, levels = c("none", "relative", "nonrelative", "centerbased", "other/multi"))
eclsk$care3<-recode(eclsk$p5primnw, recodes="0='none'; 1:2='relative'; 3:4='nonrelative'; 5='centerbased'; 6:7='other/multi';9=NA; else=NA", as.factor.result = T, levels = c("none", "relative", "nonrelative", "centerbased", "other/multi"))

#mother's education-time varying > hs ==1
eclsk$momedu1<-recode(eclsk$wkmomed, recodes = "1:3=0; 4:9=1; else =NA") 
eclsk$momedu2<-recode(eclsk$w3momed, recodes = "1:3=0; 4:9=1; else =NA") 
eclsk$momedu3<-recode(eclsk$w5momed, recodes = "1:3=0; 4:9=1; else =NA") 

#mother's work status-time varying
eclsk$momwork1<-recode(eclsk$p1hmemp, recodes = "1='fulltime';2='parttime'; 3='unemployed'; 4='nilf';else =NA", as.factor.result = T) 
eclsk$momwork2<-recode(eclsk$p1hmemp, recodes = "1='fulltime';2='parttime'; 3='unemployed'; 4='nilf';else =NA", as.factor.result = T) 
eclsk$momwork3<-recode(eclsk$p5hmemp, recodes = "1='fulltime';2='parttime'; 3='unemployed'; 4='nilf';else =NA", as.factor.result = T) 

#both parents in hh - time varying
eclsk$bothpar1<-ifelse(eclsk$p1hmemp==-1|eclsk$p1hdemp==-1, 1, ifelse(eclsk$p1hmemp==-9|eclsk$p1hdemp==-9, NA, 0))
eclsk$bothpar2<-ifelse(eclsk$p1hmemp==-1|eclsk$p1hdemp==-1, 1, ifelse(eclsk$p1hmemp==-9|eclsk$p1hdemp==-9, NA, 0))
eclsk$bothpar3<-ifelse(eclsk$p5hmemp==-1|eclsk$p5hdemp==-1, 1, ifelse(eclsk$p5hmemp==-9|eclsk$p5hdemp==-9, NA, 0))

#HH type - parents - time varying
eclsk$parenttype1<-recode(eclsk$p1hparnt, recodes="c(1,2,3)='twoparent'; 4='biomomonly'; 5='biodadonly';6:9='othernonparent'; else=NA", as.factor.result = T)
eclsk$parenttype2<-recode(eclsk$p2hparnt, recodes="c(1,2,3)='twoparent'; 4='biomomonly'; 5='biodadonly';6:9='othernonparent'; else=NA", as.factor.result = T)
eclsk$parenttype3<-recode(eclsk$p5hparnt, recodes="c(1,2,3)='twoparent'; 4='biomomonly'; 5='biodadonly';6:9='othernonparent'; else=NA", as.factor.result = T)

#Recode race with white, non Hispanic as reference using dummy vars
eclsk$hisp<-recode (eclsk$race, recodes="3:4=1;-9=NA; else=0")
eclsk$black<-recode (eclsk$race, recodes="2=1;-9=NA; else=0")
eclsk$asian<-recode (eclsk$race, recodes="5=1;-9=NA; else=0")
eclsk$nahn<-recode (eclsk$race, recodes="6:7=1;-9=NA; else=0")
eclsk$other<-recode (eclsk$race, recodes="8=1;-9=NA; else=0")
eclsk$male<-recode(eclsk$gender, recodes="1=1; 2=0; -9=NA")

Now, I need to form the transition variable, this is my event variable, and in this case it will be 1 if a child enters poverty between the first wave of the data and the third grade wave, and 0 otherwise. NOTE I need to remove any children who are already in poverty age wave 1, because they are not at risk of experiencing this particular transition.

eclsk<-subset(eclsk, is.na(care1)==F&is.na(care2)==F&is.na(care3)==F&is.na(age1)==F&is.na(age2)==F&is.na(age3)==F&is.na(eclsk$c15fpstr)==F)
head(eclsk)
##     childid gender race r1_kage r4age r5age r6age r7age c1r4mtsc c4r4mtsc
## 2  0001002C      2    1   77.20 94.73     6     4     4   68.324   58.738
## 4  0001004C      1    1   65.93 83.93     2    NA    NA   53.064   59.853
## 6  0001006C      1    1   63.00 81.00     1    NA    NA   49.382   42.676
## 7  0001007C      1    1   63.60 81.90     2     2     2   43.836   37.160
## 9  0001009C      1    1   67.73 86.00     3    NA    NA   45.589   42.724
## 10 0001010C      2    1   70.63 88.63     4     3     3   53.324   52.256
##    c5r4mtsc c6r4mtsc c7r4mtsc w1povrty w3povrty w5povrty w8povrty wkmomed
## 2    60.390   62.585   66.481        2        2        2        2       3
## 4    56.419       NA       NA        2        2       NA       NA       6
## 6    55.857       NA       NA        2        2       NA       NA       8
## 7    40.946   44.830   44.122        2        2        2        2       6
## 9    37.497       NA       NA        2        2       NA       NA       5
## 10   54.986   54.759   50.096        2        2        2        2       7
##    w3momed w5momed p1hmemp p5hmemp p1hdemp p5hdemp p1hparnt p2hparnt
## 2        5       5       2       2       1       1        1        1
## 4        7      NA       1       2       1       1        1        1
## 6        8      NA       2       2       1       4        1        1
## 7        6       6       1       1       1       2        1        1
## 9        5      NA       1       3       1       1        2        2
## 10       7       7       1       1       1       1        1        1
##    p5hparnt s2_id c1_5fp0 c15fpstr c15fppsu p1primpk p1agefrs p1primnw
## 2         1  0001  352.35       63        1        6       38        0
## 4         1  0001  358.99       63        1        7        3        5
## 6         1  0001  358.99       63        1        6        7        0
## 7         1  0001  358.99       63        1        6        5        5
## 9         2  0001  441.52       63        1        6       18        5
## 10        1  0001  352.35       63        1        6        8        5
##    p4primnw p5primnw     age1     age2     age3 pov1 pov2 pov3       care1
## 2         4        1 6.433333 7.894167 9.750000    0    0    0        none
## 4         5        5 5.494167 6.994167 8.916667    0    0   NA centerbased
## 6         0        0 5.250000 6.750000 8.750000    0    0   NA        none
## 7         4        4 5.300000 6.825000 8.916667    0    0    0 centerbased
## 9         0        0 5.644167 7.166667 9.083333    0    0   NA centerbased
## 10        5        5 5.885833 7.385833 9.333333    0    0    0 centerbased
##          care2       care3 momedu1 momedu2 momedu3 momwork1 momwork2
## 2  nonrelative    relative       0       1       1 parttime parttime
## 4  centerbased centerbased       1       1      NA fulltime fulltime
## 6         none        none       1       1      NA parttime parttime
## 7  nonrelative nonrelative       1       1       1 fulltime fulltime
## 9         none        none       1       1      NA fulltime fulltime
## 10 centerbased centerbased       1       1       1 fulltime fulltime
##      momwork3 bothpar1 bothpar2 bothpar3 parenttype1 parenttype2
## 2    parttime        0        0        0   twoparent   twoparent
## 4    parttime        0        0        0   twoparent   twoparent
## 6    parttime        0        0        0   twoparent   twoparent
## 7    fulltime        0        0        0   twoparent   twoparent
## 9  unemployed        0        0        0   twoparent   twoparent
## 10   fulltime        0        0        0   twoparent   twoparent
##    parenttype3 hisp black asian nahn other male
## 2    twoparent    0     0     0    0     0    0
## 4    twoparent    0     0     0    0     0    1
## 6    twoparent    0     0     0    0     0    1
## 7    twoparent    0     0     0    0     0    1
## 9    twoparent    0     0     0    0     0    1
## 10   twoparent    0     0     0    0     0    0
names(eclsk)
##  [1] "childid"     "gender"      "race"        "r1_kage"     "r4age"      
##  [6] "r5age"       "r6age"       "r7age"       "c1r4mtsc"    "c4r4mtsc"   
## [11] "c5r4mtsc"    "c6r4mtsc"    "c7r4mtsc"    "w1povrty"    "w3povrty"   
## [16] "w5povrty"    "w8povrty"    "wkmomed"     "w3momed"     "w5momed"    
## [21] "p1hmemp"     "p5hmemp"     "p1hdemp"     "p5hdemp"     "p1hparnt"   
## [26] "p2hparnt"    "p5hparnt"    "s2_id"       "c1_5fp0"     "c15fpstr"   
## [31] "c15fppsu"    "p1primpk"    "p1agefrs"    "p1primnw"    "p4primnw"   
## [36] "p5primnw"    "age1"        "age2"        "age3"        "pov1"       
## [41] "pov2"        "pov3"        "care1"       "care2"       "care3"      
## [46] "momedu1"     "momedu2"     "momedu3"     "momwork1"    "momwork2"   
## [51] "momwork3"    "bothpar1"    "bothpar2"    "bothpar3"    "parenttype1"
## [56] "parenttype2" "parenttype3" "hisp"        "black"       "asian"      
## [61] "nahn"        "other"       "male"
e.long<-reshape(eclsk, idvar="childid", varying=list( momedu=c("momedu1", "momedu2"), momwork=c("momwork1","momwork2"), partype=c("parenttype1", "parenttype2")), times=1:2, direction="long" , v.names = c("momedu", "momwork", "partype"))
e.long<-e.long[order(e.long$childid, e.long$time),]
head(e.long)
##             childid gender race r1_kage r4age r5age r6age r7age c1r4mtsc
## 0001002C.1 0001002C      2    1   77.20 94.73     6     4     4   68.324
## 0001002C.2 0001002C      2    1   77.20 94.73     6     4     4   68.324
## 0001004C.1 0001004C      1    1   65.93 83.93     2    NA    NA   53.064
## 0001004C.2 0001004C      1    1   65.93 83.93     2    NA    NA   53.064
## 0001006C.1 0001006C      1    1   63.00 81.00     1    NA    NA   49.382
## 0001006C.2 0001006C      1    1   63.00 81.00     1    NA    NA   49.382
##            c4r4mtsc c5r4mtsc c6r4mtsc c7r4mtsc w1povrty w3povrty w5povrty
## 0001002C.1   58.738   60.390   62.585   66.481        2        2        2
## 0001002C.2   58.738   60.390   62.585   66.481        2        2        2
## 0001004C.1   59.853   56.419       NA       NA        2        2       NA
## 0001004C.2   59.853   56.419       NA       NA        2        2       NA
## 0001006C.1   42.676   55.857       NA       NA        2        2       NA
## 0001006C.2   42.676   55.857       NA       NA        2        2       NA
##            w8povrty wkmomed w3momed w5momed p1hmemp p5hmemp p1hdemp
## 0001002C.1        2       3       5       5       2       2       1
## 0001002C.2        2       3       5       5       2       2       1
## 0001004C.1       NA       6       7      NA       1       2       1
## 0001004C.2       NA       6       7      NA       1       2       1
## 0001006C.1       NA       8       8      NA       2       2       1
## 0001006C.2       NA       8       8      NA       2       2       1
##            p5hdemp p1hparnt p2hparnt p5hparnt s2_id c1_5fp0 c15fpstr
## 0001002C.1       1        1        1        1  0001  352.35       63
## 0001002C.2       1        1        1        1  0001  352.35       63
## 0001004C.1       1        1        1        1  0001  358.99       63
## 0001004C.2       1        1        1        1  0001  358.99       63
## 0001006C.1       4        1        1        1  0001  358.99       63
## 0001006C.2       4        1        1        1  0001  358.99       63
##            c15fppsu p1primpk p1agefrs p1primnw p4primnw p5primnw     age1
## 0001002C.1        1        6       38        0        4        1 6.433333
## 0001002C.2        1        6       38        0        4        1 6.433333
## 0001004C.1        1        7        3        5        5        5 5.494167
## 0001004C.2        1        7        3        5        5        5 5.494167
## 0001006C.1        1        6        7        0        0        0 5.250000
## 0001006C.2        1        6        7        0        0        0 5.250000
##                age2     age3 pov1 pov2 pov3       care1       care2
## 0001002C.1 7.894167 9.750000    0    0    0        none nonrelative
## 0001002C.2 7.894167 9.750000    0    0    0        none nonrelative
## 0001004C.1 6.994167 8.916667    0    0   NA centerbased centerbased
## 0001004C.2 6.994167 8.916667    0    0   NA centerbased centerbased
## 0001006C.1 6.750000 8.750000    0    0   NA        none        none
## 0001006C.2 6.750000 8.750000    0    0   NA        none        none
##                  care3 momedu3 momwork3 bothpar1 bothpar2 bothpar3
## 0001002C.1    relative       1 parttime        0        0        0
## 0001002C.2    relative       1 parttime        0        0        0
## 0001004C.1 centerbased      NA parttime        0        0        0
## 0001004C.2 centerbased      NA parttime        0        0        0
## 0001006C.1        none      NA parttime        0        0        0
## 0001006C.2        none      NA parttime        0        0        0
##            parenttype3 hisp black asian nahn other male time momedu
## 0001002C.1   twoparent    0     0     0    0     0    0    1      0
## 0001002C.2   twoparent    0     0     0    0     0    0    2      1
## 0001004C.1   twoparent    0     0     0    0     0    1    1      1
## 0001004C.2   twoparent    0     0     0    0     0    1    2      1
## 0001006C.1   twoparent    0     0     0    0     0    1    1      1
## 0001006C.2   twoparent    0     0     0    0     0    1    2      1
##             momwork   partype
## 0001002C.1 parttime twoparent
## 0001002C.2 parttime twoparent
## 0001004C.1 fulltime twoparent
## 0001004C.2 fulltime twoparent
## 0001006C.1 parttime twoparent
## 0001006C.2 parttime twoparent
#find which kids failed in the first time period and remove them from the second risk period risk set
e.long$povtran<-NA
e.long$povtran[e.long$pov1==0&e.long$pov2==1&e.long$time==1]<-1
e.long$povtran[e.long$pov2==0&e.long$pov3==1&e.long$time==2]<-1
e.long$povtran[e.long$pov1==0&e.long$pov2==0&e.long$time==1]<-0
e.long$povtran[e.long$pov2==0&e.long$pov3==0&e.long$time==2]<-0

e.long$povtran2<-NA
e.long$povtran2[e.long$pov1==1&e.long$pov2==1&e.long$time==1]<-0
e.long$povtran2[e.long$pov2==1&e.long$pov3==1&e.long$time==2]<-0
e.long$povtran2[e.long$pov1==1&e.long$pov2==0&e.long$time==1]<-1
e.long$povtran2[e.long$pov2==1&e.long$pov3==0&e.long$time==2]<-1

head(e.long[, c("childid", "pov1", "pov2", "pov3","time", "povtran", "povtran2")], n=20)
##             childid pov1 pov2 pov3 time povtran povtran2
## 0001002C.1 0001002C    0    0    0    1       0       NA
## 0001002C.2 0001002C    0    0    0    2       0       NA
## 0001004C.1 0001004C    0    0   NA    1       0       NA
## 0001004C.2 0001004C    0    0   NA    2      NA       NA
## 0001006C.1 0001006C    0    0   NA    1       0       NA
## 0001006C.2 0001006C    0    0   NA    2      NA       NA
## 0001007C.1 0001007C    0    0    0    1       0       NA
## 0001007C.2 0001007C    0    0    0    2       0       NA
## 0001009C.1 0001009C    0    0   NA    1       0       NA
## 0001009C.2 0001009C    0    0   NA    2      NA       NA
## 0001010C.1 0001010C    0    0    0    1       0       NA
## 0001010C.2 0001010C    0    0    0    2       0       NA
## 0001011C.1 0001011C    0    0   NA    1       0       NA
## 0001011C.2 0001011C    0    0   NA    2      NA       NA
## 0001012C.1 0001012C    0    0   NA    1       0       NA
## 0001012C.2 0001012C    0    0   NA    2      NA       NA
## 0001013C.1 0001013C    0    0   NA    1       0       NA
## 0001013C.2 0001013C    0    0   NA    2      NA       NA
## 0002003C.1 0002003C    1    1    0    1      NA        0
## 0002003C.2 0002003C    1    1    0    2      NA        1
table(e.long$time, e.long$povtran)
##    
##        0    1
##   1 7818  512
##   2 6200  327
table(e.long$time, e.long$povtran2)
##    
##        0    1
##   1 1069  385
##   2  822  370
#[1] "none"        "relative"    "nonrelative" "centerbased" "other/multi"
e.long$caretrn<-NA
e.long$caretrn[e.long$care1=="none"&e.long$care2=="none"&e.long$time==1]<-0
e.long$caretrn[e.long$care1=="none"&e.long$care2=="relative"&e.long$time==1]<-1
e.long$caretrn[e.long$care1=="none"&e.long$care2=="nonrelative"&e.long$time==1]<-2
e.long$caretrn[e.long$care1=="none"&e.long$care2=="centerbased"&e.long$time==1]<-3
e.long$caretrn[e.long$care1=="none"&e.long$care2=="other/multi"&e.long$time==1]<-4

e.long$caretrn[e.long$care1=="relative"&e.long$care2=="relative"&e.long$time==1]<-0
e.long$caretrn[e.long$care1=="relative"&e.long$care2=="none"&e.long$time==1]<-5
e.long$caretrn[e.long$care1=="relative"&e.long$care2=="nonrelative"&e.long$time==1]<-6
e.long$caretrn[e.long$care1=="relative"&e.long$care2=="centerbased"&e.long$time==1]<-7
e.long$caretrn[e.long$care1=="relative"&e.long$care2=="other/multi"&e.long$time==1]<-8

e.long$caretrn[e.long$care1=="nonrelative"&e.long$care2=="none"&e.long$time==1]<-9
e.long$caretrn[e.long$care1=="nonrelative"&e.long$care2=="relative"&e.long$time==1]<-10
e.long$caretrn[e.long$care1=="nonrelative"&e.long$care2=="nonrelative"&e.long$time==1]<-0
e.long$caretrn[e.long$care1=="nonrelative"&e.long$care2=="centerbased"&e.long$time==1]<-11
e.long$caretrn[e.long$care1=="nonrelative"&e.long$care2=="other/multi"&e.long$time==1]<-12

e.long$caretrn[e.long$care1=="centerbased"&e.long$care2=="none"&e.long$time==1]<-13
e.long$caretrn[e.long$care1=="centerbased"&e.long$care2=="relative"&e.long$time==1]<-14
e.long$caretrn[e.long$care1=="centerbased"&e.long$care2=="nonrelative"&e.long$time==1]<-15
e.long$caretrn[e.long$care1=="centerbased"&e.long$care2=="centerbased"&e.long$time==1]<-0
e.long$caretrn[e.long$care1=="centerbased"&e.long$care2=="other/multi"&e.long$time==1]<-16

e.long$caretrn[e.long$care1=="other/multi"&e.long$care2=="none"&e.long$time==1]<-17
e.long$caretrn[e.long$care1=="other/multi"&e.long$care2=="relative"&e.long$time==1]<-18
e.long$caretrn[e.long$care1=="other/multi"&e.long$care2=="nonrelative"&e.long$time==1]<-19
e.long$caretrn[e.long$care1=="other/multi"&e.long$care2=="centerbased"&e.long$time==1]<-20
e.long$caretrn[e.long$care1=="other/multi"&e.long$care2=="other/multi"&e.long$time==1]<-0

###Second transition####
e.long$caretrn[e.long$care2=="none"&e.long$care3=="none"&e.long$time==2]<-0
e.long$caretrn[e.long$care2=="none"&e.long$care3=="relative"&e.long$time==2]<-1
e.long$caretrn[e.long$care2=="none"&e.long$care3=="nonrelative"&e.long$time==2]<-2
e.long$caretrn[e.long$care2=="none"&e.long$care3=="centerbased"&e.long$time==2]<-3
e.long$caretrn[e.long$care2=="none"&e.long$care3=="other/multi"&e.long$time==2]<-4

e.long$caretrn[e.long$care2=="relative"&e.long$care3=="relative"&e.long$time==2]<-0
e.long$caretrn[e.long$care2=="relative"&e.long$care3=="none"&e.long$time==2]<-5
e.long$caretrn[e.long$care2=="relative"&e.long$care3=="nonrelative"&e.long$time==2]<-6
e.long$caretrn[e.long$care2=="relative"&e.long$care3=="centerbased"&e.long$time==2]<-7
e.long$caretrn[e.long$care2=="relative"&e.long$care3=="other/multi"&e.long$time==2]<-8

e.long$caretrn[e.long$care2=="nonrelative"&e.long$care3=="none"&e.long$time==2]<-9
e.long$caretrn[e.long$care2=="nonrelative"&e.long$care3=="relative"&e.long$time==2]<-10
e.long$caretrn[e.long$care2=="nonrelative"&e.long$care3=="nonrelative"&e.long$time==2]<-0
e.long$caretrn[e.long$care2=="nonrelative"&e.long$care3=="centerbased"&e.long$time==2]<-11
e.long$caretrn[e.long$care2=="nonrelative"&e.long$care3=="other/multi"&e.long$time==2]<-12

e.long$caretrn[e.long$care2=="centerbased"&e.long$care3=="none"&e.long$time==2]<-13
e.long$caretrn[e.long$care2=="centerbased"&e.long$care3=="relative"&e.long$time==2]<-14
e.long$caretrn[e.long$care2=="centerbased"&e.long$care3=="nonrelative"&e.long$time==2]<-15
e.long$caretrn[e.long$care2=="centerbased"&e.long$care3=="centerbased"&e.long$time==2]<-0
e.long$caretrn[e.long$care2=="centerbased"&e.long$care3=="other/multi"&e.long$time==2]<-16

e.long$caretrn[e.long$care2=="other/multi"&e.long$care3=="none"&e.long$time==2]<-17
e.long$caretrn[e.long$care2=="other/multi"&e.long$care3=="relative"&e.long$time==2]<-18
e.long$caretrn[e.long$care2=="other/multi"&e.long$care3=="nonrelative"&e.long$time==2]<-19
e.long$caretrn[e.long$care2=="other/multi"&e.long$care3=="centerbased"&e.long$time==2]<-20
e.long$caretrn[e.long$care2=="other/multi"&e.long$care3=="other/multi"&e.long$time==2]<-0

e.long$cens[e.long$time==1&e.long$caretrn!=0]<-1
e.long$cens[e.long$time==1&e.long$caretrn==0]<-0
e.long$cens[e.long$time==2&e.long$caretrn!=0]<-1
e.long$cens[e.long$time==2&e.long$caretrn==0]<-0

head(e.long[,c("childid","care1", "care2", "care3", "caretrn", "time")], n=10)
##             childid       care1       care2       care3 caretrn time
## 0001002C.1 0001002C        none nonrelative    relative       2    1
## 0001002C.2 0001002C        none nonrelative    relative      10    2
## 0001004C.1 0001004C centerbased centerbased centerbased       0    1
## 0001004C.2 0001004C centerbased centerbased centerbased       0    2
## 0001006C.1 0001006C        none        none        none       0    1
## 0001006C.2 0001006C        none        none        none       0    2
## 0001007C.1 0001007C centerbased nonrelative nonrelative      15    1
## 0001007C.2 0001007C centerbased nonrelative nonrelative       0    2
## 0001009C.1 0001009C centerbased        none        none      13    1
## 0001009C.2 0001009C centerbased        none        none       0    2
table(e.long$caretrn, e.long$time)
##     
##         1    2
##   0  6397 6467
##   1   458  484
##   2   180  140
##   3   271  295
##   4    51   43
##   5   556  662
##   6    69   56
##   7   120  144
##   8    65   45
##   9   377  302
##   10  132  119
##   11  117  114
##   12   23   19
##   13  480  482
##   14  147  136
##   15  125   86
##   16   22   30
##   17   58   59
##   18   67   61
##   19   25   15
##   20   44   25
prop.table(table(e.long$care1,e.long$care2), margin = 1)
##              
##                     none   relative nonrelative centerbased other/multi
##   none        0.81035164 0.09047807  0.03555907  0.05353615  0.01007507
##   relative    0.31990794 0.53394707  0.03970081  0.06904488  0.03739931
##   nonrelative 0.35432331 0.12406015  0.39003759  0.10996241  0.02161654
##   centerbased 0.28070175 0.08596491  0.07309942  0.54736842  0.01286550
##   other/multi 0.27619048 0.31904762  0.11904762  0.20952381  0.07619048
table(e.long$care2,e.long$care3, e.long$black)
## , ,  = 0
## 
##              
##               none relative nonrelative centerbased other/multi
##   none        8600      816         262         494          78
##   relative    1120     1344          94         216          78
##   nonrelative  562      218         498         220          38
##   centerbased  852      242         152        1292          54
##   other/multi  104       92          30          44          28
## 
## , ,  = 1
## 
##              
##               none relative nonrelative centerbased other/multi
##   none         620      152          18          94           8
##   relative     204      306          18          72          12
##   nonrelative   42       18          22           8           0
##   centerbased  112       30          20         214           6
##   other/multi   14       30           0           6           6
library(nnet)
multnm<-multinom(factor(caretrn)~time+momedu+momwork+black+hisp+asian+nahn+other, data=e.long, weights=c1_5fp0/mean(c1_5fp0, na.rm=T), maxit=500)
## # weights:  252 (220 variable)
## initial  value 57956.896790 
## iter  10 value 31232.630138
## iter  20 value 30495.459700
## iter  30 value 29318.212680
## iter  40 value 29026.711767
## iter  50 value 28673.030070
## iter  60 value 28514.110639
## iter  70 value 28482.975658
## iter  80 value 28474.263560
## iter  90 value 28469.995547
## iter 100 value 28468.721613
## iter 110 value 28468.351130
## iter 120 value 28468.228056
## iter 130 value 28468.182375
## iter 140 value 28468.160077
## iter 150 value 28468.131257
## iter 160 value 28468.114114
## iter 170 value 28468.105813
## final  value 28468.104999 
## converged
sums<-summary(multnm)
sums$z<-sums$coefficients/sums$standard.errors
sums$pval<-2*pnorm(abs(sums$z), lower.tail=FALSE)
round(exp(sums$coefficients),3)
##    (Intercept)  time momedu momworknilf momworkparttime momworkunemployed
## 1        0.067 1.064  0.698       0.647           1.218             1.631
## 2        0.031 0.745  1.307       0.942           1.279             1.652
## 3        0.034 1.099  1.075       0.745           1.129             1.321
## 4        0.017 0.767  0.424       0.382           1.259             1.738
## 5        0.124 1.072  0.691       0.242           0.904             0.759
## 6        0.041 0.696  0.720       0.119           0.366             0.162
## 7        0.017 1.129  1.404       0.255           0.697             0.272
## 8        0.022 0.651  0.787       0.174           0.362             0.408
## 9        0.086 0.857  1.409       0.254           0.904             0.410
## 10       0.050 0.884  0.851       0.076           0.469             0.061
## 11       0.031 0.947  1.410       0.075           0.474             0.207
## 12       0.005 0.794  1.388       0.153           0.702             2.613
## 13       0.069 1.143  1.634       0.247           0.759             0.498
## 14       0.047 0.721  1.449       0.088           0.400             0.593
## 15       0.051 0.807  1.233       0.111           0.241             0.054
## 16       0.004 1.300  1.443       0.193           0.552             0.014
## 17       0.012 1.064  0.795       0.212           1.270             0.940
## 18       0.015 0.936  1.045       0.084           0.440             0.115
## 19       0.016 0.590  0.672       0.048           0.479             1.340
## 20       0.018 0.492  1.067       0.279           0.465             0.216
##    black  hisp asian  nahn other
## 1  2.058 1.586 1.776 2.986 1.724
## 2  0.875 1.076 1.263 0.428 1.881
## 3  2.108 1.342 0.848 2.367 1.861
## 4  1.593 1.318 2.715 4.700 0.796
## 5  2.073 1.610 2.176 2.402 2.483
## 6  0.733 1.884 0.887 0.503 0.913
## 7  3.003 1.349 1.323 0.625 1.212
## 8  1.899 2.004 2.240 5.070 1.335
## 9  0.506 0.956 0.436 0.590 1.025
## 10 1.105 1.213 0.693 0.892 0.741
## 11 0.900 0.985 0.881 0.466 1.045
## 12 0.704 2.624 3.206 0.640 0.028
## 13 1.039 0.839 0.963 0.738 1.992
## 14 1.622 1.350 1.410 0.898 1.935
## 15 1.008 0.544 0.615 0.276 0.374
## 16 0.755 0.058 0.397 2.347 0.636
## 17 1.945 1.209 1.208 5.523 1.145
## 18 2.548 0.797 1.186 2.554 0.400
## 19 0.697 1.254 0.032 0.534 0.035
## 20 2.348 1.248 1.116 2.592 0.256
round(sums$pval,3)
##    (Intercept)  time momedu momworknilf momworkparttime momworkunemployed
## 1            0 0.364  0.000       0.000           0.025             0.001
## 2            0 0.011  0.030       0.668           0.095             0.074
## 3            0 0.267  0.420       0.005           0.270             0.167
## 4            0 0.195  0.000       0.001           0.350             0.154
## 5            0 0.237  0.000       0.000           0.161             0.066
## 6            0 0.036  0.060       0.000           0.000             0.028
## 7            0 0.324  0.010       0.000           0.021             0.008
## 8            0 0.037  0.245       0.000           0.001             0.141
## 9            0 0.057  0.000       0.000           0.282             0.005
## 10           0 0.336  0.222       0.000           0.000             0.006
## 11           0 0.685  0.019       0.000           0.000             0.011
## 12           0 0.454  0.316       0.001           0.361             0.049
## 13           0 0.051  0.000       0.000           0.001             0.002
## 14           0 0.007  0.004       0.000           0.000             0.104
## 15           0 0.115  0.147       0.000           0.000             0.012
## 16           0 0.395  0.286       0.001           0.124             0.473
## 17           0 0.721  0.199       0.000           0.226             0.883
## 18           0 0.729  0.824       0.000           0.002             0.042
## 19           0 0.117  0.231       0.001           0.086             0.631
## 20           0 0.009  0.805       0.000           0.037             0.166
##    black  hisp asian  nahn other
## 1  0.000 0.000 0.003 0.000 0.011
## 2  0.477 0.634 0.468 0.184 0.036
## 3  0.000 0.011 0.604 0.000 0.013
## 4  0.096 0.306 0.033 0.000 0.794
## 5  0.000 0.000 0.000 0.000 0.000
## 6  0.275 0.001 0.839 0.446 0.886
## 7  0.000 0.095 0.475 0.474 0.671
## 8  0.018 0.007 0.121 0.000 0.687
## 9  0.000 0.684 0.021 0.151 0.925
## 10 0.578 0.253 0.455 0.821 0.562
## 11 0.595 0.936 0.774 0.274 0.922
## 12 0.529 0.006 0.064 0.773 0.610
## 13 0.699 0.089 0.862 0.293 0.000
## 14 0.002 0.072 0.325 0.829 0.049
## 15 0.963 0.008 0.330 0.131 0.164
## 16 0.534 0.047 0.495 0.208 0.705
## 17 0.003 0.443 0.742 0.000 0.831
## 18 0.000 0.486 0.782 0.062 0.424
## 19 0.471 0.568 0.529 0.684 0.547
## 20 0.005 0.542 0.897 0.150 0.478
dat1<-data.frame(time=1:2, momedu=1, momwork=levels(e.long$momwork)[1], black=0, hisp=0, asian=0, nahn=0, other=0)

test1<-predict(multnm, newdata=dat1, type="probs")
test1
##           0          1          2          3           4          5
## 1 0.5742358 0.02854778 0.01756713 0.02274307 0.003125732 0.05278837
## 2 0.5859232 0.03099568 0.01335048 0.02550288 0.002445317 0.05776182
##             6          7           8          9         10         11
## 1 0.011702923 0.01567970 0.006472057 0.05955997 0.02142809 0.02408080
## 2 0.008312266 0.01806438 0.004299419 0.05206791 0.01932458 0.02326736
##            12         13         14         15          16          17
## 1 0.002985629 0.07446911 0.02833848 0.02899451 0.004216369 0.005761973
## 2 0.002420345 0.08684824 0.02084453 0.02387554 0.005591943 0.006254698
##            18          19          20
## 1 0.008397824 0.003602394 0.005302312
## 2 0.008016815 0.002170410 0.002662185
barplot(test1[1,], names.arg = rownames(test1[1,]), ylim=c(0, .6))

barplot(test1[2,], names.arg = rownames(test1[2,]), col=2,ylim=c(0, .6))

dat<-expand.grid(time=c(1,2), momedu=c(0,1), momwork=levels(e.long$momwork), black=c(0,1), hisp=c(0,1),asian=c(0,1),nahn=c(0,1),other=c(0,1))
rem<-which(apply(dat[,4:8],1,sum)>1);
dat<-dat[-rem]

dat$fitted<-round(predict(multnm, newdata=dat, type="probs"), 3)
head(dat[dat$momwork=="fulltime"&dat$time==1,], n=20)
##     time momedu  momwork black hisp asian nahn other fitted.0 fitted.1
## 1      1      0 fulltime     0    0     0    0     0    0.590    0.042
## 3      1      1 fulltime     0    0     0    0     0    0.574    0.029
## 17     1      0 fulltime     1    0     0    0     0    0.492    0.072
## 19     1      1 fulltime     1    0     0    0     0    0.495    0.051
## 33     1      0 fulltime     0    1     0    0     0    0.534    0.060
## 35     1      1 fulltime     0    1     0    0     0    0.537    0.042
## 49     1      0 fulltime     1    1     0    0     0    0.419    0.097
## 51     1      1 fulltime     1    1     0    0     0    0.438    0.071
## 65     1      0 fulltime     0    0     1    0     0    0.530    0.067
## 67     1      1 fulltime     0    0     1    0     0    0.541    0.048
## 81     1      0 fulltime     1    0     1    0     0    0.396    0.103
## 83     1      1 fulltime     1    0     1    0     0    0.424    0.077
## 97     1      0 fulltime     0    1     1    0     0    0.446    0.089
## 99     1      1 fulltime     0    1     1    0     0    0.473    0.066
## 113    1      0 fulltime     1    1     1    0     0    0.311    0.128
## 115    1      1 fulltime     1    1     1    0     0    0.347    0.100
## 129    1      0 fulltime     0    0     0    1     0    0.469    0.100
## 131    1      1 fulltime     0    0     0    1     0    0.501    0.074
## 145    1      0 fulltime     1    0     0    1     0    0.326    0.143
## 147    1      1 fulltime     1    0     0    1     0    0.367    0.112
##     fitted.2 fitted.3 fitted.4 fitted.5 fitted.6 fitted.7 fitted.8
## 1      0.014    0.022    0.008    0.078    0.017    0.011    0.008
## 3      0.018    0.023    0.003    0.053    0.012    0.016    0.006
## 17     0.010    0.038    0.010    0.136    0.010    0.029    0.013
## 19     0.013    0.041    0.004    0.094    0.007    0.041    0.011
## 33     0.013    0.026    0.009    0.114    0.028    0.014    0.015
## 35     0.018    0.029    0.004    0.079    0.021    0.020    0.012
## 49     0.009    0.044    0.011    0.186    0.016    0.033    0.023
## 51     0.013    0.049    0.005    0.134    0.012    0.048    0.019
## 65     0.016    0.017    0.018    0.153    0.013    0.014    0.017
## 67     0.021    0.018    0.008    0.108    0.010    0.020    0.014
## 81     0.010    0.026    0.022    0.238    0.007    0.031    0.024
## 83     0.014    0.030    0.010    0.176    0.006    0.046    0.020
## 97     0.014    0.019    0.020    0.208    0.021    0.015    0.029
## 99     0.020    0.021    0.009    0.152    0.016    0.023    0.024
## 113    0.009    0.027    0.023    0.300    0.011    0.032    0.038
## 115    0.013    0.033    0.011    0.232    0.009    0.051    0.033
## 129    0.005    0.041    0.028    0.150    0.007    0.006    0.034
## 131    0.007    0.047    0.013    0.111    0.005    0.009    0.029
## 145    0.003    0.060    0.031    0.216    0.003    0.012    0.045
## 147    0.004    0.072    0.015    0.168    0.003    0.019    0.040
##     fitted.9 fitted.10 fitted.11 fitted.12 fitted.13 fitted.14 fitted.15
## 1      0.043     0.026     0.018     0.002     0.047     0.020     0.024
## 3      0.060     0.021     0.024     0.003     0.074     0.028     0.029
## 17     0.018     0.024     0.013     0.001     0.041     0.027     0.020
## 19     0.026     0.020     0.019     0.002     0.067     0.040     0.025
## 33     0.038     0.028     0.016     0.005     0.036     0.025     0.012
## 35     0.053     0.024     0.022     0.007     0.058     0.036     0.015
## 49     0.015     0.025     0.011     0.003     0.029     0.031     0.009
## 51     0.022     0.022     0.016     0.004     0.049     0.047     0.012
## 65     0.017     0.016     0.014     0.006     0.040     0.025     0.013
## 67     0.024     0.014     0.020     0.009     0.068     0.038     0.017
## 81     0.006     0.013     0.009     0.003     0.031     0.031     0.010
## 83     0.010     0.012     0.014     0.005     0.055     0.048     0.013
## 97     0.014     0.016     0.012     0.014     0.029     0.029     0.006
## 99     0.020     0.015     0.017     0.021     0.050     0.044     0.008
## 113    0.005     0.013     0.007     0.007     0.021     0.033     0.004
## 115    0.008     0.012     0.011     0.011     0.038     0.053     0.006
## 129    0.020     0.018     0.006     0.001     0.027     0.014     0.005
## 131    0.031     0.017     0.010     0.002     0.048     0.022     0.007
## 145    0.007     0.014     0.004     0.001     0.020     0.016     0.004
## 147    0.011     0.013     0.006     0.001     0.036     0.026     0.005
##     fitted.16 fitted.17 fitted.18 fitted.19 fitted.20
## 1       0.003     0.007     0.008     0.006     0.005
## 3       0.004     0.006     0.008     0.004     0.005
## 17      0.002     0.012     0.018     0.003     0.010
## 19      0.003     0.010     0.018     0.002     0.011
## 33      0.000     0.008     0.006     0.006     0.006
## 35      0.000     0.007     0.006     0.004     0.006
## 49      0.000     0.012     0.012     0.003     0.011
## 51      0.000     0.010     0.013     0.002     0.012
## 65      0.001     0.008     0.009     0.000     0.005
## 67      0.002     0.007     0.009     0.000     0.006
## 81      0.001     0.012     0.017     0.000     0.009
## 83      0.001     0.010     0.019     0.000     0.010
## 97      0.000     0.008     0.006     0.000     0.005
## 99      0.000     0.007     0.007     0.000     0.006
## 113     0.000     0.011     0.010     0.000     0.009
## 115     0.000     0.010     0.012     0.000     0.010
## 129     0.006     0.033     0.017     0.002     0.011
## 131     0.009     0.028     0.019     0.002     0.012
## 145     0.003     0.044     0.030     0.001     0.017
## 147     0.005     0.040     0.035     0.001     0.021

We can also use the Cox model to analyze the poverty transition outcome:

des2<-svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights=~c1_5fp0, data=e.long, nest=T)

#cox model
coxfit1<-svycoxph(Surv( time, povtran==1)~black+hisp+asian+other+nahn+momwork+momedu, des2)
coxfit2<-svycoxph(Surv( time, povtran2==1)~black+hisp+asian+other+nahn+momwork+momedu, des2)

summary(coxfit1)
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (520) clusters.
## svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights = ~c1_5fp0, 
##     data = e.long, nest = T)
## Call:
## svycoxph(formula = Surv(time, povtran == 1) ~ black + hisp + 
##     asian + other + nahn + momwork + momedu, design = des2)
## 
##   n= 14549, number of events= 813 
##    (5019 observations deleted due to missingness)
## 
##                       coef exp(coef) se(coef)       z Pr(>|z|)    
## black              1.44750   4.25249  0.14305  10.119  < 2e-16 ***
## hisp               0.92646   2.52556  0.10146   9.131  < 2e-16 ***
## asian              0.60560   1.83235  0.22205   2.727  0.00639 ** 
## other              1.10479   3.01860  0.23874   4.628 3.70e-06 ***
## nahn               0.93622   2.55033  0.39566   2.366  0.01797 *  
## momworknilf        0.26962   1.30946  0.10514   2.564  0.01033 *  
## momworkparttime    0.03108   1.03157  0.13213   0.235  0.81402    
## momworkunemployed  0.84241   2.32195  0.17177   4.904 9.38e-07 ***
## momedu            -1.32371   0.26615  0.10067 -13.149  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## black                4.2525     0.2352    3.2128    5.6287
## hisp                 2.5256     0.3960    2.0701    3.0812
## asian                1.8323     0.5457    1.1858    2.8315
## other                3.0186     0.3313    1.8906    4.8197
## nahn                 2.5503     0.3921    1.1744    5.5384
## momworknilf          1.3095     0.7637    1.0656    1.6091
## momworkparttime      1.0316     0.9694    0.7962    1.3365
## momworkunemployed    2.3220     0.4307    1.6582    3.2514
## momedu               0.2661     3.7573    0.2185    0.3242
## 
## Concordance= 0.759  (se = 0.01 )
## Rsquare= NA   (max possible= NA )
## Likelihood ratio test= NA  on 9 df,   p=NA
## Wald test            = 649.1  on 9 df,   p=0
## Score (logrank) test = NA  on 9 df,   p=NA
summary(coxfit2)
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (520) clusters.
## svydesign(ids = ~c15fppsu, strata = ~c15fpstr, weights = ~c1_5fp0, 
##     data = e.long, nest = T)
## Call:
## svycoxph(formula = Surv(time, povtran2 == 1) ~ black + hisp + 
##     asian + other + nahn + momwork + momedu, design = des2)
## 
##   n= 2564, number of events= 728 
##    (17004 observations deleted due to missingness)
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)    
## black             -0.77536   0.46054  0.16705 -4.641 3.46e-06 ***
## hisp              -0.45212   0.63628  0.15615 -2.895  0.00379 ** 
## asian             -0.32756   0.72068  0.20943 -1.564  0.11780    
## other             -1.03897   0.35382  0.36328 -2.860  0.00424 ** 
## nahn              -0.52053   0.59420  0.37804 -1.377  0.16854    
## momworknilf       -0.47783   0.62013  0.11962 -3.995 6.48e-05 ***
## momworkparttime   -0.02645   0.97389  0.12556 -0.211  0.83313    
## momworkunemployed -0.36233   0.69605  0.26539 -1.365  0.17216    
## momedu             0.52969   1.69840  0.10773  4.917 8.79e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## black                0.4605     2.1714    0.3319    0.6389
## hisp                 0.6363     1.5716    0.4685    0.8641
## asian                0.7207     1.3876    0.4780    1.0865
## other                0.3538     2.8263    0.1736    0.7211
## nahn                 0.5942     1.6829    0.2832    1.2466
## momworknilf          0.6201     1.6126    0.4905    0.7840
## momworkparttime      0.9739     1.0268    0.7614    1.2456
## momworkunemployed    0.6961     1.4367    0.4138    1.1710
## momedu               1.6984     0.5888    1.3751    2.0977
## 
## Concordance= 0.646  (se = 0.013 )
## Rsquare= NA   (max possible= NA )
## Likelihood ratio test= NA  on 9 df,   p=NA
## Wald test            = 85.27  on 9 df,   p=1.443e-14
## Score (logrank) test = NA  on 9 df,   p=NA
plot(survfit(coxfit1))

plot(survfit(coxfit2))

plot(survfit(coxfit1, newdata = data.frame(momedu=0, momwork="fulltime", black=1,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="red", ylim=c(.6, 1))
lines(survfit(coxfit1, newdata = data.frame(momedu=1, momwork="fulltime", black=1,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="green")
lines(survfit(coxfit1, newdata = data.frame(momedu=0, momwork="fulltime", black=0,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="blue")
title(main=c("Survival function for transition into Poverty between","Kindergarten and 5th Grade"))
legend("bottom", legend=c("Black, <=HS", "Black, > HS", "NH White, > HS"), lty=1, col=c("red", "green", "blue"), cex=.8)

plot(survfit(coxfit2, newdata = data.frame(momedu=0, momwork="fulltime", black=1,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="red", ylim=c(.4, 1))
lines(survfit(coxfit2, newdata = data.frame(momedu=1, momwork="fulltime", black=1,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="green")
lines(survfit(coxfit2, newdata = data.frame(momedu=0, momwork="fulltime", black=0,asian=0, hisp=0, other=0, nahn=0), conf.int=F) ,col="blue")
title(main=c("Survival function for transition out of poverty between","Kindergarten and 5th Grade"))
legend("bottom", legend=c("Black, <=HS", "Black, > HS", "NH White, > HS"), lty=1, col=c("red", "green", "blue"), cex=.8)