Theory-based Model Validation in the Generalized Multifactor Dimensionality Reduction Algorithm for Ordinal Phenotypes  

     Clinical studies indicate a close relationship between some diseases and the presence of specific interactions between genetic factors. As is the case in many studies, revealing genetic interactions that have a significant impact on the emergence of genetic diseases requires extensive statistical analyses. Because of the enormous volume of genetic data in the human race, it was necessary to develop statistical methods adapted to deal with high-dimensional data. Multifactor Dimensionality Reduction (MDR) is one of the leading nonparametric algorithms in this field. The algorithm reduces the dimensions of genetic data to obtain the most important interaction that has a direct impact on increasing the likelihood of genetic diseases appearing. In its composition, the algorithm relies on a set of nonparametric procedures to diagnose genetic interference with the highest impact exclusively on binary response variables. Like any statistical method, this algorithm is not devoid of weaknesses and application limitations, so the algorithm had to be developed to overcome the obstacles. One of the weaknesses of this algorithm is that the algorithm cannot handle data sets with ordinal response variable. Some researchers have developed a generalization of the multifactor dimensionality reduction algorithm to enable it to work with ordinal data. However, the generalized algorithm is more complex than the original algorithm. Therefore, we proposed developing the original algorithm in a simple way by employing ordinal logistic regression to classify individuals in the sample, while keeping all steps of the original algorithm unchanged. On the other hand, the MDR algorithm adopts a non-parametric method to verify the significance of the interferences nominated in the algorithm. This nonparametric procedure is based on the idea of permutational tests, and it consumes a very long time compared to parametric procedures that relies on theoretical approaches. Some researchers have suggested using the generalized extreme value distribution to verify the statistical significance of candidate interactions, but this method has only been used with continuous and binary dependent variables. In this research, the theoretical method based on the generalized extreme value distribution was employed instead of the permutational tests adopted in the algorithm when the response variable is of the ordinal type.

مقدمة

تعد الأنماط الظاهرية phenotypes انعكاساً للأنماط الجينية genotypes التي يحملها الانسان، مثل لون العينين، لون الشعر، شكل القدم، وغيرها من الصفات التي يمكن ملاحظتها بشكل ظاهري. ونعد الامراض الوراثية أحد الأنماط الظاهرية للتعابير الجينية التي يحملها الانسان، خصوصاً على مستوى التداخلات الجينية. حيث تعتبر القابلية على الإصابة ببعض الأمراض مرتبطة بشكلٍ كبير بالعوامل الجينية متعددة المواقع multilocus genetic factors على مستوى التأثيرات الرئيسة و/أو تأثيرات التداخلات (Pattin, White et al. 2009). تم استخدام العديد من الأساليب الإحصائية المعلمية لنمذجة العلاقة بين القابلية للإصابة بالمرض والعوامل الجينية. تنبني غالبية هذه الطرق على مفهوم النمذجة الخطية المعممة Generalized linear modeling (Gola, Mahachie John et al. 2016). ومع ذلك، نظراً للأبعاد العالية للبيانات الجينية و/أو حجم العينة الصغير نسبياً، حيث يُقدَّر عدد العوامل الجينية ثنائية الأليل في الجنس البشري بـ 84.7 مليون عامل جيني (Chauhan et. al, 2022) ، فقد لا تكون هذه الطرق فعالة في العمل في ظل هذه الظروف. ولرؤية ذلك، فإن مقدرات المربعات الصغرى الاعتيادية ordinary least squares (OLS) لمتجه معاملات الانحدار الخطي المتعدد  يمكن حسابه وفقًا للمعادلة رقم 1:

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

عندما نواجه مشكلة كون عدد المعلمات في أنموذج الانحدار أكبر من عدد المشاهدات ، فإنَّ رتبة المصفوفة تكون مساوية لـ  على الأكثر. مما يعني أنَّ هناك مشكلة تعدد العلاقة الخطية multicollinearity في البيانات. في مثل هذه الحالة، لا يوجد معكوس للمصفوفة ، وهذا يعني أنَّ طريقة OLS غير قابلة للتطبيق. حتى في حالة استخدام معكوس معمم generalized inverse للمصفوفة  لحساب   وفقاً للمعادلة 1، فإنه لن يكون هناك مقدر فريد بسبب كون أنَّ المعكوس المعمم غير فريد. علماً أنه لو تم استخدام معكوس معمم، فإنَّ تفسير معاملات الانحدار للمتغيرات التوضيحية المترابطة لن يكون دقيقاً (Sofroniou and Hutcheson 1999). في الدراسات الجينية، يمكن أن يصبح عدد المعلمات في الانموذج  كبيراً جداً، خصوصاً عندما يتم تضمين التداخلات من الرتب العليا بين العوامل الجينية في النمذجة. على سبيل المثال، فإنَّ الانموذج متعدد الحدود من الدرجة الثانية لبيانات فيها عشرة عوامل جينية ممكن أن تحتوي على 56 معلمة. فبالإضافة إلى معلمة المقطع، يكون لدينا المعلمات الخاصة بالتأثيرات الرئيسة والبالغ 10 معلمات، وجميع معلمات حدود التداخلات الثنائية في الأنموذج والبالغ 45 (توافيق 2 من 10). ويزداد عدد المعلمات في الأنموذج بشكل متسارع جداً في حال زيادة عدد العوامل المؤثرة في الأنموذج، وكذا الحال عند تضمين تداخلات من الدرجات العليا. علماً أن الدراسات الجينية غالبا ما تحتوي على عدد كبير من المتغيرات التوضيحية (العوامل الجينية)، ويكون التركيز فيها على التداخلات بين العوامل الجينية. عليه فإن ظهور خروقات في فروض التحليل للطرق المعلمية يكون وارداً جداً. لذلك تم تطوير بدائل لا معلمية للتغلب على صعوبات استخدام الطرق المعلمية.

خوارزمية تخفيض الأبعاد متعددة العوامل (Multifactor-Dimensionality Reduction (MDR

إحدى الخوارزميات اللامعلمية الرائدة في مجال تحليل ونمذجة البيانات الجينية هي خوارزمية تخفيض الأبعاد متعددة العوامل (Multifactor-Dimensionality Reduction (MDR، التي قدمها في الأصل ريتشي وآخرون (Ritchie, Hahn et al. 2001) ، حيث تم استخدامها على نطاق واسع لوصف العلاقة بين القابلية على الإصابة بالأمراض وتداخلات العوامل الجينية متعددة المواقع في دراسات المرضى والأصحاء case-control studies. علماً أن خوارزمية MDR تم تطويرها بالاعتماد على مفهوم طريقة التجزئة التوافقية combinatorial partitioning التي وصفها نيلسون وآخرون (Nelson, Kardia et al. 2001).

الهدف الرئيس لخوارزمية MDR هو تحديد التفاعل الجيني متعدد المواقع multilocus genetic interaction ذي المعنوية الأعلى إحصائياً مقارنةً مع باقي التداخلات. ويتم اختيار التداخل الأهم عن طريق تخفيض عدد المؤثرات الجينية إلى عامل واحد فقط، وذلك من خلال تصنيف كل التداخلات متعددة المواقع المتوفرة لدينا إلى تداخل ذي خطورة عالية أو خطورة منخفضة وفقاً لمعيار معين. كذلك يتم استخدام التحقيق التقاطعي cross-validation لتقييم صحة التداخل المقترح من أي درجة  لـ ، حيث أن  هو عدد العوامل الكلية في البيانات. علاوةً على ذلك، يتم التحقق من أهمية التفاعل النهائي المرشح باستخدام اختبار التباديل permutation testing والذي يعتمد على توليد التوزيع التجريبي لإحصاءة المختبر الاحصائي، ومن ثم استخدام التوزيع التجريبي للوقوف على مستوى المعنوية الإحصائية للأنموذج المقترح.

يتمثل أحد أوجه القصور الشائعة في خوارزمية MDR في أنها قابلة للتطبيق على البيانات المتزنة فقط، أي أن عدد المرضى والأصحاء يكون متساوياً في مجموعة البيانات، وذلك بسبب اعتماد الخوارزمية لعتبة مساوية لواحد. ضعف كبير آخر في خوارزمية MDR هو استخدام عتبة ثابتة لتصنيف الأفراد إلى مجموعات عالية المخاطر وأخرى منخفضة المخاطر. بغض النظر عن فوائد استخدام عتبة ثابتة، لأنها تقلل من العبء الحسابي، فقد يؤدي ذلك إلى خسارة هائلة في قوة الاختبار (Hua, Zhang et al. 2010). 

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

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

ويمكن تلخيص طريقة MDR في الخطوات الآتية:

1. تحديد العوامل الجينية و/أو العوامل البيئية في الدراسة.
2. استعراض التوزيع التكراري frequency distribution للبيانات في فضاء ذات بعد لكل تفاعل مُفترض من الدرجة ، . أي أن بيانات أي تفاعل ثنائي الدرجة يتم تمثيلها باستخدام جداول تكرارية ثنائية الأبعاد 2-way contingency table. وبالمثل، يتم استخدام جداول تكرارية ثلاثية الأبعاد 3-way contingency cube (أو ثلاثة جداول تكرارية ثنائية الاتجاه) لتمثيل البيانات لأي تفاعل من الدرجة الثالثة، وهكذا. يتم تحديد أبعاد هذه الجداول المتقاطعة بعدد المستويات في كل عامل. على سبيل المثال، يتم تمثيل التوزيع التكراري للتفاعل cross tabulations بين عاملين لكل منهما ثلاثة مستويات باستخدام جدول تكراري ذي ابعاد 3 × 3. تحتوي كل خلية في الجدول التكراري على أعداد المصابين وغير المصابين بالمرض، حيث تمثل كل خلية في الجدول أحد التداخلات متعددة المواقع بين العوامل الجينية المدروسة. يتم حساب نسبة المرضى إلى الأصحاء لكل خلية، ثم تقارن هذه النسب مع العتبة threshold المحددة مسبقاً لتحديد فيما إذا كانت الخلية (التداخل الجيني) ذات خطورة عالية أم منخفضة على ظهور المرض. حيث يتم تقرير فيما إذا كان الأفراد في كل خلية معرضين لخطر الإصابة بالمرض عندما تكون نسبة المرضى إلى الأصحاء تتجاوز أو تساوي العتبة المحددة. بخلاف ذلك، فإن الأفراد الحاملين للتداخل الجيني المتعدد المواقع يُصنفون على أنهم أقل عرضة للإصابة بالمرض عندما تكون نسبة المرضى إلى الأصحاء أقل من العتبة. علما أن العتبة المستخدمة في خوارزمية MDR هي واحد (Ritchie, Hahn et al. 2001). الهدف من عملية التصنيف هو تقليل أبعاد البيانات إلى متغير تنبؤ ثنائي واحد فقط.
3. يتم اختيار نموذج مقترح (تفاعل) للرتبة باعتباره الانموذج الذي يحتوي على أصغر خطأ تصنيف classification error (CE) لكل انموذج من الدرجة  ممكن بناؤه. للحصول على CE لكل انموذج، يتم تسجيل العدد الإجمالي للأفراد الذين تم تصنيفهم بشكل خاطئ (المرضى الذين تم تصنيفهم على أنهم منخفضو المخاطر، الأصحاء الذين تم تصنيفهم على أنهم مرتفعو المخاطر) لكل انموذج. عادةً ما يُطلق على التصنيف الخاطئ للمرضى بـ false negative error (FN)، في حين يُطلق على عدد الأصحاء المصنفين خطأً على أنهم مصابون بـ false positive error (FP).
4. بعد ذلك، ومن أجل تقييم صحة الأنموذج الذي تم اختياره في الخطوة الثالثة، يتم استخدام التحقيق التقاطعي بعشرة طويات (10-fold cross-validation CV) لكل تداخل من الدرجة . لتنفيذ إجراء CV، يتم تقسيم البيانات عشوائياً إلى عشر مجاميع متساوية الحجم تقريباً، بحيث يكون لكل مجموعة العدد نفسه المحدد من المصابين وغير المصابين لتحقيق حالة التوازن بين عدد المرضى والأصحاء في العينة. في كل مرة، يتم استبعاد مجموعة واحدة كمجموعة بيانات اختبار، بينما تعتبر المجموعات التسعة المتبقية مجموعة بيانات تدريب. في وقت لاحق، يتم إجراء تصنيف البيانات وإجراءات اختيار الأنموذج الموضحة في الخطوتين الثانية والثالثة على مجموعة بيانات التدريب. بعد ذلك، يتم تصنيف الأفراد الذين ينتمون إلى مجموعة بيانات الاختبار (المجموعة المستبعدة) وحسب انتماءاتهم الجينية وفقاً للتنبؤ الثنائي الذي تم الحصول عليه من تنفيذ الخطوتين الثانية والثالثة على مجموعة بيانات التدريب. يتم تنفيذ إجراء الاستبعاد على كل مجموعة من المجموعات العشر، ويتم حساب أخطاء التصنيف CEs لكل انموذج من الدرجة ممكن بناؤه باستخدام بيانات مجموعات التدريب. وبطريقة مماثلة لحساب CE، يتم حساب خطأ التنبؤ prediction error (PE) لكل مجموعة مستبعدة في الطيات العشرة. بمعنى آخر، فإن خطأ التنبؤ PE هو عدد الأفراد المصنفين بشكل خاطئ في بيانات مجموعة الاختبار. وللتخلص من التأثيرات المحتملة للتقسيم العشوائي للبيانات الى عشرة مجاميع، يتم تكرار جميع إجراءات الـ CV بالكامل عدة مرات (على سبيل المثال خمس مرات). حيث يتم إجراء تقسيم عشوائي جديد للبيانات الى عشر مجموعات متساوية الحجم في كل تكرار. ثم يتم حساب متوسط أخطاء التصنيف   لكل انموذج من الدرجة . بعد ذلك، يتم تحديد النماذج التي تصغر قيمة   لكل درجة من درجات التداخلات الممكنة . أخيراً، الأنموذج الذي يمثل بشكل أفضل العلاقة بين التفاعل الوراثي متعدد المواقع وقابلية الإصابة بالمرض بين جميع النماذج المُختارة هو الأنموذج الذي يحتوي على الحد الأدنى لمتوسط أخطاء التنبؤ ، حيث يتم حساب بطريقة مماثلة لـ  لكن باستخدام بيانات مجموعة الاختبار. يتم استخدام اتساق التحقيق التقاطعي Cross-validation consistency (CVC) لتقييم صلاحية الأنموذج المحدد. أي أن خوارزمية MDR تحسب عدد المرات التي يتم فيها تحديد كل أنموذج معين من جميع الطيات العشرة. يتم حساب متوسط ​​ اتساق التحقيق التقاطعي النهائي  لكل أنموذج مقترح بناءً على نتائج جميع التكرارات.
5. أخيرًا، للتحقق من معنوية التداخل المرشح، يتم استخدام اختبار التباديل مع 1000 مجموعة من البيانات المبدلة permuted data sets. حيث يتم خلط مشاهدات متغير الاستجابة (مصاب وغير مصاب) بشكل عشوائي بينما يتم الاحتفاظ بمعلومات المتغيرات التوضيحية كما هي. ثم يتم تطبيق خوارزمية MDR بالكامل على كل مجموعة من مجموعات البيانات المبدلة للحصول على انموذج جديد مع الـ  الخاص به. وبذلك سيكون لدينا 1000 قيمة من قيم  تستخدم لبناء توزيع العدم التجريبي لـ   الناتجة من البيانات الأصلية. لفحص الأهمية الإحصائية للأنموذج المقترح، تتم مقارنة  المحسوبة من مجموعة البيانات الأصلية بالتوزيع التجريبي لـ  الناتج من 1000 اختبار التباديل. حيث أن:

حيث أن  تمثل دالة الدليل indicator function. يُعدُّ الانموذج المقترح ذا دلالة إحصائية إذا كانت .

التعديلات المقترحة على خوارزمية MDR

نظراً لعدم تعامل خوارزمية MDR مع البيانات التي يكون فيها متغير الاستجابة من النوع الترتيبي، فقد تم اقتراح تعديل الخوارزمية لتتعامل مع هكذا نوع من المتغيرات وذلك من خلال تعديل الجزء الخاص بتصنيف خلايا الجداول التوافقية للتداخلات من الدرجات  إلى عوامل الخطورة. حيث أن التصنيف في خوارزمية MDR ينحصر بين (مصاب، غير مصاب). وتحسب بعدها نسبة المرضى إلى الاصحاء لتقرير فيما إذا كان التداخل الجيني ذا خطورة عالية أم لا، ومنها يتم حساب خطأي التصنيف والتقدير CE and PE. التعديل المقترح سيعتمد على توظيف الانحدار اللوجستي الترتيبي لتصنيف الخلايا إلى فئات تماثل فئات المتغير المعتمد. فلو فرضنا أن المتغير المعتمد له ثلاث فئات مختلفة وقابلة للترتيب، مثل مستويات الكولسترول في الدم (طبيعي، فوق الطبيعي، مرتفع)، فعندئذٍ يمكن استخدام الانحدار اللوجستي الترتيبي لتصنيف الخلايا وفقا لإحدى هذه الفئات الثلاثة. وبعد ذلك يتم حساب خطأ التصنيف وخطأ التقدير بنفس الطريقة السابقة. وتبقى آلية اختيار التداخلات ذات التأثير العالي على تطور المرض كما هي عليه في الخوارزمية الأصلية. وتكون صيغة الانحدار اللوجستي الترتيبي كالآتي:

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

والتداخلات الأليلية هي التداخلات التي تحدث بين عاملين جينيين ثنائيي الألائل. ويتم تمثيل بيانات التداخل من الدرجة الثانية أعلاه بمتغير وهمي وحسب التداخلات المتعلقة بتطور المرض. ثم يحسب خطأ التصنيف CE وخطأ التقدير PE ومعامل الاتساق CVC كما هو الحال عليه في خوارزمية تخفيض الأبعاد الأصلية.

وتم أيضاً اقتراح تعديل اسلوب التحقق من معنوية التداخلات المرشحة في خوارزمية MDR ليعتمد على التوزيعات النظرية بدلا من اختبار التباديل. حيث تكمن الفكرة في توظيف توزيع القيمة العظمى المعمم GEVD لحساب قيمة الـ  للاحصاءة المحسوبة التي هي . ويعود السبب في اختيار توزيع القيمة العظمى المعمم إلى أن النماذج المرشحة من خوارزمية MDR يتم اختيارها بحيث تعظم قيمة  . حيث يتم توليد توزيع العدم null distribution للإحصاءة عن طريق عدد محدود جداً من اختبارات التباديل (50 تباديل على الاكثر). والذي بدوره سيؤدي إلى تخفيض زمن التنفيذ بشكل كبير جدا مقارنة بأسلوب التباديل الذي يعتمد على توليد 1000 عينة تباديلية. ويتم تقدير معلمات توزيع القيمة العظمى المعمم من خلال العينة ذات الحجم 50، وتوظف لحساب قيمة الـ  للتداخل المرشح. حيث يتم اشتقاق التوزيع النظري بالاعتماد على قيم  المولدة من اختبار التباديل في خوارزمية MDR المعدلة. حيث يتم تقدير معلمات توزيع GEV وتطبيق الدالة التراكمية لحساب الـ  لمعامل الاتساق للتداخل المرشح، وكالآتي:

حيث أن  تمثل المتغير الذي يتبع توزيع GEV والذي يمثل معدل معامل الاتساق  بالنسبة لنا، وأن  تمثل القيمة المشاهدة لـ  والخاصة بالتداخل المرشح من خوارزمية MDR، وأن  يمثلون مقدرات الإمكان الأعظم لمعلمات توزيع GEV، حيث أن  تمثل معلمة الموقع،  تمثل معلمة القياس، وأن  تمثل معلمة الشكل. ويعدُّ التداخل الجيني المقترح من خوارزمية MDR معنويا إذا كانت  .

وبذلك يمكن سرد خطوات الخوارزمية المعدلة كالآتي:

1. تحديد العوامل الجينية و/أو العوامل البيئية في الدراسة.
2. استعراض التوزيع التكراري frequency distribution للبيانات في جداول توافقية ذات بعد لكل تفاعل مُفترض من الدرجة ، . تحتوي كل خلية في الجدول التوافقي على أعداد الأشخاص الحاملين لكل فئة من فئات متغير الاستجابة، حيث تمثل كل خلية في الجدول أحد التداخلات متعددة المواقع بين العوامل الجينية المدروسة. يتم استخدام الانحدار اللوجستي الترتيبي لإعادة تصنيف الخلايا ضمن كل تداخل من الدرجة ، وذلك لغرض حساب خطأ التصنيف CE لكل تداخل.
3. يتم اختيار أفضل انموذج (تداخل) من كل رتبة بوصفه الانموذج الذي يحتوي على أصغر خطأ تصنيف classification error (CE). للحصول على CE لكل انموذج، يتم تسجيل العدد الإجمالي للأفراد الذين تم تصنيفهم بشكل خاطئ لكل انموذج، حيث يعتبر التصنيف غير صحيح إذا تم تصنيف الأشخاص بشكل مخالف للفئة التي ينتمون اليها في الواقع. ونلاحظ أنه لا يوجد أخطاء من نوع False Positive وFalse Negative في حالتنا هذه، كون أن الفئات قد تمثل جميعها أشخاصا مصابين، مثل درجات الإصابة بأحد أنواع الأورام الخبيثة.
4. بعد ذلك، ومن أجل تقييم صحة الأنموذج الذي تم اختياره في الخطوة الثالثة، يتم استخدام التحقيق التقاطعي بعشرة طويات (10-fold cross-validation CV) لكل تداخل من الدرجة . لتنفيذ إجراء CV، يتم تقسيم البيانات عشوائياً إلى عشر مجاميع متساوية الحجم تقريباً. في كل مرة، يتم استبعاد مجموعة واحدة كمجموعة بيانات اختبار، بينما تعتبر المجموعات التسعة المتبقية مجموعة بيانات تدريب. في وقت لاحق، يتم إجراء تصنيف البيانات وإجراءات اختيار الأنموذج الموضحة في الخطوتين الثانية والثالثة على مجموعة بيانات التدريب. بعد ذلك، يتم تصنيف الأفراد الذين ينتمون إلى مجموعة بيانات الاختبار (المجموعة المستبعدة) وحسب انتماءاتهم الجينية وفقاً للتصنيف الذي تم الحصول عليه من تنفيذ الخطوتين الثانية والثالثة على مجموعة بيانات التدريب. يتم تنفيذ إجراء الاستبعاد على كل مجموعة من المجموعات العشر، ويتم حساب أخطاء التصنيف CEs لكل انموذج من الدرجة ممكن بناؤه باستخدام بيانات مجموعات التدريب. وبطريقة مماثلة لحساب CE، يتم حساب خطأ التنبؤ prediction error (PE) لكل مجموعة مستبعدة في الطيات العشرة. بمعنى آخر، فإن خطأ التنبؤ PE هو عدد الأفراد المصنفين بشكل خاطئ في بيانات مجموعة الاختبار. وللتخلص من التأثيرات المحتملة للتقسيم العشوائي للبيانات الى عشر مجاميع، يتم تكرار جميع إجراءات الـ CV بالكامل خمس مرات. حيث يتم إجراء تقسيم عشوائي جديد للبيانات إلى عشر مجموعات متساوية الحجم في كل تكرار. ثم يتم حساب متوسط أخطاء التصنيف   لكل انموذج من الدرجة . بعد ذلك، يتم تحديد النماذج التي تصغر قيمة   لكل درجة من درجات التداخلات الممكنة . أخيراً، الأنموذج الذي يمثل بشكل أفضل العلاقة بين التفاعل الوراثي متعدد المواقع وقابلية الإصابة بالمرض بين جميع النماذج المُختارة هو الأنموذج الذي يحتوي على الحد الأدنى لمتوسط أخطاء التنبؤ ، حيث يتم حساب بطريقة مماثلة لـ  لكن باستخدام بيانات مجموعة الاختبار. يتم استخدام اتساق التحقيق التقاطعي Cross-validation consistency (CVC) لتقييم صلاحية الأنموذج المحدد. أي أن خوارزمية MDR تحسب عدد المرات التي يتم فيها تحديد كل أنموذج معين من جميع الطيات العشرة. يتم حساب متوسط ​​اتساق التحقيق التقاطعي النهائي  لكل أنموذج مقترح بناءً على نتائج جميع التكرارات.
5. أخيرًا، للتحقق من معنوية الأنموذج المرشح، يتم استخدام 50 عينة من البيانات المبدلة permuted data sets. حيث يتم خلط مشاهدات متغير الاستجابة (فئات متغير الاستجابة) بشكل عشوائي بينما يتم الاحتفاظ بمعلومات المتغيرات التوضيحية كما هي. ثم يتم تطبيق الخطوات السابقة بالكامل على كل مجموعة من مجموعات البيانات المبدلة للحصول على انموذج جديد مع الـ  الخاص به. وبذلك سيكون لدينا 50 قيمة من قيم  تستخدم لتقدير معلمات توزيع العدم لإحصاءة  الناتجة من البيانات الأصلية وذلك وفقا لتوزيع GEV. وللوقوف على المعنوية الإحصائية للأنموذج المقترح، يتم حساب قيمة الـ p-value كالآتي:

ويعدُّ الانموذج المقترح ذا دلالة إحصائية إذا كانت .

دراسة محاكاة Simulation Study

لغرض اختبار فاعلية الخوارزمية وفقا للتعديلات والمقترحة في الجانب النظري في تشخيصها للتداخلات ذات التأثير العالي على الإصابة بالمرض، تم تطبيق الخوارزمية المعدلة على بيانات جينية مولدة وفقا لسيناريوهات مختلفة سيرد ذكرها لاحقاً. ولضيق الوقت وقلة الإمكانات الحاسوبية المتاحة، فقد تم الاقتصار في جميع المجموعات المولدة على توليد المعلومات الجينية لستة عوامل جينية فقط . وبذلك يكون عدد التداخلات من الدرجة الثانية الممكنة هو 15 تداخلاً. بالمقابل، يكون عدد تداخلات العوامل الجينية من الدرجة الثالثة هو 20 تداخلاً، وكما موضح في الجدول رقم 2:

الجدول رقم 2: التداخلات الجينية الممكنة في دراسة المحاكاة

وعُدَّ وجود أليلين فقط لكل عامل جيني، أليل سائد يرمز له بحرف الكبير، وأليل متنحي يرمز له بحرف صغير. وعُدَّ التوزيع النظري للألائل المولدة خاضعا لقانون هاردي- وينبيرغ وفقا للاحتمالات

وبذلك تكون الاحتمالية النظرية المعتمدة في التوليد لظهور التداخل الأليلي السائد (Major Allele Combination) هي:

بينما تكون احتمالية ظهور التداخل الأليلي المتنحي (Minor Allele Combination) هي:

حيث تم تطبيق الاحتمالات أعلاه على جميع العوامل الجينية المولدة.

وتم أولا توليد مجتمع بحجم مئة ألف مشاهدة، ثم قمنا بسحب عينة عشوائية واحدة فقط بحجم 500 مشاهدة. وتم تكرار هذا الإجراء عشرين مرة، أي أن مجموع العينات بحجم 500 التي تم سحبها عشوائياً هو 20 عينة، كل عينة مسحوبة من مجتمع مولد بشكل منفصل عن باقي المجتمعات. وتم توليد عشرة من هذه المجتمعات بحيث يكون التداخل الجيني المسبب لزيادة احتمالية تطور المرض من الدرجة الثانية، أي . فيما تم توليد المجتمعات العشرة المتبقية بحيث يكون التداخل الجيني المرتبط بتفاقم حالة المرض من الدرجة الثالثة، أي . وتم إعادة عملية التوليد اعلاه لتوليد عشرة مجموعات جديدة من البيانات بحجم 1000 مشاهدة مرتبطة بتداخل محدد من الدرجة الثانية، وعشر مجموعات أخرى بحجم 1000 مرتبطة بتداخل محدد من الدرجة الثالثة. وبالطريقة نفسها تم توليد 20 عينة بحجم 2000 مشاهدة. أي أن مجموع العينات المولدة بلغ 60 عينة في المجمل. علماً أن الدراسات الجينية تحتاج إلى عدد كبير نسبياً من المشاهدات لضمان ظهور النمط الحقيقي للأمراض الوراثية السارية في المجتمع.

ولربط مستويات الإصابة بالمرض (المتغير المعتمد) بأحد التداخلات الجينية، تم توليد متغير الاستجابة الترتيبي بثلاث فئات مرتبة (1,2,3) وفقاً لأنموذج انحدار لوجستي مفترض بحيث يكون فيه أحد التداخلات الجينية المدرجة في الجدول رقم 2 معنوياً بمعامل انحدار محدد مسبقاً، بينما تهمل بقية المتغيرات التوضيحية في الانموذج، أي اعتبار أن معاملات الانحدار لباقي المتغيرات والتداخلات تساوي صفراً. مثلاً في التداخل من الدرجة الثانية بين العاملين  نفترض أن  وباقي المعاملات أصفار. أو في التداخل من الدرجة الثالثة بين العوامل  نفترض أن  وباقي المعاملات أصفار، وهكذا. وتمت إضافة حد عشوائي يتوزع وفقاً للتوزيع الطبيعي بمتوسط وانحراف معياري محددين وبمعلمة انحدار تساوي واحد، وذلك لإضفاء السمة العشوائية على عملية توليد المتغير المعتمد. وتم أيضاً تحديد قيم لمعلمتي المقطع . وأخيراً، تم تحديد التداخلات متعددة المواقع multilocus genetic interaction المرتبطة في تطور المرض وحسب درجة التداخل. ولتوضيح ذلك، افرض أن التداخل المستهدف هو بين العوامل الثلاثة . وبما أن التداخل بين العوامل الثلاثة يمكن تمثيله بجدول تقاطعي ثلاثي الأبعاد (ثلاثة جداول توافقية كما موضح في الجدول رقم 3). فقد تم ربط التداخلات متعددة المواقع المؤشرة باللون الأحمر مع تطور المرض وذلك لضمان معنوية التداخل من الدرجة الثالثة بين العوامل الثلاثة.

بعد ذلك نطبق المعادلة الآتية:

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

وبعد إتمام توليد البيانات وسحب العينة المطلوبة، يتم تنفيذ جميع خطوات الخوارزمية المعدلة والموضحة في الفصل الثاني من هذه الرسالة. حيث يتم استخدام بيانات كل عينة لبناء انموذج انحدار لوجستي يشمل جميع التداخلات الممكنة والموضحة في الجدول رقم 2 أعلاه بغض النظر عن معنويتها. وتم بناء الانموذج باستخدام الدالة polr من مكتبة MASS في برنامج R (Venables & Ripley, 2002). ثم يتم إعادة تصنيف الأفراد في كل عينة باستخدام الانموذج المقدر ومنها يتم حساب خطأ التصنيف  وخطأ التقدير  واتساق التحقق المتقاطع . وأخيراً يتم التحقق من معنوية التداخل المرشح بإعادة تنفيذ الخوارزمية على بيانات تبادلية بعدد 50 لكل عينة، ومنها نحسب  لكل تداخل مرشح من كل عينة تباديلية. حيث تم حساب مقدرات الإمكان الأعظم لمعلمات توزيع القيمة العظمى المعمم GEVD للإحصاءة  من كل عينة من خلال العينة التباديلية وذلك باستخدام الدالة egevd من المكتبة EnvStats في برنامج R (Millard SP, 2013). ومن التوزيع المقدر يتم حساب الـ  للتداخل المرشح وفقا للمعادلة رقم 4 وذلك من خلال الدالة pgevd من المكتبة EnvStats في برنامج R.

يذكر أنه تم إجراء جميع إجراءات المحاكاة الموضحة في أعلاه باستخدام برنامج  (R Core Team (2023)) في بيئة نظام Windows 11 باستخدام حاسبة مدعمة بمعالج Intel Core i7-8550U.

نتائج المحاكاة

عندما  والتداخل من الدرجة الثانية

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث أن  و  تمثلان معلمتي المقطع لفئتي المتغير المعتمد الأولى والثانية، وأن  تمثل معامل الانحدار للمتغير الذي يمثل التداخل بين العاملين الجينيين ، وأن  تمثل معامل الحد العشوائي الذي تمت اضافته لضمان عشوائية عملية التوليد.

وتم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 4 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 500 والتداخل المسبب للمرض من الدرجة الثانية.

الجدول رقم 4: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 500 والتداخل المسبب للمرض من الدرجة الثانية

نلاحظ أن الخوارزمية المعدلة تمكنت من تشخيص التداخل الحقيقي من الدرجة الثانية في 8 عينات وفشلت في حالتين فقط. كما نلاحظ أنه بالرغم من فشل الخوارزمية في تحديد التداخل الحقيقي في الحالتين السادسة والتاسعة، مع ذلك فإن الخوارزمية تمكنت من تحديد أحد العوامل الجينية الداخلة في تكوين التداخل المعتمد في توليد البينات. حيث تم تشخيص العامل  في العينة السادسة، وكذلك العامل  في العينة التاسع. ويلاحظ أيضاً أن التداخلين المذكورين غير معنويين، ما يدل على نجاعة الأسلوب المستخدم في التحقق من معنوية التداخلات المرشحة. ونلاحظ أيضا استقرار عمل الخوارزمية من ناحية الزمن المستهلك في التنفيذ، وهو وقت قصير جداً مقارنة فيما لو تم اعتماد اسلوب الاختبارات التباديلية التقليدي والذي يعتمد على توليد التوزيع التجريبي لـ  باستخدام 1000 عينة تباديلية.

عندما  والتداخل من الدرجة الثالثة

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث أن  و  تمثلان معلمتي المقطع لفئتي المتغير المعتمد الأولى والثانية، وأن  تمثل معامل الانحدار للمتغير الذي يمثل التداخل بين العوامل الجينية ، وأن  تمثل معامل الحد العشوائي الذي تمت اضافته لضمان عشوائية عملية التوليد عملية التوليد.

وتم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 5 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 500 والتداخل المسبب للمرض من الدرجة الثالثة.

الجدول رقم 5: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 500 والتداخل المسبب للمرض من الدرجة الثالثة

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

ويلاحظ بشكل عام أن قيم  منخفضة نسبياً في الحالات التي تكون فيها حجم العينة 500 بغض النظر عن معنويتها، خصوصاً للنماذج من الدرجة الثالثة. وذلك يعود إلى صغر حجم العينة مما يقلل من فرصة تحديد التداخل الحقيقي. يذكر أن القيمة العظمى الممكنة للـ  هي 10 وذلك لأننا نستخدم عشر طويات في عملية التحقق المتقاطع.

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

عندما  والتداخل من الدرجة الثانية

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث تم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 6 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 1000 والتداخل المسبب للمرض من الدرجة الثانية.

الجدول رقم 6: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 1000 والتداخل المسبب للمرض من الدرجة الثانية

عندما  والتداخل من الدرجة الثالثة

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث تم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 7 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 1000 والتداخل المسبب للمرض من الدرجة الثالثة.

الجدول رقم 7: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 1000 والتداخل المسبب للمرض من الدرجة الثالثة

يلاحظ من النتائج في الجدول رقم 6 ورقم 7 عند مقارنتها بالنتائج الموضحة في الجدولين رقم 4 و5 تحسن أداء الخوارزمية عندما ازداد حجم العينة إلى 1000. حيث تمكنت الخوارزمية من تحديد التداخلات الحقيقية من الدرجة الثانية، وثمانية تداخلات حقيقية من الدرجة الثالثة. كما يلاحظ تحسن ملحوظ في قيم الـ  في كلتا الحالتين، حيث اقتربت القيم من 10 أكثر. وفيما يخص معنوية التداخلات المقترحة، فإن جميع التداخلات أظهرت معنوية إحصائية باستثناء التداخلين في العينتين الرابعة والتاسعة للتداخلات من الدرجة الثالثة، وذلك لكون أن الخوارزمية فشلت في تحديد التداخل الصحيح في هاتين العينتين. وأخيراً فإن زمن التنفيذ قد ارتفع بشكل طفيف وذلك لازدياد العبء الحسابي بسبب ازدياد حجم العينة.

عندما  والتداخل من الدرجة الثانية

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث تم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 8 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 2000 والتداخل المسبب للمرض من الدرجة الثانية.

الجدول رقم 8: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 2000 والتداخل المسبب للمرض من الدرجة الثانية

عندما  والتداخل من الدرجة الثالثة

تم استخدام متجه المعلمات الآتي لتوليد مشاهدات المتغير المعتمد في جميع العينات العشرة:

حيث تم توليد الحد العشوائي من التوزيع الطبيعي بمتوسط صفري وانحراف معياري يساوي 3. ويوضح الجدول رقم 9 أدناه نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 2000 والتداخل المسبب للمرض من الدرجة الثالثة.

الجدول رقم 9: نتائج المحاكاة في الحالات التي يكون فيها حجم العينة 2000 والتداخل المسبب للمرض من الدرجة الثالثة

يلاحظ من النتائج في الجداول رقم 8 ورقم 9 عند مقارنتها بالنتائج الموضحة في الجداول رقم 4، 5، 6 و7 تحسن أداء الخوارزمية عندما ازداد حجم العينة الى 2000. حيث تمكنت الخوارزمية من تحديد التداخلات الحقيقية من الدرجة الثانية، والتداخلات الحقيقية من الدرجة الثالثة في جميع العينات المولدة. كما يلاحظ تحسن ملحوظ وكبير في قيم الـ  في كلتا الحالتين، حيث اقتربت القيم من 10 أكثر من السابق. وفيما يخص معنوية التداخلات المقترحة، فإن جميع التداخلات من الدرجتين الثانية والثالثة أظهرت معنوية إحصائية. وأخيراً فإن زمن التنفيذ قد ارتفع بشكل طفيف وذلك لازدياد العبء الحسابي بسبب ازدياد حجم العينة.

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