النماذج متعددة المستويات والبايزية: lme4 وbrms وبيانات البانل في R

بياناتك عن طلبة، لكن الطلبة موزعون على مدارس، والمدارس على مديريات. أو عن أسر تعيش في مناطق تختلف في خدماتها وأسواقها. أو عن قياسات متكررة للأشخاص أنفسهم عبر الزمن. في كل هذه الحالات تخرق البيانات افتراض الاستقلال الذي يقوم عليه الانحدار العادي، فتأتي الأخطاء المعيارية أصغر مما ينبغي، وتظهر علاقات دالة لا تصمد. النماذج متعددة المستويات (Multilevel Models)، وتسمى أيضًا الهرمية أو المختلطة، صممت لهذا النوع من البيانات، ونسختها البايزية تضيف طريقة مباشرة للتعبير عن عدم اليقين. يشرح هذا المقال متى تحتاج إليها، وكيف تقرأ معامل الارتباط داخل الصنف، وكيف تطبق المقاطع والميول العشوائية في lme4، ثم النسخة البايزية في brms مع فحوص التقارب، وأخيرًا التأثيرات الثابتة والعشوائية في بيانات البانل بحزمة plm، ومتى يناسب كل خيار.

لماذا لا يكفي الانحدار العادي مع البيانات المتداخلة؟

الطلبة في المدرسة نفسها يتشاركون المعلمين والموارد والبيئة، فتتشابه نتائجهم أكثر من طلبة مدارس مختلفة. هذا التشابه يعني أن ألف طالب من عشرين مدرسة لا يحملون معلومات ألف مشاهدة مستقلة. وتجاهل البنية يؤدي إلى ثلاث مشكلات:

• أخطاء معيارية متحيزة إلى الأسفل: خاصة لمتغيرات المستوى الأعلى، مثل نوع المدرسة، فتكثر الدلالة الزائفة.
• خلط المستويات: العلاقة بين متغيرين داخل المدارس قد تختلف عن العلاقة بين متوسطات المدارس، بل قد تعاكسها. وتفسير علاقة على مستوى المجموعات كأنها علاقة فردية هو المغالطة البيئية (Ecological Fallacy).
• إهدار المعلومات: تقدير نموذج منفصل لكل مدرسة يعطي تقديرات غير مستقرة للمدارس الصغيرة، ودمج الجميع في نموذج واحد يتجاهل الاختلاف بينها.

يقع النموذج متعدد المستويات بين الطرفين: يقدر لكل مجموعة انحرافها عن المتوسط العام، لكنه يسحب تقديرات المجموعات الصغيرة أو المتذبذبة نحو هذا المتوسط. يسمى ذلك التجميع الجزئي (Partial Pooling) أو الانكماش (Shrinkage)، وهو مصدر قوة هذه النماذج في التقدير والتنبؤ معًا.

معامل الارتباط داخل الصنف (ICC)

أول نموذج تقدره هو النموذج الصفري بلا متغيرات مستقلة، لتقسيم التباين الكلي إلى جزء بين المجموعات وجزء داخلها. معامل الارتباط داخل الصنف هو نسبة التباين بين المجموعات إلى التباين الكلي: ICC = τ₀₀ / (τ₀₀ + σ²)، وهو أيضًا الارتباط المتوقع بين وحدتين من المجموعة نفسها.

نستخدم بيانات sleepstudy المرفقة بحزمة lme4: متوسط زمن الاستجابة بالمللي ثانية لـ18 مشاركًا قيس يوميًا خلال عشرة أيام في دراسة عن تقييد النوم، فالقياسات متداخلة داخل الأشخاص.

library(lme4)
data("sleepstudy", package = "lme4")    # 180 rows: 18 subjects x 10 days

# Null (intercept-only) model with a random intercept per subject
m0 <- lmer(Reaction ~ 1 + (1 | Subject), data = sleepstudy)

vc  <- as.data.frame(VarCorr(m0))
icc <- vc$vcov[vc$grp == "Subject"] / sum(vc$vcov)
icc                                      # same as performance::icc(m0)

# Design effect for clusters of equal size m: 1 + (m - 1) * ICC
1 + (10 - 1) * icc

أثر التصميم (Design Effect) يترجم ICC إلى لغة حجم العينة. مثال حسابي افتراضي: إذا كان ICC يساوي 0.05 وفي كل مدرسة 30 طالبًا، فأثر التصميم 1 + 29 × 0.05 = 2.45، أي أن حجم العينة الفعلي أقل من نصف الحجم الاسمي. لذلك لا تحكم بأن ICC الصغير يعني الاستغناء عن النموذج متعدد المستويات، خاصة مع مجموعات كبيرة، وأدخل أثر التصميم في حساب حجم العينة حين تخطط لعينة عنقودية.

المقاطع العشوائية والميول العشوائية في lme4

في نموذج المقطع العشوائي (Random Intercept) لكل شخص مستوى أساسي خاص به، لكن أثر الأيام واحد للجميع. وفي نموذج الميل العشوائي (Random Slope) يسمح لأثر الأيام نفسه بأن يختلف من شخص لآخر، مع تقدير الارتباط بين المقطع والميل.

library(lmerTest)    # adds Satterthwaite df and p-values to lmer() output

m1 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
m2 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)

anova(m1, m2)        # LRT (refits with ML); conservative for variance terms
summary(m2)          # fixed effects, variance components, correlation
confint(m2, parm = "beta_", method = "profile")   # CIs for fixed effects
head(coef(m2)$Subject)   # subject-specific (shrunken) intercepts and slopes
isSingular(m2)       # TRUE would signal an over-complex random structure

ثلاث نقاط تفصل التحليل المتقن عن الآلي. الأولى أن lme4 يقدر بطريقة REML افتراضيًا، وهي مناسبة لمكونات التباين، لكن مقارنة نماذج تختلف في التأثيرات الثابتة تحتاج إلى ML، وanova() يعيد التقدير بها تلقائيًا. واختبار نسبة الاحتمال لمكون تباين محافظ، لأن القيمة الصفرية تقع على حد فضاء المعالم، فقيمة p الناتجة أكبر من الصحيحة تقريبًا.

الثانية تمركز المتغيرات. في مثال الطلبة، يختلط أثر المستوى الاجتماعي الاقتصادي للطالب بأثر متوسط مدرسته. الحل أن تدخل المتغير متمركزًا حول متوسط المجموعة مع متوسط المجموعة نفسه، فتحصل على أثر داخل المدارس وأثر بينها منفصلين، بصيغة مثل score ~ ses_within + ses_school_mean + school_type + (1 + ses_within | school). وإذا أضفت تفاعلًا بين ses_within وschool_type صار لديك تفاعل عبر المستويات (Cross-level Interaction) يجيب عن سؤال: هل تختلف قوة العلاقة باختلاف نوع المدرسة؟

الثالثة تحذير التقدير الشاذ (Singular Fit): يظهر حين يقدر تباين قريب من الصفر أو ارتباط قريب من ±1، ويعني غالبًا أن البيانات لا تدعم البنية العشوائية المطلوبة. يدور نقاش منهجي بين من يوصي بأوسع بنية عشوائية يسمح بها التصميم ومن يفضل البنية التي تدعمها البيانات؛ والمهم أن تبرر اختيارك ولا تتجاهل التحذير. وإذا كانت أساسيات النمذجة تحتاج مراجعة، فابدأ بمقال تحليل الانحدار.

النسخة البايزية بحزمة brms

تكتب حزمة brms النموذج بصيغة lme4 نفسها تقريبًا، ثم تقدره عبر Stan بخوارزمية مونت كارلو الهاميلتونية (HMC). الفرق المفاهيمي أن المعالم هنا متغيرات عشوائية لها توزيع قبلي (Prior) يعبر عما نعرفه قبل البيانات، وتحدثه البيانات إلى توزيع بعدي (Posterior). ومن التوزيع البعدي تستخرج ما تحتاجه: المتوسط، وفترة المصداقية (Credible Interval)، واحتمال أن يتجاوز الأثر قيمة لها معنى عملي.

library(brms)

# Weakly informative priors on the scale of the data (reaction time in ms).
# brms centres predictors, so "Intercept" is the mean response at mean Days.
priors <- c(
  prior(normal(300, 100), class = "Intercept"),
  prior(normal(0, 20), class = "b", coef = "Days"),
  prior(student_t(3, 0, 50), class = "sd"),
  prior(student_t(3, 0, 50), class = "sigma"),
  prior(lkj(2), class = "cor")
)

b2 <- brm(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy,
          prior = priors, chains = 4, iter = 2000, warmup = 1000,
          cores = 4, seed = 2026)

summary(b2)                  # estimates, 95% CrI, Rhat, Bulk_ESS, Tail_ESS
plot(b2)                     # trace and density plots per parameter
pp_check(b2, ndraws = 100)   # posterior predictive check
hypothesis(b2, "Days > 5")  # posterior probability that the slope exceeds 5 ms/day

فحوص التقارب قبل قراءة أي رقم

النتائج البايزية لا تقرأ قبل التأكد من أن سلاسل المحاكاة تقاربت. افحص ثلاثة أمور. الأول R-hat، الذي يقارن التباين بين السلاسل بالتباين داخلها؛ والتوصية الشائعة حاليًا ألا يتجاوز 1.01 لأي معلمة. الثاني حجم العينة الفعال (ESS) بنوعيه: Bulk_ESS لموثوقية تقدير المركز، وTail_ESS لموثوقية حدود الفترات؛ ويقترح كثيرون ألا يقل كل منهما عن 400 تقريبًا مع أربع سلاسل. الثالث رسائل الانتقالات المتباعدة (Divergent Transitions)، فهي تعني أن الخوارزمية لم تستكشف جزءًا من التوزيع البعدي جيدًا؛ والعلاج رفع control = list(adapt_delta = 0.95)، أو توزيعات قبلية أدق، أو إعادة صياغة النموذج. وتبدو رسوم الأثر (Trace Plots) السليمة سلاسل متداخلة بلا اتجاه.

أما التوزيعات القبلية فيجب أن تكون قابلة للدفاع: القيم أعلاه ضعيفة المعلومات، تستبعد الأرقام المستحيلة كأزمنة استجابة بالثواني الكاملة دون أن تفرض نتيجة. افحص ما تتنبأ به قبل رؤية البيانات بإعادة التقدير مع sample_prior = "only"، واعرض تحليل حساسية بتوزيعات بديلة. وللمقارنة بين نماذج بايزية استخدم التحقق المتقاطع التقريبي بترك مشاهدة واحدة عبر loo() وloo_compare().

تفسير فترة المصداقية 95% مباشر: بشرط النموذج والتوزيع القبلي، احتمال وقوع المعلمة داخل الفترة 95%. وهذا ما يظنه كثيرون خطأ عن فترة الثقة التكرارية، التي تعني أن 95% من الفترات المبنية بالطريقة نفسها عبر عينات متكررة ستحتوي القيمة الحقيقية.

بيانات البانل: التأثيرات الثابتة أم العشوائية؟

في الاقتصاد القياسي تسمى البيانات التي تتابع الوحدات نفسها عبر الزمن بيانات البانل (Panel Data)، والسؤال المركزي فيها: هل الخصائص غير المقيسة الثابتة لكل وحدة مرتبطة بالمتغيرات المستقلة؟ نموذج التأثيرات الثابتة (Fixed Effects) يزيلها بالتحويل الداخلي، أي طرح متوسط كل وحدة من قيمها، فيسمح بارتباطها بالمتغيرات، لكنه يعجز عن تقدير أثر أي متغير ثابت عبر الزمن. أما نموذج التأثيرات العشوائية (Random Effects)، الذي يقابل في منطقه نموذج المقطع العشوائي، فيستخدم التباين داخل الوحدات وبينها فيكون أكفأ، لكنه يفترض أن التأثير الفردي غير مرتبط بالمتغيرات المستقلة.

library(plm)
data("Grunfeld", package = "plm")    # 10 US firms observed 1935-1954

pg <- pdata.frame(Grunfeld, index = c("firm", "year"))

pool <- plm(inv ~ value + capital, data = pg, model = "pooling")
fe   <- plm(inv ~ value + capital, data = pg, model = "within")
re   <- plm(inv ~ value + capital, data = pg, model = "random")

pFtest(fe, pool)               # individual effects vs pooled OLS
plmtest(pool, type = "bp")     # Breusch-Pagan LM test for random effects
phtest(fe, re)                 # Hausman test, H0: RE is consistent

# Cluster-robust (by firm) standard errors for the FE model
summary(fe, vcov = vcovHC(fe, method = "arellano", cluster = "group", type = "HC1"))

اختبار هاوسمان (Hausman) يقارن تقديرات النموذجين: إذا اختلفت اختلافًا منهجيًا فالتأثيرات العشوائية غير متسقة، ويفضل النموذج الثابت. لكن صيغته التقليدية تفترض أن نموذج التأثيرات العشوائية كفء تمامًا، وهو افتراض يسقط مع عدم تجانس الأخطاء أو ارتباطها التسلسلي. البديل المتين نهج موندلاك (Mundlak): أضف متوسطات الوحدات للمتغيرات المتغيرة عبر الزمن إلى نموذج التأثيرات العشوائية، واختبر دلالتها المشتركة بأخطاء معيارية متينة؛ وتقدم phtest() صيغة انحدارية قريبة عبر method = "aux". ولاحظ أن Grunfeld فيها عشر شركات فقط، والأخطاء المعيارية المجمعة غير موثوقة مع عدد قليل من المجموعات، فالمثال هنا لعرض الصيغة لا للتوصية بها في هذا الحجم. وإذا كانت الصدمات الزمنية المشتركة مهمة فاستخدم effect = "twoways".

متى تستخدم كل نهج؟

• التأثيرات الثابتة: حين يكون هدفك أثرًا سببيًا لمتغير يتغير داخل الوحدات، وتخشى خصائص ثابتة غير مقيسة تربكه. هذا منطق قريب من الفروق في الفروق الذي يشرحه مقال تقييم أثر البرامج والسياسات.
• النموذج متعدد المستويات أو التأثيرات العشوائية: حين تهتم بمتغيرات على مستوى المجموعة، كنوع المدرسة أو خصائص المنطقة، أو بتباين الأثر بين المجموعات، أو بالتنبؤ لمجموعات جديدة. أضف متوسطات المجموعات لتفصل الأثر داخلها عن الأثر بينها.
• النسخة البايزية: حين يكون عدد المجموعات قليلًا، أو البنية العشوائية معقدة وتتكرر تحذيرات التقدير الشاذ، أو تحتاج عبارات احتمالية مباشرة لصانع القرار، أو لديك معرفة سابقة موثقة تريد إدخالها بشفافية.
• متغير تابع غير متصل: للنتائج الثنائية أو العددية استخدم glmer() في lme4 أو عائلة التوزيع المناسبة في brms، والمنطق نفسه ينطبق.

أخطاء شائعة

• تجاهل التداخل والاكتفاء بانحدار عادي على بيانات عنقودية.
• إسقاط النموذج متعدد المستويات لأن ICC صغير، رغم أن المجموعات الكبيرة تجعل أثر التصميم ملموسًا.
• مقارنة نماذج تختلف في التأثيرات الثابتة وهي مقدرة بطريقة REML.
• تجاهل تحذيرات التقدير الشاذ أو عدم التقارب.
• خلط الأثر داخل المجموعات بالأثر بينها دون تمركز مناسب.
• في brms: قبول التوزيعات القبلية الافتراضية دون فحص، أو تجاهل الانتقالات المتباعدة وقيم R-hat المرتفعة.
• الاختيار بين FE وRE بنتيجة اختبار هاوسمان وحده دون منطق نظري، أو تطبيق صيغته التقليدية مع أخطاء غير متجانسة.
• الثقة بأخطاء معيارية مجمعة مع عدد قليل جدًا من المجموعات.

أسئلة شائعة

كم مجموعة أحتاج لنموذج متعدد المستويات؟

لا يوجد حد متفق عليه، وتتداول الأدبيات قواعد تقريبية مختلفة. المؤكد أن تقدير تباين المستوى الأعلى يحتاج عددًا معقولًا من المجموعات، وأن متغيرات هذا المستوى تقدر بعدد المجموعات لا بعدد الأفراد. مع مجموعات قليلة تفيد النسخة البايزية بتوزيعات قبلية معقولة، أو تكفي التأثيرات الثابتة للمجموعات إذا لم يكن تباينها موضع اهتمام.

هل يجعل التوزيع القبلي النتائج ذاتية؟

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

ما الفرق بين التأثيرات العشوائية في plm والنموذج متعدد المستويات في lme4؟

المنطق واحد في حالة المقطع العشوائي: تأثير خاص بكل وحدة يعامل متغيرًا عشوائيًا. يختلفان في طريقة التقدير والتقاليد: plm يخدم أسئلة الاقتصاد القياسي حول الاتساق والاختيار بين FE وRE، وlme4 أوسع في البنى العشوائية، كالميول العشوائية والمستويات المتعددة والتصنيف المتقاطع.

هل تصلح هذه النماذج لمتغير تابع ثنائي؟

نعم، بالدالة glmer() مع family = binomial في lme4، أو بـfamily = bernoulli() في brms. ويحسب ICC حينها على المقياس الكامن بافتراض أن تباين المستوى الأول في النموذج اللوجستي يساوي π²/3.

اختيار البنية المناسبة قرار منهجي يسبق أي كود: ما مستويات التداخل، وما السؤال، وهل يهمك الأثر داخل المجموعات أم بينها. يعمل فريق مركز المستشار بلغتي R وPython على النماذج متعددة المستويات والبايزية ونماذج البانل للباحثين والمؤسسات؛ تفاصيل ذلك في صفحة خدمة تحليل البيانات المتقدم، والاستشارة الأولى مجانية ونرد على الاستفسارات خلال ساعة.