الانحدار اللوجستي في R: ‏glm() مع family = binomial

حين يكون المتغير التابع نعم/لا

يتنبأ الانحدار الخطي برقم. لكن كثيرًا من الأسئلة الجديرة بالنمذجة ثنائية: هل يغادر العميل، وهل يتعافى المريض، وهل يُنقر على البريد الإلكتروني. ومعايرة خط مستقيم لمتغير تابع قيمته 0/1 تنهار فورًا، إذ يتنبأ الخط بمرح باحتمالات مثل −0.3 أو 1.4، وهي هراء.

ويحل الانحدار اللوجستي ذلك بنمذجة احتمال الحدث عبر تحويل اللوغاريتم الأرجحي (logit): ‏log(p / (1 − p)) = intercept + slope × x. ويمتد مقياس اللوغاريتم الأرجحي على خط الأعداد كله، فتلائمه معادلة خطية ملاءمة طبيعية، ويضغط التحويل العكسي كل تنبؤ داخل المجال (0, 1) على امتداد المنحنى المألوف بشكل حرف S. وثمن الحيلة أن المعاملات تعيش على مقياس اللوغاريتم الأرجحي، وأن لعبة قراءة الانحدار اللوجستي كلها هي ترجمتها إلى شيء يفهمه البشر.

المعايرة: ‏glm() مع family = binomial

‏glm() (النموذج الخطي المعمم) هي الشقيقة الكبرى لـ lm()، ويختار family = binomial الانحدار اللوجستي. وفي mtcars يسجّل am نوع ناقل الحركة (1 = يدوي، 0 = أوتوماتيكي)، فهل تميل السيارات الموفرة للوقود إلى أن تكون يدوية؟

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

وأمران قبل قراءة المخرجات. أولًا، family = binomial ليست اختيارية، فإذا حذفتها عايرت glm() بصمت انحدار المربعات الصغرى العادي. ثانيًا، يجب أن يكون المتغير التابع ثنائيًا: ‏0/1 أو منطقيًا أو عاملًا بمستويين (وينمذج R احتمال المستوى الثاني).

والآن الملخص كتلة كتلة:

• ‏Coefficients: قيمة Estimate لـ mpg تبلغ نحو 0.31، وهي ميل باللوغاريتم الأرجحي: فكل ميل إضافي لكل غالون يضيف 0.31 إلى اللوغاريتم الأرجحي لكون ناقل الحركة يدويًا. والموجب يعني "يرفع الاحتمال" والسالب يعني "يخفضه"، وبعد الإشارة لا يعمل حدس أحد على هذا المقياس، ولهذا وُجد القسم التالي.
• قيمة z وPr(>|z|): المنطق نفسه المتبع في اختبارات t في الانحدار (Estimate ÷ Std. Error ثم قيمة p لسؤال "هل يمكن أن يكون هذا صفرًا؟")، لكن باستخدام تقريب طبيعي، ومن هنا z بدل t. وهنا p ≈ 0.011: فمن غير المرجّح أن يكون الارتباط بين اقتصاد الوقود ونوع ناقل الحركة ضوضاء.
• ‏Null deviance مقابل Residual deviance: الانحراف هو مقياس سوء الملاءمة في عالم glm (والأصغر أفضل). وNull deviance (43.2 عند 31 درجة حرية) هو نموذج الحد الثابت وحده، وResidual deviance (29.7 عند 30 درجة حرية) هو نموذجك. والانخفاض البالغ نحو 13.6 مقابل درجة حرية واحدة لمتنبئ واحد هو نظير "ارتفع معامل التحديد" في عالم glm.
• ‏AIC: درجة لمقارنة النماذج توازن بين الملاءمة والتعقيد، والأدنى يفوز. وهي بلا معنى منفردة، ومفيدة بين نماذج مرشحة على البيانات نفسها.

من اللوغاريتم الأرجحي إلى نسب الأرجحية: ‏exp(coef())

يَنقل الرفع الأسّي المعاملات من مقياس اللوغاريتم الأرجحي الجمعي إلى مقياس الأرجحية الضربي:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

‏exp(0.307) ≈ 1.36، وهذا هو نمط الجملة الصادق الجدير بالحفظ: "كل ميل إضافي لكل غالون يضرب أرجحية أن يكون ناقل الحركة يدويًا بنحو 1.36". ونسبة الأرجحية فوق 1 ترفع الأرجحية، وتحتها تخفضها، والقيمة 1 بالضبط تعني لا أثر، ولهذا يكون حكم فترة الثقة لنسب الأرجحية هو "هل تستبعد الفترة القيمة 1؟" (لا الصفر، فالصفر كان الحد على مقياس اللوغاريتم الأرجحي).

وانتبه إلى اللغة: الأرجحية ليست احتمالًا. فالأرجحية = p / (1 − p)، لذا فإن الاحتمال 0.75 يقابل أرجحية 3. وضرب الأرجحية في 1.36 ليس كضرب الاحتمال في 1.36، والفجوة كبيرة حين يكون الحدث شائعًا. فنسبة أرجحية تساوي 2 لحدث نادر تتصرف كـ "مضاعفة الخطر تقريبًا"، أما لحدث معدله 50% فلا تفعل ذلك إطلاقًا. فلا تُبلغ أبدًا عن نسبة أرجحية بصياغة نسبة الخطر ("أرجح بمقدار 1.36 مرة") ما لم يكن الحدث نادرًا.

الاحتمالات المتنبأ بها: مصيدة type = "response"

أشيع خلل في الانحدار اللوجستي على أرض الواقع:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

يُرجع الاستدعاء الأول القيمة الافتراضية type = "link"، أي تنبؤات على مقياس اللوغاريتم الأرجحي بقيم سالبة وغيرها. ويُرجع الثاني احتمالات حقيقية. فإذا خرجت "احتمالاتك" يومًا سالبة أو أكبر من 1، فهذا هو السبب. شغّل الكتلة: فالسيارة التي تحقق 15 ميلًا لكل غالون لا فرصة لها عمليًا لأن تكون يدوية، والسيارة التي تحقق 30 يرجَّح جدًا أن تكون يدوية، وينثني منحنى S عبر المنتصف.

التصنيف: العتبة وجدول الالتباس

تصبح الاحتمالات فئات متنبأ بها باختيار عتبة قطع، وقيمة 0.5 هي الخيار الافتراضي، وبطاقة النتائج الصادقة جدول بالمتنبأ به مقابل الفعلي:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

خلايا القطر هي القرارات الصحيحة، والخليتان خارج القطر هما الخطآن المختلفان (التنبؤ بيدوي لسيارة أوتوماتيكية والعكس). والدقة الإجمالية وحدها قد تجمّل نموذجًا تجميلًا سيئًا، فإذا كان 95% من العملاء لا يغادرون فإن "تنبأ بألا أحد يغادر" يسجّل 95% دون التقاط أي مغادر، لذا انظر دائمًا إلى نوعَي الخطأ. والقيمة 0.5 عرف لا قانون: فحين تختلف تكلفة الخطأين، حرّك العتبة تبعًا لذلك.

وثمة تحفظ من باب الصدق: هذا الجدول يقيّم النموذج على البيانات نفسها التي عُوير عليها، وهو ما يجمّله. والتقييم الحقيقي يحجب بيانات لم يرها النموذج قط.

متنبئات متعددة

تمامًا كما في lm(): أضف الحدود بـ +، ويكتسب كل تفسير القيد "مع تثبيت البقية":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

وكل معامل مرفوع أسّيًا صار الآن مضاعِف الأرجحية لزيادة وحدة واحدة في ذلك المتنبئ بين السيارات المتشابهة في المتنبئات الأخرى. وتتوسع الآلية، لكن تتوسع معها تحفظات الانحدار الخطي: فالمتنبئات المترابطة تعيد خلط معاملات بعضها بعضًا.

تحذيرات

• الفصل التام. إذا قسّم متنبئ ما المتغير التابع تقسيمًا مثاليًا (كأن تكون كل سيارة فوق حد معين من اقتصاد الوقود يدوية وكل ما دونه أوتوماتيكيًا)، فإن معامل الإمكان الأعظم يريد أن يكون لانهائيًا. ويحذّر R بالرسالة glm.fit: fitted probabilities numerically 0 or 1 occurred ويُبلغ عن معاملات ضخمة بأخطاء معيارية سخيفة. فلا تسلّم تلك الأرقام، بل بسّط النموذج، أو احصل على بيانات أكثر، أو استخدم طريقة معاقَبة (كحزمتَي brglm2 أو logistf).
• أحداث كافية. القيد الملزم هو عدد الحالات في الفئة الأندر لا مجموع الصفوف. وتطلب قاعدة عملية قديمة نحو 10 إلى 15 حدثًا لكل متنبئ، أما أمثلة السيارات الاثنتين والثلاثين هنا فهي لتعليم الآلية لا قالبًا لأحجام عينات قابلة للنشر.
• نسب الأرجحية ليست نسب خطر حين يكون الحدث شائعًا، وقد غُطّي ذلك أعلاه، ويُكرَّر لأن المحكّمين سيلتقطونه حتى لو لم تلتقطه أنت.

الخلاصة

• متغير تابع ثنائي ← glm(y ~ x, data = df, family = binomial)، ولا تنسَ family أبدًا.
• المعاملات الخام لوغاريتمات أرجحية، وتعطي exp(coef(fit)) نسب الأرجحية، والقيمة العدمية لفتراتها هي 1.
• نمط الجملة: "كل زيادة بوحدة واحدة في x تضرب أرجحية وقوع الحدث بمقدار exp(b)".
• استخدم predict(..., type = "response") للاحتمالات، فالقيمة الافتراضية تُرجع لوغاريتمات أرجحية، وهي مصدر الالتباس الأول.
• صنّف باستخدام عتبة واحكم باستخدام جدول التباس، فالدقة وحدها قد تكذب.
• انتبه إلى تحذيرات الفصل، وعُدّ أحداثك، ولا تُلبس نسب الأرجحية ثوب نسب الخطر.

في المستند التالي: الآلية الكامنة خلف كل فترة رأيتها حتى الآن، أي فترات الثقة عبر t.test() وconfint() وprop.test().