婚外情数据 R 二元logistic
install.packages(c("AER", "robust", "qcc")) #下载安装包
data(Affairs, package="AER") #打开 AER包中的,数据集Affairs
table(Affairs$affairs) #Affairs$affairs,
summary(Affairs)
str(Affairs)
affairs变量是指 婚外情的次数,是结果变量
int 整数
将Affairs$ affairs变量 变成二分类变量 Affairs$ ynaffair
Affairs$ynaffair[Affairs$affairs > 0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair,
levels=c(0,1),
labels=c("No","Yes"))
factor函数
Affairs $ ynaffair <- factor(Affairs$ ynaffair,
levels=c(0,1),
labels=c(“No”,“Yes”))
labels=c(“No”,“Yes”)
设置标签,1=有婚外情,0=无婚外情
table(Affairs$ynaffair)
fit full model,该模型纳入了所有的自变量
#glm 广义线性模型
fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +
religiousness + education + occupation +rating,
data=Affairs,family=binomial())
summary(fit.full)
coef(summary(fit.full))
fit reduced model
fit.reduced <- glm(ynaffair ~ age + yearsmarried +
religiousness +
rating,
data=Affairs, family=binomial())
summary(fit.reduced)
Fisher计分迭代次数:4
#Number of Fisher Scoring iterations: 4
compare models 使用方差 比较两个模型####
anova(fit.reduced, fit.full, test=“Chisq”)
两个模型没有差别,p=0.2108
summary(fit.reduced)
求每个自变量的系数####
coef(summary(fit.reduced))
coef(fit.reduced) #求每个自变量的系数
confint(fit.reduced) #得到系数的95%区间
计算OR值#####
exp(coef(fit.reduced)) #计算OR值
exp(confint(fit.reduced)) #计算OR值的置信区间
age OR值的95%区间
age 的系数 -0.07006400 -0.001854759
age OR 0.9323341 0.998147
一 calculate probability of extramariatal affair by marital ratings#####
#构建一个新的数据集testdata
#查看rating 对婚姻的自我评分这个变量 ,对结局婚外情次数的影响#####
如何操作呢
#1 把其他的变量取均数,得到一个新的数据框是testdata
#2 用predict函数,模型是fit.reduced,
新的数据框是testdata,
#类型为response,数据框是testdata中 产生一个新的变量 testdata$
prob
testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
age = mean(Affairs
a
g
e
)
,
y
e
a
r
s
m
a
r
r
i
e
d
=
m
e
a
n
(
A
f
f
a
i
r
s
age), yearsmarried = mean(Affairs
age),yearsmarried=mean(Affairsyearsmarried),
religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata=testdata,
type=“response”)
testdata
输出结果可以看到,rating打分越高,prob结局婚外情 几率会减少。
#二 calculate probabilites of extramariatal affair by age#####
#查看age年龄 这个变量 ,对结局婚外情次数的影响#####
testdata1 <- data.frame(rating = mean(Affairs
r
a
t
i
n
g
)
,
a
g
e
=
s
e
q
(
17
,
57
,
10
)
,
y
e
a
r
s
m
a
r
r
i
e
d
=
m
e
a
n
(
A
f
f
a
i
r
s
rating), age = seq(17, 57, 10), yearsmarried = mean(Affairs
rating),age=seq(17,57,10),yearsmarried=mean(Affairsyearsmarried),
religiousness = mean(Affairs
r
e
l
i
g
i
o
u
s
n
e
s
s
)
)
t
e
s
t
d
a
t
a
1
religiousness)) testdata1
religiousness))testdata1prob <- predict(fit.reduced, newdata=testdata1, type=“response”)
testdata1
age = seq(17, 57, 10) 10年为一个等级
年龄Min. :17.50。 Max. :57.00