## 10.3节图



################## 图10.3.1 ##################
curve( 1/(1+ exp(-x)), -10,10,lty=1, add=T, col = 'blue',lwd=2, xlab="   ", ylab="  ")  
abline( v = 0,   lwd = 3, lty = 2,col = 'red' )
abline( h = 1,   lwd = 2, lty = 3,col = 'gray60' )
abline( h = 0,   lwd = 2, lty = 3,col = 'gray60' )






################## 图10.3.2 ##################
# 创建绘图区域
plot(0, 0, type = "n", xlim = c(-6, 6), ylim = c(0, 1),
     ann = FALSE, xaxt = "n", yaxt = "n")
# 绘制曲线
curve(pnorm(x), -6, 6, lty = 1, add = TRUE, col = 'black', lwd = 3)
curve(1/(1 + exp(-x)), -6, 6, lty = 2, add = TRUE, col = 'red', lwd = 3)
curve(1 - exp(-exp(x)), -6, 6, lty = 3, add = TRUE, col = 'blue', lwd = 3)
legend(x = 2.1, y = 0.35,  # 向右上方移动，远离边缘
       legend = c("Normal", "Sigmoid", "Gumbel 生存"),
       col = c("black", "red", "blue"),
       lwd = 3,
       lty = c(1, 2, 3))  # 移除图例边框
# 添加坐标轴
axis(1)  # x轴
axis(2)  # y轴
# 添加斜体坐标轴标题
title(
  xlab = expression(italic(x)),  # x轴标题斜体
  ylab = expression(italic(y)),  # y轴标题斜体
  cex.lab = 1.2,  # 适当增大标题字体
  line=2.3)







################## 例10.3.1 ##################
data<- read.table("/Users/macbook/Desktop/线性模型例题/Chapter10/data/例10.3.1.txt",header=TRUE)
xx <- with(data, seq(min(x), max(x), len = 200))

fit <- glm(y ~ x, family = binomial(link = "logit"), data) ## Logistic-MLE
fit1 <- glm(y ~ 1, family = binomial(link = "logit"), data) ## Logistic-MLE
summary(fit)

# 图10.3.3
par(mar = c(5, 5, 4, 4) + 0.3)
plot(y ~ x, data, pch = 19, col = "gray40", xlab = expression(italic(X)), ylab = expression(hat(italic(Y))))
lines(xx, predict(loess( y ~  x,data ), data.frame(x = xx)),  lwd = 3, lty = 2, col = 'red')
lines(xx, predict(fit, data.frame(x = xx), type = "resp"), lwd = 3,lty = 1)
legend(list(x=22,y=0.3),    
       legend=c("Logistic ","LOWESS" ),  lwd = 2,      
       col=c("gray40","red"),                   
       lty=c(1,2))  
#title("Scatter Plot with Loess (blue) and Logistic Mean Response Functions")

## 拟合值和残差
fit <- glm(y ~ x, family = binomial(link = "logit"), data) ## Logistic-MLE
datars<-data.frame(Experience = data$x, TaskSuc= data$y, Fitted = fitted(fit), Residual=resid(fit))
datars[,3]<-round(datars[,3],3)
datars[,4]<-round(datars[,4],3) 
datars 

## 图10.3.4
xx <- with(data, seq(min(x), max(x), len = 200))
plot(y ~ x, data, pch = 19, col = "gray40", xlab = expression(italic(X)), ylab = expression(hat(italic(Y))))
#title("Logistic (Gray), Probit (Blue), and \nComplementary Log-Log (Red) Fitted Models")
fit <- glm(y ~ x, data = data, family = binomial(link = "probit"))
lines(xx, predict(fit, data.frame(x =xx), type = "resp"), col = 'red', lwd = 3, lty = 2)

fit <- glm(y ~ x, data = data, family = binomial(link = "cloglog"))
lines(xx, predict(fit, data.frame(x = xx), type = "resp"), col = 'blue', lwd = 3, lty = 3)

fit <- glm(y ~ x, family = binomial(link = "logit"), data)
lines(xx, predict(fit, data.frame(x = xx), type = "resp"), col = 'gray40', lwd = 3,lty = 1)

legend(list(x=24,y=0.3),    
       legend=c("Logistic "," probit", " cloglog"),        
       col=c("gray40","red","blue"), lwd = 3,                   
       lty=c(1,2,3))  







################## 例10.3.2 ##################
n1<-c(4,28,38,51,32)
n2<-c(70,147,207,202,92)

####列联表chi-square 独立性检验
N<-sum(n1)+sum(n2)
n<-n1+n2
m1<-n*sum(n1)/N
m2<-n*sum(n2)/N
M1<-diag(m1); M2<-diag(m2)
chisq4<- t(n1-m1)%*% solve(M1)%*% (n1-m1)+t(n2-m2)%*% solve(M2)%*% (n2-m2) 
round( chisq4,4) 
####列联表Wilks G^2 独立性检验
Gsq<-2*sum(n1*log(n1/m1))+2*sum(n2*log(n2/m2))
round(Gsq,4)

n<-n1+n2
pi<-n1/n
y<-log(n1/n2)
v<-(n/(n1*n2))
i<-c(1,2,3,4,5)
round(cbind(i,pi,y,v),4)    ###Table 10.3.3

C<-matrix(0, 4,5)
for (i  in 1:4 ){
  C[i,i]=1
  C[i,i+1]=-1}
V=diag(v) 
chi2<-t(y)%*%t(C)%*% solve(C %*%  V  %*%t(C))%*% C %*% y
c(chi2,qchisq(0.95,4), 1-pchisq(chi2,4)) ## chisquare test
c(y[2]-y[1]-qnorm(0.975)*(v[1]+v[2]), y[2]-y[1]+qnorm(0.975)*(v[1]+v[2]))

## WLSE
x<-c(35,45,55,65,75)
v<-round(v,4)
fit<-lm(y~x, weight=v)
summar<-summary(fit) 
b<- round(summar$coefficients[,1],4)
yhat<- b[1]+b[2]*x
pihat<- round(exp(yhat)/(1+exp(yhat)),3)
y<-round(y,4)
pi<-round(n1/n,3)
expected_n1<-n*pihat
results<-cbind(x, y, yhat, pi, pihat,  expected_n1, n1 )
results

## MLE
fit <- glm(y~x, family = binomial(link = "logit") )








################## 例10.3.7 ##################
data<- read.table("/Users/macbook/Desktop/第十章/data/10.3.7案例分析.txt",header=TRUE)

reg <- glm(Y~X1+X2+X3+X4,  data, family = binomial(link = "logit") )
summary(reg)

X=model.matrix(reg)  ## 设计阵
pi=predict(reg, type="response") ## pi估计
w<-pi*(1-pi)
Cov<-solve(t(X)%*% diag(w)%*% X)  ##beta MLE 的协方差
H<- diag(sqrt(w))%*%X%*% Cov%*%t(X)%*%diag(sqrt(w))
round(Cov, 4)
#round(H, 4)
round(exp(reg$coef[-1]),4) ## odds ratio


logit<-log(pi/(1-pi))  ###y*=logit(pi)
reg1 <- glm(Y~X1+X4,  data, family = binomial(link = "logit") )

##########tests
## z统计量的计算
s<-c(sqrt(Cov[1,1]), sqrt(Cov[2,2]), sqrt(Cov[3,3]),sqrt(Cov[4,4]),sqrt(Cov[5,5]))
coefi<-reg$coefficients
z<-coefi/s   ###z统计量的计算
2*(1-pnorm(abs(z))) ## p-value 

## 多个回归系数似然比检验
reg0 <- glm(Y~X1+X4,  data, family = binomial(link = "logit") )
pi0=predict(reg0, type="response")

l1<-sum(data$Y*log(pi)+(1-data$Y)*log(1-pi))  ###对数似然函数
l0<-sum(data$Y*log(pi0)+(1-data$Y)*log(1-pi0))
G2<--2*(l0-l1)
p=1-pchisq(G2,2)
c(l1,l0, G2,p)  #results
#-50.5270751 -51.1296645   1.2051787   0.5473924

## 逐步选择
step(reg,  direction='both'  )   ## 默认为AIC  k=2

step(reg,  direction='both', k=log(nrow(data)))  ##BIC

########################残差分析
H<- diag(sqrt(w))%*%X%*% Cov%*%t(X)%*%diag(sqrt(w))##帽子矩阵
h<-diag(H)
e<-data$Y- pi0      ##ordinal residual  
rp<-e/sqrt(w)      #Pearson residual
rsp<-rp/sqrt(1-h)  # standard Pearson residual
d<- -2*(data$Y*log(pi0)+ (1-data$Y )*log(1-pi0)) 
dev<-sign(e)*sqrt(d) ## deviance residual

############# glm自带的残差函数 
e<-residuals(reg0, type = "response") 
p<-  residuals(reg0, type = "pearson") #Pearson residual
dev<- residuals(reg0, type = "deviance") 



## 绘图
## 图10.3.5
par(
  mfrow = c(2, 2),  # 2x2布局
  ann = TRUE,       # 允许自动添加标签
  mar = c(3, 3, 2, 1),   # 减小每个子图的边距（下,左,上,右）
  oma = c(2, 2, 1, 1) ,  # 减小外边界
  mgp = c(1.8, 0.8, 0)        # 调整标签与坐标轴的距离（标题,标签,轴线）
)
plot(e~pi0, lwd = 2, pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(italic(hat(e)[i])))
plot(rp~pi0, lwd = 2, pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(italic(r["Pi"])))
plot(rsp~pi0,lwd = 2,  pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(italic(r[sPi])))
plot(dev~pi0, lwd = 2, pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(dev[i]))


logit<-log(pi0/(1-pi0)) ##
## 图10.3.6
par(
  mfrow = c(2, 2),  # 2x2布局
  ann = TRUE,       # 允许自动添加标签
  mar = c(3, 3, 2, 1),   # 减小每个子图的边距（下,左,上,右）
  oma = c(2, 2, 1, 1) ,  # 减小外边界
  mgp = c(1.8, 0.8, 0)        # 调整标签与坐标轴的距离（标题,标签,轴线）
)   
plot(rsp~pi0,  pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(italic(r[sPi])))
lines(lowess(pi0,rsp), col="red", lwd=2)
#abline(h = 0,  ity=3 , col = 'gray', lwd = 2)
plot(rsp~logit,  pch = "o", xlab = expression(logit(hat(pi)[i])), ylab = expression(dev[i]))
lines(lowess(logit,rsp), col="red", lwd=2)

plot(dev~pi0,  pch = "o", xlab = expression(italic(hat(pi)[i])), ylab = expression(italic(r[sPi])))
lines(lowess(pi0,dev), col="red", lwd=2)

plot(dev~logit,  pch = "o", xlab = expression(logit(hat(pi)[i])), ylab = expression(dev[i]))
lines(lowess(logit,dev), col="red", lwd=2)


## 图10.3.7
install.packages("hnp")
library(MASS)
library(hnp)
par(mfrow=c(1,1),ann=TRUE, mar = c(6, 6, 4, 2) + 0.1, mgp = c(3, 1, 0))
my.hnp<-hnp(reg0, sim = 99, conf = 0.95, resid.type="deviance",
            halfnormal = T, scale = F, how.many.out = T ) 
plot(my.hnp,  xlab = expression(Phi^-1 * group("(", (k + n - 1/8)/(2 * n + 1/2), ")")),
     ylab=expression(abs(dev)["s(k)"]), legpos="bottomright", lwd = 3) 



########################影响分析
reg0 <- glm(Y~X1+X4,  data, family = binomial(link = "logit") )
pi0=predict(reg0, type="response")
X <- model.matrix(reg0)
y <- model.response(model.frame(reg0))

##### Cook distance
w<-pi0*(1-pi0)
Cov<-solve(t(X)%*% diag(w) %*% X)   # 协方差
H<- diag(sqrt(w))%*%X%*% Cov%*%t(X)%*%diag(sqrt(w))
h<- diag(H) #杠杆值
rp<-residuals(reg0, type = "pearson")
cook.D<-h*rp^2/(ncol(X))

########Crossvalidation (CV) errors
#####function for chi^2
chisq<-function(y,pi){
  c<-sum((y-pi)^2/(pi*(1-pi)))
  return(c)
}
#####functions for DEV
DEV<-function(y,pi){   
  d<- - 2*(y*log(pi)+(1-y)*log(1-pi)) 
  return(sum(d))
}

#######Leave-one-out CV
DEVcv<-rep(0, nrow(X))
chisqcv<-rep(0, nrow(X))
for (i in 1:nrow(X))           ## Leave-one-out CV
{
  regi<-glm(y[-i]~X[-i,1]+X[-i,2], family = binomial(link = "logit"))
  pi_i=predict(regi, type="response")
  chisqcv[i]<-chisq(y[-i],pi_i)
  DEVcv[i]<-DEV(y[-i],pi_i)
}
deltachisqcv=chisq(y,pi0)-chisqcv
delta_dev_cv=DEV(y,pi0)-DEVcv
id<-seq(1:nrow(X))
result<-cbind(deltachisqcv, delta_dev_cv, cook.D, h) ###影响三个统计量的值
round(result,3)

c(which(deltachisqcv==max(deltachisqcv)), ###各自最大值对应的id
  which(delta_dev_cv==max(delta_dev_cv)),
  which(cook.D==max(cook.D)),
  which(h==max(h)))


## 图10.3.8
par(mfrow = c(2, 2), mar = c(5, 5, 2, 2) )  
plot(deltachisqcv~id, lwd = 2, type = "o", xlab =  expression(italic(i)), ylab = expression(Delta * chi[i]^2))
text(32, max(deltachisqcv) - 0.03, paste("点", 19))
plot(delta_dev_cv~id,lwd = 2, type = "o",  xlab = expression(italic(i)), ylab = expression(Delta * dev[i]))
text(32, max(delta_dev_cv) - 0.03, paste("点", 19))
plot(cook.D~id,  lwd = 2,type = "o" , xlab = expression(italic(i)), ylab = expression(italic(D[i])))
text(32, max(cook.D)-0.003, paste("点", 48))
plot(h~id, lwd = 2, type = "o", xlab = expression(italic(i)), ylab = expression(italic(h[(ii)])))
text(32, max(h)-0.003, paste("点", 83))



###########列线图
install.packages("rms",type="binary")
install.packages("DynNom")
install.packages("regplot")
install.packages("htmlTable")
install.packages("Hmisc",type="binary")
rm(list=ls())
library(rms)
library(Hmisc)
library(lattice)
library(survival)
library(ggplot2)

data<- read.table("/Users/macbook/Desktop/第十章/data/10.3.7案例分析.txt",header=TRUE)
data1<-data.frame(data)#生成数据框

dim(data1)#查看数据（98行6列）
str(data1)#查看变量 #打包数据
dd <- datadist(data1)
options(datadist="dd")#显示有效位数，默认7位
fit1 <- lrm(Y ~ X1+X4,data = data1, x=T,y=T)#逻辑回归
par(ann=FALSE)
nom1 <- nomogram(fit1, fun=plogis,
                 fun.at=c(0.001,0.1,0.25,0.5,0.75,0.9,0.99),#刻度
                 lp=T, # 显示线性概率
                 funlabel="diseaseoutbreak") #命名  
# nom1$X1$label <- "指标1（斜体可结合expression）"
# nom1$X4$label <- "指标4（斜体可结合expression）"
par(mfrow=c(1,1))
plot(nom1)                
      

help(nomogram)

##计算(X1=20, X4=1)处的响应概率
new_data <- data.frame(X1 = 20, X4 = 1)
predict(reg0, newdata = new_data, type = "response")


p201<-1/(1+exp(2.33515-0.02929*20-1.67345*1))
logit201<-log(p201/(1-p201))
c(p201, logit201)
L<-logit201-qnorm(0.975)* sqrt(t(c(1,20,1))%*% Cov %*% c(1,20,1)) 
U<-logit201+qnorm(0.975)* sqrt(t(c(1,20,1))%*% Cov %*% c(1,20,1))  

Lpi<-1/(1+exp(-L))
Upi<-1/(1+exp(-U))

c(L,U,Lpi,Upi)


#####计算ROC 曲线
library(ggplot2)
install.packages("ROCR")
install.packages("gplots")
install.packages("pROC")
library(ROCR) 
library(ggplot2)
data<- read.table("/Users/macbook/Desktop/第十章/data/10.3.7案例分析.txt",header=TRUE)
data1<-data.frame(data) #生成数据框
#随机分组，其中70%作为训练集，30%作为测试集
set.seed(1)
train_data1<-sample(98,69) 
train<-data1[train_data1,]
test<-data1[-train_data1,]


## 图10.3.1
library(pROC)  
fit0<-glm(Y~X1+X4,data=train,family = binomial)#family = binomial()表示连接函数为二项分布
prob0<-predict(fit0,newdata=test,type="response")#使用测试集进行预测
p1 <- roc(test$Y, prob0)
plot(p1,
     legacy.axes = TRUE,
     main="ROC 曲线",
     xlab = "1-特异度", 
     ylab = "敏感度",
     thresholds="best", # 基于youden指数选择roc曲线最佳阈值点
     print.thres="best", lwd=3, col = "red")

p1$auc ##AUC 的值
ci.auc(p1)##AUC 的95%置信区间

