تمرین R (طرح آزمایش ها)

تمرین درس R نام دانشجوها : منیر اسدی – زینب لاغری نام استاد : اقای رامین خاورزاده طرح آزمایش ها طرح آزمایش  تک عاملی در این نوع آزمایش ها متغیر اندازه گیری شده پاسخ واحد های آزمایشی به سطوح مختلف یک عامل است . هدف از انجام این گونه آزمایش ها تشخیص اثر سطوح مختلف یک عامل روی متغیر پاسخ است . وارد کردن داده ها برای داده هایی زیر که مربوط به اثر چهار نوع کاتالیزور بر غلظت یک جز از سه جز مواد ترکیب دهنده ی یک مایع است ، بررسی می کنیم که آیا چهار کاتالیزور اثر یکسانی دارند .

		کاتالیزور
D	C	B	A
52.9	52.0	55.3	55.2
51.0	54.2	54.5	56.3
52.1	55.4	56.0	55.0
51.7		55.3	55.8
		54.0

ابتدا داده ها را به صورت زیر در R وارد کنید :

()kata<-scan<

1: 55.2 56.3 55.0 55.8 54.0

6: 55.3 54.5 56.0 55.3

10: 52.0 54.2 55.4

13: 52.9 51.0 52.1 51.7

17:

Read 16 items

حال برای تفکیک این بردار به چهار ستون ابتدا بردار کاراکتری زیر را بسازید :

(( dens<-rep(LETTERS[1:4],c(5,4,3,4<

سپس ب استفاده از تابع  factor() ، بردار dens را به یک ساختار عاملی تبدیل کنید :

( dens<-factor(dens<

اکنون می توان یک چارچوب اطلاعاتی از بردارهای kata و dens به صورت زیر ساخت :

( kata.dens<-data.frame(dens,kata<

 kata.dens<

   dens kata

1     A 55.2

2     A 56.3

3     A 55.0

4     A 55.8

5     A 54.0

6     B 55.3

7     B 54.5

8     B 56.0

9     B 55.3

10    C 52.0

11    C 54.2

12    C 55.4

13    D 52.9

14    D 51.0

15    D 52.1

16    D 51.7

بررسی داده ها با استفاده از تابع plot.design()  می توان نمودار هایی برای متغیر پاسخ در هر یک از سطوح عامل رسم کرد:

(( par(mfrow=c(1,2<

( plot.design(kata.dens<

( plot.design(kata.dens,fun=median<

میانگین ها و میانه های تیمار

در شکل در نمودار سمت چپی (راست) ، علامت های روی خط عمودی نمایانگر میانگین های (میانه ای ) سطوح مختلف تیمار و خط افقی ، میانگین (میانه) کل داده ها است . در این مثال چون بین میانه و میانگین داده ها اختلاف زیادی نیست ، نقطه دور افتاده وجود ندارد . با استفاده از تابع plot()  می توان نمودارهای جعبه ای متغیر پاسخ را در هر سطح از عامل رسم کرد .

( plot(kata.dens<

نمودارهای جعبه ای برای هر تیمار

تحلیل آماری طرح : تابع aov() محاسبات تحلیل واریانس را برای این طرح انجام می دهد :

(aov.kata<-aov(kata~dens,data=kata.dens<

 aov.kata<

:Call

  ( aov(formula = kata ~ dens, data = kata.dens

:Terms

                    dens Residuals

Sum of Squares  31.26071  11.99367

Deg. of Freedom        3        12

Residual standard error: 0.9997361

Estimated effects may be unbalanced

برای نمایش جدول تحلیل واریانس (ANOVA) از تابع summary() استفاده می کنیم :

( summary(aov.kata<

            (Df Sum Sq Mean Sq F value  Pr(>F  

**dens         3  31.26  10.420   10.43 0.00116 

Residuals   12  11.99   0.999                  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

مقادیر برازش داده شده و باقیمانده ها را می توان به ترتیب با به کارگیری توابع fitted.values() (یا مخفف آن fitted) و residual() (یا مخفف آن resid) روی تابع aov() به دست آورد ، مثلا برای مقادیر برازش شده داریم :

( fitted.values(aov.kata<

       1        2        3        4        5        6        7        8

55.26000 55.26000 55.26000 55.26000 55.26000 55.27500 55.27500 55.27500

       9       10       11       12       13       14       15       16

55.27500 53.86667 53.86667 53.86667 51.92500 51.92500 51.92500 51.92500

برای بررسی توزیع باقیمانده ها ، از بافت نگار و نمودار احتمال نرمال برای تابع residuals() استفاده می کنیم:

((hist(resid(aov.kata

((qqnorm(resid(aov.kata<

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

(( plot(fitted(aov.kata),resid(aov.kata<

با توجه به شکل این نمودار ساختار غیر معمولی را نشان نمی دهد . طرح کاملا تصادفی : در طرح کاملا تصادفی (CRD) اثر تعدادی تیمار هر یک در چند تکرار بر تغیر (صفت مورد نظر ) مورد بررسی قرار می گیرد. در این طرح یک منبع تغییرات قابل کنترل (تیمار) وجود دارد . برای آشنایی با نحوه  تجزیه واریانس مثال زیر بررسی می شود . مثال : در آزمایشی اثر هفت تیمار شامل سه مقدار کود گوگرد در فصل بهار و پاییز به همراه عدم مصرف کود (شاهد) هر یک در چهار تکرار بر میزان بیماری گال پوستی سیب زمینی مورد بررسی قرار گرفته است تیمارها در آزمایش بیماری گال پوستی سیب زمینی

G	F	E	D	C	B	A
17	18	30	10	16	9	30
7	24	7	4	10	9	18
16	12	21	4	18	16	32
17	19	9	5	18	4	26

برای خواندن داده ها (28 عدد) از دو روش وارد نمودن مستقیم و خواندن از طریق فایل های اکسل می توان استفاده نمود . ایتدا داده ها در یک فایل متنی (tab.delimited) با پسوند txt در فضای نرم افزار excel در یک آدرس فایل دلخواه (d://my data/…)  ذخیره و در محیط نرم افزار R با تابع read.table() به صورت زیر بازخوانی می گردد و در آبجکت data.potato  ذخیره می شود . برای مشاهده داده ها  می توانید آبجکت حاوی داده ها (data.potato)  را بازخوانی کنید :

> data.potato<-read.table("d://mydata/potato.txt",header=T)

> data.potato

Fertilizer Disease

1          A         30

2          A         18

3          A         32

4          A         26

5          B          9

.           .           .

.           .           .

.           .           .

28       G         17

28 عددی حاصل از آزمایش (متغیر Disease) و نوع تیمار هر یک (Fertilizer) نشان داده شده است . پس از حصول اطمینان از صحت خوانده شدن داده ها  تجزیه ی واریانس با تابع aov() صورت می گیرد. نتایج تجزیه واریانس در آبجکت anova.potato  نگهداری می شود . برای مشاهده نتایج تجزیه واریانس در قالب جدول از تابع summary استفاده می شود :

>anova.potato<-aov(disease~Fertillizer,data=data.potato)

Summary(anova.potato)

            Df Sum Sq Mean Sq F value  Pr(>F)

Fertillizer  6   1049  174.79  4.907   0.00279**

Residuals   21     748     35.62

--

Signif.codes:  0’***’ 0.001’**’ 0.01’*’  0.05’.’  0.1’ ‘ 1

همانطور که مشاهده می شود مدل تجزیه واریانس طرح CRD به صورت تیمار~ متغیر مورد ارزیابی نوشته شده است که شبیه به مدل y~x در تجزیه رگرسیون تابع lm() می باشد .

Data=data.potato به برنامه می گوید که متغیرهای Disease و Fertillizer در کجا قرار دارند . چنانچه از این برهان استفاده نشود ، باید قبل از نوشتن مدل داده های مورد نظر پیوست شوند

(attach(data.potato))

.

آزمون مقایسه میانگین LSD برای انجام مقایسه میانگین تیمارها با آزمون حداقل اختلافمعنی دار (LSD) در محیط نرم افزارR  می توانیم از تابع  pairwise.t.test() برای مدل تجزیه واریانس (یا آبجکت حاوی آن ) استفاده نماییم . برای آشنایی بیشتر با این تابع داده های مثال قبلی (اثر کود گوگرده بر بیماری گال پوستی در سیب زمینی ) را به کار می گیریم .

>Pairwise.t.test(disease,Fertilizer,p.adj=”none”)

          Pairwise comparisons using t test with pooled DS

Data:Disease and Fertilizer

  A       B          C          E          F

B 0.00061-   -           -           -           -

C 0.01648 0.16978 -         -           -           -

D 7.3e-05 0.38429 0.03112 -           -                -          

E 0.03112 0.10052 0.76999 0.01648   -                    -

F 0.06403   0.05063   0.52170   0.00744     0.72581   -

G 0.00851   0.27305   0.76999   0.05698    0.55991    0.35399

P value adjustment method:none

همانگونه که در خروجی فوق مشاهده می شود مقایسه دوبه دوی میانگین همه تیمار ها بر اساس LSD صورت گرفته و میزان احتمال خطا (α) برای هر مقایسه  ارائه شده است . به عنوان مثال تفاوت میانگین دو تیمار Aو B در سطح احتمال 00061/0   معنی دار است در حالی که میانگین دو تیمار A و F فاقد تفاوت معنی دار است . آزمون LDS با استفاده از بسته نرم افزاری agricolae ، برای انجام آزمون LSD باید از تابع  LSD.test() از بسته نرم افزاری  agricolae استفاده نماییم . چنانچه این بسته روی سیستم شما نصب نشده باشد ، لبتدا باید آن را لود و نصب نمایید .در صورتی که این بسته قبلا نصب شده باشد تنها کافی است آن را با تابع library() لود نمایید . تابع LSD.test() و برهان هایش به کاربر این امکان را می دهد تا میانگین صفت  مورد ارزیابی (Disease) را برای تیمارهای مورد بررسی (fertilizer) در سطح آماری دلخواه (alpha=) برای درجه آزادی (Dferror=) و واریانس خطای آزمایش (Mserror=) مورد مقایسه قرار دهند قبل از این که تابع را درج نماییم برای اختصار در برنامه نویسی بهتر است ابتدا داده ها را با تابع attach() بازخوانی نماییم.

>attach(data.potato)

>LSD.test(Disease,fertilizer,Dferror=21,Mserror=35.62,alpha=0.05)

Study:

LSD t Test for Disease

Mean Square Error:35.62

Fertilizer, means and individual (95%)Cl

  Disease        std.err   r       LCL      UCL Min. Max.

A   26.50       3.095696   4   20.062148   32.937852   18   32

B   9.50          2.466441   4   4.370754   14.629246   4   16

C  15.50        1.892969   4   11.563355   19.436645   10   18

D   5.75         1.436141   4   2.763382   8.736618   4   10

E   16.75        5.390965   4   5.538875   27.961125   7   30

F   18.25        2.462214   4   13.129545   23.370455   12   24

G   14.25       2.428134   4   9.200420    19.299580   8   17

Alpha:0.05 ; Df Error: 21

Critical Value of t: 2.079614

Least Significant Difference 8.776365

Means with the same letter are not significantly different.

Groups, Treatments and means

a          A         26.5

ab       F          18.25

b          E          16.75

b          C          15.5

bc        G         14.25

bc        B          9.5

c          D         5.75

همانطور که در خروجی برنامه مشاهده می شود بیشترین میانگین مربوط به تیمار (حرف a) و کمترین آن در D (حرفc) مشاهده می گردد . همچینین مقدار LSD در سطح 5 درصد 776365/8  شده است . به همین ترتیب با استفاده از توابع بسته نرم افزاری agricolae می توانید مقایسات میانگین باروش دانکن (تابعduncan.test()) و  SNK (تابع SNK.test()) را نیز انجام داد

آزمون duncan.test()

برای این آزمون مثالی در زیر داده شده است:

Examples

library(agricolae)

data(sweetpotato)

model<-aov(yield~virus,data=sweetpotato)

comparison <- duncan.test(model,”virus”,

main=”Yield of sweetpotato. Dealt with different virus”)

duncan.test(model,”virus”,alpha=0.01,console=TRUE)

# version old duncan.test()

attach(sweetpotato)

df<-df.residual(model)

Mserror<-deviance(model)/df

comparison <- duncan.test(yield,virus,df,Mserror, group=TRUE)

print(comparison$groups)

detach(sweetpotato)

آزمون LSD.test برای این آزمون مثالی در زیر داده شده است:

Examples

library(agricolae)

data(sweetpotato)

model<-aov(yield~virus, data=sweetpotato)

out <- LSD.test(model,”virus”, p.adj=”bonferroni”)

#stargraph

bar.group(out$groups,ylim=c(0,45),density=4,border=”blue”)

#endgraph

# Old version LSD.test()

attach(sweetpotato)

df<-df.residual(model)

Mserror<-deviance(model)/df

out <- LSD.test(yield,virus,df,Mserror, p.adj=”bonferroni”, group=FALSE)

detach(sweetpotato)

#stargraph

# Variation range: Max – Min

bar.err(out$means,variation=”range”,ylim=c(0,44),bar=FALSE,col=0)

#endgraph

LSD.test(model,”virus”,p.adj=”bon”,console=TRUE)

آزمون waller.test برای این آزمون مثالی در زیر داده شده است:

Examples

library(agricolae)data(sweetpotato)model<-aov(yield~virus, data=sweetpotato)out <- waller.test(model,”virus”, group=TRUE)#startgraphpar(mfrow=c(2,2))# variation: SE is error standard# variation: range is Max – Minbar.err(out$means,variation=”SD”,horiz=TRUE,xlim=c(0,45),bar=FALSE,col=colors()[25],space=2, main=”Standard deviation”,las=1)bar.err(out$means,variation=”SE”,horiz=FALSE,ylim=c(0,45),bar=FALSE,col=colors()[15],space=2,main=”SE”,las=1)bar.err(out$means,variation=”range”,ylim=c(0,45),bar=FALSE,col=”green”,space=3,main=”Range = Max – Min”,las=1)bar.group(out$groups,horiz=FALSE,ylim=c(0,45),density=8,col=”red”, main=”Groups”,las=1)#endgraph# Old version HSD.test()attach(sweetpotato)df<-df.residual(model)Mserror<-deviance(model)/dfFc<-anova(model)[“virus”,4]out <- waller.test(yield, virus, df, Mserror, Fc, group=TRUE)print(out$mean)detach(sweetpotato)

آزمون HSD.test برای این آزمون  مثالی در زیر داده شده است:

Examples

library(agricolae)

data(sweetpotato)

model<-aov(yield~virus, data=sweetpotato)

out <- HSD.test(model,”virus”, group=TRUE,console=TRUE,

main=”Yield of sweetpotato\nDealt with different virus”)

#stargraph

bar.group(out$groups,ylim=c(0,45),density=4,border=”blue”)

#endgraph

out<-HSD.test(model,”virus”, group=FALSE)

means<-out$means

# Old version HSD.test()

attach(sweetpotato)

df<-df.residual(model)

Mserror<-deviance(model)/df

HSD.test(yield,virus,df,Mserror, group=TRUE,console=TRUE,

main=”Yield of sweetpotato. Dealt with different virus”)

detach(sweetpotato)

آزمون مقایسه میانگین توکی برای انجام این آزمون اختلاف معنی دار قابل اعتماد توکی (HSD) در محیط نرم افزار R از تابع TukeyHSD() برای مدل تجزیه واریانس (یا آبجکت حاوی آن) استفاده می شود .

>TurkeyHSD(anova.potato)

 Tukey multiple comparisons of mean

  95% family-wise confidence level

Fit: aov(formula=Disease~Fertillizer,data=data.potato)

$Fertillizer

            Diff     lwe     upr     p          adj

B-A     -17      -30.718738 -3.281262    0.0092170

C-A     -11      -24.718738 2.718738     0.1738769

D-A     -20.75           -34.468738 -7.031262    0.0012167

E-A      -9.75  -23.468738 3.968738     0.2847369

.

.

.

G-E     -2.50  -16.218738 11.218738   09963747

G-F     -4.00  -17.718738 9.718738     0.9597399

همانگونه که مشاهده می کنید در مقایسه دو به دوی تیمارها همراه با اختلاف میانگین احتمال معنی دار بودن این اختلافات مشاهده می گردد . تنها تفاوت دو تیمار A و B در سطح 5درصد و تفاوت دو تیمار A و D در سطح 1 دزصد معنی دار ده است . ضریب تغیرات خطا (CV) برای محاسبه ضریب تغییرات خطا (CV%) در یک طرح آزمایشی می توانید با استفاده از اطلاعات موجود و 2 تابع mean() و sqrt() به صورت CV را محاسبه می نماییم .

>attach(data.potato)

>cv<-sqrt(35.62)-100/mean(Disease)

>cv

[1]39.22793

همانطور که مشاهده می شود CV خطا حدود 2/39 درصد است . البته چنانچه بسته نرم افزاری  agricolae  را نصب نمودید  می توانید تابع cv.model() را به صورت زیر به کار ببرید :

>library(agricolae)

>cv.model(anova.potato)

[1] 39.22741

طرح بلوک کامل تجزیه واریانس داده های حاصل از آزمایشاتی که در قالب طرح بلوک کامل (RB) اجرا شده اند نیز معمولا با استفاده از تابع aov() صورت می گیرد . به مثال زیر توجه کنید در آزمایشی عملکرد علوفه 5 رقم یونجه (رنجر A ،همدانیB،قره یونجهC ،محلی تبریز D، محلی خویE) در قالب یک طرح بلوک کامل با 4 تکرار مورد مقایسه قرار گرفته است .

داده های آزمایش مقایسه عملکرد علوفه 5 رقم یونجه
E	D	C	B	A	بلوک
محلی خوی	محلی تبریز	قره یونجه	همدانی	رنجر
78/5	74/5	82/6	17/8	15/7	1
62/6	63/4	58/8	78/8	74/7	2
81/5	56/5	06/6	0/8	39/6	3
22/4	94/4	78/6	75/6	63/4	4

قبل از شروع تجزیه تحلیل در R داده ها که شامل سه متغیر block  و variety و yield می باشند . ایتدا داده ها را با استفاه از تابع read.tabel() در محیط R خوانده شده و در آبجکت data.alfalfa نگهداری می شود  برای اطمینان از صحت خوانده شدن داده ها آنها را فرا خوانی می نماییم  سپس متغیر بلوک که سطوح آن عددی است با استفاده از تابع as.factor() به فاکتور تبدیل می شود . در مرحله بعد داده ها را جهت دسترسی مستقیم به آن ها پیوست می گردد .  داده ها با استفاده از تابع aov() مورد تجزیه واریانس قرار می گیرند  و برای نمایش نتایج در قالب جدول  تجزیه واریانس تابع summary() به کار برده می شود .

> data.alfalfa<-read.table("C:\\Users\\1401\\Desktop\\111.txt")

> data.alfalfa

   block variety yield

1      1       A  7.15

2      2       A  7.74

3      3       A  6.39

4      4       A  4.63

5      1       B  8.17

6      2       B  8.78

7      3       B  8.00

8      4       B  6.75

9      1       C  6.82

10     2       C  8.58

11     3       C  7.06

12     4       C  6.78

13     1       D  5.74

14     2       D  4.63

15     3       D  5.56

16     4       D  4.94

17     1       E  5.78

18     2       E  6.62

19     3       E  5.81

20     4       E  4.22

> data.alfalfa$block<-as.factor(data.alfalfa$block)

> attach(data.alfalfa)

> anova.alfalfa<-aov(yield~block+variety)

> summary(anova.alfalfa)

            Df Sum Sq Mean Sq F value   Pr(>F)            

block        3  8.619   2.873   6.821 0.006182 **

variety      4 20.513   5.128  12.174 0.000347 ***

Residuals   12  5.055   0.421                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

همانگونه که مشاهده می شود میانگین مربعات بلوک و تیمار (رقم) به ترتیب 887/2 و088/5 بود که هر دو در سطح احتمال یک درصد معنی دار شده اند . مدل تجزیه واریانس  طرح RB به صورت زیر نوشته شده است :

تیمار+بلوک ~ متغیر مورد ارزیابی

برای مقایسه میانگین تیمار ها با استفاده از آزمون LSD و یا دانکن می توانید از امکانات نرم افزاری agricolae استفاده می شود . چنانچه این بسته از قبل نصب شده باشد تنها کافی است با استفاده تز تابع library() آن را لود نایید :

Study:

LSD tTest for yield

Mean Square Error:0.421

Variety, means and individual(99%) Cl

          Yield std.err        r        LCL    UCL  Min. Max.

A       6.4775       0.6749738    4    4.415766  8.539234   4.63  7.74

B       7.9250       0.4259597    4    6.623889  9.226111   6.75  8.78

C       7.3100       0.4278240    4    6.003195  8.616805   6.78  8.58

D       5.2175       0.2602042    4    4.422696  6.012304    4.63  5.74

E        5.6300       0.50748267   4   4.078987  7.172013    4.22  6.62

Alpha:0.01 ; Df Error:12

Critical value of t:3.05454

Least Significant Difference 1.401431

Means with the same letter are not significantly different.

Groups, Treatments and means

a        B       7.925

 ab    C       7.31

bc      A       6.478

c        E        5.63

c        D       5.218

همانگونه که مشاهده می شود واریته B (یونجه همدانی) با عملکرد 925/7 بهترین (حرفa) و D واریته (محلی تبریز) و E (محلی خوی) ضعیف ترین (حرف c) بودند . تجزیه طرح بلوک با داده های گمشده چنانچه بنا به دلایلی یک یا چند واحد آزمایشی دز یک یا چند بلوک از دست برود تابع aov() توانایی تجزیه طرح بلوک به صورت نا متعادل را نیز دارد . جدول زیر داده های آزمایشی را نشان می دهد که در آن ارتفاع بوته چهار ژنوتیپ ذرت در شرایط مزرعه در قالب یک طرح بلوک کامل با 4 تکرار مورد بررسی قرار گرفته است . در حین آزمایش دو واحد آزمایشی که در جدول با علامت × مشخص شده اند از دست رفته اند . همانگونه که مشاهده می شود تابع aov() به خوبی این داده ها را که نامتعادل شده اند در قالب طرح بلوک مورد تجزیه قرار داد است .

ارقام ذرت
4	3	2	1	بلوک
125	155	135	161	1
130	167	141	175	2
128	190	×	149	3
×	160	127	167	4

14 داده موجود به صورت دستی با استفاده از تابع data.frame()  در آبجکت data.corn در نرم افزار R وارد شده اند . اگر چه این روش برای داده های آزمایشات واقعی که معمولا حجیم هستند چندان مناسب نمی باشد .

> data.corn<-data.frame(block=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4),variety=c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4),height=c(161,175,149,167,135,141,NA,127,155,167,190,160,125,130,125,NA))

برای مشاهده داده های وارد شده در آبجکت data.corn را بازخوانی می کنیم . همانطور که مشاهده می شود دو عدد از دست رفته را به صورت NA وارد نموده ایم .

> data.corn

   block variety height

1      1       1    161

2      2       1    175

3      3       1    149

4      4       1    167

5      1       2    135

6      2       2    141

7      3       2     NA

8      4       2    127

9      1       3    155

10     2       3    167

11     3       3    190

12     4       3    160

13     1       4    125

14     2       4    130

15     3       4    125

16     4       4     NA

با استفاده از دستور زیر داده ها را مورد تجزیه قرار می دهیم :

> anova.corn<-aov(height~as.factor(block)+as.factor(variety))

> summary(anova.corn)

                   Df Sum Sq Mean Sq F value  Pr(>F)  

as.factor(block)    3    253    84.5   0.618 0.62491  

as.factor(variety)  3   4322  1440.6  10.546 0.00546 **

Residuals           7    956   136.6                                       

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2 observations deleted due to missingness

همانگونه که مشاهده شد به دلیل اینکه سطوح variety و block به صورت عدد تعریف شده بود از تابع as.factor() استفاده شده است تا تبدیل به فاکتور گروه بندی شوند . بر خلاف مواد قبلی تابع as.factor() مستقیما درون تابع aov() استفاده شده است . مشاهده می شود که بین بلوک ها تفاوتی وجود ندارد ولی چهار رقم مورد بررسی از نظر ارتفاع بوته در سطح یک درصد تفاوت آماری دارند . در پایان خروجی نیز مشاهده می نمایید که پیغام وجود دو عدد گم شده ارائه شده است . آزمایشات فاکتوریل آزمایشات فاکتوریل آزمایشاتی هستند که در آنها اثرات بیش از یک فاکتور (تیمار) مورد بررسی قرار می گیرد . هدف اصلی این آزمایشات بررسی اثرات اصلی و اثر متقابل فاکتورها می باشد . نحوه ی تجزیه داده های حاصل از آزمایشات فاکتوریل بسته به طرح پایه ای که آزمایش در قالب آن انجام شده متفاوت خواهد بود . آزمایش فاکتوریل در قالب طرح کاملا تصادفی برای آشنایی با نحوه تجزیه واریانس آزمایشات فاکتوریل (در قالب طرح تصادفی) با استفاده از تابع aov() مثال زیر را مورد تجزیه واریانس قرار می دهیم . در یک آزمایش فاکتوریل افزایش وزن بدن نوعی موش صحرایی بر حسب گرم (weight)  با یک طرح پایه کاملا تصادفی تحت پوشش شش تیمار رژیم غذایی متشکل از دو فاکتور نوع پروتئین (دارای p1,p2,p3  ) و میزان پروتئین (دارای دو سطح low و high ) در چهار تکرار مورد بررسی قرار گرفته است و داده ها در جدول زیر تنظیم شده اند .

high			low
P3	P2	P1	P3	P2	P1
94	98	100	49	86	90
79	74	102	82	95	76
96	95	118	73	97	90
102	111	104	86	67	64

ابتدا با استفاده از توابع data.frame() و rep() داده ها را به آبجکت data.protein وارد می کنیم . برای آزمون صحت وارد نمودن و خوانده شدن داده ها توسط برنامه ، داده ها را فراخوانی نموده ایم . برای دسترسی مستقیم به متغییر های موجود در داده ها از تابع attach() استفاده شده است .سپس با استفاده از تابع aov() مدل تجزیه واریانس را تعریف و تجزیه را انجام می دهیم و در آبجکت anova.protein نگهداری می نماییم. برای ارائه نتایج در قالب جدول تجزیه واریانس از تابع summary() استفاده شده است :

> data.protein<-data.frame(amount=rep(c("low","high"),ech=12),protein=rep(c("p1","p2","p3"),2,each=4),weight=c(90,76,90,64,86,95,97,67,49,82,73,86,100,102,118,104,98,74,95,111,94,79,96,102))

> data.protein

   amount protein weight

1     low      p1     90

2    high      p1     76

3     low      p1     90

4    high      p1     64

5     low      p2     86

6    high      p2     95

7     low      p2     97

8    high      p2     67

9     low      p3     49

10   high      p3     82

11    low      p3     73

12   high      p3     86

13    low      p1    100

14   high      p1    102

15    low      p1    118

16   high      p1    104

17    low      p2     98

18   high      p2     74

19    low      p2     95

20   high      p2    111

21    low      p3     94

22   high      p3     79

23    low      p3     96

24   high      p3    102

> anova.protein<-aov(weight~protein*amount)

>summary(anova.protein)

          Df Sum Sq Mean Sq F value Pr(>F)

Protein       2   465.6    232.8      1.373   0.27860

Amount     1       1980.2   1980.2  11.680   0.00307**

Protein:amount  2   328.1   164.0    0.968   0.39889

Residuals             18   3051.5    169.5           

---

Signif.codes: 0’***’ 0.001’**’  0.01’*’ 0.05’.’ 0.1’ ‘ 1

همانگونه که مشاهده می شود تنها منبع متغیر amount (میزان پروتئین) معنی دار است (α=0.00307) همانطور که در بالا مشاهده می شود مدل تجزیه واریانس آزمایش فاکتوریل (با دو فاکتور) در قالب طرح کاملا تصادفی زیر نوشته شده است : فاکتور دوم *فاکتور اول ~ متغیر مورد بازیابی همانگونه که در مدل مشاهده می شود برای نوشتن مدل ها ی فاکتوریل از علامت "*" بین فاکتورها استفاده می کنیم . نوشتن A*B در مدل کلیه منابع تغیر ساده و فاکتوریل (A،B وA*B) لازم برای درج در جدول تجزیه واریانس را ایجاد می کند . اثرات ساده فاکتوریل ها در مثال شاملprotein  و amount و اثر متقابل دو فاکتور A و B در نر افزار R با عبارات A:B بیان می شود .