משטחי תגובה

מחשבון

מדריך

מה זה משטח תגובה, ומתי משתמשים בו

מתודולוגיית משטחי התגובה (Response Surface Methodology, RSM) היא "אוסף של טכניקות מתמטיות וסטטיסטיות" לבניית מודל של תגובה שתלויה בכמה משתנים, כש"המטרה היא לייעל את התגובה" (Montgomery, 2020, p. 408). התגובה הממוצעת היא פונקציה של הגורמים, η = f(x₁, x₂), והגרף שלה הוא "משטח התגובה"; מציגים אותו במפת שווה-תגובה – "קווים של תגובה קבועה" במישור של שני גורמים (עמ' 408-409).

ההבדל מניתוח שונות רגיל: שם הגורמים קטגוריאליים, ושואלים אם הממוצעים שונים. כאן הגורמים כמותיים ורציפים – אורך שיעור, מספר תרגילים, מינון, טמפרטורה – ושואלים איפה בתחום שלהם התגובה הכי טובה, ובאיזה כיוון כדאי להתקדם. את f האמיתית אין יודעים, ולכן מקרבים אותה בפולינום: מסדר ראשון (משוואה 11.1) או מסדר שני (11.2; עמ' 409). זה קירוב מקומי בלבד – "לא סביר שפולינום יהיה קירוב סביר על פני כל המרחב", אבל "באזור קטן יחסית הוא עובד בדרך כלל היטב" (עמ' 409).

הדוגמאות של Montgomery באות מהנדסה כימית – זמן תגובה וטמפרטורה מול תפוקה. אותה לוגיקה מתאימה לכל מחקר שבו החוקרת קובעת רמות של משתנים רציפים: אורך ותדירות של התערבות, מינון, מספר פריטים, זמן חשיפה. קובץ הדוגמה בעמוד – נתונים בדויים על שיעור תרגול – בנוי כך.

עבודה בשלבים

"RSM היא הליך רציף" (Montgomery, 2020, p. 410). כשנמצאים רחוק מהמיטב, למשטח כמעט אין עקמומיות ומודל מסדר ראשון מספיק; מטפסים איתו לעבר המיטב. קרוב למיטב העקמומיות מופיעה, ומתאימים מודל מסדר שני כדי לאתר אותו – כמו "טיפוס על גבעה" (עמ' 410). המחשבון בנוי לשני השלבים: מודל ליניארי עם מסלול העלייה התלולה ביותר, ומודל מסדר שני עם נקודה סטציונרית.

בשלב הראשון מספיק ניסוי 2^k (או שבר שלו) עם כמה נקודות מרכז. נקודות המרכז מאפשרות שלושה דברים: לאמוד את הטעות, לבדוק אינטראקציות, ולבדוק עקמומיות (עמ' 412). בדיקת העקמומיות היא בדרגת חופש אחת: אם המשטח מישורי, ממוצע הנקודות הפקטוריאליות ȳ_F שווה לממוצע נקודות המרכז ȳ_C; ההפרש אומד את Σβ_jj, ו-SS = n_Fn_C(ȳ_F − ȳ_C)²/(n_F + n_C) נבחן מול הטעות הטהורה (עמ' 413). העמוד מחשב אותה כשהנתונים הם בדיוק 2^k עם נקודות מרכז, ומפרק את שאריות המודל הליניארי לאינטראקציות, לריבועי הטהור ולטעות הטהורה – כמו לוחות 11.2 ו-11.5 (עמ' 413, 415).

מסלול העלייה התלולה ביותר

"שיטת העלייה התלולה ביותר היא הליך לתנועה רציפה בכיוון העלייה המרבית בתגובה"; למזעור – הירידה התלולה ביותר (Montgomery, 2020, p. 411). הכיוון מאונך לקווי שווה-התגובה של המודל הליניארי, ו"הצעדים לאורך המסלול פרופורציוניים למקדמי הרגרסיה" (עמ' 411). האלגוריתם (עמ' 415): בוחרים גודל צעד Δx_j בגורם אחד – "בדרך כלל את הגורם שאנחנו יודעים עליו הכי הרבה, או את זה עם המקדם הגדול ביותר בערכו המוחלט"; בכל גורם אחר Δx_i = b_i/(b_j/Δx_j); וממירים ליחידות הטבעיות. בדוגמה 11.1 הצעד הוא 5 דקות בזמן ו-2 מעלות בטמפרטורה (עמ' 416).

"את גודל הצעד בפועל קובע הנסיין, לפי הידע על התהליך" (עמ' 411). מריצים ניסויים לאורך המסלול "עד שלא נראית עוד עלייה בתגובה", ושם מתאימים מודל ליניארי חדש; ההגעה לסביבת המיטב "מסומנת בדרך כלל על ידי חוסר התאמה של מודל מסדר ראשון" (עמ' 411). בדוגמה 11.1 התגובה עלתה עד הצעד העשירי וירדה אחריו (עמ' 413-414), והמודל החדש סביב הצעד העשירי הראה עקמומיות מובהקת – "המודל מסדר ראשון אינו קירוב הולם" (עמ' 414). בעמוד – בוחרים גורם בסיס, גודל צעד ומספר צעדים, והלוח נותן את המסלול ביחידות מקודדות ובטבעיות, עם תחזית רק בתוך טווח הנתונים.

המודל מסדר שני: נקודה סטציונרית וניתוח קנוני

בכתיב מטריצות, ŷ = b₀ + x'b + x'Bx (משוואה 11.5): b – וקטור המקדמים הליניאריים, ו-B – מטריצה סימטרית עם המקדמים הריבועיים באלכסון וחצי מקדמי האינטראקציה מחוצה לו (Montgomery, 2020, pp. 417-418). הנקודה הסטציונרית – שבה כל הנגזרות מתאפסות – היא x_s = −½B⁻¹b, והתגובה בה ŷ_s = b₀ + ½x_s'b (משוואות 11.7-11.8, עמ' 418). היא "יכולה להיות נקודת מקסימום, מינימום או אוכף" (עמ' 416).

כדי לדעת מה היא, עוברים לצורה הקנונית: מזיזים את הראשית לנקודה הסטציונרית ומסובבים את הצירים לצירים הראשיים של המשטח, וכך ŷ = ŷ_s + λ₁w₁² + … + λ_kw_k² (משוואה 11.9), כש-λ הם הערכים העצמיים של B (עמ' 418). כולם שליליים – מקסימום; כולם חיוביים – מינימום; בסימנים מעורבים – אוכף; ו"המשטח תלול ביותר בכיוון w_i שבו |λ_i| הגדול ביותר" (עמ' 418). הצירים החדשים הם w = M'(x − x_s), כש-M – הווקטורים העצמיים המנורמלים (משוואה 11.10, עמ' 423). בדוגמה 11.2 של Montgomery: x_s = (0.389, 0.306), כלומר 87 דקות ו-176.5 מעלות, ŷ_s = 80.21, ושני הערכים העצמיים שליליים (−0.963, −1.414) – מקסימום (עמ' 420); העמוד משחזר את כל המספרים האלה.

רכסים. כשאחד הערכים העצמיים קרוב לאפס, התגובה "כמעט אינה רגישה" לתנועה לאורך הציר שלו (עמ' 424). אם הנקודה הסטציונרית בתוך אזור הניסוי – זה "רכס סטציונרי": קו של נקודות כמעט שוות, וגמישות לבחור ביניהן לפי עלות או נוחות. אם היא רחוקה מחוץ לאזור – רכס עולה (או יורד): "אי-אפשר להסיק על המשטח האמיתי או על הנקודה הסטציונרית", אבל כדאי לחקור הלאה בכיוון הציר השטוח (עמ' 424-425). העמוד מסמן ערך עצמי קטן מעשירית הגדול בערכו המוחלט – הסף של canonical בחבילת rsm (מוסכמת האתר) – ומודיע כשהנקודה מחוץ לטווח הנתונים.

כמה תגובות. כשיש כמה תגובות, מתאימים מודל לכל אחת ומחפשים צירוף שמספק את כולן – בחפיפה של מפות שווה-תגובה, באופטימיזציה עם אילוצים, או בפונקציית רצייה (desirability) שמצרפת את כל התגובות למדד אחד (עמ' 425-429). העמוד מנתח תגובה אחת בכל הרצה; לכמה תגובות – מריצים לכל אחת ומשווים את המפות.

איזה מערך?

למודל מסדר ראשון – 2^k או שבר שבו האפקטים הראשיים אינם מכונים זה בזה; נקודות מרכז אינן משנות את המקדמים הליניאריים ואת האורתוגונליות, ומוסיפות אומדן לטעות (Montgomery, 2020, p. 430). בעמוד "תכנון ניסויים" – 2^k ושברים.

מערך מרכזי מורכב (CCD) – "המחלקה הפופולרית ביותר" למודל מסדר שני (עמ' 430): 2^k (או שבר ברזולוציה V) עם n_F הרצות, 2k נקודות ציריות ב-±α על כל ציר, ו-n_C נקודות מרכז (עמ' 431). הוא מתאים לבנייה רציפה: מריצים קודם את 2^k עם נקודות המרכז, ואם יש חוסר התאמה – מוסיפים את הנקודות הציריות (עמ' 431). שתי החלטות: α ומספר נקודות המרכז.

• סיבובי – α = n_F^1/4 (1.414 ל-k = 2, 1.682 ל-k = 3): שונות החיזוי שווה בכל הנקודות שבאותו מרחק מהמרכז – "שונות התגובה החזויה קבועה על פני כדורים" (עמ' 431-432).
• כדורי – α = √k: כל הנקודות על כדור אחד; לאזור כדורי זו הבחירה הטובה מבחינת שונות החיזוי, ו"אין חשיבות לסיבוביות מדויקת כדי שהמערך יהיה טוב" (עמ' 432). מומלצות שלוש עד חמש נקודות מרכז (עמ' 432).
• ממורכז-פאות – α = 1: שלוש רמות בלבד, לאזור קובייתי, כשאי-אפשר לצאת מהטווח; אינו סיבובי, ו"שתיים או שלוש" נקודות מרכז מספיקות (עמ' 433).
• חסימה אורתוגונלית – כשהחלק הפקטוריאלי והצירי רצים בזמנים שונים, הם שני בלוקים. המערך "נחסם אורתוגונלית" כשאפקט הבלוקים אינו משפיע על המקדמים (עמ' 437); ב-CCD זה קורה כש-α = [n_F(n_A + n_CA)/(2(n_F + n_CF))]^1/2 (משוואה 11.18, עמ' 439). אפשר "לוותר מעט על הסיבוביות או על החסימה האורתוגונלית בלי אובדן מידע גדול" (עמ' 439). רק נקודות מרכז שרצו באותו בלוק הן חזרות (עמ' 440).

Box-Behnken – מערך בשלוש רמות, "צירוף של 2^k עם מערכי בלוקים לא-שלמים": כל זוג גורמים ב-±1 והשאר במרכז, ועוד נקודות מרכז. הוא "יעיל מאוד במספר ההרצות", "סיבובי או כמעט סיבובי", כל נקודותיו על כדור ברדיוס √2, ואין בו פינות – יתרון כשהצירוף הקיצוני של כל הגורמים יקר או בלתי אפשרי (עמ' 432; לוח 11.8 ל-k = 3, עמ' 433). בעמוד – k = 3 עד 5, באותה בנייה.

בתכנון העמוד מציג את שונות החיזוי המתוקננת N·x'(X'X)⁻¹x (משוואה 11.16, עמ' 436) במרכז, בנקודה פקטוריאלית ובנקודה צירית – במערך סיבובי השתיים האחרונות שוות כשהן באותו מרחק. מערכים אופטימליים (D, I), מערכי תערובות ומערכים לניסויי מחשב (עמ' 440-472) – אינם בעמוד.

האם המודל מתאים?

חוסר התאמה מול טעות טהורה. שאריות המודל מתפרקות לשניים: טעות טהורה – השונות בין הרצות חוזרות באותה נקודה, שאינה תלויה במודל; וחוסר התאמה – מה שנשאר. אם חוסר ההתאמה גדול מהטעות הטהורה, המודל מחמיץ משהו (Montgomery, 2020, Tables 11.2, 11.5 and 11.7). לכן נקודות המרכז חשובות: בלי חזרות אין טעות טהורה ואין המבחן הזה.

השוואת מודלים. לוח 1 בפלט משווה את שלושת המודלים – ליניארי, עם אינטראקציות, ומסדר שני – כמו "Model Summary Statistics" ו"Lack of Fit Tests" אצל Montgomery (לוח 11.7, עמ' 421): סטיית התקן של השאריות, R², R² מתוקן, R² לחיזוי ו-PRESS, וחוסר ההתאמה. בדוגמה 11.2 הליניארי וזה עם האינטראקציות נכשלים (R² לחיזוי שלילי, חוסר התאמה מובהק), והמודל מסדר שני מתאים (R² לחיזוי .918, חוסר התאמה לא מובהק, p = .290). R² לחיזוי נשען על PRESS – הסכום של ריבועי השאריות כשכל הרצה נחזית בלי עצמה.

אבחון. כמו בכל רגרסיה – שאריות, מנוף והשפעה (לוח 11.7, עמ' 422). בעמוד – לוח האבחון המקופל בסוף הפלט.

איך קוראים את הפלט?

בתכנון: לוח 1 – המערך ותכונותיו: החלקים, α, מספר ההרצות והמקדמים, דרגות החופש לחוסר התאמה ולטעות טהורה, החסימה, ושונות החיזוי. לוח 2 – דף ההרצה בסדר אקראי, עם הערכים המקודדים והטבעיים. תרשים 1 – נקודות המערך. ההורדה לאקסל כוללת מקרא עם הקידוד, והעמוד קורא אותו בחזרה בניתוח.

בניתוח: לפני ההרצה – לוח הקידוד: לכל גורם מרכז וחצי-טווח (x = (ξ − מרכז)/חצי-טווח), מהמקרא או לפי הערכים. לוח 1 – השוואת המודלים. לוח 2 – ניתוח השונות של המודל שנבחר: כל איבר עם סכום ריבועים חלקי, ושאריות שמתפרקות לחוסר התאמה וטעות טהורה. לוח 3 – המקדמים ביחידות מקודדות, עם רווחי סמך ו-VIF, ומתחתיו המשוואה ביחידות המקור. במודל מסדר שני – לוח 4, הנקודה הסטציונרית, ולוח 5, הניתוח הקנוני; במודל ליניארי – מסלול העלייה התלולה ביותר. תרשים – מפת שווה-תגובה לכל זוג גורמים, כשהשאר במרכז או בנקודה הסטציונרית. אחריו – חיזוי בנקודה עם רווח סמך ורווח ניבוי, לוח אבחון, פסקת ממצאים ותיבת "קריאת התוצאה".

הנחות ואזהרות

הקירוב מקומי. הפולינום מתאים לאזור שנבדק, לא מעבר לו (Montgomery, 2020, p. 409). כל תחזית מחוץ לנקודות המערך היא אקסטרפולציה; העמוד לא מחשב תחזיות במסלול מחוץ לטווח, ומזהיר כשהנקודה הסטציונרית מחוץ לו. ניתוח הנקודה הסטציונרית מניח שהיא "בתוך אזור החקירה" (עמ' 418).

אקראיות. כמו בכל ניסוי, סדר ההרצה אקראי – דף ההרצה שהעמוד בונה. ההנחות על הטעות הן של רגרסיה רגילה: בלתי תלויות, שונות שווה, התפלגות קרובה לנורמלית – בודקים בשאריות.

קידוד. המקדמים, הנקודה הסטציונרית והניתוח הקנוני – ביחידות מקודדות, כמו אצל Montgomery. המשוואה ביחידות המקור שקולה, אבל המקדמים שלה תלויים ביחידות ואין להשוות ביניהם (עמ' 422). הערכים העצמיים תלויים בקידוד; קידוד שבו ±1 הם הרמות הפקטוריאליות שומר על המשמעות שלהם.

מה העמוד אינו מחשב. מערכים אופטימליים (D, I, G), מערכי תערובות, ניסויי מחשב, EVOP (Montgomery, 2020, Sections 11.4.4-11.7), פונקציית רצייה לכמה תגובות, ניסויי 3^k (פרק 9) ומערכים חסינים (פרק 12). בעמוד "תכנון ניסויים" – גורמים קטגוריאליים, 2^k ושברים; ובעמוד "גורמים אקראיים, מקוננים וספליט-פלוט" – מערכים עם הגבלה על האקראיות.

ואיך עושים זאת ב-SPSS?

ב-SPSS Statistics אין פרוצדורה למשטחי תגובה. מחשבים את המשתנים המקודדים, את מכפלותיהם ואת ריבועיהם ב-Transform › Compute Variable, ומריצים Analyze › General Linear Model › Univariate עם כולם כ-Covariates, עם סכומי ריבועים מסוג III – הם ה-SS החלקיים של לוח 2 – ועם Lack of fit ב-Options. הסינטקס שהעמוד מייצר עושה בדיוק את זה (COMPUTE ואז UNIANOVA … WITH … /PRINT=PARAMETER LOF). את הנקודה הסטציונרית, הניתוח הקנוני והמפות SPSS אינו מחשב – בשבילם R או העמוד. בלי בלוקים, הסינטקס מוסיף גם REGRESSION עם אותם איברים – אותם מקדמים. הנתיבים הם מידע כללי על התוכנה, והסינטקס לא הורץ ב-SPSS. ב-PSPP, התוכנה החופשית שקוראת סינטקס של SPSS, אין UNIANOVA; את החלק של REGRESSION הרצנו בו על קובץ הדוגמה, והמקדמים זהים לעמוד.

ואיך עושים זאת ב-R?

חבילת rsm: coded.data לקידוד, rsm עם FO, TWI, PQ או SO למודל, ו-summary נותן את המקדמים ואת ניתוח השונות לפי קבוצות עם חוסר התאמה וטעות טהורה; canonical – הנקודה הסטציונרית והערכים העצמיים; contour – המפות; ccd ו-bbd – המערכים. canonical מאפס ערכים עצמיים קטנים מעשירית הגדול לפני חישוב הנקודה, ולכן במשטח עם רכס הנקודה שלו שונה מ-−½B⁻¹b של Montgomery שבעמוד. מסלול העלייה של steepest נמדד במרחק אחיד ולא לפי צעד בגורם אחד. הסקריפט הורץ כלשונו על קובץ הדוגמה האנגלי, והערכים בהערות זהים לעמוד. כלל האתר לעברית ב-R: קובץ UTF-8, כל שורת הערה מתחילה במילה עברית, ושמות משתנים באנגלית.

# משטחי תגובה – פרופסור סטאט. קובץ זה נשמר ב-UTF-8; כל שורת הערה מתחילה במילה עברית
library(readxl)   # קריאת אקסל
library(rsm)      # משטחי תגובה (install.packages("rsm"))

# טעינה: התבנית התקנית – שורה = הרצה; הגיליון "Data"
d <- as.data.frame(read_excel("C:/Rdata/DATA_rsm_example_eng.xlsx", sheet = "Data"))

# קידוד: x = (ξ − מרכז)/חצי-טווח – הרמות הפקטוריאליות הופכות ל-1− ול-1+
cd <- coded.data(d, x1 ~ (lesson_length - 45) / 10, x2 ~ (practice_items - 10) / 3, x3 ~ (review_weeks - 3) / 1)

# מודל מסדר שני: הטבלה של rsm – FO, TWI, PQ, שאריות, חוסר התאמה וטעות טהורה
fit <- rsm(retention ~ SO(x1, x2, x3), data = cd)
summary(fit)       # מקדמים: b0 71.74; x1 4.26, x2 −2.68, x3 1.69; x1:x2 1.10, x1:x3 0.23, x2:x3 −1.68; x1² −3.26, x2² −2.57, x3² −1.68; חוסר התאמה: F(5, 5) = 0.89, p = .551

# SS חלקי לכל איבר (לוח 2 בעמוד): עמודות מפורשות, ואז drop1 – כל איבר אחרי כל האחרים
e <- transform(as.data.frame(cd), x12 = x1 * x2, x13 = x1 * x3, x23 = x2 * x3, x1q = x1^2, x2q = x2^2, x3q = x3^2)
drop1(lm(retention ~ x1 + x2 + x3 + x12 + x13 + x23 + x1q + x2q + x3q, data = e), test = "F")   # x1 F = 123.97 … x3q F = 21.07

# הנקודה הסטציונרית והערכים העצמיים (לוח 4 בעמוד)
canonical(fit)     # xs = 0.568, −0.687, 0.883; ביחידות המקור 50.7 דקות, 7.9 תרגילים, 3.9 שבועות; ערכים עצמיים −1.16, −2.64, −3.71 – מקסימום

# מפת שווה-תגובה לשני הגורמים הראשונים, כשהשלישי בנקודה הסטציונרית
contour(fit, ~ x1 + x2, image = TRUE, at = canonical(fit)$xs)

# השוואת המודלים (לוח 1 בעמוד): R² מתוקן ו-PRESS
for (f in list(retention ~ FO(x1, x2, x3), retention ~ FO(x1, x2, x3) + TWI(x1, x2, x3), retention ~ SO(x1, x2, x3))) {
  m <- rsm(f, data = cd)
  cat(sprintf("%-45s adjR2 = %.3f  PRESS = %.2f\n", deparse(f), summary(m)$adj.r.squared, sum((resid(m) / (1 - hatvalues(m)))^2)))
}   # R² מתוקן .482 / .432 / .944; PRESS 423.72 / 625.07 / 88.27

דיווח התוצאות

הדיווח צריך לאפשר לשחזר את המערך ואת הניתוח: סוג המערך (CCD, ו-α שלו; Box-Behnken), הגורמים וטווחיהם, הקידוד, מספר נקודות המרכז והבלוקים, ומספר ההרצות; המודל שנבחר ולמה (חוסר ההתאמה, R² מתוקן ו-R² לחיזוי); לוח מקדמים; ואם יש – הנקודה הסטציונרית ביחידות המקור, התגובה החזויה בה ואופייה (הערכים העצמיים). APA (2020, סעיף 6.43) מבקש "די מידע כדי שהקורא יוכל לאשש את הניתוחים". כללי הצורה (סעיף 6.36): R² ו-p בלי אפס לפני הנקודה; מקדמים, ממוצעים וערכים עצמיים – עם. פסקת הממצאים בפלט בנויה כך, בעברית ובאנגלית; ולפני שמאמצים את הנקודה – מריצים בה הרצות אישור.

טעויות נפוצות

"גורם אחד בכל פעם". מזיזים כל גורם בנפרד ומחמיצים אינטראקציות ואת כיוון השיפור האמיתי – המסלול התלול ביותר משנה את כל הגורמים יחד, בפרופורציה למקדמים (Montgomery, 2020, p. 415).

מודל מסדר שני ממערך שאינו תומך בו. ב-2^k עם נקודות מרכז כל הריבועיים זהים, ואי-אפשר לאמוד אותם בנפרד – רק את העקמומיות הכוללת (עמ' 413). מוסיפים נקודות ציריות.

לפרש נקודה סטציונרית מחוץ לאזור הניסוי כאילו היא המיטב. מחוץ לאזור המודל אינו תקף (עמ' 425).

לבחור את המודל לפי R² בלבד. R² עולה תמיד עם עוד איברים; R² מתוקן, R² לחיזוי וחוסר ההתאמה אומרים אם התוספת שווה (לוח 11.7, עמ' 421).

להשוות מקדמים במשוואה הטבעית. הם תלויים ביחידות; משווים במקודדת (עמ' 422).

לוותר על נקודות מרכז. בלעדיהן אין טעות טהורה, אין מבחן חוסר התאמה, ואין בדיקת עקמומיות (עמ' 412-413).

לעגל את α. בדוגמה 11.2 הלוח מציג α = 1.414, והפלט חושב עם √2 – והתוצאות שונות בספרה השלישית (ממצא באימות העמוד). בתכנון שמים את α המדויק; אם רמה מעוגלת בפועל, מקלידים את הערך שהורץ.

מקורות

American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000

Montgomery, D. C. (2020). Design and analysis of experiments (10th ed.). Wiley.

Box ו-Hunter (1957) ו-Box ו-Behnken (1960) – כמצוטטים אצל Montgomery (2020). לפי כלל האתר, ברשימת המקורות מופיעים רק מקורות שנקראו במקור; הציטוטים לפי העמוד המודפס. פרק 11 של Montgomery נקרא במלואו מתמונות העמודים (שכבת הטקסט של הסריקה משמיטה את רוב המשוואות). החישוב בעמוד מאומת מול R 4.3.3 (lm, drop1, predict, hatvalues, rstudent, cooks.distance) וחבילת rsm (rsm, canonical, ccd, bbd), ומול כל המספרים המודפסים בדוגמאות 11.1 ו-11.2 ובלוחות 11.8 ו-11.11.