R语言logistic回归分析

##logistic回归分析########################################################

##二分类logistic回归################################################

##logistic回归模型
p <- seq(from=0,to=1,by=0.01)

odds <- p/(1-p)
plot(log(odds),p,type ='l',col='blue',ylab = 'probability',las=1)
abline(h=0.5,lty='dashed')
abline(v=0,lty='dashed')

abline(v=-2,lty='dashed')
abline(v=2,lty='dashed')

##logistic回归实例
library(MASS)
data("birthwt")
str(birthwt)  

library(epiDisplay)
tab1(birthwt$ptl)
tab1(birthwt$ftv)

library(dplyr)
birthwtweight <- birthwt%>%
  mutate(race=factor(race,labels = c('white','black','other')),
         smoke=factor(smoke,labels = c('no','yes')),
         ptl=ifelse(ptl>0,'1+',ptl),
         ptl=factor(ptl),
         ht=factor(ht,labels = c('no','yes')),
         ui=factor(ui,labels = c('no','yes')),
         ftv=ifelse(ftv>1,'2+',ftv),
         ftv=factor(ftv))
str(birthwtweight)

glm1 <- glm(low~age+lwt+race+smoke+ptl+ht+ui+ftv,
            family = binomial,
            data=birthwtweight)
summary(glm1)

drop1(glm1)

glm11 <- glm(low~age+lwt+race+smoke+ptl+ht+ui,
             family = binomial,
             data=birthwtweight)
drop1(glm11)

glm111 <- glm(low~lwt+race+smoke+ptl+ht+ui,
              family = binomial,
              data=birthwtweight)
drop1(glm111)

glm2 <- step(glm1,trace = FALSE)
summary(glm2)

anova(glm1,glm2,test = 'Chisq')
anova(glm1,glm2,test = 'LRT')

AIC(glm1,glm2)

coef(glm2)

exp(coef(glm2))

confint(glm2)

exp(confint(glm2))

library(epiDisplay)
logistic.display(glm2)

##无序多分类logistic回归################################################
library(epiDisplay)
data("Ectopic")
str(Ectopic)

library(nnet)
multi1 <- multinom(outc~hia,data = Ectopic)
summary(multi1)

st <- summary(multi1)$standard.errors
z <- coef(multi1)/st;z

p.values <- pnorm(abs(z),lower.tail = FALSE)*2

confint(multi1)

exp(coef(multi1))
exp(confint(multi1))

mlogit.display(multi1)

##有序logistic回归
dat <- array(c(10,7,19,6,0,2,7,5,1,5,6,16),
             dim = c(2,2,3),
             dimnames = list(method=c('old','new'),
                             sex=c('male','female'),
                             outcome=c('effectless','effective','recover')))
dat <- as.table(dat)
data1 <- as.data.frame(dat)
data1
str(data1)
data1$outcome <- ordered(data1$outcome)
data1$outcome

library(MASS)
polr1 <- polr(outcome~sex+method,weights = Freq,data = data1)
summary(polr1)

exp(coef(polr1))
exp(confint(polr1))

x <- seq(-5,5,length.out=100)
y <- dnorm(x,0,1)
plot(x,y,type = 'l',las=1)
abline(v=-2,lty='dashed')
abline(v=2,lty='dashed')