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
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchartAs 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/dem7903_App_Hier/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 0names(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" , drop = names(eclsk)[4:27])
e.long<-e.long[order(e.long$childid, e.long$time),]
head(e.long)## childid gender race s2_id c1_5fp0 c15fpstr c15fppsu p1primpk
## 0001002C.1 0001002C 2 1 0001 352.35 63 1 6
## 0001002C.2 0001002C 2 1 0001 352.35 63 1 6
## 0001004C.1 0001004C 1 1 0001 358.99 63 1 7
## 0001004C.2 0001004C 1 1 0001 358.99 63 1 7
## 0001006C.1 0001006C 1 1 0001 358.99 63 1 6
## 0001006C.2 0001006C 1 1 0001 358.99 63 1 6
## p1agefrs p1primnw p4primnw p5primnw age1 age2 age3
## 0001002C.1 38 0 4 1 6.433333 7.894167 9.750000
## 0001002C.2 38 0 4 1 6.433333 7.894167 9.750000
## 0001004C.1 3 5 5 5 5.494167 6.994167 8.916667
## 0001004C.2 3 5 5 5 5.494167 6.994167 8.916667
## 0001006C.1 7 0 0 0 5.250000 6.750000 8.750000
## 0001006C.2 7 0 0 0 5.250000 6.750000 8.750000
## pov1 pov2 pov3 care1 care2 care3 momedu3
## 0001002C.1 0 0 0 none nonrelative relative 1
## 0001002C.2 0 0 0 none nonrelative relative 1
## 0001004C.1 0 0 NA centerbased centerbased centerbased NA
## 0001004C.2 0 0 NA centerbased centerbased centerbased NA
## 0001006C.1 0 0 NA none none none NA
## 0001006C.2 0 0 NA none none none NA
## momwork3 bothpar1 bothpar2 bothpar3 parenttype3 hisp black
## 0001002C.1 parttime 0 0 0 twoparent 0 0
## 0001002C.2 parttime 0 0 0 twoparent 0 0
## 0001004C.1 parttime 0 0 0 twoparent 0 0
## 0001004C.2 parttime 0 0 0 twoparent 0 0
## 0001006C.1 parttime 0 0 0 twoparent 0 0
## 0001006C.2 parttime 0 0 0 twoparent 0 0
## asian nahn other male time momedu1 momwork1 parenttype1
## 0001002C.1 0 0 0 0 1 0 parttime twoparent
## 0001002C.2 0 0 0 0 2 1 parttime twoparent
## 0001004C.1 0 0 0 1 1 1 fulltime twoparent
## 0001004C.2 0 0 0 1 2 1 fulltime twoparent
## 0001006C.1 0 0 0 1 1 1 parttime twoparent
## 0001006C.2 0 0 0 1 2 1 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 1table(e.long$time, e.long$povtran)##
## 0 1
## 1 7818 512
## 2 6200 327table(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 2table(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 25prop.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.07619048table(e.long$care2,e.long$care3)##
## none relative nonrelative centerbased other/multi
## none 9222 968 280 590 86
## relative 1324 1650 112 288 90
## nonrelative 604 238 520 228 38
## centerbased 964 272 172 1508 60
## other/multi 118 122 30 50 34library(nnet)
multnm<-multinom(factor(caretrn)~time+momedu1+momwork1+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
## convergedsums<-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 momedu1 momwork1nilf momwork1parttime
## 1 0.067 1.064 0.698 0.647 1.218
## 2 0.031 0.745 1.307 0.942 1.279
## 3 0.034 1.099 1.075 0.745 1.129
## 4 0.017 0.767 0.424 0.382 1.259
## 5 0.124 1.072 0.691 0.242 0.904
## 6 0.041 0.696 0.720 0.119 0.366
## 7 0.017 1.129 1.404 0.255 0.697
## 8 0.022 0.651 0.787 0.174 0.362
## 9 0.086 0.857 1.409 0.254 0.904
## 10 0.050 0.884 0.851 0.076 0.469
## 11 0.031 0.947 1.410 0.075 0.474
## 12 0.005 0.794 1.388 0.153 0.702
## 13 0.069 1.143 1.634 0.247 0.759
## 14 0.047 0.721 1.449 0.088 0.400
## 15 0.051 0.807 1.233 0.111 0.241
## 16 0.004 1.300 1.443 0.193 0.552
## 17 0.012 1.064 0.795 0.212 1.270
## 18 0.015 0.936 1.045 0.084 0.440
## 19 0.016 0.590 0.672 0.048 0.479
## 20 0.018 0.492 1.067 0.279 0.465
## momwork1unemployed black hisp asian nahn other
## 1 1.631 2.058 1.586 1.776 2.986 1.724
## 2 1.652 0.875 1.076 1.263 0.428 1.881
## 3 1.321 2.108 1.342 0.848 2.367 1.861
## 4 1.738 1.593 1.318 2.715 4.700 0.796
## 5 0.759 2.073 1.610 2.176 2.402 2.483
## 6 0.162 0.733 1.884 0.887 0.503 0.913
## 7 0.272 3.003 1.349 1.323 0.625 1.212
## 8 0.408 1.899 2.004 2.240 5.070 1.335
## 9 0.410 0.506 0.956 0.436 0.590 1.025
## 10 0.061 1.105 1.213 0.693 0.892 0.741
## 11 0.207 0.900 0.985 0.881 0.466 1.045
## 12 2.613 0.704 2.624 3.206 0.640 0.028
## 13 0.498 1.039 0.839 0.963 0.738 1.992
## 14 0.593 1.622 1.350 1.410 0.898 1.935
## 15 0.054 1.008 0.544 0.615 0.276 0.374
## 16 0.014 0.755 0.058 0.397 2.347 0.636
## 17 0.940 1.945 1.209 1.208 5.523 1.145
## 18 0.115 2.548 0.797 1.186 2.554 0.400
## 19 1.340 0.697 1.254 0.032 0.534 0.035
## 20 0.216 2.348 1.248 1.116 2.592 0.256round(sums$pval,3)## (Intercept) time momedu1 momwork1nilf momwork1parttime
## 1 0 0.364 0.000 0.000 0.025
## 2 0 0.011 0.030 0.668 0.095
## 3 0 0.267 0.420 0.005 0.270
## 4 0 0.195 0.000 0.001 0.350
## 5 0 0.237 0.000 0.000 0.161
## 6 0 0.036 0.060 0.000 0.000
## 7 0 0.324 0.010 0.000 0.021
## 8 0 0.037 0.245 0.000 0.001
## 9 0 0.057 0.000 0.000 0.282
## 10 0 0.336 0.222 0.000 0.000
## 11 0 0.685 0.019 0.000 0.000
## 12 0 0.454 0.316 0.001 0.361
## 13 0 0.051 0.000 0.000 0.001
## 14 0 0.007 0.004 0.000 0.000
## 15 0 0.115 0.147 0.000 0.000
## 16 0 0.395 0.286 0.001 0.124
## 17 0 0.721 0.199 0.000 0.226
## 18 0 0.729 0.824 0.000 0.002
## 19 0 0.117 0.231 0.001 0.086
## 20 0 0.009 0.805 0.000 0.037
## momwork1unemployed black hisp asian nahn other
## 1 0.001 0.000 0.000 0.003 0.000 0.011
## 2 0.074 0.477 0.634 0.468 0.184 0.036
## 3 0.167 0.000 0.011 0.604 0.000 0.013
## 4 0.154 0.096 0.306 0.033 0.000 0.794
## 5 0.066 0.000 0.000 0.000 0.000 0.000
## 6 0.028 0.275 0.001 0.839 0.446 0.886
## 7 0.008 0.000 0.095 0.475 0.474 0.671
## 8 0.141 0.018 0.007 0.121 0.000 0.687
## 9 0.005 0.000 0.684 0.021 0.151 0.925
## 10 0.006 0.578 0.253 0.455 0.821 0.562
## 11 0.011 0.595 0.936 0.774 0.274 0.922
## 12 0.049 0.529 0.006 0.064 0.773 0.610
## 13 0.002 0.699 0.089 0.862 0.293 0.000
## 14 0.104 0.002 0.072 0.325 0.829 0.049
## 15 0.012 0.963 0.008 0.330 0.131 0.164
## 16 0.473 0.534 0.047 0.495 0.208 0.705
## 17 0.883 0.003 0.443 0.742 0.000 0.831
## 18 0.042 0.000 0.486 0.782 0.062 0.424
## 19 0.631 0.471 0.568 0.529 0.684 0.547
## 20 0.166 0.005 0.542 0.897 0.150 0.478dat<-expand.grid(time=c(1,2), momedu1=c(0,1), momwork1=levels(e.long$momwork1), 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$momwork1=="fulltime",], n=10)## time momedu1 momwork1 black hisp asian nahn other fitted.0 fitted.1
## 1 1 0 fulltime 0 0 0 0 0 0.590 0.042
## 2 2 0 fulltime 0 0 0 0 0 0.602 0.046
## 3 1 1 fulltime 0 0 0 0 0 0.574 0.029
## 4 2 1 fulltime 0 0 0 0 0 0.586 0.031
## 17 1 0 fulltime 1 0 0 0 0 0.492 0.072
## 18 2 0 fulltime 1 0 0 0 0 0.497 0.077
## 19 1 1 fulltime 1 0 0 0 0 0.495 0.051
## 20 2 1 fulltime 1 0 0 0 0 0.501 0.055
## 33 1 0 fulltime 0 1 0 0 0 0.534 0.060
## 34 2 0 fulltime 0 1 0 0 0 0.546 0.066
## fitted.2 fitted.3 fitted.4 fitted.5 fitted.6 fitted.7 fitted.8 fitted.9
## 1 0.014 0.022 0.008 0.078 0.017 0.011 0.008 0.043
## 2 0.010 0.024 0.006 0.086 0.012 0.013 0.006 0.038
## 3 0.018 0.023 0.003 0.053 0.012 0.016 0.006 0.060
## 4 0.013 0.026 0.002 0.058 0.008 0.018 0.004 0.052
## 17 0.010 0.038 0.010 0.136 0.010 0.029 0.013 0.018
## 18 0.008 0.042 0.008 0.147 0.007 0.033 0.009 0.016
## 19 0.013 0.041 0.004 0.094 0.007 0.041 0.011 0.026
## 20 0.010 0.046 0.003 0.102 0.005 0.046 0.007 0.022
## 33 0.013 0.026 0.009 0.114 0.028 0.014 0.015 0.038
## 34 0.010 0.030 0.007 0.125 0.020 0.016 0.010 0.033
## fitted.10 fitted.11 fitted.12 fitted.13 fitted.14 fitted.15 fitted.16
## 1 0.026 0.018 0.002 0.047 0.020 0.024 0.003
## 2 0.023 0.017 0.002 0.055 0.015 0.020 0.004
## 3 0.021 0.024 0.003 0.074 0.028 0.029 0.004
## 4 0.019 0.023 0.002 0.087 0.021 0.024 0.006
## 17 0.024 0.013 0.001 0.041 0.027 0.020 0.002
## 18 0.021 0.013 0.001 0.047 0.020 0.017 0.002
## 19 0.020 0.019 0.002 0.067 0.040 0.025 0.003
## 20 0.018 0.018 0.001 0.077 0.029 0.021 0.004
## 33 0.028 0.016 0.005 0.036 0.025 0.012 0.000
## 34 0.026 0.015 0.004 0.042 0.018 0.010 0.000
## fitted.17 fitted.18 fitted.19 fitted.20
## 1 0.007 0.008 0.006 0.005
## 2 0.008 0.008 0.003 0.003
## 3 0.006 0.008 0.004 0.005
## 4 0.006 0.008 0.002 0.003
## 17 0.012 0.018 0.003 0.010
## 18 0.013 0.017 0.002 0.005
## 19 0.010 0.018 0.002 0.011
## 20 0.010 0.017 0.001 0.005
## 33 0.008 0.006 0.006 0.006
## 34 0.009 0.006 0.004 0.003