library(car)## Loading required package: carDatalibrary(VGAM)## Warning: package 'VGAM' was built under R version 3.6.3## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:car':
##
## logitlibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerlibrary(survey)## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:VGAM':
##
## calibrate## The following object is masked from 'package:graphics':
##
## dotchartload(url("https://github.com/coreysparks/data/blob/master/brfss_2017.Rdata?raw=true"))
#The names in the data are very ugly, so I make them less ugly
nams<-names(brfss_17)
head(nams, n=10)## [1] "dispcode" "statere1" "safetime" "hhadult" "genhlth" "physhlth"
## [7] "menthlth" "poorhlth" "hlthpln1" "persdoc2"#we see some names are lower case, some are upper and some have a little _ in the first position. This is a nightmare.
newnames<-tolower(gsub(pattern = "_",replacement = "",x = nams))
names(brfss_17)<-newnames
#Poor or fair self rated health
brfss_17$badhealth<-Recode(brfss_17$genhlth, recodes="4:5=1; 1:3=0; else=NA")
#sex
brfss_17$male<-as.factor(ifelse(brfss_17$sex==1, "Male", "Female"))
#Age cut into intervals
brfss_17$agec<-cut(brfss_17$age80, breaks=c(0,24,39,59,79,99))
#race/ethnicity
brfss_17$black<-Recode(brfss_17$racegr3, recodes="2=1; 9=NA; else=0")
brfss_17$white<-Recode(brfss_17$racegr3, recodes="1=1; 9=NA; else=0")
brfss_17$other<-Recode(brfss_17$racegr3, recodes="3:4=1; 9=NA; else=0")
brfss_17$hispanic<-Recode(brfss_17$racegr3, recodes="5=1; 9=NA; else=0")
brfss_17$race_eth<-Recode(brfss_17$racegr3, recodes="1='nhwhite'; 2='nh black'; 3='nh other';4='nh multirace'; 5='hispanic'; else=NA", as.factor = T)
#insurance
brfss_17$ins<-Recode(brfss_17$hlthpln1, recodes ="7:9=NA; 1=1;2=0")
#income grouping
brfss_17$inc<-ifelse(brfss_17$incomg==9, NA, brfss_17$incomg)
#education level
brfss_17$educ<-Recode(brfss_17$educa, recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA", as.factor=T)
brfss_17$educ<-relevel(brfss_17$educ, ref='2hsgrad')
#employment
brfss_17$employ<-Recode(brfss_17$employ1, recodes="1:2='Employed'; 2:6='nilf'; 7='retired'; 8='unable'; else=NA", as.factor=T)
brfss_17$employ<-relevel(brfss_17$employ, ref='Employed')
#marital status
brfss_17$marst<-Recode(brfss_17$marital, recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA", as.factor=T)
brfss_17$marst<-relevel(brfss_17$marst, ref='married')
#Age cut into intervals
brfss_17$agec<-cut(brfss_17$age80, breaks=c(0,24,39,59,79,99))
#BMI, in the brfss_17a the bmi variable has 2 implied decimal places, so we must divide by 100 to get real bmi's
brfss_17$bmi<-brfss_17$bmi5/100
brfss_17$obese<-ifelse(brfss_17$bmi>=30, 1, 0)
#smoking currently
brfss_17$smoke<-Recode(brfss_17$smoker3, recodes="1:2='Current'; 3='Former';4='NeverSmoked'; else=NA", as.factor=T)
brfss_17$smoke<-relevel(brfss_17$smoke, ref = "NeverSmoked")names(brfss_17)## [1] "dispcode" "statere1" "safetime" "hhadult" "genhlth" "physhlth"
## [7] "menthlth" "poorhlth" "hlthpln1" "persdoc2" "medcost" "checkup1"
## [13] "bphigh4" "bpmeds" "cholchk1" "toldhi2" "cholmed1" "cvdinfr4"
## [19] "cvdcrhd4" "cvdstrk3" "asthma3" "asthnow" "chcscncr" "chcocncr"
## [25] "chccopd1" "havarth3" "addepev2" "chckidny" "diabete3" "diabage2"
## [31] "lmtjoin3" "arthdis2" "arthsocl" "joinpai1" "sex" "marital"
## [37] "educa" "renthom1" "numhhol2" "numphon2" "cpdemo1a" "veteran3"
## [43] "employ1" "children" "income2" "internet" "weight2" "height3"
## [49] "pregnant" "deaf" "blind" "decide" "diffwalk" "diffdres"
## [55] "diffalon" "smoke100" "smokday2" "stopsmk2" "lastsmk2" "usenow3"
## [61] "ecigaret" "ecignow" "alcday5" "avedrnk2" "drnk3ge5" "maxdrnks"
## [67] "fruit2" "fruitju2" "fvgreen1" "frenchf1" "potatoe1" "vegetab2"
## [73] "exerany2" "exract11" "exeroft1" "exerhmm1" "exract21" "exeroft2"
## [79] "exerhmm2" "strength" "seatbelt" "flushot6" "flshtmy2" "pneuvac3"
## [85] "shingle2" "hivtst6" "hivtstd3" "hivrisk5" "casthdx2" "casthno2"
## [91] "callbckz" "wdusenow" "wdinftrk" "wdhowoft" "wdshare" "namtribe"
## [97] "namothr" "urbnrrl" "ststr" "impsex" "rfhlth" "phys14d"
## [103] "ment14d" "hcvu651" "rfhype5" "cholch1" "rfchol1" "michd"
## [109] "ltasth1" "casthm1" "asthms1" "drdxar1" "lmtact1" "lmtwrk1"
## [115] "lmtscl1" "prace1" "mrace1" "hispanc" "race" "raceg21"
## [121] "racegr3" "ageg5yr" "age65yr" "age80" "ageg" "wtkg3"
## [127] "bmi5" "bmi5cat" "rfbmi5" "educag" "incomg" "smoker3"
## [133] "rfsmok3" "ecigsts" "curecig" "drnkany5" "rfbing5" "drnkwek"
## [139] "rfdrhv5" "ftjuda2" "frutda2" "grenda1" "frnchda" "potada1"
## [145] "vegeda2" "misfrt1" "misveg1" "frtres1" "vegres1" "frutsu1"
## [151] "vegesu1" "frtlt1a" "veglt1a" "frt16a" "veg23a" "fruite1"
## [157] "vegete1" "totinda" "minac11" "minac21" "pacat1" "paindx1"
## [163] "pa150r2" "pa300r2" "pa30021" "pastrng" "parec1" "pastae1"
## [169] "rfseat2" "rfseat3" "flshot6" "pneumo2" "aidtst3" "mmsa"
## [175] "mmsawt" "seqno" "mmsaname" "badhealth" "male" "agec"
## [181] "black" "white" "other" "hispanic" "race_eth" "ins"
## [187] "inc" "educ" "employ" "marst" "bmi" "obese"
## [193] "smoke"library(car)
brfss_17$mentalhealth<-Recode(brfss_17$menthlth,
recodes="1:10=1;11:20=2;21:30=3; else=NA", as.factor = T)
brfss_17$mentalnum<-car::Recode(brfss_17$menthlth,
recodes="1:10=1;11:20=2;21:30=3; else=NA", as.factor = F)
library(dplyr)##
## 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, unionlibrary(survey)
sub<-brfss_17%>%
select(badhealth,mentalnum, mentalhealth, mmsaname, bmi,
agec,race_eth, marst, educ,white, black, hispanic,
other, smoke, ins, mmsawt, ststr) %>%
filter( complete.cases(.))
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~ststr, weights=~mmsawt, data =sub )fit.solr1<-svyolr(mentalhealth~race_eth+educ+agec,des)
summary(fit.solr1)## Call:
## svyolr(mentalhealth ~ race_eth + educ + agec, des)
##
## Coefficients:
## Value Std. Error t value
## race_ethnh black 0.05569935 0.06824710 0.8161423
## race_ethnh multirace 0.47828269 0.09367709 5.1056525
## race_ethnh other 0.02280236 0.10170216 0.2242073
## race_ethnhwhite 0.08666030 0.05560065 1.5586203
## educ0Prim 0.31686947 0.11774042 2.6912549
## educ1somehs 0.23693022 0.07322870 3.2354832
## educ3somecol -0.24522916 0.04344020 -5.6452120
## educ4colgrad -0.85704756 0.04151284 -20.6453622
## agec(24,39] 0.11589007 0.05466254 2.1201001
## agec(39,59] 0.29631728 0.05318091 5.5718733
## agec(59,79] 0.30443874 0.05835449 5.2170574
## agec(79,99] 0.19178270 0.11030487 1.7386603
##
## Intercepts:
## Value Std. Error t value
## 1|2 0.6429 0.0710 9.0542
## 2|3 1.4902 0.0707 21.0679fit.solr1$deviance+2*length(fit.solr1$coefficients)## [1] 113193.3Which gives us the results from the proportional odds model. In this case, we see that multiracial have greater odds of reporting fewer mental health days, compared to Hispanics, also NH Other respondents are less likely to report fewer mental health days, compared to Hispanics.
We also see that as education increases, the odds of reporting fewer mental health days increases.
As age increases, the odds of reporting fewer mental health days increases until the oldest age group. The oldest age group is more likely to report fewer meantl health days tan the other age groups.
ex1<-svyglm(I(mentalnum>1)~race_eth+educ+agec,des, family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!ex2<-svyglm(I(mentalnum>2)~race_eth+educ+agec,des, family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!library(VGAM)
#Proportional odds
fit.vgam<-vglm(as.ordered(mentalhealth)~race_eth+educ+agec,
brfss_17, weights =mmsawt/mean(mmsawt, na.rm=T),
family=cumulative(parallel = T, reverse = T)) #<-parallel = T == proportional odds
summary(fit.vgam)##
## Call:
## vglm(formula = as.ordered(mentalhealth) ~ race_eth + educ + agec,
## family = cumulative(parallel = T, reverse = T), data = brfss_17,
## weights = mmsawt/mean(mmsawt, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -6.635 -0.5365 -0.2636 0.39732 15.08
## logitlink(P[Y>=3]) -8.091 -0.3432 -0.1717 -0.06898 12.10
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.68699 0.02644 -25.982 < 2e-16 ***
## (Intercept):2 -1.52877 0.02699 -56.644 < 2e-16 ***
## race_ethnh black 0.07899 0.02734 2.890 0.00386 **
## race_ethnh multirace 0.57686 0.05551 10.392 < 2e-16 ***
## race_ethnh other 0.01607 0.03491 0.460 0.64535
## race_ethnhwhite 0.11794 0.02150 5.486 4.10e-08 ***
## educ0Prim 0.27290 0.03932 6.940 3.91e-12 ***
## educ1somehs 0.27767 0.02741 10.132 < 2e-16 ***
## educ3somecol -0.22187 0.01891 -11.731 < 2e-16 ***
## educ4colgrad -0.84339 0.02154 -39.153 < 2e-16 ***
## agec(24,39] 0.12183 0.02298 5.302 1.15e-07 ***
## agec(39,59] 0.29705 0.02242 13.247 < 2e-16 ***
## agec(59,79] 0.30533 0.02563 11.915 < 2e-16 ***
## agec(79,99] 0.24343 0.05188 4.692 2.71e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 136179.7 on 146000 degrees of freedom
##
## Log-likelihood: -68089.86 on 146000 degrees of freedom
##
## Number of Fisher scoring iterations: 4
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## race_ethnh black race_ethnh multirace race_ethnh other
## 1.0821988 1.7804368 1.0161948
## race_ethnhwhite educ0Prim educ1somehs
## 1.1251737 1.3137709 1.3200556
## educ3somecol educ4colgrad agec(24,39]
## 0.8010156 0.4302509 1.1295570
## agec(39,59] agec(59,79] agec(79,99]
## 1.3458862 1.3570762 1.2756107#Non-proportional odds
fit.vgam2<-vglm(as.ordered(mentalhealth)~race_eth+educ+agec,brfss_17,
weights =mmsawt/mean(mmsawt, na.rm=T),
family=cumulative(parallel = F, reverse = T)) #<-parallel = F == Nonproportional odds
summary(fit.vgam2)##
## Call:
## vglm(formula = as.ordered(mentalhealth) ~ race_eth + educ + agec,
## family = cumulative(parallel = F, reverse = T), data = brfss_17,
## weights = mmsawt/mean(mmsawt, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logitlink(P[Y>=2]) -6.576 -0.5335 -0.2611 0.40590 14.98
## logitlink(P[Y>=3]) -8.390 -0.3443 -0.1718 -0.06958 12.47
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.639957 0.026910 -23.781 < 2e-16 ***
## (Intercept):2 -1.718692 0.034045 -50.483 < 2e-16 ***
## race_ethnh black:1 0.086059 0.028028 3.071 0.00214 **
## race_ethnh black:2 0.050762 0.033979 1.494 0.13519
## race_ethnh multirace:1 0.619293 0.057769 10.720 < 2e-16 ***
## race_ethnh multirace:2 0.491846 0.067482 7.289 3.13e-13 ***
## race_ethnh other:1 0.007571 0.035543 0.213 0.83131
## race_ethnh other:2 0.022383 0.044702 0.501 0.61657
## race_ethnhwhite:1 0.111544 0.022016 5.066 4.05e-07 ***
## race_ethnhwhite:2 0.135368 0.026754 5.060 4.20e-07 ***
## educ0Prim:1 0.258441 0.041078 6.291 3.14e-10 ***
## educ0Prim:2 0.277231 0.046184 6.003 1.94e-09 ***
## educ1somehs:1 0.285296 0.028607 9.973 < 2e-16 ***
## educ1somehs:2 0.256917 0.032443 7.919 2.39e-15 ***
## educ3somecol:1 -0.223620 0.019449 -11.498 < 2e-16 ***
## educ3somecol:2 -0.218560 0.023162 -9.436 < 2e-16 ***
## educ4colgrad:1 -0.826545 0.021924 -37.700 < 2e-16 ***
## educ4colgrad:2 -0.914026 0.027689 -33.010 < 2e-16 ***
## agec(24,39]:1 0.071359 0.023336 3.058 0.00223 **
## agec(24,39]:2 0.324652 0.030159 10.765 < 2e-16 ***
## agec(39,59]:1 0.238156 0.022837 10.428 < 2e-16 ***
## agec(39,59]:2 0.524203 0.029185 17.962 < 2e-16 ***
## agec(59,79]:1 0.224362 0.026223 8.556 < 2e-16 ***
## agec(59,79]:2 0.591832 0.032564 18.174 < 2e-16 ***
## agec(79,99]:1 0.161944 0.053618 3.020 0.00253 **
## agec(79,99]:2 0.528144 0.062564 8.442 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 135899.8 on 145988 degrees of freedom
##
## Log-likelihood: -67949.88 on 145988 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## race_ethnh black:1 race_ethnh black:2 race_ethnh multirace:1
## 1.0898709 1.0520729 1.8576151
## race_ethnh multirace:2 race_ethnh other:1 race_ethnh other:2
## 1.6353323 1.0076002 1.0226350
## race_ethnhwhite:1 race_ethnhwhite:2 educ0Prim:1
## 1.1180028 1.1449577 1.2949094
## educ0Prim:2 educ1somehs:1 educ1somehs:2
## 1.3194707 1.3301550 1.2929380
## educ3somecol:1 educ3somecol:2 educ4colgrad:1
## 0.7996186 0.8036756 0.4375585
## educ4colgrad:2 agec(24,39]:1 agec(24,39]:2
## 0.4009069 1.0739672 1.3835495
## agec(39,59]:1 agec(39,59]:2 agec(59,79]:1
## 1.2689071 1.6891123 1.2515237
## agec(59,79]:2 agec(79,99]:1 agec(79,99]:2
## 1.8072964 1.1757938 1.6957816AIC(fit.vgam)## [1] 136207.7AIC(fit.vgam2)## [1] 135951.8mfit<-vglm(mentalhealth~race_eth+educ+agec,
family=multinomial(refLevel = 1),
data=brfss_17,
weights=mmsawt/mean(mmsawt, na.rm=T))
summary(mfit)##
## Call:
## vglm(formula = mentalhealth ~ race_eth + educ + agec, family = multinomial(refLevel = 1),
## data = brfss_17, weights = mmsawt/mean(mmsawt, na.rm = T))
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## log(mu[,2]/mu[,1]) -4.498 -0.4097 -0.2265 -0.09901 16.57
## log(mu[,3]/mu[,1]) -5.511 -0.4375 -0.2411 -0.09478 12.97
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.22192 0.03439 -35.534 < 2e-16 ***
## (Intercept):2 -1.46838 0.03500 -41.949 < 2e-16 ***
## race_ethnh black:1 0.08845 0.03674 2.407 0.0161 *
## race_ethnh black:2 0.08550 0.03505 2.439 0.0147 *
## race_ethnh multirace:1 0.59319 0.07297 8.129 4.34e-16 ***
## race_ethnh multirace:2 0.65484 0.07082 9.246 < 2e-16 ***
## race_ethnh other:1 -0.03293 0.04649 -0.708 0.4788
## race_ethnh other:2 0.04419 0.04570 0.967 0.3336
## race_ethnhwhite:1 0.05605 0.02893 1.937 0.0527 .
## race_ethnhwhite:2 0.15477 0.02762 5.603 2.11e-08 ***
## educ0Prim:1 0.16435 0.05593 2.938 0.0033 **
## educ0Prim:2 0.31997 0.04815 6.645 3.03e-11 ***
## educ1somehs:1 0.24401 0.03788 6.441 1.19e-10 ***
## educ1somehs:2 0.31753 0.03395 9.352 < 2e-16 ***
## educ3somecol:1 -0.17369 0.02588 -6.710 1.94e-11 ***
## educ3somecol:2 -0.26076 0.02391 -10.908 < 2e-16 ***
## educ4colgrad:1 -0.60058 0.02890 -20.778 < 2e-16 ***
## educ4colgrad:2 -1.03017 0.02831 -36.389 < 2e-16 ***
## agec(24,39]:1 -0.14759 0.02961 -4.984 6.22e-07 ***
## agec(24,39]:2 0.29726 0.03086 9.632 < 2e-16 ***
## agec(39,59]:1 -0.04444 0.02908 -1.528 0.1265
## agec(39,59]:2 0.51970 0.02993 17.365 < 2e-16 ***
## agec(59,79]:1 -0.13503 0.03450 -3.914 9.09e-05 ***
## agec(59,79]:2 0.56516 0.03343 16.905 < 2e-16 ***
## agec(79,99]:1 -0.18893 0.07413 -2.549 0.0108 *
## agec(79,99]:2 0.48866 0.06447 7.580 3.45e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
##
## Residual deviance: 135891.7 on 145988 degrees of freedom
##
## Log-likelihood: -67945.83 on 145988 degrees of freedom
##
## Number of Fisher scoring iterations: 4
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Reference group is level 1 of the responseround(exp(coef(mfit)), 3)## (Intercept):1 (Intercept):2 race_ethnh black:1
## 0.295 0.230 1.092
## race_ethnh black:2 race_ethnh multirace:1 race_ethnh multirace:2
## 1.089 1.810 1.925
## race_ethnh other:1 race_ethnh other:2 race_ethnhwhite:1
## 0.968 1.045 1.058
## race_ethnhwhite:2 educ0Prim:1 educ0Prim:2
## 1.167 1.179 1.377
## educ1somehs:1 educ1somehs:2 educ3somecol:1
## 1.276 1.374 0.841
## educ3somecol:2 educ4colgrad:1 educ4colgrad:2
## 0.770 0.548 0.357
## agec(24,39]:1 agec(24,39]:2 agec(39,59]:1
## 0.863 1.346 0.957
## agec(39,59]:2 agec(59,79]:1 agec(59,79]:2
## 1.682 0.874 1.760
## agec(79,99]:1 agec(79,99]:2
## 0.828 1.630round(exp(confint(mfit)), 3)## 2.5 % 97.5 %
## (Intercept):1 0.275 0.315
## (Intercept):2 0.215 0.247
## race_ethnh black:1 1.017 1.174
## race_ethnh black:2 1.017 1.167
## race_ethnh multirace:1 1.569 2.088
## race_ethnh multirace:2 1.675 2.211
## race_ethnh other:1 0.883 1.060
## race_ethnh other:2 0.956 1.143
## race_ethnhwhite:1 0.999 1.119
## race_ethnhwhite:2 1.106 1.232
## educ0Prim:1 1.056 1.315
## educ0Prim:2 1.253 1.513
## educ1somehs:1 1.185 1.375
## educ1somehs:2 1.285 1.468
## educ3somecol:1 0.799 0.884
## educ3somecol:2 0.735 0.807
## educ4colgrad:1 0.518 0.580
## educ4colgrad:2 0.338 0.377
## agec(24,39]:1 0.814 0.914
## agec(24,39]:2 1.267 1.430
## agec(39,59]:1 0.904 1.013
## agec(39,59]:2 1.586 1.783
## agec(59,79]:1 0.817 0.935
## agec(59,79]:2 1.648 1.879
## agec(79,99]:1 0.716 0.957
## agec(79,99]:2 1.437 1.850dat<-expand.grid(race_eth=levels(brfss_17$race_eth),
educ=levels(brfss_17$educ), agec=levels(brfss_17$agec))
fitm<-predict(mfit, newdat=dat,type="response")
dat<-cbind(dat, round(fitm,3))
head(dat, n=20)## race_eth educ agec 1 2 3
## 1 hispanic 2hsgrad (0,24] 0.656 0.193 0.151
## 2 nh black 2hsgrad (0,24] 0.636 0.205 0.159
## 3 nh multirace 2hsgrad (0,24] 0.506 0.270 0.224
## 4 nh other 2hsgrad (0,24] 0.655 0.187 0.158
## 5 nhwhite 2hsgrad (0,24] 0.633 0.197 0.170
## 6 hispanic 0Prim (0,24] 0.601 0.209 0.191
## 7 nh black 0Prim (0,24] 0.580 0.220 0.200
## 8 nh multirace 0Prim (0,24] 0.447 0.281 0.273
## 9 nh other 0Prim (0,24] 0.600 0.202 0.199
## 10 nhwhite 0Prim (0,24] 0.576 0.211 0.213
## 11 hispanic 1somehs (0,24] 0.591 0.222 0.187
## 12 nh black 1somehs (0,24] 0.570 0.234 0.196
## 13 nh multirace 1somehs (0,24] 0.437 0.297 0.266
## 14 nh other 1somehs (0,24] 0.590 0.215 0.195
## 15 nhwhite 1somehs (0,24] 0.566 0.225 0.209
## 16 hispanic 3somecol (0,24] 0.702 0.174 0.125
## 17 nh black 3somecol (0,24] 0.683 0.185 0.132
## 18 nh multirace 3somecol (0,24] 0.559 0.250 0.191
## 19 nh other 3somecol (0,24] 0.702 0.168 0.130
## 20 nhwhite 3somecol (0,24] 0.681 0.178 0.141fitted.ord<-round(predict(fit.vgam, newdat=dat[,1:3], type="response"), 3)
dat<-cbind(dat, fitted.ord)
names(dat)<-c(names(dat)[1:3], "mp1", "mp2", "mp3", "op1", "op2", "op3")
head(dat, n=20)## race_eth educ agec mp1 mp2 mp3 op1 op2 op3
## 1 hispanic 2hsgrad (0,24] 0.656 0.193 0.151 0.665 0.157 0.178
## 2 nh black 2hsgrad (0,24] 0.636 0.205 0.159 0.647 0.162 0.190
## 3 nh multirace 2hsgrad (0,24] 0.506 0.270 0.224 0.528 0.194 0.278
## 4 nh other 2hsgrad (0,24] 0.655 0.187 0.158 0.662 0.158 0.181
## 5 nhwhite 2hsgrad (0,24] 0.633 0.197 0.170 0.639 0.165 0.196
## 6 hispanic 0Prim (0,24] 0.601 0.209 0.191 0.602 0.176 0.222
## 7 nh black 0Prim (0,24] 0.580 0.220 0.200 0.583 0.181 0.236
## 8 nh multirace 0Prim (0,24] 0.447 0.281 0.273 0.459 0.204 0.336
## 9 nh other 0Prim (0,24] 0.600 0.202 0.199 0.598 0.177 0.224
## 10 nhwhite 0Prim (0,24] 0.576 0.211 0.213 0.574 0.184 0.243
## 11 hispanic 1somehs (0,24] 0.591 0.222 0.187 0.601 0.177 0.223
## 12 nh black 1somehs (0,24] 0.570 0.234 0.196 0.582 0.182 0.236
## 13 nh multirace 1somehs (0,24] 0.437 0.297 0.266 0.458 0.204 0.338
## 14 nh other 1somehs (0,24] 0.590 0.215 0.195 0.597 0.178 0.225
## 15 nhwhite 1somehs (0,24] 0.566 0.225 0.209 0.572 0.184 0.244
## 16 hispanic 3somecol (0,24] 0.702 0.174 0.125 0.713 0.139 0.148
## 17 nh black 3somecol (0,24] 0.683 0.185 0.132 0.696 0.145 0.158
## 18 nh multirace 3somecol (0,24] 0.559 0.250 0.191 0.582 0.182 0.236
## 19 nh other 3somecol (0,24] 0.702 0.168 0.130 0.709 0.141 0.150
## 20 nhwhite 3somecol (0,24] 0.681 0.178 0.141 0.688 0.149 0.163AIC(fit.vgam) #proportional odds## [1] 136207.7AIC(fit.vgam2) #cumulative logit, non proportional## [1] 135951.8AIC(mfit) #multinomial## [1] 135943.7