| UP |
Firstly, we need to create a data.frame. On this dataframe we execute a random effect model.
Steps to follow:
| 1) |
Create the data-frame on the data (textfile ch24pr17 ; 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 "\\".
>Ch24pr17<-read.table("c:\\temp\\Ch24pr17.txt");Ch24pr17
V1 V2 V3 V4
1 72.0
1 1 1
2 74.6
1 1 2
3 67.4
1 1 3
(rest of dataset)
Or directly from the webserver:
Ch24pr17<-read.table("http://www.biw.kuleuven.be/VAKKEN/statisticsbyR/datasetsTXT/CH24PR17.txt")
We give names to the variables.
names(Ch24pr17)=c("Value","Coats","Batch","Replic")
We attach the data.frame for direct access to the data.
attach(Ch24pr17)
Important is to transfrom the qualitative variables expressed by numbers into factors.
Coats<-factor(Coats)
Batch<-factor(Batch)
An useful visualisation is by calling the interaction plot function with (Factor A, Factor B , Response Var) as arguments
interaction.plot(Coats,Batch,Value)
Parallel lines or not? Look and judge carefully. What is your opinion?
We see differences depending on the Beer type (Bier)
More meaningfull is
boxplot(Value~Batch)
boxplot(Value~Coats)
What do you see?
Firstly check on the replications and the number of levels.
tapply(Value,Coats:Batch,length)
1:1 1:2 1:3 1:4 2:1 2:2 2:3 2:4 3:1 3:2 3:3
3:4
4 4 4
4 4 4 4 4
4 4 4 4
nA=nlevels(Coats);nB= nlevels(Batch);nA;nB
[1] 3
[1] 4
Now we can execute a two-way factor with interaction by lm
(the data are balanced so
lm can be used; however the interpretation in the case of a mixed effects
is specific; see chapter 24 ).
Let us fit the ANOVA model like we would do for a fixed effects; however the interpretation will be different.
ResLml<-lm(Value~Coats+Batch+Coats:Batch)
anova(ResLml)
Analysis of Variance Table
Response: Value
Df Sum Sq
Mean Sq F value Pr(>F)
Coats 2
150.388 75.194
15.591 1.327e-05 ***
Batch 3
152.852 50.951
10.564 3.984e-05 ***
Coats:Batch 6
1.852
0.309 0.064
0.9988
Residuals 36 173.625
4.823
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Compare F under H0 at alfa 5%: Do we accept Ho or Ha for interaction, fixed effect and the random effect?? What does it mean?
Values of F :
qf(0.95,6,36);qf(0.95,2,6);qf(0.95,3,36)
[1] 2.363751
[1] 5.143253
[1] 2.866266
Test factor A and factor B at alfa 5%. See tabel 24.5 and 24.6 in Neter et al (1996); pag 981-982.
Solve the exercise 24.17 c.
tapply(Value,Coats,mean)
1
2
3
73.10625 76.79375 76.92500
Bonferroni for two contrasts at alfa 10%
qt(0.975,6)
[1] 2.446912
Remark (niet voor Proeftechniek):
This analysis can also be done by specific functions for random and mixed models:
Firstly load the package nlme (standard but not in the base); and use lme. This function one has to specify
Data can be structured in a groupedData ( is data.frame with the grouping formula):
Ch24pr17<-read.table("http://www.biw.kuleuven.be/VAKKEN/statisticsbyR/datasetsTXT/CH24PR17.txt")
names(Ch24pr17)=c("Value","Coats","Batch","Replic")
Ch24pr17GD=groupedData(formula=Value~Coats|Batch,data=Ch24pr17)
Ch24pr17GD$Batch<-factor(Ch24pr17GD$Batch)
Ch24pr17GD$Coats<-factor(Ch24pr17GD$Coats)
Remark that the groupedData does not appear to take over the factors from the originating data.frame.
ResLme=lme(Value~Coats,random=~1|Batch,data=Ch24pr17GD);summary(ResLme)
Linear mixed-effects model fit by REML
Data: Ch24pr17GD
AIC
BIC logLik
217.8680 226.9013 -103.9340
Random effects:
Formula: ~1 | Batch
(Intercept) Residual
StdDev: 1.974262 2.044022
Fixed effects: Value ~ Coats
Value Std.Error DF
t-value p-value
(Intercept) 73.10625 1.1115549
42 65.76936 <.0001
Coats2 3.68750
0.7226709 42 5.10260
<.0001
Coats3 3.81875
0.7226709 42 5.28422
<.0001
Correlation:
(Intr) Coats2
Coats2 -0.325
Coats3 -0.325 0.500
Standardized Within-Group Residuals:
Min
Q1 Med
Q3 Max
-2.1990458 -0.6281479 0.1033012 0.6554653
1.3637802
Number of Observations: 48
Number of Groups: 4
This is a mixed model procedure. In this cas not a LS criterium is used
but the REML ( REstricted Maximum Likelihood) is used for fitting. This procedure is
appropiate for unbalanced data.
plot(ResLme)
intervals(ResLme,0.95)
Approximate 95% confidence intervals
Fixed effects:
lower
est. upper
(Intercept) 70.863041 73.10625
75.349459
Coats2 2.229091
3.68750 5.145909
Coats3 2.360341
3.81875 5.277159
Random Effects:
Level: Batch
lower
est. upper
sd((Intercept)) 0.8249685 1.974262 4.724678
Within-group standard error:
lower
est. upper
1.650472 2.044022
2.531413
Somewhat more easy for interpretation is when we remove the intercept and see the fixed effects for each Coating separately:
ResLme2=lme(Value~Coats-1,random=~1|Batch,data=Ch24pr17GD)
intervals(ResLme2,0.95)
Approximate 95% confidence intervals
Fixed effects:
lower
est. upper
Coats1 70.86304 73.10625 75.34946
Coats2 74.55054 76.79375 79.03696
Coats3 74.68179 76.92500 79.16821
Random Effects:
Level: Batch
lower
est. upper
sd((Intercept)) 0.8257672 1.974262
4.720109
Within-group standard error:
lower
est. upper
1.650480 2.044022 2.531401
| UP |
23 April, 2003 by Guido Wyseure