广义线性模型
1 广义线性模型
1.1 Logisitic回归：妊娠糖尿病分析
因变量：糖尿病(Diabetes){阳性：pos，阴性：neg}，设阳性为1、阴性为0
自变量：年龄(Age)、体重指数(BMI，kg/m2 )、血糖浓度(Glucose )、舒张压(Diastolic blood pressure，(mm)Hg )、怀孕次数(Number of times pregnant )
数据文件：diabetes.csv，共 724个观察值
1.1.1 划分训练集和测试集
- 前450条个案为训练集，用于估计Logist模型
- 后274条个案为测试集，用于评价模型的估计效果
- 训练集糖尿病率36.44%，测试集糖尿病率为31.02%，两者大致相等。
1.1.2 训练集估计回归方程
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -7.950 | 0.97 | -8.21 | 0.00 |
| Age | 0.012 | 0.01 | 1.00 | 0.32 |
| BMI | 0.089 | 0.02 | 4.81 | 0.00 |
| Glucose | 0.032 | 0.00 | 7.39 | 0.00 |
| Pressure | -0.005 | 0.01 | -0.50 | 0.61 |
| Pregnant | 0.098 | 0.04 | 2.45 | 0.01 |
\[ log(\frac{p}{1-p}) =-7.95+0.012\times Age+0.089\times BMI+0.032\times Glucose-0.005\times Pressure+0.098\times Pregnant\]
1.1.3 测试集预测效果评价
| pos_pred | neg_pred | |
|---|---|---|
| pos | 53 | 32 |
| neg | 21 | 168 |
由训练集预测混淆矩阵可知
- 准确率(accuracy)：80.66%
- 精确率(precision)：62.35%
- 召回率(recall)：71.62%
- \(F_1\)得分(\(F_1\) score)：66.67%
1.1.4 回归模型边际效应
Age BMI Glucose Pressure Pregnant
0.26 1.94 0.70 -0.11 2.14 年龄每增加一岁患病风险提高0.26%；体重指数每增加1患病风险提高1.94%；血糖浓度每增加1患病风险提高0.7%；舒张压每增加1患病风险降低0.11%；怀孕次数每增加一次患病风险提高2.14%。
体重指数、血糖浓度、怀孕次数对患病风险呈正向影响，符合预期。
年龄和舒张压对患病风险影响小且不显著，考虑逐步回归选择更合适的模型
1.2 Logisitic回归逐步回归
1.2.1 逐步回归的回归方程
1.2.2 逐步回归的预测效果
1.2.3 逐步回归的边际效应
2 判别分析
2.1 线性贝叶斯判别
2.1.1 判别函数
由于只有目标变量只有两类，线性贝叶斯判别等价于Fisher判别。以下为Fisher判别的判别函数：
\[ w =5.981+0.008\times Age+0.068\times BMI+0.028\times Glucose-0.003\times Pressure+0.082\times Pregnant\]
该判别函数和Logsitic回归的方程近似等价，各系数存在近似的倍数关系。注意，这里\(w\)大于0判为neg，小于0判为pos。
2.1.2 测试集合的预测结果
测试集合预测的后验概率:
两类后验概率差异越大代表判别越有把握，错判的概率越小。
测试集合预测的混淆矩阵
| pos_pred | neg_pred | |
|---|---|---|
| pos | 53 | 32 |
| neg | 22 | 167 |
预测效果与Logistic模型基本一致。
2.2 二次贝叶斯判别
2.2.1 测试集合的预测
2.2.2 测试集合的混淆矩阵
3 聚类分析
利用例子7.2中2007年城镇居民消费数据作聚类分析，并比较不同聚类的效果
3.1 系统聚类
3.1.1 类平均法
从树状图看两类聚合为一类时聚类距离明显突变，分为两类比较合理。其中第一类包括：北京、上海、浙江、广东，其余为一类。
从各类的类中心看，第一类为沿海发达地区，各类消费水平明显高于第二类。
| 类别 | 食品 | 衣着 | 设备 | 医疗 | 交通 | 教育 | 居住 | 杂项 |
|---|---|---|---|---|---|---|---|---|
| 1 | 5252.19 | 1265.92 | 864.95 | 940.69 | 2730.43 | 2297.70 | 1317.82 | 583.77 |
| 2 | 3252.73 | 983.08 | 521.21 | 635.12 | 1012.19 | 1072.89 | 880.07 | 316.79 |
3.1.2 离差平方和法
plot(hclust(D,'ward.D'),hang=-1) #ward.D法
plot(hclust(D,'ward.D2'),hang=-1) #ward.D2法
3.2 K-means聚类
km=kmeans(X,2) #kmeans聚类
km$cluster #分类结果 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1 2 1 2 1 2 2 1 2 1 1 1 1 2 2 2 2 2 2 1
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
2 1 2 1 2 1 2 2 1 2 2 2 2 1 2 2 2 1 1 2
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
1 1 1 1 2 2 2 2 2 1 1 2 1 2 1 2 1 2 1 2
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
2 2 2 2 1 2 1 1 2 2 2 2 2 2 2 2 1 2 2 2
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
1 2 2 2 2 1 1 1 2 2 2 2 1 1 2 2 2 2 2 1
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
2 2 1 1 2 2 1 1 2 2 2 2 1 2 2 1 2 2 2 2
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
2 2 1 2 1 2 2 2 2 2 2 2 1 2 2 2 1 2 2 1
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
2 1 2 1 1 1 1 2 2 2 1 1 2 2 2 1 2 1 2 2
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
2 2 1 2 2 1 2 1 1 1 2 2 2 1 1 1 1 2 1 2
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
2 1 2 1 2 2 2 1 2 1 2 2 2 2 1 1 2 1 2 1
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
1 1 2 1 2 2 2 2 1 1 1 2 2 2 1 1 2 1 1 2
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
2 2 1 1 1 1 2 2 2 1 2 1 1 2 1 2 2 2 2 2
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
2 2 2 2 2 1 1 1 2 1 2 2 2 2 2 2 2 2 2 2
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
2 2 2 2 1 1 1 1 2 1 1 2 2 2 2 2 2 2 1 1
281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
1 2 2 2 1 2 2 1 2 1 1 2 2 2 2 1 2 2 2 2
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
1 2 1 2 2 2 1 2 1 2 1 2 2 2 2 2 2 1 2 1
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
1 2 2 2 2 1 1 2 2 1 2 2 2 1 2 2 1 1 1 2
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
1 1 2 2 2 2 1 1 2 2 2 1 2 2 1 2 2 2 2 2
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380
2 2 2 2 1 2 2 1 2 2 1 2 2 2 2 1 2 1 1 2
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400
1 2 2 2 1 1 2 2 1 1 1 1 2 1 2 2 2 2 2 2
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
1 1 1 1 2 2 2 1 2 1 1 2 2 1 2 2 2 2 1 2
421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
2 2 2 2 1 2 2 1 1 2 1 1 2 2 2 2 2 2 2 2
441 442 443 444 445 446 447 448 449 450
1 1 1 2 1 2 1 2 2 2 plot(X,pch=km$cluster)
set.seed(123)
x1=matrix(rnorm(10000,0,0.25),ncol=10)
x2=matrix(rnorm(10000,1,0.25),ncol=10)
X=rbind(x1,x2)
km=kmeans(X,2) #kmeans聚类
kc=km$cluster;kc [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1037] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1074] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1111] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1148] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1185] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1222] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1259] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1296] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1333] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1370] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1407] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1444] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1481] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1518] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1555] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1592] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1629] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1666] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1703] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1740] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1777] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1814] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1851] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1888] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1925] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1962] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[1999] 2 2plot(X,pch=kc,col=kc)
3.2.1 分为三类
第一类包括：河北、山西、内蒙古、吉林、黑龙江、江西、河南、贵州、陕西、甘肃、青海、宁夏、新疆 ，多为中部和沿海地区
第二类包括：天津、辽宁、江苏、安徽、福建、山东、湖北、湖南、广西、海南、重庆、四川、云南、西藏，多为东北和西部地区 ，多为东北和西部地区。
第三类包括：北京、上海、浙江、广东 ，为沿海发达地区。
从各类的类中心看，第三类地区的各项消费水平最高，尤其在交通、教育等项目；第一类地区的的各项消费水平均低于第三类但高于第二类地区。
| 类别 | 食品 | 衣着 | 设备 | 医疗 | 交通 | 教育 | 居住 | 杂项 |
|---|---|---|---|---|---|---|---|---|
| 1 | 2840.50 | 1030.10 | 492.95 | 634.22 | 887.03 | 1006.11 | 833.05 | 322.25 |
| 2 | 3635.52 | 939.41 | 547.45 | 635.95 | 1128.40 | 1134.90 | 923.74 | 311.73 |
| 3 | 5252.19 | 1265.92 | 864.95 | 940.69 | 2730.43 | 2297.70 | 1317.82 | 583.77 |