تمرین هشتم درس R

به نام خدا تمرین درس R نام دانشجو : سپیده جاوید نام استاد: آقای خاورزاده دستور lm و دستور aov

()Lm

برای برازش یک مدل خطی ساده به داده‌ها از تابع ()lm ، که به فرم کلی lm(formula, data) است، استفاده می‌کنیم که در آن شناسه‌ی formula ، فرم تابع خطی را که بصورت کلی

( ... + متغیر مستقل 2 + متغیر مستقل1) ~ متغیر وابسته

است مشخص می‌کند. در شناسه‌ی data ، داده‌هایی که مدل خطی به آن برازش داده می‌شود قرار می‌گیرد. از این دستور برای انجام رگرسیون، تحلیل قشر تک واریانس و تحلیل کوواریانس استفاده می‌شود. با اجرای این دستور، این تابع فقط برآورد ضرایب را محاسبه می‌کند. برای درک بهتر به ذکر مثالی می‌پردازیم که از قسمت help دستور آورده‌ایم. مثال:

Examples require(graphics)

برای رسم نمودارهای گرافیکی استفاده می‌شود.

## Annette Dobson (1990) "An Introduction to Generalized Linear Models". ## Page 9: Plant Weight Data

در این نرم‌افزار، با قرار دادن # می‌توان توضیحاتی راجع به دستور مورد نظر نوشت و تأثیری در خروجی نرم‌افزار ندارد. اینجا به همین صورت عمل شده است: ## آنت دابسون (1990) "مقدمه‌ای بر مدل‌های خطی تعمیم یافته". ## صفحه 9: کارخانه داده‌های وزن.

ctl <- c (4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) ctl [1] 4.17 5.58 5.18 6.11 4.50 4.61 5.17 4.53 5.33 5.14

trt <- c (4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) trt [1] 4.81 4.17 4.41 3.59 5.87 3.83 6.03 4.89 4.32 4.69

در قسمت بالا داده‌هایی بصورت برداری تعریف کرده و آنها را متغیر ctl و trt نامیده‌ است. خروجی‌ها، داده‌ها را از راست به چپ معرفی کرده است.

group <- gl (2, 10, 20, labels = c("Ctl","Trt")) group [1] Ctl Ctl Ctl Ctl Ctl Ctl Ctl Ctl Ctl Ctl Trt Trt Trt Trt Trt Trt Trt Trt Trt [20] Trt Levels: Ctl Trt

دستور ()gl ، به معنی این است که 2 متغیر ctl و trt را هر کدام 10 بار، به تعداد 20 بار تکرار می‌کند. در واقع داده‌های گروه را اینگونه ساخته و تعریف می‌کند.

weight <- c (ctl, trt) weight [1] 4.17 5.58 5.18 6.11 4.50 4.61 5.17 4.53 5.33 5.14 4.81 4.17 4.41 3.59 5.87 3.83 6.03 4.89 4.32 4.69

داده‌های دو متغیر ctl و trt را بصورت برداری ترکیب کرده و متغیر وزن نامیده است.

lm.D9 <- lm(weight ~ group) lm.D9

Call: lm(formula = weight ~ group)

Coefficients: (Intercept) groupTrt 5.032 -0.371

مدل رگرسیونی ساختیم که متغیر وزن، "متغیر وابسته" و متغیر گروه، "متغیر مستقل" می‌باشد. خروجی آن مقدار ضرایب را مشخص کرده است. می‌دانیم که خط رگرسیونی به صورت y= α+βx است؛ در واقع intercept ، مقدار α و groupTrt ، مقدار β است.

lm.D90 <- lm(weight ~ group - 1) # omitting intercept lm.D90

Call: lm(formula = weight ~ group - 1)

Coefficients: groupCtl groupTrt 5.032 4.661

وقتی (‌‌1-) آخر دستور lm بیاید، به معنی این است که ضریب ثابت ندارد (α ندارد). در قسمت توضیحات هم به حذف ضریب ثابت اشاره شده است.

anova(lm.D9)

Analysis of Variance Table

Response: weight Df Sum Sq Mean Sq F value Pr(>F) group 1 0.6882 0.68820 1.4191 0.249 Residuals 18 8.7293 0.48496

دستور ()anova ، جدول آنالیز واریانس را مشخص می‌کند و مقادیر درجه آزادی، مجموع مربعات، میانگین مربعات، آماره F و p-value را بدست می‌دهد. اگر p-value از 0.05 بیشتر باشد برازش مدل خوب نیست. در این مثال p-value از 0.05 بیشتر شده است پس مدل، مدل خوبی نیست.

summary(lm.D90)

Call: lm(formula = weight ~ group - 1)

Residuals: Min 1Q Median 3Q Max -1.0710 -0.4938 0.0685 0.2462 1.3690

Coefficients: Estimate Std. Error t value Pr(>|t|) groupCtl 5.0320 0.2202 22.85 9.55e-15 *** groupTrt 4.6610 0.2202 21.16 3.62e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6964 on 18 degrees of freedom Multiple R-squared: 0.9818, Adjusted R-squared: 0.9798 F-statistic: 485.1 on 2 and 18 DF, p-value: < 2.2e-16

دستور ()summary ، خلاصه‌ی کلی از مدل را بدست می‌دهد. به عبارت دیگر نتایج کامل‌تری از مشاهدات مدل را در اختیار ما می‌گذارد. علاوه بر برآورد ضرایب، اطلاعات و آماره‌های مفید دیگری را نیز بدست می‌دهد. خروجی بالا قسمت باقی‌مانده‌ها (خطاها)، مقادیر ماکزیمم و مینیمم، چارک‌های اول و سوم و میانه نمایش داده شده است. مجموع خطاها همیشه صفر است در این صورت میانگین هم صفر است. قسمت ضرایب، مقادیر برآورد، انحراف استاندارد و p-value مشخص شده است. اگر p-value بزرگتر از 0.05 بیشتر باشد فرض صفر بودن ضریب رد می‌شود و بالعکس. سطر بعد سطوح معنی‌دار را نشان داده است.

opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0)) plot(lm.D9, las = 1) # Residuals, Fitted

دستور ()par ، صفحه نمایش نمودارها را با 2 سطر و 2 ستون تقسیم‌بندی می‌کند.خروجی آن، ابتدا صفحه خالی است؛ با وارد کردن دستور ()plot صفحه‌ی نمودارها نمایش داده می‌شود.

دستور aov

()Aov

از این دستور برای انجام تجزیه و تحلیل مدل واریانس استفاده می‌شود. هم‌چنین در تحلیل آماری طرح آزمایش‌ها کاربرد دارد. تابع ()aov محاسبات تحلیل واریانس را برای طرح آزمایش‌ها نیز انجام می‌دهد. در واقع تجزیه و تحلیل واریانس (ANOVA) یک روش آماری است که معمولا برای بررسی مقایسه‌ی زیر مجموعه‌ای از داده‌ها استفاده می‌شود. در مورد پایه یک طرفه ANOVA، از آزمون t دو نمونه برای گروه های مستقل پوشش شرایط که در آن بیش از دو گروه مورد مقایسه قرار می‌گیرد، استفاده می‌شود. در آنالیز واریانس یک طرفه، گروه‌ها بر اساس یک عامل طبقه بندی شده‌اند و اصطلاحات استاندارد مورد استفاده، برای توصیف مجموعه‌ای از سطوح عامل کاربرد خاص دارند. تغییر در اندازه‌گیری‌های گرفته شده بر روی اجزاء منحصر به فرد از مجموعه داده‌ها تأثیر دارند و ANOVA، اینکه آیا این تغییر را می‌توان با گروه بندی‌های معرفی شده توسط عامل طبقه بندی بررسی کرد، توضیح داده است. فرم کلی این دستور بصورت زیر است: AOV (فرمول، داده = NULL، پیش بینی = FALSE، QR = TRUE، تضاد = NULL، ...) برای درک بهتر به ذکر مثالی می‌پردازیم که از قسمت help دستور آورده‌ایم. مثال 1:

Examples ## From Venables and Ripley (2002) p.165. ## Set orthogonal contrast

همانطور که قبلا گفته شد، با قرار دادن # می‌توان توضیحاتی راجع به دستور مورد نظر نوشت و تأثیری در خروجی نرم‌افزار ندارد. ## از ونبلز و ریپلی (2002) p.165. ## تنظیم ارتوگونال.

op <- options(contrasts = c("contr.helmert", "contr.poly")) op

$contrasts unordered ordered "contr.treatment" "contr.poly"

دستور ()options ، به ما اجازه می‌دهد برای تنظیم و بررسی از گزینه‌های مختلف جهانی استفاده کنیم که روشی است که در R محاسبه می‌شود و در نمایش نتایج آن تاثیر می‌گذارد. خروجی بالا، اثر تیمار را بصورت مرتب نشده و اثرهای متقابل را مرتب شده نشان داده است.

npk.aov <- aov(yield ~ block + N*P*K, npk)

Call: aov(formula = yield ~ block + N * P * K, data = npk)

Terms: block N P K N:P N:K P:K Sum of Squares 343.2950 189.2817 8.4017 95.2017 21.2817 33.1350 0.4817 Deg. of Freedom 5 1 1 1 1 1 1 Residuals Sum of Squares 185.2867 Deg. of Freedom 12

Residual standard error: 3.929447 1 out of 13 effects not estimable Estimated effects are balanced

آزمون aov را برای داده های عملکرد و بلوک با در نظر گرفتن N * P * K انجام داده است. خروجی آن شامل مقادیر مجموع مربعات و درجات آزادی برای بلوک، K ،P ،N و اثرات متقابل N:P ،N:K ،P:K است؛ هم‌چنین برای باقی‌مانده‌ها نیز مجموع مربعات و درجه آزادی محاسبه شده است.

summary(npk.aov)

Df Sum Sq Mean Sq F value Pr(>F) block 5 343.3 68.66 4.447 0.01594 * N 1 189.3 189.28 12.259 0.00437 ** P 1 8.4 8.40 0.544 0.47490 K 1 95.2 95.20 6.166 0.02880 * N:P 1 21.3 21.28 1.378 0.26317 N:K 1 33.1 33.14 2.146 0.16865 P:K 1 0.5 0.48 0.031 0.86275 Residuals 12 185.3 15.44 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

همانطور که قبلا گفته شد، این دستور خلاصه‌ی کلی از آزمون انجام شده را بدست می‌دهد. به عبارت دیگر نتایج کامل‌تری از مشاهدات آزمون را در اختیار ما می‌گذارد. برای نمایش جدول تحلیل واریانس aov از تابع ()summary استفاده می‌کنیم. خروجی مقادیر درجه آزادی، مجموع مربعات، میانگین مربعات، آماره F و p-value را برای بلوک، K ،P ،N و اثرات متقابل آنها مشخص کرده است. اگر p-value بزرگتر از 0.05 بیشتر باشد فرض صفر بودن میانگین‌ها رد می‌شود و بالعکس. سطر بعد سطوح معنی‌دار را نشان داده است.

coefficients(npk.aov)

(Intercept) block1 block2 block3 block4 block5 N1 P1 K1 N1:P1 N1:K1 P1:K1 54.8750000 1.7125000 1.6791667 -1.8229167 -1.0137500 0.2950000 2.8083333 -0.5916667 -1.9916667 -0.9416667 -1.1750000 0.1416667

دستور ()coefficients : ضریب، یک تابع عمومی که عصاره ضرایب مدل از اشیاء بازگردانده شده توسط توابع مدل سازی است. در واقع تمام مقادیر ضرایب مربوط به آزمون را معین کرده است.

## to show the effects of re-ordering terms contrast the two fits

برای نشان دادن اثراتی از اصطلاحات مرتب شده مقابل دو تناسب

aov(yield ~ block + N * P + K, npk)

Call: aov(formula = yield ~ block + N * P + K, data = npk)

Terms: block N P K N:P Residuals Sum of Squares 343.2950 189.2817 8.4017 95.2017 21.2817 218.9033 Deg. of Freedom 5 1 1 1 1 14

Residual standard error: 3.954232 Estimated effects are balanced

آزمون aov را برای داده های عملکرد و بلوک با در نظر گرفتن N * P + K انجام داده است. خروجی آن شامل مقادیر مجموع مربعات و درجات آزادی برای بلوک، K ،P ،N و اثر متقابل N:P و باقی‌مانده‌ها است.

aov(terms(yield ~ block + N * P + K, keep.order = TRUE), npk)

Call: aov(formula = terms(yield ~ block + N * P + K, keep.order = TRUE), data = npk)

Terms: block N P N:P K Residuals Sum of Squares 343.2950 189.2817 8.4017 21.2817 95.2017 218.9033 Deg. of Freedom 5 1 1 1 1 14

Residual standard error: 3.954232 Estimated effects are balanced

دستور ()terms : شرایط تابع، یک تابع عمومی که می‌تواند مورد استفاده برای استخراج شرایط اشیاء از انواع مختلف اشیاء داده R است. منظور از keep.order= TRUE ، یعنی اینکه حفظ ترتیب مد نظر است. ادامه عملیات همانند قبل صورت گرفته است. آزمون aov با خروجی شامل مقادیر مجموع مربعات و درجات آزادی می‌باشد.

## as a test, not particularly sensible statistically

به عنوان یک آزمون، به خصوص معقول از نظر آماری

npk.aovE <- aov(yield ~ N*P*K + Error(block), npk) npk.aovE

Call: aov(formula = yield ~ N * P * K + Error(block), data = npk) Grand Mean: 54.875 Stratum 1: block Terms: N:P:K Residuals Sum of Squares 37.00167 306.29333 Deg. of Freedom 1 4

Residual standard error: 8.750619 Estimated effects are balanced

Stratum 2: Within Terms: N P K N:P N:K P:K Residuals Sum of Squares 189.28167 8.40167 95.20167 21.28167 33.13500 0.48167 185.28667 Deg. of Freedom 1 1 1 1 1 1 12

Residual standard error: 3.929447 Estimated effects are balanced

آزمون aov را برای متغیر عملکرد و N * P * K با در نظر گرفتن خطای بلوک انجام داده است. خروجی شامل میانگین کل و نیز دو طبقه است. طبقه‌ی اول، مقادیر مجموع مربعات و نیز درجات آزادی را برای N:P:K و باقی‌مانده‌ها مشخص کرده است و طبقه‌ی دوم، مقادیر مجموع مربعات و نیز درجات آزادی K ،P ،N و اثرات متقابل N:P ،N:K ،P:K و باقی‌مانده‌ها است.

summary(npk.aovE)

Error: block Df Sum Sq Mean Sq F value Pr(>F) N:P:K 1 37.0 37.00 0.483 0.525 Residuals 4 306.3 76.57

Error: Within Df Sum Sq Mean Sq F value Pr(>F) N 1 189.28 189.28 12.259 0.00437 ** P 1 8.40 8.40 0.544 0.47490 K 1 95.20 95.20 6.166 0.02880 * N:P 1 21.28 21.28 1.378 0.26317 N:K 1 33.14 33.14 2.146 0.16865 P:K 1 0.48 0.48 0.031 0.86275 Residuals 12 185.29 15.44 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

مطابق قبل، این دستور خلاصه‌ی کلی از آزمون انجام شده را بدست می‌دهد. به عبارت دیگر نتایج کامل‌تری از مشاهدات آزمون را در اختیار ما می‌گذارد. برای نمایش جدول تحلیل واریانس aov از تابع() summary استفاده می‌کنیم. خروجی نرم‌افزار برای خطای بلوک، مقادیر درجه آزادی، مجموع مربعات، میانگین مربعات، آماره F و p-value را برای N:P:K و باقی‌مانده‌ها مشخص کرده است و قسمت بعد مقادیر درجه آزادی، مجموع مربعات، میانگین مربعات، آماره F و p-value را برای K ،P ،N و اثرات متقابل N:P ،N:K ،P:K و باقی‌مانده‌ها مشخص کرده است. اگر p-value بزرگتر از 0.05 بیشتر باشد فرض صفر بودن میانگین‌ها رد می‌شود و بالعکس. سطر بعد سطوح معنی‌دار را نشان داده است.

مثال 2: ما یکی از مجموعه‌ داده‌های موجود در R مربوط به یک آزمایش رشد گیاه را در نظر می‌گیریم. هدف از این آزمایش برای مقایسه عملکرد در گیاهان برای گروه کنترل و توجه دو تیمار است. متغیر واکنش یک اندازه‌گیری بر روی وزن خشک از گیاهان می‌باشد.

plant.df <- antGrowth plant.df$group <- factor(plant.df$group, labels = c("Control", "Treatment 1", "Treatment 2")) plant.mod1 <- lm(weight ~ group, data = plant.df) plant.mod1

Call: lm(formula = weight ~ group, data = plant.df) Coefficients: (Intercept) groupTreatment 1 groupTreatment 2 5.032 -0.371 0.494

مدل بالا به یک سری از داده‌ها برازش داده شده است در این صورت می‌توانیم اقدامات لازم برای مطالعه‌ی بهتر در برازش داده‌ها و دیگر فرضیات مدل انجام داد. خروجی ضرایب دو تیمار را مشخص کرده است.

summary(plant.mod1)

Call: lm(formula = weight ~ group, data = plant.df) Residuals: Min 1Q Median 3Q Max -1.0710 -0.4180 -0.0060 0.2627 1.3690 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.0320 0.1971 25.527 <2e-16 *** groupTreatment 1 -0.3710 0.2788 -1.331 0.1944 groupTreatment 2 0.4940 0.2788 1.772 0.0877 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6234 on 27 degrees of freedom Multiple R-squared: 0.2641, Adjusted R-squared: 0.2096 F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591

خروجی مدل شواهدی از تفاوت در رشد متوسط برای دو تیمار نسبت به گروه کنترل نشان می‌دهد. برای تجزیه و تحلیل جدول واریانس برای این مدل می توان از دستور ANOVA استفاده کرد.

anova(plant.mod1) Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) group 2 3.7663 1.8832 4.8461 0.01591 * Residuals 27 10.4921 0.3886 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

این جدول تایید می‌کند که تفاوت بین دو گروه معنی‌دار است؛ بطور خلاصه مدل برجسته و خوبی وجود دارد.

confint(plant.mod1) 2.5 % 97.5 % (Intercept) 4.62752600 5.4364740 groupTreatment 1 -0.94301261 0.2010126

groupTreatment 2 -0.07801261 1.0660126

تابع confint برای محاسبه فاصله اطمینان بر پارامترهای تیمار، با فواصل اطمینان پیش فرض 95٪ استفاده می‌شود.

پایان