الانحدار اللوجستي (logistic regression) باستخدام لغة البرمجة آر R

الانحدار اللوجستي (logistic regression) هو أسلوب إحصائي يستخدم لوصف العلاقة بين متغير تابع (dependent variable) ثنائي (binary) ومتغير مستقل (independent variable) واحد او أكثر وفق الصيغة التالية:

\[
logit(p)=\alpha+\beta X
\]

حيث

\[
logit(p)=\ln{\left(\frac{p}{1-p}\right)}
\]

حيث $p=p(y=1|X)$. وهناك عدد من الفروض (assumptions) الهامة التي ينبغي التنبه لها عند استخدام الانحدار اللوجستي وهي:

الآن سوف نوضح كيف يمكن تطبيق نموذج الانحدار اللوجستي على البيانات (PimaIndiansDiabetes) بإستخدام لغة البرمجة آر R وفق الخطوات التالية:

1. تحميل البيانات

In:
library(mlbench)
data(PimaIndiansDiabetes)
head(PimaIndiansDiabetes)
Out:
  pregnant glucose pressure triceps insulin mass pedigree age diabetes
1        6     148       72      35       0 33.6    0.627  50      pos
2        1      85       66      29       0 26.6    0.351  31      neg
3        8     183       64       0       0 23.3    0.672  32      pos
4        1      89       66      23      94 28.1    0.167  21      neg
5        0     137       40      35     168 43.1    2.288  33      pos
6        5     116       74       0       0 25.6    0.201  30      neg

في هذه الحالة سوف يكون المتغير التابع هو diabetes وهو ثنائي كما يتضح، وجميع المتغيرات المتبقية هي متغيرات مستقلة.

2. تقسيم البيانات الى مجموعة تدريب (train) ومجموعة إختبار (test)

In:
m=nrow(PimaIndiansDiabetes)
train_i= sample(seq_len(m),size =0.75*m)  
train_data =PimaIndiansDiabetes[train_i,]
test_data =PimaIndiansDiabetes[-train_i,]

 3. تطبيق نموذج الانحدار اللوجستي

In:
Model = glm(diabetes ~., data = train_data,family = binomial) 
summary(Model)
Out:
Call:
glm(formula = diabetes ~ ., family = binomial, data = train_data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4879  -0.7398  -0.4308   0.7455   2.8940  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.103335   0.813837  -9.957  < 2e-16 ***
pregnant     0.114101   0.037850   3.015  0.00257 ** 
glucose      0.037800   0.004364   8.662  < 2e-16 ***
pressure    -0.012371   0.005913  -2.092  0.03642 *  
triceps      0.003894   0.008233   0.473  0.63622    
insulin     -0.001816   0.001065  -1.705  0.08826 .  
mass         0.074861   0.017213   4.349 1.37e-05 ***
pedigree     0.831963   0.337920   2.462  0.01382 *  
age          0.009841   0.010634   0.925  0.35473    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 742.59  on 575  degrees of freedom
Residual deviance: 548.64  on 567  degrees of freedom
AIC: 566.64

Number of Fisher Scoring iterations: 5

حيث يجب تحديد المتغيرات المستقله الهامة التي يجب ان تكون في النموذج وحذف الأخرى للوصول لأفضل نموذج يجب تطبيقة.

 4. تشخيص النموذج (Model Diagnostics) كما يلي:

• التأكد من الخطية (linearity) بين المتغيرات المستقله والمتغير التابع على logit scale.
  

In:
library(tidyverse)
probability = predict(Model, train_data, type="response")
logit=log(probability/(1-probability))
newdata=cbind(logit,train_data[,-9])
newdata = newdata %>%
  gather(key = "predictors", value = "predicted_values", -logit)
ggplot(newdata, aes(predicted_values,logit))+
  geom_point(size = 0.6, alpha = 0.6) +
  theme_bw() + 
  xlab("predicted values")+
  facet_wrap(~predictors, scales = "free_x")

• التأكد من عدم وجود المشاهدات المتطرفة من خلال شكل الإنتشار للقيم التنبؤية (predictive values) ضد البواقي المعيارية (Standardized residuals).
  

In:
library(broom)
Model.out = augment(Model)
ggplot(Model.out, aes(.fitted, .std.resid)) + 
  geom_point(aes(color = diabetes), alpha = .5) +
  xlab("Predicted values")+
  ylab("Standardized residuals")+
  theme_bw()

وذلك من خلال التأكد من عدم وجود مشاهدات لها بواقي معيارية خارج النطاق 3±.

• التأكد من عدم وجود خطية متعددة (multicollinearity).
  

In:
library(car)
vif(Model)
Out:
pregnant  glucose pressure  triceps  insulin     mass pedigree      age 
1.428726 1.248356 1.222852 1.482516 1.472432 1.306988 1.043258 1.505727 

وذلك من خلال حساب التضخم في التباين (variance-inflation)، حيث أن تجاوز 10 يشير إلى وجود علاقة خطية متعددة.

• الأرتباط الذاتي (autocorrelation).

In:
acf(Model.out$.resid,main='Autocorrelation')

• مصفوفة الخطأ (confiusion matrix)

In:
predicted_probability = predict(Model, test_data, type="response") 
predicted_values=ifelse(predicted_probability>0.5,1,0)
actual_values=ifelse(test_data$diabetes=="pos",1,0) 
confiusion_matrix=table(predicted_values,actual_values) 
confiusion_matrix
Out:
                actual_values
predicted_values   0   1
               0 112  34
               1  11  35

• الحساسية (sensitivity ) والنوعية (specificity).

In:
library(caret)
sensitivity(confiusion_matrix)
Out:
[1] 0.9105691
In:
specificity(confiusion_matrix)
Out:
[1] 0.5072464