婚外情数据 R 二元logistic

婚外情数据 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$age),yearsmarried=mean(Affairs$yearsmarried),
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$rating),age=seq(17,57,10),yearsmarried=mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness))testdata1$prob <- predict(fit.reduced, newdata=testdata1, type=“response”)
testdata1

age = seq(17, 57, 10) 10年为一个等级

年龄Min. :17.50。 Max. :57.00