Married example
Szu-Yu Chen
29 Oct 2020
1 Data Management
Load data
dtaM<-read.table("C:/Users/ASUS/Desktop/data/Married.txt", header= T)
head(dtaM)## 教育 性別 婚姻 年齡
## 1 國小以下 女 未婚 32.5
## 2 專科 男 有偶 57.5
## 3 國中 男 有偶 42.5
## 4 國小以下 男 有偶 42.5
## 5 高中 女 喪偶 62.5
## 6 高中 女 有偶 52.5dtaM$教育<-factor(dtaM$教育, levels= c('國小以下','國中','高中','專科','大學以上'))
dtaM$婚姻<-factor(dtaM$婚姻, levels= c('有偶','未婚','離婚','喪偶'))
dtaM$年次<-(109-dtaM$年齡)
dtaM$結過婚<-ifelse(dtaM$婚姻=='未婚', 0, 1)
apply(dtaM, 2, function(x) table(x)/dim(dtaM)[1])## $教育
## x
## 大學以上 高中 國小以下 國中 專科
## 0.1070 0.3297 0.2406 0.1360 0.1867
##
## $性別
## x
## 女 男
## 0.5073 0.4927
##
## $婚姻
## x
## 未婚 有偶 喪偶 離婚
## 0.1966 0.6550 0.0341 0.1143
##
## $年齡
## x
## 32.5 37.5 42.5 47.5 52.5 57.5 62.5
## 0.1619 0.1516 0.1435 0.1434 0.1558 0.1306 0.1132
##
## $年次
## x
## 46.5 51.5 56.5 61.5 66.5 71.5 76.5
## 0.1132 0.1306 0.1558 0.1434 0.1435 0.1516 0.1619
##
## $結過婚
## x
## 0 1
## 0.1966 0.8034str(dtaM)## 'data.frame': 10000 obs. of 6 variables:
## $ 教育 : Factor w/ 5 levels "國小以下","國中",..: 1 4 2 1 3 3 1 2 3 1 ...
## $ 性別 : chr "女" "男" "男" "男" ...
## $ 婚姻 : Factor w/ 4 levels "有偶","未婚",..: 2 1 1 1 4 1 2 2 2 1 ...
## $ 年齡 : num 32.5 57.5 42.5 42.5 62.5 52.5 32.5 37.5 42.5 62.5 ...
## $ 年次 : num 76.5 51.5 66.5 66.5 46.5 56.5 76.5 71.5 66.5 46.5 ...
## $ 結過婚: num 0 1 1 1 1 1 0 0 0 1 ...2 Analysis
2.1 Logistic regression
m0 (simple model)
summary(m0<-glm(結過婚~年次,data = dtaM, family = binomial))##
## Call:
## glm(formula = 結過婚 ~ 年次, family = binomial, data = dtaM)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5499 0.2811 0.4498 0.7040 1.0540
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.72820 0.21685 35.64 <2e-16 ***
## 年次 -0.09713 0.00319 -30.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9913.1 on 9999 degrees of freedom
## Residual deviance: 8766.2 on 9998 degrees of freedom
## AIC: 8770.2
##
## Number of Fisher Scoring iterations: 5exp(cbind(coef=coef(m0), confint(m0)))## Waiting for profiling to be done...## coef 2.5 % 97.5 %
## (Intercept) 2271.5014815 1491.0940567 3489.2110673
## 年次 0.9074337 0.9017331 0.9130795m_full (full model)
summary (m_full <- glm (結過婚~年次*性別*教育, data = dtaM, family = binomial))##
## Call:
## glm(formula = 結過婚 ~ 年次 * 性別 * 教育, family = binomial,
## data = dtaM)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0831 0.2328 0.4607 0.6584 1.3921
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.242274 0.569092 9.212 < 2e-16 ***
## 年次 -0.066837 0.008136 -8.215 < 2e-16 ***
## 性別男 7.167431 1.059813 6.763 1.35e-11 ***
## 教育國中 -0.631569 0.993666 -0.636 0.525041
## 教育高中 0.522242 0.802506 0.651 0.515199
## 教育專科 0.574094 1.012367 0.567 0.570659
## 教育大學以上 0.993446 1.476495 0.673 0.501049
## 年次:性別男 -0.098016 0.014898 -6.579 4.73e-11 ***
## 年次:教育國中 0.015964 0.014465 1.104 0.269759
## 年次:教育高中 0.006055 0.011738 0.516 0.605950
## 年次:教育專科 0.007865 0.015370 0.512 0.608822
## 年次:教育大學以上 0.015564 0.025722 0.605 0.545132
## 性別男:教育國中 -2.238388 1.687470 -1.326 0.184682
## 性別男:教育高中 -5.033484 1.319284 -3.815 0.000136 ***
## 性別男:教育專科 -3.991911 1.508908 -2.646 0.008156 **
## 性別男:教育大學以上 -5.842571 1.995969 -2.927 0.003420 **
## 年次:性別男:教育國中 0.027590 0.024293 1.136 0.256089
## 年次:性別男:教育高中 0.056880 0.018937 3.004 0.002667 **
## 年次:性別男:教育專科 0.035498 0.022259 1.595 0.110763
## 年次:性別男:教育大學以上 0.044030 0.033610 1.310 0.190183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9913.1 on 9999 degrees of freedom
## Residual deviance: 8489.7 on 9980 degrees of freedom
## AIC: 8529.7
##
## Number of Fisher Scoring iterations: 6Drop 3way interaction
drop1(m_full, test ='Chisq')## Single term deletions
##
## Model:
## 結過婚 ~ 年次 * 性別 * 教育
## Df Deviance AIC LRT Pr(>Chi)
## <none> 8489.7 8529.7
## 年次:性別:教育 4 8499.1 8531.1 9.4573 0.05063 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m_drop1 <- update (m_full, ~ . -年次:性別:教育))##
## Call:
## glm(formula = 結過婚 ~ 年次 + 性別 + 教育 + 年次:性別 + 年次:教育 +
## 性別:教育, family = binomial, data = dtaM)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8957 0.2352 0.4540 0.6443 1.4225
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.024854 0.511003 11.790 < 2e-16 ***
## 年次 -0.078070 0.007260 -10.753 < 2e-16 ***
## 性別男 4.750472 0.511680 9.284 < 2e-16 ***
## 教育國中 -1.237958 0.799252 -1.549 0.12141
## 教育高中 -1.019981 0.621028 -1.642 0.10051
## 教育專科 -0.259934 0.716219 -0.363 0.71666
## 教育大學以上 -0.145389 0.955842 -0.152 0.87910
## 年次:性別男 -0.063821 0.007125 -8.958 < 2e-16 ***
## 年次:教育國中 0.024595 0.011556 2.128 0.03331 *
## 年次:教育高中 0.028792 0.009027 3.190 0.00143 **
## 年次:教育專科 0.019908 0.010677 1.865 0.06223 .
## 年次:教育大學以上 0.033298 0.016165 2.060 0.03941 *
## 性別男:教育國中 -0.270610 0.172757 -1.566 0.11725
## 性別男:教育高中 -1.104423 0.139475 -7.918 2.41e-15 ***
## 性別男:教育專科 -1.491088 0.178288 -8.363 < 2e-16 ***
## 性別男:教育大學以上 -2.882467 0.314888 -9.154 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9913.1 on 9999 degrees of freedom
## Residual deviance: 8499.1 on 9984 degrees of freedom
## AIC: 8531.1
##
## Number of Fisher Scoring iterations: 62.2 Visdualization01
library(ggplot2)
ggplot(data=dtaM, aes(x=年次, y=結過婚, shape=性別, color=性別)) + stat_summary(fun.data = mean_se, size = 1) + stat_summary (aes(y = fitted(m_drop1)), fun.y = mean, geom = 'line', size = 1) + facet_grid(.~教育)+ scale_x_continuous(breaks = seq(30,70, by = 5)) + labs(x='年次', y= '結過婚百分比') + theme (legend.position = c(.9,.1))## Warning: `fun.y` is deprecated. Use `fun` instead.
2.3 Dispersion parameter
library(dispmod)## Warning: package 'dispmod' was built under R version 4.0.3#將資料整理成 dispmod 要的形式
dta_f <- data.frame(ftable(dtaM, row.vars = c(1, 2, 5), col.vars = '結過婚'))
dta_f1 <- subset(dta_f , 結過婚 == 1)
dta_f1$Tot <- dta_f[1:70, 5] + dta_f[71:140, 5]
dta_f1 <- dta_f1[,-4]
#將年次設定為連續變項
dta_f1$年次 <- as.numeric(as.vector(dta_f1$年次))
#看看資料
#程式報表11.9
head(dta_f1)## 教育 性別 年次 Freq Tot
## 71 國小以下 女 46.5 47 53
## 72 國中 女 46.5 32 37
## 73 高中 女 46.5 117 132
## 74 專科 女 46.5 78 89
## 75 大學以上 女 46.5 302 308
## 76 國小以下 男 46.5 76 77#以整理過資料重新分析，確認與前面相同
summary(m_last <- glm(cbind(Freq, Tot-Freq) ~ 年次 * 性別 * 教育 - 年次:性別:教育,
data = dta_f1, family = binomial))##
## Call:
## glm(formula = cbind(Freq, Tot - Freq) ~ 年次 * 性別 * 教育 -
## 年次:性別:教育, family = binomial, data = dta_f1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7255 -0.8572 0.0203 0.5868 3.1982
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.024854 0.511003 11.790 < 2e-16 ***
## 年次 -0.078070 0.007260 -10.753 < 2e-16 ***
## 性別男 4.750472 0.511680 9.284 < 2e-16 ***
## 教育國中 -1.237958 0.799252 -1.549 0.12141
## 教育高中 -1.019981 0.621028 -1.642 0.10051
## 教育專科 -0.259934 0.716218 -0.363 0.71666
## 教育大學以上 -0.145389 0.955842 -0.152 0.87910
## 年次:性別男 -0.063821 0.007125 -8.958 < 2e-16 ***
## 年次:教育國中 0.024595 0.011556 2.128 0.03331 *
## 年次:教育高中 0.028792 0.009027 3.190 0.00143 **
## 年次:教育專科 0.019908 0.010677 1.865 0.06223 .
## 年次:教育大學以上 0.033298 0.016165 2.060 0.03941 *
## 性別男:教育國中 -0.270610 0.172757 -1.566 0.11725
## 性別男:教育高中 -1.104423 0.139475 -7.918 2.41e-15 ***
## 性別男:教育專科 -1.491088 0.178288 -8.363 < 2e-16 ***
## 性別男:教育大學以上 -2.882467 0.314888 -9.154 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1530.13 on 69 degrees of freedom
## Residual deviance: 116.19 on 54 degrees of freedom
## AIC: 446.31
##
## Number of Fisher Scoring iterations: 4#加入分散參數
#程式報表11.10, 11.11
summary(m_lastd <- glm.binomial.disp(m_last))##
## Binomial overdispersed logit model fitting...
## Iter. 1 phi: 0.009144343
## Iter. 2 phi: 0.0068854
## Iter. 3 phi: 0.007239847
## Iter. 4 phi: 0.007179674
## Iter. 5 phi: 0.007189756
## Iter. 6 phi: 0.007188063
## Iter. 7 phi: 0.007188347
## Iter. 8 phi: 0.0071883
## Iter. 9 phi: 0.007188308
## Iter. 10 phi: 0.007188306
## Iter. 11 phi: 0.007188306
## Converged after 11 iterations.
## Estimated dispersion parameter: 0.007188306
##
## Call:
## glm(formula = cbind(Freq, Tot - Freq) ~ 年次 * 性別 * 教育 -
## 年次:性別:教育, family = binomial, data = dta_f1, weights = disp.weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.34206 -0.58236 0.01665 0.49197 1.91197
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.58425 0.70791 7.888 3.06e-15 ***
## 年次 -0.07189 0.01043 -6.895 5.38e-12 ***
## 性別男 4.87049 0.71603 6.802 1.03e-11 ***
## 教育國中 -0.84564 1.03835 -0.814 0.4154
## 教育高中 -0.88963 0.90997 -0.978 0.3283
## 教育專科 0.05569 0.99733 0.056 0.9555
## 教育大學以上 0.03362 1.26618 0.027 0.9788
## 年次:性別男 -0.06503 0.01015 -6.407 1.49e-10 ***
## 年次:教育國中 0.01878 0.01533 1.225 0.2204
## 年次:教育高中 0.02688 0.01355 1.984 0.0473 *
## 年次:教育專科 0.01487 0.01495 0.995 0.3200
## 年次:教育大學以上 0.03072 0.02041 1.505 0.1322
## 性別男:教育國中 -0.28931 0.24870 -1.163 0.2447
## 性別男:教育高中 -1.13265 0.22766 -4.975 6.52e-07 ***
## 性別男:教育專科 -1.48291 0.25459 -5.825 5.72e-09 ***
## 性別男:教育大學以上 -2.87409 0.41432 -6.937 4.01e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 633.454 on 69 degrees of freedom
## Residual deviance: 52.675 on 54 degrees of freedom
## AIC: 236.78
##
## Number of Fisher Scoring iterations: 4##
## Call:
## glm(formula = cbind(Freq, Tot - Freq) ~ 年次 * 性別 * 教育 -
## 年次:性別:教育, family = binomial, data = dta_f1, weights = disp.weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.34206 -0.58236 0.01665 0.49197 1.91197
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.58425 0.70791 7.888 3.06e-15 ***
## 年次 -0.07189 0.01043 -6.895 5.38e-12 ***
## 性別男 4.87049 0.71603 6.802 1.03e-11 ***
## 教育國中 -0.84564 1.03835 -0.814 0.4154
## 教育高中 -0.88963 0.90997 -0.978 0.3283
## 教育專科 0.05569 0.99733 0.056 0.9555
## 教育大學以上 0.03362 1.26618 0.027 0.9788
## 年次:性別男 -0.06503 0.01015 -6.407 1.49e-10 ***
## 年次:教育國中 0.01878 0.01533 1.225 0.2204
## 年次:教育高中 0.02688 0.01355 1.984 0.0473 *
## 年次:教育專科 0.01487 0.01495 0.995 0.3200
## 年次:教育大學以上 0.03072 0.02041 1.505 0.1322
## 性別男:教育國中 -0.28931 0.24870 -1.163 0.2447
## 性別男:教育高中 -1.13265 0.22766 -4.975 6.52e-07 ***
## 性別男:教育專科 -1.48291 0.25459 -5.825 5.72e-09 ***
## 性別男:教育大學以上 -2.87409 0.41432 -6.937 4.01e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 633.454 on 69 degrees of freedom
## Residual deviance: 52.675 on 54 degrees of freedom
## AIC: 236.78
##
## Number of Fisher Scoring iterations: 4#看看分散參數估計值
m_lastd$dispersion## [1] 0.007188306#看看模型是否改善
anova(m_last, m_lastd)## Analysis of Deviance Table
##
## Model 1: cbind(Freq, Tot - Freq) ~ 年次 * 性別 * 教育 - 年次:性別:教育
## Model 2: cbind(Freq, Tot - Freq) ~ 年次 * 性別 * 教育 - 年次:性別:教育
## Resid. Df Resid. Dev Df Deviance
## 1 54 116.189
## 2 54 52.675 0 63.5152.4 Visdualization02
#準備再一次模型診斷，收集資料，畫圖
#圖11.3
dta_mlastd <- data.frame(dta_f1, phat = fitted(m_lastd))
ggplot(data = dta_mlastd, aes(x = 年次, y = Freq/Tot, color = 性別, shape = 性別)) +
geom_point(size = 3) +
stat_smooth(aes(x = 年次, y = phat), method = 'loess', se = F, size = 1) +
facet_grid(. ~ 教育) +
scale_x_continuous(breaks = seq(30, 70, by = 5)) +
labs(x= '年次', y = '結過婚百分比') +
theme(legend.position = 'NONE' )## `geom_smooth()` using formula 'y ~ x'
#準備繪圖。為了繪製方便，模型中去掉截距
summary(m_lastd0 <- update(m_lastd, . ~ . - 1))##
## Call:
## glm(formula = cbind(Freq, Tot - Freq) ~ 年次 + 性別 + 教育 +
## 年次:性別 + 年次:教育 + 性別:教育 - 1, family = binomial,
## data = dta_f1, weights = disp.weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.34206 -0.58236 0.01665 0.49197 1.91197
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## 年次 -0.07189 0.01043 -6.895 5.38e-12 ***
## 性別女 5.58425 0.70791 7.888 3.06e-15 ***
## 性別男 10.45474 0.83601 12.506 < 2e-16 ***
## 教育國中 -0.84564 1.03835 -0.814 0.4154
## 教育高中 -0.88963 0.90997 -0.978 0.3283
## 教育專科 0.05569 0.99733 0.056 0.9555
## 教育大學以上 0.03362 1.26618 0.027 0.9788
## 年次:性別男 -0.06503 0.01015 -6.407 1.49e-10 ***
## 年次:教育國中 0.01878 0.01533 1.225 0.2204
## 年次:教育高中 0.02688 0.01355 1.984 0.0473 *
## 年次:教育專科 0.01487 0.01495 0.995 0.3200
## 年次:教育大學以上 0.03072 0.02041 1.505 0.1322
## 性別男:教育國中 -0.28931 0.24870 -1.163 0.2447
## 性別男:教育高中 -1.13265 0.22766 -4.975 6.52e-07 ***
## 性別男:教育專科 -1.48291 0.25459 -5.825 5.72e-09 ***
## 性別男:教育大學以上 -2.87409 0.41432 -6.937 4.01e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2504.909 on 70 degrees of freedom
## Residual deviance: 52.675 on 54 degrees of freedom
## AIC: 236.78
##
## Number of Fisher Scoring iterations: 4require(coefplot)## Loading required package: coefplot## Warning: package 'coefplot' was built under R version 4.0.3#繪製估計值
#圖11.4
coefplot(m_lastd0) + labs(x = '估計值', y = '變項', title = '')
2.5 多分類邏輯迴歸
##多分類邏輯迴歸
#先看看年代與婚姻列聯表
#程式報表11.12
with(dtaM, prop.table(table(年次, 婚姻), 1))## 婚姻
## 年次 有偶 未婚 離婚 喪偶
## 46.5 0.719964664 0.056537102 0.113957597 0.109540636
## 51.5 0.728177642 0.080398162 0.129402757 0.062021440
## 56.5 0.713735558 0.087291399 0.155327343 0.043645700
## 61.5 0.716875872 0.140167364 0.120641562 0.022315202
## 66.5 0.671080139 0.183275261 0.135888502 0.009756098
## 71.5 0.622031662 0.278364116 0.093007916 0.006596306
## 76.5 0.455836936 0.478690550 0.058060531 0.007411983#載入 mlogit 套件，準備執行多分類邏輯迴歸
require(mlogit)## Loading required package: mlogit## Warning: package 'mlogit' was built under R version 4.0.3## Loading required package: dfidx## Warning: package 'dfidx' was built under R version 4.0.3##
## Attaching package: 'dfidx'## The following object is masked from 'package:stats':
##
## filter#將資料轉成 mlogit 需要格式
dta_w <- mlogit.data(dtaM, shape = 'wide', choice = '婚姻')
#看看資料
#程式報表11.13
head(dta_w)## ~~~~~~~
## first 10 observations out of 40000
## ~~~~~~~
## 教育 性別 婚姻 年齡 年次 結過婚 chid alt idx
## 1 國小以下 女 TRUE 32.5 76.5 0 1 未婚 1:未婚
## 2 國小以下 女 FALSE 32.5 76.5 0 1 有偶 1:有偶
## 3 國小以下 女 FALSE 32.5 76.5 0 1 喪偶 1:喪偶
## 4 國小以下 女 FALSE 32.5 76.5 0 1 離婚 1:離婚
## 5 專科 男 FALSE 57.5 51.5 1 2 未婚 2:未婚
## 6 專科 男 TRUE 57.5 51.5 1 2 有偶 2:有偶
## 7 專科 男 FALSE 57.5 51.5 1 2 喪偶 2:喪偶
## 8 專科 男 FALSE 57.5 51.5 1 2 離婚 2:離婚
## 9 國中 男 FALSE 42.5 66.5 1 3 未婚 3:未婚
## 10 國中 男 TRUE 42.5 66.5 1 3 有偶 3:有偶
##
## ~~~ indexes ~~~~
## chid alt
## 1 1 未婚
## 2 1 有偶
## 3 1 喪偶
## 4 1 離婚
## 5 2 未婚
## 6 2 有偶
## 7 2 喪偶
## 8 2 離婚
## 9 3 未婚
## 10 3 有偶
## indexes: 1, 2#只展示年代效果
#程式報表11.14
summary(ml0 <- mlogit(婚姻 ~ 0 | 年次, data = dta_w ) )##
## Call:
## mlogit(formula = 婚姻 ~ 0 | 年次, data = dta_w, method = "nr")
##
## Frequencies of alternatives:choice
## 未婚 有偶 喪偶 離婚
## 0.1966 0.6550 0.0341 0.1143
##
## nr method
## 7 iterations, 0h:0m:2s
## g'(-H)^-1g = 0.000103
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):有偶 7.3104740 0.2200718 33.219 < 2.2e-16 ***
## (Intercept):喪偶 9.8996465 0.4503860 21.980 < 2.2e-16 ***
## (Intercept):離婚 5.9703619 0.2805370 21.282 < 2.2e-16 ***
## 年次:有偶 -0.0935790 0.0032415 -28.869 < 2.2e-16 ***
## 年次:喪偶 -0.1905932 0.0079392 -24.007 < 2.2e-16 ***
## 年次:離婚 -0.1002414 0.0043105 -23.255 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -8918.3
## McFadden R^2: 0.07105
## Likelihood ratio test : chisq = 1364.2 (p.value = < 2.22e-16)#勝算比信賴區間
cbind(b_hat = exp(coef(ml0)), b_ci = exp(confint(ml0)))## b_hat 2.5 % 97.5 %
## (Intercept):有偶 1.495886e+03 971.7913780 2.302629e+03
## (Intercept):喪偶 1.992333e+04 8241.2398186 4.816495e+04
## (Intercept):離婚 3.916474e+02 225.9964875 6.787171e+02
## 年次:有偶 9.106660e-01 0.9048986 9.164702e-01
## 年次:喪偶 8.264687e-01 0.8137080 8.394295e-01
## 年次:離婚 9.046190e-01 0.8970085 9.122941e-01#準備畫圖，先得到預測值
newdta <- with(dtaM, data.frame(
年次 = rep(seq(from = 39.5, to = 69.5, length = 7), 4),
婚姻 = factor(rep(c('有偶', '未婚', '離婚', '喪偶'), each = 7) ) ) )
newdta_w <- mlogit.data(newdta, shape = 'wide', choice = '婚姻')
ml0_phat <- as.data.frame(predict(ml0, newdta_w)[1:7, ])
#需轉換資料形式，載入 reshape 套件
require(reshape2)## Loading required package: reshape2#製造年次婚姻列聯表、將預測值轉成需要格式、彙整在一起
obs_p <- data.frame(with(dtaM, prop.table(table(年次, 婚姻), 1)))
dta_op <- data.frame(obs_p, melt(ml0_phat))## No id variables; using all as measure variables2.6 Visdualization03
#繪圖
#圖11.5
ggplot(data = dta_op, aes(x = 年次, y = Freq, shape = 婚姻)) +
geom_point(size = 3) +
stat_smooth(aes(x = 年次, y = value, group = 婚姻), method = 'loess', se = F, size = 1) +
labs(x = '年次', y = '分類百分比') +
theme(legend.position = c(.92, .88))## `geom_smooth()` using formula 'y ~ x'
The End