第2章
image.png
image.png
代码案例
数据形式
image.png
#第 2章 代码开始
> ## 利用 nnet 包中的函数multinom ，建立多元logistic回归模型：
> library(foreign)
> library(nnet)
> library(ggplot2)
> library(reshape2)
> ml <- read.dta("hsbdemo.dta")#读取.dta格式数据
> #View(ml)#查看数据的具体形式
> with(ml, table(ses, prog))#ses行和prog列的表格
prog
ses general academic vocation
low 16 19 12
middle 20 44 31
high 9 42 7
> #with在数据框上的操作
> #do.call(操作，数据)
> with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
M SD
general 51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754
> ml$prog2 <- relevel(ml$prog, ref = "academic")#以academic作为参考
> test <- multinom(prog2 ~ ses + write, data = ml)#回归
# weights: 15 (8 variable)
initial value 219.722458
iter 10 value 179.982880
final value 179.981726
converged
> summary(test)
Call:
multinom(formula = prog2 ~ ses + write, data = ml)
Coefficients:
(Intercept) sesmiddle seshigh write
general 2.852198 -0.5332810 -1.1628226 -0.0579287
vocation 5.218260 0.2913859 -0.9826649 -0.1136037
#行general 和vocation是指相对于academic的值，列sesmiddle seshigh是指相对于seslow的值
Std. Errors:
(Intercept) sesmiddle seshigh write
general 1.166441 0.4437323 0.5142196 0.02141097
vocation 1.163552 0.4763739 0.5955665 0.02221996
Residual Deviance: 359.9635
AIC: 375.9635
> # 2-tailed z test
> z <- summary(test)$coefficients/summary(test)$standard.errors #Z值的计算
> z
(Intercept) sesmiddle seshigh write
general 2.445214 -1.2018081 -2.261334 -2.705562
vocation 4.484769 0.6116747 -1.649967 -5.112689
> p <- (1 - pnorm(abs(z), 0, 1)) * 2##P值的计算
> p
(Intercept) sesmiddle seshigh write
general 0.0144766100 0.2294379 0.02373856 6.818902e-03
vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07
#general 相对于academic，因素seshigh和write是有意义的
#vocation 相对于academic，因素write是有意义的
> # extract the coefficients from the model and exponentiate
> exp(coef(test)) #OR值
(Intercept) sesmiddle seshigh write
general 17.32582 0.5866769 0.3126026 0.9437172
vocation 184.61262 1.3382809 0.3743123 0.8926116
#解释OR值，write每提高一个单位，产生general是产生academic的0.9437172倍
#write每提高一个单位，产生vocation是产生academic的0.8926116倍
#
> head(pp <- fitted(test))
academic general vocation
1 0.1482764 0.3382454 0.5134781
2 0.1202017 0.1806283 0.6991700
3 0.4186747 0.2368082 0.3445171
4 0.1726885 0.3508384 0.4764731
5 0.1001231 0.1689374 0.7309395
6 0.3533566 0.2377976 0.4088458
> dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
#这里把写作write取均值，即，每个人的写作贡献固定了，只有社会经济地位的变化了
> predict(test, newdata = dses, "probs")#预测
academic general vocation
1 0.4396845 0.3581917 0.2021238
2 0.4777488 0.2283353 0.2939159
3 0.7009007 0.1784939 0.1206054
> dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),3))
> # store the predicted probabilities for each value of ses and write
> pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))
> # calculate the mean probabilities within each level of ses
> by(pp.write[, 3:5], pp.write$ses, colMeans)
pp.write$ses: high
academic general vocation
0.6164315 0.1808037 0.2027648
--------------------------------------------------------
pp.write$ses: low
academic general vocation
0.3972977 0.3278174 0.2748849
--------------------------------------------------------
pp.write$ses: middle
academic general vocation
0.4256198 0.2010864 0.3732938
> #melt data set to long for ggplot2
> lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
> head(lpp) # view first few rows
ses write variable probability
1 low 30 academic 0.09843588
2 low 31 academic 0.10716868
3 low 32 academic 0.11650390
4 low 33 academic 0.12645834
5 low 34 academic 0.13704576
6 low 35 academic 0.14827643
> # plot predicted probabilities across write values for each level of ses
> # facetted by program type
> ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~., scales = "free")
> install.packages("mlogit")
> library(Formula)
> library(maxLik)
> library(miscTools)
> library(mlogit)
> data("Fishing", package = "mlogit")
> Fish <- mlogit.data(Fishing,shape = "wide",choice = "mode")
> summary(mlogit(mode ~ 0|income, data = Fish))
Call:
mlogit(formula = mode ~ 0 | income, data = Fish, method = "nr")
Frequencies of alternatives:choice
beach boat charter pier
0.11337 0.35364 0.38240 0.15059
nr method
4 iterations, 0h:0m:0s
g'(-H)^-1g = 8.32E-07
gradient close to zero
Coefficients :
Estimate Std. Error z-value Pr(>|z|)
(Intercept):boat 7.3892e-01 1.9673e-01 3.7560 0.0001727 ***
(Intercept):charter 1.3413e+00 1.9452e-01 6.8955 5.367e-12 ***
(Intercept):pier 8.1415e-01 2.2863e-01 3.5610 0.0003695 ***
income:boat 9.1906e-05 4.0664e-05 2.2602 0.0238116 *
income:charter -3.1640e-05 4.1846e-05 -0.7561 0.4495908
income:pier -1.4340e-04 5.3288e-05 -2.6911 0.0071223 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Log-Likelihood: -1477.2
McFadden R^2: 0.013736
Likelihood ratio test : chisq = 41.145 (p.value = 6.0931e-09)
2.2章 有序分类logistic回归 即 等级资料
有序分类logistic回归
案例
有可能把等级资料转换成二分类资料就转换，否则只能当作等级资料或多分类资料；
image.png
是否从学校毕业：不可能，有可能，极有可能；父母是否有学位：1表示有，0表示无；公立还是私立学校：1公里，0私立
## 程序包MASS提供polr()函数可以进行ordered logit或probit回归
require(foreign)
require(ggplot2)
require(MASS)
require(Hmisc)
require(reshape2)
dat <- read.dta("ologit.dta")
#View(dat)
head(dat)
## one at a time, table apply, pared, and public
lapply(dat[, c("apply", "pared", "public")], table)
## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ public + apply + pared, data = dat))
summary(dat$gpa)
sd(dat$gpa)
ggplot(dat, aes(x = apply, y = gpa)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
facet_grid(pared ~ public, margins = TRUE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
## fit ordered logit model and store results 'm'
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)
## view a summary of the model
summary(m)
## store table
(ctable <- coef(summary(m)))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable <- cbind(ctable, "p value" = p))
(ci <- confint(m)) # default method gives profiled CIs
confint.default(m) # CIs assuming normality
## odds ratios
exp(coef(m))
## OR and CI
exp(cbind(OR = coef(m), ci))
sf <- function(y) {
c('Y>=1' = qlogis(mean(y >= 1)),
'Y>=2' = qlogis(mean(y >= 2)),
'Y>=3' = qlogis(mean(y >= 3)))
}
(s <- with(dat, summary(as.numeric(apply) ~ pared + public + gpa, fun=sf)))
glm(I(as.numeric(apply) >= 2) ~ pared, family="binomial", data = dat)
glm(I(as.numeric(apply) >= 3) ~ pared, family="binomial", data = dat)
s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]
s # print
plot(s, which=1:3, pch=1:3, xlab='logit', main=' ', xlim=range(s[,3:4]))
newdat <- data.frame(
pared = rep(0:1, 200),
public = rep(0:1, each = 200),
gpa = rep(seq(from = 1.9, to = 4, length.out = 100), 4))
newdat <- cbind(newdat, predict(m, newdat, type = "probs"))
##show first few rows
head(newdat)
lnewdat <- melt(newdat, id.vars = c("pared", "public", "gpa"),
variable.name = "Level", value.name="Probability")
## view first few rows
head(lnewdat)
ggplot(lnewdat, aes(x = gpa, y = Probability, colour = Level)) +
geom_line() + facet_grid(pared ~ public, labeller="label_both")
