### https://www.kaggle.com/datasets/pentakrishnakishore/diabetes-csv

rm(list=ls())


Logistic_Residuals = function(X,Y,beta){
     n = length(Y)
     X = as.matrix(X)
     beta = matrix(beta)
     one = rep(1,n)
     X = cbind(one,X)
     L = X %*% matrix(beta)
     p = exp(L)/(1+exp(L))

     a = rep(0,n)
     b = rep(1,n)
     b[Y == 0] = 1-p[Y == 0]
     a[Y == 1] = 1-p[Y == 1]

     U = runif(n,a,b)      
     r = qnorm(U)
     return(r)
     }



D = read.csv("diabetes.csv",header=TRUE)
attach(D)
n = nrow(D)

##  names(D)
## [1] "Pregnancies"              "Glucose"                 
## [3] "BloodPressure"            "SkinThickness"           
## [5] "Insulin"                  "BMI"                     
## [7] "DiabetesPedigreeFunction" "Age"                     
## [9] "Outcome"      

out = glm(Outcome ~ Glucose,data=D,family="binomial")
## summary(out)


### deviance residuals first
pdf("LogisticResid_1.pdf")
dr = resid(out,type="deviance")
plot(Glucose,dr,pch=21,bg="pink",cex=2.2,ylab="",xlab="",cex.axis=1.4)  ## 
title(ylab="Glucose", line=2.1, cex.lab=2)
title(xlab="Residuals",line=2.6, cex.lab=2)
abline(h=0,col="red",lwd=2)
dev.off()


### ranomized quantile residuals
pdf("LogisticResid_2.pdf")

beta = out$coef
X = D[,-9]
Y = D[,9]

r = Logistic_Residuals(Glucose,Outcome,beta)
plot(Glucose,r,pch=21,bg="pink",cex=2,ylab="",xlab="",cex.axis=1.4)  ## 
title(ylab="Glucose", line=2.1, cex.lab=2)
title(xlab="Residuals",line=2.6, cex.lab=2)
abline(h=0,col="red",lwd=2)


dev.off()