## [1] "dispcode" "statere1" "safetime" "hhadult" "genhlth" "physhlth"
## [7] "menthlth" "poorhlth" "hlthpln1" "persdoc2"#education level
brfss_17$educ<-Recode(brfss_17$educa, recodes="1:2='1Prim'; 3='2Somehs'; 4='3Hsgrad'; 5='4Somecol'; 6='5Colgrad';9=NA",
as.factor=T)
# brfss_17$educ<-relevel(brfss_17$educ, ref=3Hsgrad)
brfss_17$educnum<- car::recode(brfss_17$educa, recodes="1:2=1;3=2;4=3;5=4;6=5;else=NA", as.factor= F)
#First we tell R our survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~ststr, weights=~mmsawt, data =brfss_17 )Research Question: Does age, race/ethnicity, income level, and occupation affect educational outcomes?
To fit an ordinal logit to survey data in R, we use the svyolr function in the survey library.
#Here I fit three nested models for the education outcome
fit.solr1<-svyolr(educ~race_eth+agec+inc,des)
summary(fit.solr1)## Call:
## svyolr(educ ~ race_eth + agec + inc, des)
##
## Coefficients:
## Value Std. Error t value
## race_ethnh black 0.9180126 0.037849616 24.254211
## race_ethnh multirace 1.0278742 0.060479319 16.995465
## race_ethnh other 1.7677114 0.063066669 28.029250
## race_ethnhwhite 1.0547407 0.032568201 32.385599
## agec(24,39] 0.4681750 0.037647992 12.435590
## agec(39,59] 0.1391676 0.036550078 3.807587
## agec(59,79] 0.1414724 0.036823518 3.841904
## agec(79,99] -0.1147814 0.052531250 -2.185012
## inc 0.5686744 0.007432816 76.508608
##
## Intercepts:
## Value Std. Error t value
## 1Prim|2Somehs -0.3603 0.0515 -7.0011
## 2Somehs|3Hsgrad 0.8442 0.0461 18.2938
## 3Hsgrad|4Somecol 2.6278 0.0465 56.5646
## 4Somecol|5Colgrad 4.1966 0.0479 87.6931
## (42457 observations deleted due to missingness)The results from the proportional odds models shows that NH whites have lower odds of reporting worst education, compared to Hispanics, also NH Other respondents are less likely to report worse education, compared to Hispanics.
I also see that as age increases, the odds of reporting worse education decreases, compared to those < 25 years of age.
ex1<-svyglm(I(educnum>1)~race_eth+inc+agec,des, family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!ex2<-svyglm(I(educnum>2)~race_eth+inc+agec,des, family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!#Just print the odds ratios
round(exp(rbind(coef(ex1)[-1],
coef(ex2)[-1])),3)## race_ethnh black race_ethnh multirace race_ethnh other race_ethnhwhite
## [1,] 13.733 13.245 5.736 13.697
## [2,] 4.282 5.233 5.411 5.125
## inc agec(24,39] agec(39,59] agec(59,79] agec(79,99]
## [1,] 1.820 0.302 0.17 0.141 0.067
## [2,] 1.886 0.619 0.41 0.435 0.323I will say the proportional odds assumption are met as most of the relationship between each pair of outcome groups are the same across transitions.
I evaluated the proportional odd assumption by fitting a binomial logistic regressions and examining the relationship between each pair of outcome groups to if they are the same across transitions.
#Proportional odds
fit.vgam<-vglm(as.ordered(educ)~race_eth+inc+agec,
brfss_17, weights =mmsawt/mean(mmsawt, na.rm=T),
family=cumulative(parallel = T, reverse = T), subset = is.na(educ)==F) #<-parallel = T == proportional odds
summary(fit.vgam)##
## Call:
## vglm(formula = as.ordered(educ) ~ race_eth + inc + agec, family = cumulative(parallel = T,
## reverse = T), data = brfss_17, weights = mmsawt/mean(mmsawt,
## na.rm = T), subset = is.na(educ) == F)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.360200 0.019151 18.808 < 2e-16 ***
## (Intercept):2 -0.844170 0.017724 -47.629 < 2e-16 ***
## (Intercept):3 -2.627778 0.018675 -140.714 < 2e-16 ***
## (Intercept):4 -4.196631 0.020081 -208.980 < 2e-16 ***
## race_ethnh black 0.918006 0.015203 60.382 < 2e-16 ***
## race_ethnh multirace 1.027893 0.037785 27.204 < 2e-16 ***
## race_ethnh other 1.767692 0.019616 90.116 < 2e-16 ***
## race_ethnhwhite 1.054732 0.012154 86.781 < 2e-16 ***
## inc 0.568680 0.003344 170.045 < 2e-16 ***
## agec(24,39] 0.468175 0.015324 30.553 < 2e-16 ***
## agec(39,59] 0.139173 0.014841 9.378 < 2e-16 ***
## agec(59,79] 0.141475 0.015760 8.977 < 2e-16 ***
## agec(79,99] -0.114773 0.025976 -4.418 9.94e-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]),
## logitlink(P[Y>=4]), logitlink(P[Y>=5])
##
## Residual deviance: 475287.5 on 753659 degrees of freedom
##
## Log-likelihood: -237643.8 on 753659 degrees of freedom
##
## Number of Fisher scoring iterations: 6
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## race_ethnh black race_ethnh multirace race_ethnh other
## 2.5042916 2.7951696 5.8573184
## race_ethnhwhite inc agec(24,39]
## 2.8712052 1.7659342 1.5970775
## agec(39,59] agec(59,79] agec(79,99]
## 1.1493229 1.1519721 0.8915687#Non-proportional odds
fit.vgam2<-vglm(as.ordered(educ)~race_eth+inc+agec,brfss_17,
weights =mmsawt/mean(mmsawt, na.rm=T),
family=cumulative(parallel = F, reverse = T), subset = is.na(educ)==F) #<-parallel = F == Nonproportional odds## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1
## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y = y, extra = extra): It seems
## that the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals =
## residuals, : fitted values close to 0 or 1## Warning in slot(family, "validparams")(eta, y, extra = extra): It seems that
## the nonparallelism assumption has resulted in intersecting linear/additive
## predictors. Try propodds() or fitting a partial nonproportional odds model or
## choosing some other link function, etc.## Warning in vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, :
## iterations terminated because half-step sizes are very small## Warning in vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2, : some
## quantities such as z, residuals, SEs may be inaccurate due to convergence at a
## half-step## Warning in log(prob): NaNs producedsummary(fit.vgam2)## Warning in matrix.power(wz, M = M, power = 0.5, fast = TRUE): Some weight
## matrices have negative eigenvalues. They will be assigned NAs##
## Call:
## vglm(formula = as.ordered(educ) ~ race_eth + inc + agec, family = cumulative(parallel = F,
## reverse = T), data = brfss_17, weights = mmsawt/mean(mmsawt,
## na.rm = T), subset = is.na(educ) == F)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -8.362e-01 2.579e-02 -3.243e+01 < 2e-16 ***
## (Intercept):2 -1.377e+00 2.151e-02 -6.402e+01 < 2e-16 ***
## (Intercept):3 -2.019e+00 1.634e-02 -1.236e+02 < 2e-16 ***
## (Intercept):4 -3.529e+00 2.447e-02 -1.442e+02 < 2e-16 ***
## race_ethnh black:1 3.716e+00 3.378e-02 1.100e+02 < 2e-16 ***
## race_ethnh black:2 2.011e+00 1.871e-02 1.075e+02 < 2e-16 ***
## race_ethnh black:3 6.488e-01 1.454e-02 4.463e+01 < 2e-16 ***
## race_ethnh black:4 1.895e-01 1.887e-02 1.004e+01 < 2e-16 ***
## race_ethnh multirace:1 3.480e+00 1.090e-01 3.193e+01 < 2e-16 ***
## race_ethnh multirace:2 2.075e+00 6.545e-02 3.171e+01 < 2e-16 ***
## race_ethnh multirace:3 8.706e-01 4.454e-02 1.955e+01 < 2e-16 ***
## race_ethnh multirace:4 2.423e-01 4.586e-02 5.283e+00 1.27e-07 ***
## race_ethnh other:1 3.079e+00 4.262e-02 7.225e+01 < 2e-16 ***
## race_ethnh other:2 2.043e+00 2.963e-02 6.894e+01 < 2e-16 ***
## race_ethnh other:3 1.137e+00 2.129e-02 5.340e+01 < 2e-16 ***
## race_ethnh other:4 1.318e+00 2.162e-02 6.096e+01 < 2e-16 ***
## race_ethnhwhite:1 3.443e+00 1.753e-02 1.964e+02 < 2e-16 ***
## race_ethnhwhite:2 2.010e+00 1.107e-02 1.815e+02 < 2e-16 ***
## race_ethnhwhite:3 7.414e-01 9.791e-03 7.572e+01 < 2e-16 ***
## race_ethnhwhite:4 4.272e-01 1.418e-02 3.013e+01 < 2e-16 ***
## inc:1 5.777e-01 9.021e-09 6.404e+07 < 2e-16 ***
## inc:2 6.224e-01 8.246e-09 7.548e+07 < 2e-16 ***
## inc:3 4.782e-01 2.201e-03 2.173e+02 < 2e-16 ***
## inc:4 4.686e-01 3.878e-03 1.208e+02 < 2e-16 ***
## agec(24,39]:1 -8.053e-01 2.744e-02 -2.935e+01 < 2e-16 ***
## agec(24,39]:2 -4.446e-01 2.346e-02 -1.895e+01 < 2e-16 ***
## agec(24,39]:3 3.316e-01 1.801e-02 1.841e+01 < 2e-16 ***
## agec(24,39]:4 8.769e-01 2.030e-02 4.319e+01 < 2e-16 ***
## agec(39,59]:1 -1.334e+00 2.579e-02 -5.174e+01 < 2e-16 ***
## agec(39,59]:2 -7.955e-01 2.151e-02 -3.698e+01 < 2e-16 ***
## agec(39,59]:3 6.650e-02 1.606e-02 4.141e+00 3.45e-05 ***
## agec(39,59]:4 6.694e-01 1.958e-02 3.418e+01 < 2e-16 ***
## agec(59,79]:1 -1.407e+00 2.579e-02 -5.455e+01 < 2e-16 ***
## agec(59,79]:2 -7.011e-01 2.151e-02 -3.259e+01 < 2e-16 ***
## agec(59,79]:3 1.248e-01 1.606e-02 7.773e+00 7.64e-15 ***
## agec(59,79]:4 6.106e-01 2.045e-02 2.985e+01 < 2e-16 ***
## agec(79,99]:1 -2.090e+00 2.579e-02 -8.107e+01 < 2e-16 ***
## agec(79,99]:2 -9.192e-01 2.151e-02 -4.273e+01 < 2e-16 ***
## agec(79,99]:3 -1.798e-01 2.458e-02 -7.313e+00 2.60e-13 ***
## agec(79,99]:4 4.678e-01 3.262e-02 1.434e+01 < 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]),
## logitlink(P[Y>=4]), logitlink(P[Y>=5])
##
## Residual deviance: 3177180 on 753632 degrees of freedom
##
## Log-likelihood: NA on 753632 degrees of freedom
##
## Number of Fisher scoring iterations: 2
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):2', '(Intercept):3', '(Intercept):4', 'race_ethnh black:1', 'race_ethnh multirace:1', 'race_ethnh other:1', 'race_ethnhwhite:1', 'inc:1', 'inc:3', 'agec(24,39]:2', 'agec(39,59]:2', 'agec(79,99]:2'
##
##
## Exponentiated coefficients:
## race_ethnh black:1 race_ethnh black:2 race_ethnh black:3
## 41.1165639 7.4704991 1.9131880
## race_ethnh black:4 race_ethnh multirace:1 race_ethnh multirace:2
## 1.2086341 32.4501030 7.9679311
## race_ethnh multirace:3 race_ethnh multirace:4 race_ethnh other:1
## 2.3882294 1.2741329 21.7439010
## race_ethnh other:2 race_ethnh other:3 race_ethnh other:4
## 7.7121089 3.1171439 3.7363678
## race_ethnhwhite:1 race_ethnhwhite:2 race_ethnhwhite:3
## 31.2687098 7.4605636 2.0989327
## race_ethnhwhite:4 inc:1 inc:2
## 1.5329497 1.7819707 1.8633837
## inc:3 inc:4 agec(24,39]:1
## 1.6131846 1.5977073 0.4469330
## agec(24,39]:2 agec(24,39]:3 agec(24,39]:4
## 0.6410685 1.3931484 2.4034621
## agec(39,59]:1 agec(39,59]:2 agec(39,59]:3
## 0.2633511 0.4513489 1.0687644
## agec(39,59]:4 agec(59,79]:1 agec(59,79]:2
## 1.9530743 0.2449760 0.4960511
## agec(59,79]:3 agec(59,79]:4 agec(79,99]:1
## 1.1329550 1.8414756 0.1236303
## agec(79,99]:2 agec(79,99]:3 agec(79,99]:4
## 0.3988462 0.8354572 1.5964361mfit<-vglm(educ~race_eth+inc+agec,
family=multinomial(refLevel = 1),
data=brfss_17,
weights=mmsawt/mean(mmsawt, na.rm=T))
summary(mfit) ##
## Call:
## vglm(formula = educ ~ race_eth + inc + agec, family = multinomial(refLevel = 1),
## data = brfss_17, weights = mmsawt/mean(mmsawt, na.rm = T))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.84615 0.06800 12.444 < 2e-16 ***
## (Intercept):2 0.92878 0.06396 14.522 < 2e-16 ***
## (Intercept):3 0.15489 0.06466 2.396 0.0166 *
## (Intercept):4 -3.10362 0.07034 -44.125 < 2e-16 ***
## race_ethnh black:1 1.95982 0.05598 35.012 < 2e-16 ***
## race_ethnh black:2 2.66060 0.05293 50.266 < 2e-16 ***
## race_ethnh black:3 2.89517 0.05322 54.405 < 2e-16 ***
## race_ethnh black:4 2.89161 0.05537 52.227 < 2e-16 ***
## race_ethnh multirace:1 1.68952 0.19337 8.737 < 2e-16 ***
## race_ethnh multirace:2 2.49043 0.18035 13.808 < 2e-16 ***
## race_ethnh multirace:3 2.99992 0.17939 16.722 < 2e-16 ***
## race_ethnh multirace:4 2.95400 0.18240 16.195 < 2e-16 ***
## race_ethnh other:1 0.51196 0.07431 6.889 5.6e-12 ***
## race_ethnh other:2 1.41874 0.06265 22.645 < 2e-16 ***
## race_ethnh other:3 1.76311 0.06228 28.310 < 2e-16 ***
## race_ethnh other:4 2.96407 0.06270 47.271 < 2e-16 ***
## race_ethnhwhite:1 1.74351 0.03919 44.486 < 2e-16 ***
## race_ethnhwhite:2 2.60298 0.03610 72.111 < 2e-16 ***
## race_ethnhwhite:3 2.82946 0.03633 77.878 < 2e-16 ***
## race_ethnhwhite:4 3.09691 0.03790 81.705 < 2e-16 ***
## inc:1 0.10293 0.01120 9.188 < 2e-16 ***
## inc:2 0.46061 0.01009 45.657 < 2e-16 ***
## inc:3 0.67625 0.01013 66.737 < 2e-16 ***
## inc:4 1.15267 0.01081 106.646 < 2e-16 ***
## agec(24,39]:1 -0.97971 0.06927 -14.144 < 2e-16 ***
## agec(24,39]:2 -1.47733 0.06503 -22.718 < 2e-16 ***
## agec(24,39]:3 -1.33310 0.06522 -20.441 < 2e-16 ***
## agec(24,39]:4 -0.12449 0.06761 -1.841 0.0656 .
## agec(39,59]:1 -1.34157 0.06755 -19.860 < 2e-16 ***
## agec(39,59]:2 -1.98198 0.06354 -31.191 < 2e-16 ***
## agec(39,59]:3 -2.03549 0.06381 -31.901 < 2e-16 ***
## agec(39,59]:4 -0.96744 0.06619 -14.616 < 2e-16 ***
## agec(59,79]:1 -1.65810 0.07061 -23.481 < 2e-16 ***
## agec(59,79]:2 -2.17389 0.06606 -32.906 < 2e-16 ***
## agec(59,79]:3 -2.12918 0.06628 -32.125 < 2e-16 ***
## agec(59,79]:4 -1.15699 0.06874 -16.832 < 2e-16 ***
## agec(79,99]:1 -2.34442 0.09122 -25.701 < 2e-16 ***
## agec(79,99]:2 -2.74598 0.08200 -33.490 < 2e-16 ***
## agec(79,99]:3 -3.01573 0.08318 -36.254 < 2e-16 ***
## agec(79,99]:4 -2.03067 0.08701 -23.339 < 2e-16 ***
## ---
## 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]),
## log(mu[,4]/mu[,1]), log(mu[,5]/mu[,1])
##
## Residual deviance: 464429.5 on 753632 degrees of freedom
##
## Log-likelihood: -232214.8 on 753632 degrees of freedom
##
## Number of Fisher scoring iterations: 7
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Reference group is level 1 of the responseAIC(fit.vgam) #proportional odds## [1] 475313.5AIC(mfit) #multinomial## [1] 464509.5The multinomial model is the best model compared to the proportional odd model because it has the lowest AIC though the AIC’s are very close.
tab<- round(exp(cbind(OR = coef(mfit), confint(mfit))),3)
colnames(tab) <- c("Oddratio", "CI(2.5%)", "CI(97.5%)")
knitr::kable(tab)| Oddratio | CI(2.5%) | CI(97.5%) | |
|---|---|---|---|
| (Intercept):1 | 2.331 | 2.040 | 2.663 |
| (Intercept):2 | 2.531 | 2.233 | 2.869 |
| (Intercept):3 | 1.168 | 1.029 | 1.325 |
| (Intercept):4 | 0.045 | 0.039 | 0.052 |
| race_ethnh black:1 | 7.098 | 6.360 | 7.921 |
| race_ethnh black:2 | 14.305 | 12.895 | 15.869 |
| race_ethnh black:3 | 18.087 | 16.295 | 20.075 |
| race_ethnh black:4 | 18.022 | 16.169 | 20.088 |
| race_ethnh multirace:1 | 5.417 | 3.708 | 7.913 |
| race_ethnh multirace:2 | 12.066 | 8.473 | 17.183 |
| race_ethnh multirace:3 | 20.084 | 14.130 | 28.546 |
| race_ethnh multirace:4 | 19.182 | 13.417 | 27.426 |
| race_ethnh other:1 | 1.669 | 1.442 | 1.930 |
| race_ethnh other:2 | 4.132 | 3.654 | 4.672 |
| race_ethnh other:3 | 5.831 | 5.161 | 6.587 |
| race_ethnh other:4 | 19.377 | 17.136 | 21.910 |
| race_ethnhwhite:1 | 5.717 | 5.295 | 6.174 |
| race_ethnhwhite:2 | 13.504 | 12.582 | 14.494 |
| race_ethnhwhite:3 | 16.936 | 15.772 | 18.186 |
| race_ethnhwhite:4 | 22.129 | 20.545 | 23.836 |
| inc:1 | 1.108 | 1.084 | 1.133 |
| inc:2 | 1.585 | 1.554 | 1.617 |
| inc:3 | 1.966 | 1.928 | 2.006 |
| inc:4 | 3.167 | 3.100 | 3.234 |
| agec(24,39]:1 | 0.375 | 0.328 | 0.430 |
| agec(24,39]:2 | 0.228 | 0.201 | 0.259 |
| agec(24,39]:3 | 0.264 | 0.232 | 0.300 |
| agec(24,39]:4 | 0.883 | 0.773 | 1.008 |
| agec(39,59]:1 | 0.261 | 0.229 | 0.298 |
| agec(39,59]:2 | 0.138 | 0.122 | 0.156 |
| agec(39,59]:3 | 0.131 | 0.115 | 0.148 |
| agec(39,59]:4 | 0.380 | 0.334 | 0.433 |
| agec(59,79]:1 | 0.191 | 0.166 | 0.219 |
| agec(59,79]:2 | 0.114 | 0.100 | 0.129 |
| agec(59,79]:3 | 0.119 | 0.104 | 0.135 |
| agec(59,79]:4 | 0.314 | 0.275 | 0.360 |
| agec(79,99]:1 | 0.096 | 0.080 | 0.115 |
| agec(79,99]:2 | 0.064 | 0.055 | 0.075 |
| agec(79,99]:3 | 0.049 | 0.042 | 0.058 |
| agec(79,99]:4 | 0.131 | 0.111 | 0.156 |