| UP |
Firstly, we need to create a data.frame. On this dataframe we execute a two-way ANOVA. Description of the Data-set describes an experiment of tomatos grown in a greenhouse in function of a fertilizer concentration. The experiment is done in a Randomized Complete Blocked Design (RCBD).
Steps to follow:
| 1) |
Create the data-frame on the data (textfile greenhouse-experiment ; right click on the file and save on C:\Temp ). Look to the file in R by File/display file. |
| 2) | Create factors in the frame |
| 3) | Execute the ANOVA two-way without interaction;
Check the the graphically the predictions, residuals, QQ-plot, Cook distance plot...for every ANOVA and note your observations |
| 4) | Balanced and fixed factors; so SS-decomposition is useful and correct
The General Linear Test or classical F tests can be done. |
| 5) | If appropiate execute multiple comparison or test specific contrast questions. |
We create a data-frame by typing (or copy/paste) the red letters after > followed by <ENTER> The blue letters are the answer from the R-programme. We check the content by typing in the name of the data-frame. Remark systems like R originate in Unix in which they discriminate uppercase and undercase and replace the usual backslash "\" by two "\\".
Greenhouse<-read.table("c:\\temp\\serreexperiment.txt",header=TRUE);Greenhouse
Treatm Block Dose response
1
2 1 25
40.51
2
1 1 0
33.98
(rest of dataset)
Or directly from the webserver
Greenhouse<-read.table("http://www.biw.kuleuven.be/vakken/statisticsbyR/ANOVAbyRr/serreExperiment.txt",header=TRUE)
We attach the data.frame for direct access to the data.
attach(Greenhouse)
Important is to transfrom the qualitative variables expressed by numbers into factors.
Block<-factor(Block);Dose<-factor(Dose)
Let us check if it is indeed a factor
> is.factor(Dose) ; is.factor(Block)
[1] TRUE
[1] TRUE
An useful visualisation is by calling the interaction plot function with (Factor A, Factor B , Response Var) as arguments
>
interaction.plot(Dose,Block,response)
Parallel lines or not? Look and judge carefully. What is your opinion?
Anyhow we cannot test the interaction between Block-factor and Dose-factor in a ANOVA-model.
Now we can execute a two-way factor without interaction ANOVA (the data are balanced so the decompostion of SSTR in SSA SSB is allowed and correct).
As a general approach (valid for unbalanced and balanced) we execute the GLT. In R just execute the Full model and the Reduced model and then ask execute the 'anova' function.
> RCBDfull<-lm(response~Dose+Block)
> anova(RCBDfull)
Analysis of Variance Table
Response:
response
Df Sum Sq Mean Sq
F value Pr(>F)
Dose 3
1139.02 379.67
40.437 0.0002249 ***
Block 2
774.83 387.42 41.262
0.0003114 ***
Residuals 6
56.33 9.39
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
We can accept that both factors are significantly different from the H0?
As it is balanced; it is correct to use the Sum of Squares-decompostion, but we could equally use the GLT:
RCBDnoblock<-lm(response~Dose)
RCBDnoDose<-lm(response~Block)
anova(RCBDnoDose,RCBDfull)
anova(RCBDnoblock,RCBDfull)
We get the same results with the GLT-approach as with the sum-decompostion.
We accept the two-factor model. For the model we accept we have to check the residuals. We can plot them on one page. We subdivide 4 windows on one page organized 2 by 2 by the layout command:
layout(matrix(c(1,2,3,4),2,2))
plot(RCBDfull)
Are we happy with the residuals???
We know we have differences in response in function of the factor Dose and the assumptions about the error are
summary(RCBDfull)
Call:
lm(formula = response ~ Dose + Block)
Residuals:
Min 1Q Median 3Q Max
-3.68667 -1.02667 0.07333 0.96250 3.72333
Coefficients:
Estimate Std. Error
t value Pr(>|t|)
(Intercept) 33.172
2.167 15.310 4.91e-06 ***
Dose25 7.570
2.502 3.026
0.023225 *
Dose50 13.990
2.502 5.592
0.001391 **
Dose100 26.563
2.502 10.617
4.11e-05 ***
Block2 -14.425
2.167 -6.658
0.000555 ***
Block3 -18.810
2.167 -8.681
0.000129 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 3.064 on 6 degrees of freedom
Multiple R-Squared: 0.9714, Adjusted R-squared: 0.9476
F-statistic: 40.77 on 5 and 6 DF, p-value: 0.0001485
One can see that the standard error on the estimate of the parameters is constant within a factor.
Remark: if we try to test interaction we get:
> RCBDinteract<-lm(response~Dose+Block+Dose:Block)
> anova(RCBDinteract)
Analysis of Variance Table
Response:
response
Df Sum Sq Mean Sq
F v alue Pr(>F)
Dose
3 1139.02 379.67
Block
2 774.83
387.42
Dose:Block
6 56.33
9.39
Residuals
0 0.00
It is not possible to evaluate the F-test by lack of degrees of freedom ( replicates are needed within each block)
But in general terms we can test the interaction by the Tukey additivity test. Load the function and execute:
tukeys.add.test(response,Dose,Block)
Tukey's one df F test for Additivity
data: Dose and Block on response
F = 0.4221, num df = 1, denom df = 5, p-value = 0.5446
sample estimates:
D estimate
0.007721938
Alternative Mixed model approach (NIET voor Proeftechniek)
An RCBD can be seen as a mixed model: fixed factor is the Dose and random factor is the Block.
Firstly load the package nlme
standard package but not in the base
either via menu: "Package/Load package" and select "nlme")
or by the comand:
library(nlme)
Loading required package: nls
Loading required package: lattice
and then
Greenhouse<-read.table("http://www.biw.kuleuven.be/vakken/statisticsbyR/ANOVAbyRr/serreExperiment.txt",header=TRUE)
Greenhouse$Block=factor(Greenhouse$Block)
Greenhouse$Dose=factor(Greenhouse$Dose)
Remark that attach gives problems in the package nlme (so we are not using it)
The data can be stored in a special data type: groupedData by
GH.grDat=groupedData(formula=response~Dose|Block,data=Greenhouse)
GH.grDat
Grouped Data: response ~ Dose | Block
Treatm Block
Dose response
1
2 1
25 40.51
2
1 1
0 33.98
3
4 1
100 58.78
etc....
GH.lme=lme(fixed=response~Dose-1,data=GH.grDat,random=~1|Block)
summary(GH.lme)
Gives in the console:
Compare this results to the classical two way ANOVA analysis.
| UP |
14 April, 2003 by Guido Wyseure