Lesson 25 - Logistic Regression

Load Packages and Data

require(nnet)

Key to Variables

• CC - Categorical variable. 1 = Overt diabetes, 2 = Chemical diabetes, 3 - Non-diabetic
• RW - Relative weight
• IR - Insulin response
• SSPG - Steady state plasma glucose

df <- read.table("Diabetes.txt", sep="\t", header=TRUE)
summary(df)

       RW               IR             SSPG             CC       
 Min.   :0.7100   Min.   : 10.0   Min.   : 29.0   Min.   :1.000  
 1st Qu.:0.8800   1st Qu.:118.0   1st Qu.:100.0   1st Qu.:2.000  
 Median :0.9800   Median :156.0   Median :159.0   Median :3.000  
 Mean   :0.9773   Mean   :186.1   Mean   :184.2   Mean   :2.297  
 3rd Qu.:1.0800   3rd Qu.:221.0   3rd Qu.:257.0   3rd Qu.:3.000  
 Max.   :1.2000   Max.   :748.0   Max.   :480.0   Max.   :3.000  

df$CC <- factor(df$CC, levels= c("3", "2", "1"))
summary(df)

       RW               IR             SSPG       CC    
 Min.   :0.7100   Min.   : 10.0   Min.   : 29.0   3:76  
 1st Qu.:0.8800   1st Qu.:118.0   1st Qu.:100.0   2:36  
 Median :0.9800   Median :156.0   Median :159.0   1:33  
 Mean   :0.9773   Mean   :186.1   Mean   :184.2         
 3rd Qu.:1.0800   3rd Qu.:221.0   3rd Qu.:257.0         
 Max.   :1.2000   Max.   :748.0   Max.   :480.0         

Generate Boxplots

par(mfrow=c(1,3))
plot(RW ~ CC, df)
plot(IR ~ CC, df)

plot(SSPG ~ CC, df)
par(mfrow=c(1,1))

Create Model

mod <- multinom(CC ~ RW + IR + SSPG, df)

# weights:  15 (8 variable)
initial  value 159.298782 
iter  10 value 69.027793
iter  20 value 68.418245
iter  30 value 68.414665
final  value 68.414644 
converged

summary(mod)

Call:
multinom(formula = CC ~ RW + IR + SSPG, data = df)

Coefficients:
  (Intercept)        RW           IR       SSPG
2   -7.615261  3.472572  0.003586749 0.01641449
1   -1.845230 -5.867196 -0.013353688 0.04550552

Std. Errors:
  (Intercept)       RW          IR        SSPG
2    2.335615 2.446151 0.002349168 0.004981886
1    3.463507 3.866580 0.005019289 0.009241721

Residual Deviance: 136.8293 
AIC: 152.8293 

Generate Predictions

nd <- data.frame(
  RW = c(0.87, 0.91, 0.85, 1.15),
  IR = c(120, 620, 150, 150),
  SSPG = c(150, 160, 230, 230)  )
predict(mod, nd, type="probs")

          3         2            1
1 0.7351087 0.1340857 0.1308056224
2 0.4024971 0.5973904 0.0001124671
3 0.1467817 0.1034117 0.7498066056
4 0.2580303 0.5152335 0.2267361597

predict(mod, nd, type="class")

[1] 3 2 1 2
Levels: 3 2 1

Generating Predictions the Hard Way

x <- as.numeric(nd[1,])
x

[1]   0.87 120.00 150.00

cf <- summary(mod)$coefficients
cf

  (Intercept)        RW           IR       SSPG
2   -7.615261  3.472572  0.003586749 0.01641449
1   -1.845230 -5.867196 -0.013353688 0.04550552

logodds23 <- cf[1,1]  + sum(x * cf[1,2:4])
logodds13 <- cf[2,1]  + sum(x * cf[2,2:4])
c(logodds13, logodds23)

[1] -1.726306 -1.701539

odds13 <- exp(logodds13)
odds23 <- exp(logodds23)
c(odds13, odds23)

[1] 0.1779405 0.1824026

p1 <- odds13 / (odds13 + odds23 + 1)
p2 <- odds23 / (odds13 + odds23 + 1)
p3 <- 1 - p1 - p2
c(p3, p2, p1)

[1] 0.7351087 0.1340857 0.1308056