library(car)
## Loading required package: carData
library(mice)
## Warning: package 'mice' was built under R version 3.5.2
## Loading required package: lattice
##
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(ggplot2)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(foreign)
brfss2017=read.xport("~/Desktop/DataFolder/brfss.XPT ")
brfss2017$menth<-Recode(brfss2017$MENTHLTH,recodes="88=0;77=NA;99=NA",as.factor=T)
brfss2017$hlthinsur<-Recode(brfss2017$HLTHPLN1,recodes="1='Y';2='N';else=NA",as.factor=T)
brfss2017$medcost<-Recode(brfss2017$MEDCOST,recodes="1='Y';2='N';else=NA",as.factor=T)
brfss2017$sxorient<-Recode(brfss2017$SXORIENT,recodes="1='S';2='G';else=NA",as.factor=T)
brfss2017$marital<-Recode(brfss2017$MARITAL, recodes="1='M';5='S';else=NA",as.factor=T)
summary(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## menth hlthinsur medcost sxorient marital
## 0 :300134 N : 35743 N :401190 G : 3167 M :232891
## 30 : 24843 Y :412502 Y : 47528 S :190310 S : 73939
## 2 : 21661 NA's: 1771 NA's: 1298 NA's:256539 NA's:143186
## 5 : 15711
## 1 : 13973
## (Other): 66491
## NA's : 7203
#Look at patterns of missingness
md.pattern(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## medcost hlthinsur menth marital sxorient
## 128426 1 1 1 1 1 0
## 172744 1 1 1 1 0 1
## 61305 1 1 1 0 1 1
## 77596 1 1 1 0 0 2
## 1471 1 1 0 1 1 1
## 2251 1 1 0 1 0 2
## 1235 1 1 0 0 1 2
## 1976 1 1 0 0 0 3
## 415 1 0 1 1 1 1
## 739 1 0 1 1 0 2
## 149 1 0 1 0 1 2
## 293 1 0 1 0 0 3
## 25 1 0 0 1 1 2
## 42 1 0 0 1 0 3
## 14 1 0 0 0 1 3
## 37 1 0 0 0 0 4
## 225 0 1 1 1 1 1
## 404 0 1 1 1 0 2
## 156 0 1 1 0 1 2
## 316 0 1 1 0 0 3
## 16 0 1 0 1 1 2
## 44 0 1 0 1 0 3
## 21 0 1 0 0 1 3
## 59 0 1 0 0 0 4
## 9 0 0 1 1 1 2
## 16 0 0 1 1 0 3
## 6 0 0 1 0 1 3
## 14 0 0 1 0 0 4
## 2 0 0 0 1 1 3
## 1 0 0 0 1 0 4
## 2 0 0 0 0 1 4
## 7 0 0 0 0 0 5
## 1298 1771 7203 143186 256539 409997
#Let's see how pairs of variables are missing together...
md.pairs(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## $rr
## menth hlthinsur medcost sxorient marital
## menth 442813 441172 441667 190691 302978
## hlthinsur 441172 448245 447004 192855 305581
## medcost 441667 447004 448718 193040 306113
## sxorient 190691 192855 193040 193477 130589
## marital 302978 305581 306113 130589 306830
##
## $rm
## menth hlthinsur medcost sxorient marital
## menth 0 1641 1146 252122 139835
## hlthinsur 7073 0 1241 255390 142664
## medcost 7051 1714 0 255678 142605
## sxorient 2786 622 437 0 62888
## marital 3852 1249 717 176241 0
##
## $mr
## menth hlthinsur medcost sxorient marital
## menth 0 7073 7051 2786 3852
## hlthinsur 1641 0 1714 622 1249
## medcost 1146 1241 0 437 717
## sxorient 252122 255390 255678 0 176241
## marital 139835 142664 142605 62888 0
##
## $mm
## menth hlthinsur medcost sxorient marital
## menth 7203 130 152 4417 3351
## hlthinsur 130 1771 57 1149 522
## medcost 152 57 1298 861 581
## sxorient 4417 1149 861 256539 80298
## marital 3351 522 581 80298 143186
#find the most common value
mcv.sxorient<-factor(names(which.max(table(brfss2017$sxorient))),levels=levels(brfss2017$sxorient))
mcv.sxorient
## [1] S
## Levels: G S
mcv.marital<-factor(names(which.max(table(brfss2017$marital))),levels=levels(brfss2017$marital))
mcv.marital
## [1] M
## Levels: M S
mcv.medcost<-factor(names(which.max(table(brfss2017$medcost))),levels=levels(brfss2017$medcost))
mcv.medcost
## [1] N
## Levels: N Y
mcv.hlthinsur<-factor(names(which.max(table(brfss2017$hlthinsur))),levels=levels(brfss2017$hlthinsur))
mcv.hlthinsur
## [1] Y
## Levels: N Y
mcv.menth<-factor(names(which.max(table(brfss2017$menth))),levels=levels(brfss2017$menth))
mcv.menth
## [1] 0
## 31 Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ... 30
brfss2017$sxorient.imp<-as.factor(ifelse(is.na(brfss2017$sxorient)==T, mcv.sxorient, brfss2017$sxorient))
levels(brfss2017$sxorient.imp)<-levels(brfss2017$sxorient)
prop.table(table(brfss2017$sxorient))
##
## G S
## 0.01636887 0.98363113
prop.table(table(brfss2017$sxorient.imp))
##
## G S
## 0.007037528 0.992962472
barplot(prop.table(table(brfss2017$sxorient)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$sxorient.imp)), main="Imputed Data",ylim=c(0, .6))
brfss2017$marital.imp<-as.factor(ifelse(is.na(brfss2017$marital)==T, mcv.marital, brfss2017$marital))
levels(brfss2017$marital.imp)<-levels(brfss2017$marital)
prop.table(table(brfss2017$marital))
##
## M S
## 0.7590229 0.2409771
prop.table(table(brfss2017$marital.imp))
##
## M S
## 0.835697 0.164303
barplot(prop.table(table(brfss2017$marital)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$marital.imp)), main="Imputed Data",ylim=c(0, .6))
brfss2017$medcost.imp<-as.factor(ifelse(is.na(brfss2017$medcost)==T, mcv.medcost, brfss2017$medcost))
levels(brfss2017$medcost.imp)<-levels(brfss2017$medcost)
prop.table(table(brfss2017$medcost))
##
## N Y
## 0.8940805 0.1059195
prop.table(table(brfss2017$medcost.imp))
##
## N Y
## 0.894386 0.105614
barplot(prop.table(table(brfss2017$medcost)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$medcost.imp)), main="Imputed Data",ylim=c(0, .6))
brfss2017$hlthinsur.imp<-as.factor(ifelse(is.na(brfss2017$hlthinsur)==T, mcv.hlthinsur, brfss2017$hlthinsur))
levels(brfss2017$hlthinsur.imp)<-levels(brfss2017$hlthinsur)
prop.table(table(brfss2017$hlthinsur))
##
## N Y
## 0.07973987 0.92026013
prop.table(table(brfss2017$hlthinsur.imp))
##
## N Y
## 0.07942606 0.92057394
barplot(prop.table(table(brfss2017$hlthinsur)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$hlthinsur.imp)), main="Imputed Data",ylim=c(0, .6))
brfss2017$menth.imp<-as.factor(ifelse(is.na(brfss2017$menth)==T, mcv.menth, brfss2017$menth))
levels(brfss2017$menth.imp)<-levels(brfss2017$menth)
prop.table(table(brfss2017$menth))
##
## 0 1 2 3 4
## 6.777895e-01 3.155508e-02 4.891681e-02 2.904612e-02 1.454564e-02
## 5 6 7 8 9
## 3.547999e-02 4.236551e-03 1.424529e-02 2.838670e-03 4.674660e-04
## 10 11 12 13 14
## 2.595678e-02 1.987295e-04 1.987295e-03 3.251937e-04 5.417637e-03
## 15 16 17 18 19
## 2.478247e-02 4.629494e-04 3.635846e-04 5.239232e-04 6.774869e-05
## 20 21 22 23 24
## 1.495665e-02 8.829912e-04 2.913194e-04 1.648551e-04 2.235707e-04
## 25 26 27 28 29
## 5.207616e-03 1.987295e-04 4.132670e-04 1.497246e-03 8.536335e-04
## 30
## 5.610269e-02
prop.table(table(brfss2017$menth.imp))
##
## 0 1 2 3 4
## 0.6829468286 0.0310500071 0.0481338441 0.0285812060 0.0143128244
## 5 6 7 8 9
## 0.0349120920 0.0041687407 0.0140172794 0.0027932340 0.0004599836
## 10 11 12 13 14
## 0.0255413141 0.0001955486 0.0019554860 0.0003199886 0.0053309216
## 15 16 17 18 19
## 0.0243857996 0.0004555394 0.0003577651 0.0005155372 0.0000666643
## 20 21 22 23 24
## 0.0147172545 0.0008688580 0.0002866565 0.0001622165 0.0002199922
## 25 26 27 28 29
## 0.0051242622 0.0001955486 0.0004066522 0.0014732810 0.0008399701
## 30
## 0.0552047038
barplot(prop.table(table(brfss2017$menth)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$menth.imp)), main="Imputed Data",ylim=c(0, .6))
#Multiple imputation
imp<-mice(data=brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")],seed=22,m=5)
##
## iter imp variable
## 1 1 menth hlthinsur medcost sxorient marital
## 1 2 menth hlthinsur medcost sxorient marital
## 1 3 menth hlthinsur medcost sxorient marital
## 1 4 menth hlthinsur medcost sxorient marital
## 1 5 menth hlthinsur medcost sxorient marital
## 2 1 menth hlthinsur medcost sxorient marital
## 2 2 menth hlthinsur medcost sxorient marital
## 2 3 menth hlthinsur medcost sxorient marital
## 2 4 menth hlthinsur medcost sxorient marital
## 2 5 menth hlthinsur medcost sxorient marital
## 3 1 menth hlthinsur medcost sxorient marital
## 3 2 menth hlthinsur medcost sxorient marital
## 3 3 menth hlthinsur medcost sxorient marital
## 3 4 menth hlthinsur medcost sxorient marital
## 3 5 menth hlthinsur medcost sxorient marital
## 4 1 menth hlthinsur medcost sxorient marital
## 4 2 menth hlthinsur medcost sxorient marital
## 4 3 menth hlthinsur medcost sxorient marital
## 4 4 menth hlthinsur medcost sxorient marital
## 4 5 menth hlthinsur medcost sxorient marital
## 5 1 menth hlthinsur medcost sxorient marital
## 5 2 menth hlthinsur medcost sxorient marital
## 5 3 menth hlthinsur medcost sxorient marital
## 5 4 menth hlthinsur medcost sxorient marital
## 5 5 menth hlthinsur medcost sxorient marital
## Warning: Number of logged events: 83
print(imp)
## Class: mids
## Number of multiple imputations: 5
## Imputation methods:
## menth hlthinsur medcost sxorient marital
## "polyreg" "logreg" "logreg" "logreg" "logreg"
## PredictorMatrix:
## menth hlthinsur medcost sxorient marital
## menth 0 1 1 1 1
## hlthinsur 1 0 1 1 1
## medcost 1 1 0 1 1
## sxorient 1 1 1 0 1
## marital 1 1 1 1 0
## Number of logged events: 83
## it im dep meth out
## 1 1 1 hlthinsur logreg menth19
## 2 1 1 medcost logreg menth19
## 3 1 1 sxorient logreg menth19
## 4 1 1 marital logreg menth19
## 5 1 2 hlthinsur logreg menth19
## 6 1 2 medcost logreg menth19
dat.imp<-complete(imp,action=1)
head(dat.imp,n=20)
## menth hlthinsur medcost sxorient marital
## 1 0 Y N S M
## 2 0 Y N S M
## 3 0 Y N S M
## 4 0 Y N S M
## 5 0 Y N S M
## 6 0 Y N S M
## 7 0 Y N S M
## 8 0 Y N S M
## 9 25 Y N S M
## 10 1 Y N S M
## 11 0 Y N S M
## 12 0 Y N S M
## 13 0 Y N S M
## 14 0 Y N S M
## 15 0 Y N S M
## 16 0 Y N S M
## 17 0 Y N S M
## 18 0 Y N S M
## 19 30 Y N S M
## 20 0 N N S M
#compare to the original data
head(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")],n=20)
## menth hlthinsur medcost sxorient marital
## 1 0 Y N <NA> <NA>
## 2 0 Y N <NA> M
## 3 0 Y N <NA> M
## 4 0 Y N <NA> <NA>
## 5 0 Y N <NA> <NA>
## 6 0 Y N <NA> <NA>
## 7 0 Y N <NA> M
## 8 0 Y N <NA> <NA>
## 9 25 Y N <NA> M
## 10 1 Y N <NA> <NA>
## 11 0 Y N <NA> <NA>
## 12 0 Y N <NA> <NA>
## 13 0 Y N <NA> <NA>
## 14 0 Y N <NA> M
## 15 0 Y N <NA> M
## 16 0 Y N <NA> <NA>
## 17 0 Y N <NA> <NA>
## 18 <NA> Y N <NA> <NA>
## 19 30 Y N <NA> <NA>
## 20 0 N N <NA> <NA>
fit.medcost<-with(data=imp,expr=lm(medcost~sxorient+marital+hlthinsur+menth))
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
fit.medcost
## call :
## with.mids(data = imp, expr = lm(medcost ~ sxorient + marital +
## hlthinsur + menth))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 7203 1771 1298 256539 143186
##
## analyses :
## [[1]]
##
## Call:
## lm(formula = medcost ~ sxorient + marital + hlthinsur + menth)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY menth1
## 1.315824 0.007322 0.022407 -0.278471 0.033247
## menth2 menth3 menth4 menth5 menth6
## 0.044017 0.058818 0.074817 0.073632 0.102440
## menth7 menth8 menth9 menth10 menth11
## 0.110330 0.118964 0.175612 0.107438 0.100749
## menth12 menth13 menth14 menth15 menth16
## 0.164317 0.155460 0.144571 0.137274 0.221114
## menth17 menth18 menth19 menth20 menth21
## 0.182974 0.140609 0.552699 0.169977 0.176432
## menth22 menth23 menth24 menth25 menth26
## 0.130470 0.228769 0.315181 0.169919 0.219150
## menth27 menth28 menth29 menth30
## 0.191353 0.209948 0.168221 0.191851
##
##
## [[2]]
##
## Call:
## lm(formula = medcost ~ sxorient + marital + hlthinsur + menth)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY menth1
## 1.30605 0.01663 0.02362 -0.27800 0.03357
## menth2 menth3 menth4 menth5 menth6
## 0.04323 0.05910 0.07428 0.07232 0.10004
## menth7 menth8 menth9 menth10 menth11
## 0.10524 0.12393 0.18364 0.10582 0.08172
## menth12 menth13 menth14 menth15 menth16
## 0.15892 0.17399 0.14637 0.13889 0.19176
## menth17 menth18 menth19 menth20 menth21
## 0.17907 0.14139 0.18094 0.17019 0.19286
## menth22 menth23 menth24 menth25 menth26
## 0.18983 0.20212 0.29457 0.17085 0.26537
## menth27 menth28 menth29 menth30
## 0.16628 0.20945 0.16834 0.19105
##
##
## [[3]]
##
## Call:
## lm(formula = medcost ~ sxorient + marital + hlthinsur + menth)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY menth1
## 1.320329 0.002152 0.023929 -0.278106 0.034140
## menth2 menth3 menth4 menth5 menth6
## 0.043804 0.058773 0.077379 0.072838 0.098975
## menth7 menth8 menth9 menth10 menth11
## 0.104329 0.117892 0.204228 0.107136 0.074838
## menth12 menth13 menth14 menth15 menth16
## 0.165840 0.143811 0.138470 0.135658 0.194247
## menth17 menth18 menth19 menth20 menth21
## 0.169491 0.167640 0.215022 0.170305 0.184884
## menth22 menth23 menth24 menth25 menth26
## 0.118580 0.237609 0.261767 0.166496 0.198076
## menth27 menth28 menth29 menth30
## 0.171988 0.208537 0.171109 0.193082
##
##
## [[4]]
##
## Call:
## lm(formula = medcost ~ sxorient + marital + hlthinsur + menth)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY menth1
## 1.315898 0.007552 0.023063 -0.278454 0.031945
## menth2 menth3 menth4 menth5 menth6
## 0.042079 0.057917 0.072290 0.071466 0.097126
## menth7 menth8 menth9 menth10 menth11
## 0.103850 0.127157 0.169804 0.108018 0.057621
## menth12 menth13 menth14 menth15 menth16
## 0.160279 0.155768 0.138171 0.138768 0.204975
## menth17 menth18 menth19 menth20 menth21
## 0.169220 0.141436 0.150977 0.168838 0.186721
## menth22 menth23 menth24 menth25 menth26
## 0.139654 0.250025 0.250326 0.167148 0.251222
## menth27 menth28 menth29 menth30
## 0.171380 0.221416 0.185457 0.189881
##
##
## [[5]]
##
## Call:
## lm(formula = medcost ~ sxorient + marital + hlthinsur + menth)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY menth1
## 1.324529 -0.001018 0.022636 -0.278643 0.031968
## menth2 menth3 menth4 menth5 menth6
## 0.043694 0.058835 0.073432 0.074154 0.095260
## menth7 menth8 menth9 menth10 menth11
## 0.099728 0.117560 0.185639 0.109510 0.059035
## menth12 menth13 menth14 menth15 menth16
## 0.159508 0.197981 0.139320 0.138713 0.197408
## menth17 menth18 menth19 menth20 menth21
## 0.172766 0.149258 0.272098 0.169639 0.181810
## menth22 menth23 menth24 menth25 menth26
## 0.130215 0.227693 0.235637 0.167465 0.239681
## menth27 menth28 menth29 menth30
## 0.259925 0.207774 0.170602 0.188507
esp.p<-pool(fit.medcost)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
print(esp.p)
## Class: mipo m = 5
## estimate ubar b t dfcom df riv lambda fmi
## (Intercept) 1.316526733 NA 4.726706e-05 NA 449982 NA NA NA NA
## sxorientS 0.006527692 NA 4.495935e-05 NA 449982 NA NA NA NA
## maritalS 0.023130156 NA 4.112672e-07 NA 449982 NA NA NA NA
## hlthinsurY -0.278335346 NA 7.250145e-08 NA 449982 NA NA NA NA
## menth1 0.032974222 NA 9.655629e-07 NA 449982 NA NA NA NA
## menth2 0.043365554 NA 5.993608e-07 NA 449982 NA NA NA NA
## menth3 0.058687621 NA 2.014692e-07 NA 449982 NA NA NA NA
## menth4 0.074439551 NA 3.610693e-06 NA 449982 NA NA NA NA
## menth5 0.072882282 NA 1.125303e-06 NA 449982 NA NA NA NA
## menth6 0.098768680 NA 7.537712e-06 NA 449982 NA NA NA NA
## menth7 0.104694764 NA 1.439055e-05 NA 449982 NA NA NA NA
## menth8 0.121100188 NA 1.801720e-05 NA 449982 NA NA NA NA
## menth9 0.183784314 NA 1.709126e-04 NA 449982 NA NA NA NA
## menth10 0.107584565 NA 1.808701e-06 NA 449982 NA NA NA NA
## menth11 0.074793317 NA 3.162268e-04 NA 449982 NA NA NA NA
## menth12 0.161773711 NA 9.620781e-06 NA 449982 NA NA NA NA
## menth13 0.165401636 NA 4.482311e-04 NA 449982 NA NA NA NA
## menth14 0.141381362 NA 1.453403e-05 NA 449982 NA NA NA NA
## menth15 0.137860212 NA 1.950415e-06 NA 449982 NA NA NA NA
## menth16 0.201899975 NA 1.400641e-04 NA 449982 NA NA NA NA
## menth17 0.174704157 NA 3.711406e-05 NA 449982 NA NA NA NA
## menth18 0.148066144 NA 1.321813e-04 NA 449982 NA NA NA NA
## menth19 0.274348057 NA 2.623717e-02 NA 449982 NA NA NA NA
## menth20 0.169789363 NA 3.469835e-07 NA 449982 NA NA NA NA
## menth21 0.184541511 NA 3.681791e-05 NA 449982 NA NA NA NA
## menth22 0.141749310 NA 7.782614e-04 NA 449982 NA NA NA NA
## menth23 0.229242532 NA 3.100889e-04 NA 449982 NA NA NA NA
## menth24 0.271496767 NA 1.067403e-03 NA 449982 NA NA NA NA
## menth25 0.168375500 NA 3.591787e-06 NA 449982 NA NA NA NA
## menth26 0.234699694 NA 7.053979e-04 NA 449982 NA NA NA NA
## menth27 0.192184464 NA 1.525341e-03 NA 449982 NA NA NA NA
## menth28 0.211424940 NA 3.189427e-05 NA 449982 NA NA NA NA
## menth29 0.172744842 NA 5.219573e-05 NA 449982 NA NA NA NA
## menth30 0.190875128 NA 3.113271e-06 NA 449982 NA NA NA NA