Generalized linear models
requires packages AER, robust, gcc
install.packages(c("AER", "robust", "qcc"))
微信截圖_20200313210929.png
Logistic Regression
# get summary statistics,示例數據見上
data(Affairs, package="AER")
summary(Affairs)
table(Affairs$affairs)
# create binary outcome variable甸陌，婚外情次數分為（有或無）二分類變量
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"))
table(Affairs$ynaffair)
# fit full model楼誓，納入所有的變量
fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children +
religiousness + education + occupation +rating,
data=Affairs,family=binomial())
summary(fit.full)
# fit reduced model，納入顯著的變量
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, data=Affairs, family=binomial())
summary(fit.reduced)
# compare models
anova(fit.reduced, fit.full, test="Chisq")
# interpret coefficients
coef(fit.reduced)
exp(coef(fit.reduced))
# calculate probability of extramariatal affair by marital ratings
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
# calculate probabilites of extramariatal affair by age
testdata <- data.frame(rating = mean(Affairs$rating),
age = seq(17, 57, 10),
yearsmarried = mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")
testdata
# evaluate overdispersion
fit <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = binomial(), data = Affairs)
fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = quasibinomial(), data = Affairs) # robust模型
pchisq(summary(fit.od)$dispersion * fit$df.residual,
fit$df.residual, lower = F) #P值不顯著红竭，就選第一個模型
# Logistic回歸案例2
Example11_4 <- read.table ("example11_4.csv", header=TRUE, sep=",")
attach(Example11_4)
fit1 <- glm(y~ x1 + x2, family= binomial(), data=Example11_4)
summary(fit1)
coefficients(fit1)
exp(coefficients(fit1))
exp (confint(fit1))
fit2 <- glm(y~ x1 + x2 + x1:x2 , family= binomial(), data=Example11_4)
summary(fit2)
coefficients(fit2)
exp(coefficients(fit2))
exp (confint(fit2))
Example11_4$x11 <- ifelse (x1==1 & x2==1, 1, 0)
Example11_4$x10 <- ifelse (x1==1 & x2==0, 1, 0)
Example11_4$x01 <- ifelse (x1==0 & x2==1, 1, 0)
Example11_4$x00 <- ifelse (x1==0 & x2==0, 1, 0)
fit3 <- glm(y~ x11 + x10 + x01, family= binomial(), data=Example11_4)
summary(fit3)
coefficients(fit3)
exp(coefficients(fit3))
exp(confint(fit3))
detach (Example11_4)
# Logistic回歸案例3
Example11_5 <- read.table ("example11_5.csv", header=TRUE, sep=",")
attach(Example11_5)
fullfit <- glm(y~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 , family= binomial(), data=Example11_5)
summary(fullfit)
nothing <- glm(y~ 1, family= binomial(), data=Example11_5)
summary(nothing)
bothways <- step(nothing, list(lower=formula(nothing), upper=formula(fullfit)), direction="both")
fit1 <- glm(y~ x6 + x5 + x8 + x1 + x2 , family= binomial(), data=Example11_5)
summary(fit1)
fit2 <- glm(y~ x6 + x5 + x8 + x1, family= binomial(), data=Example11_5)
summary(fit2)
coefficients(fit2)
exp(coefficients(fit2))
exp (confint(fit2))
detach (Example11_5)
#條件Logistic回歸案例 4
#install.packages("survival")
library(survival)
Example11_6 <- read.table ("example11_6.csv", header=TRUE, sep=",")
attach(Example11_6)
model <- clogit(outcome~ exposure+ strata(id))
summary(model)
detach(Example11_6)
#Logistic回歸案例 5 低出生體重兒列線圖
library(foreign)
library(rms)
mydata<-read.spss("lweight.sav")
mydata<-as.data.frame(mydata)
head(mydata)
mydata$low <- ifelse(mydata$low =="低出生體重",1,0)
mydata$race1 <- ifelse(mydata$race =="白種人",1,0)
mydata$race2 <- ifelse(mydata$race =="黑種人",1,0)
mydata$race3 <- ifelse(mydata$race =="其他種族",1,0)
attach(mydata)
dd<-datadist(mydata)
options(datadist='dd')
fit1<-lrm(low~age+ftv+ht+lwt+ptl+smoke+ui+race1+race2,data=mydata,x=T,y=T)
fit1
summary(fit1)
nom1 <- nomogram(fit1, fun=plogis,fun.at=c(.001, .01, .05, seq(.1,.9, by=.1), .95, .99, .999),lp=F, funlabel="Low weight rate")
plot(nom1)
cal1 <- calibrate(fit1, method='boot', B=1000)
plot(cal1,xlim=c(0,1.0),ylim=c(0,1.0))
mydata$race <- as.factor(ifelse(mydata$race=="白種人", "白種人","黑人及其他種族"))
dd<-datadist(mydata)
options(datadist='dd')
fit2<-lrm(low~age+ftv+ht+lwt+ptl+smoke+ui+race,data=mydata,x=T,y=T)
fit2
summary(fit2)
nom2 <- nomogram(fit2, fun=plogis,fun.at=c(.001, .01, .05, seq(.1,.9, by=.1), .95, .99, .999),lp=F, funlabel="Low weight rate")
plot(nom2)
cal2 <- calibrate(fit2, method='boot', B=1000)
plot(cal2,xlim=c(0,1.0),ylim=c(0,1.0))
fit3<-lrm(low~ht+lwt+ptl+smoke+race,data=mydata,x=T,y=T)
fit3
summary(fit3)
nom3 <- nomogram(fit3, fun=plogis,fun.at=c(.001, .01, .05, seq(.1,.9, by=.1), .95, .99, .999),lp=F, funlabel="Low weight rate")
plot(nom3)
cal3 <- calibrate(fit3, method='boot', B=1000)
plot(cal3,xlim=c(0,1.0),ylim=c(0,1.0))
#C-statistics計算
library(foreign)
library(rms)
mydata<-read.spss("lweight.sav")
mydata<-as.data.frame(mydata)
head(mydata)
mydata$low <- ifelse(mydata$low =="低出生體重",1,0)
mydata$race1 <- ifelse(mydata$race =="白種人",1,0)
mydata$race2 <- ifelse(mydata$race =="黑種人",1,0)
mydata$race3 <- ifelse(mydata$race =="其他種族",1,0)
attach(mydata)
dd<-datadist(mydata)
options(datadist='dd')
fit1<-lrm(low~age+ftv+ht+lwt+ptl+smoke+ui+race1+race2,data=mydata,x=T,y=T)
fit1 #直接讀取模型中Rank Discrim.參數 C
mydata$predvalue<-predict(fit1)
library(ROCR)
pred <- prediction(mydata$predvalue, mydata$low)
perf<- performance(pred,"tpr","fpr")
plot(perf)
abline(0,1)
auc <- performance(pred,"auc")
auc #auc即是C-statistics
#library(Hmisc)
somers2(mydata$predvalue, mydata$low) #somers2 {Hmisc}
#Logistic回歸案例 6 亞組分析森林圖
library(forestplot)
rs_forest <- read.csv('rs_forest.csv',header = FALSE)
# 讀入數據的時候大家一定要把header設置成FALSE，保證第一行不被當作列名稱卸留。
# tiff('Figure 1.tiff',height = 1600,width = 2400,res= 300)
forestplot(labeltext = as.matrix(rs_forest[,1:3]),
#設置用于文本展示的列蟋滴，此處我們用數據的前三列作為文本，在圖中展示
mean = rs_forest$V4, #設置均值
lower = rs_forest$V5, #設置均值的lowlimits限
upper = rs_forest$V6, #設置均值的uplimits限
is.summary = c(T,T,T,F,F,T,F,F,T,F,F),
#該參數接受一個邏輯向量咳焚，用于定義數據中的每一行是否是匯總值，若是庞溜，則在對應位置設置為TRUE革半，若否，則設置為FALSE；設置為TRUE的行則以粗體出現
zero = 1, #設置參照值又官，此處我們展示的是OR值延刘，故參照值是1，而不是0
boxsize = 0.4, #設置點估計的方形大小
lineheight = unit(10,'mm'),#設置圖形中的行距
colgap = unit(3,'mm'),#設置圖形中的列間距
lwd.zero = 2,#設置參考線的粗細
lwd.ci = 1.5,#設置區(qū)間估計線的粗細
col=fpColors(box='#458B00', summary= "#8B008B",lines = 'black',zero = '#7AC5CD'),
#使用fpColors()函數定義圖形元素的顏色六敬，從左至右分別對應點估計方形碘赖，匯總值，區(qū)間估計線觉阅，參考線
xlab="The estimates",#設置x軸標簽
graph.pos = 3)#設置森林圖的位置崖疤，此處設置為3秘车，則出現在第三列
Poisson Regression
# look at dataset典勇，計數資料，發(fā)生次數很小
data(breslow.dat, package="robust")
names(breslow.dat)
summary(breslow.dat[c(6, 7, 8, 10)])
# plot distribution of post-treatment seizure counts
opar <- par(no.readonly=TRUE)
par(mfrow=c(1, 2))
attach(breslow.dat)
hist(sumY, breaks=20, xlab="Seizure Count",
main="Distribution of Seizures")
boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons")
par(opar)
# fit regression
fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())
summary(fit)
# interpret model parameters
coef(fit)
exp(coef(fit))
# evaluate overdispersion
deviance(fit)/df.residual(fit)
library(qcc)
qcc.overdispersion.test(breslow.dat$sumY, type="poisson")
# fit model with quasipoisson
fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,
family=quasipoisson())
summary(fit.od)
