[R] Split plot analysis problems

Mark Difford mark_difford at yahoo.co.uk
Thu Jul 23 21:11:13 CEST 2009


Hi Jean-Paul,

>> However, I've tried both solutions on my model, and I got different
>> residuals :...
>> What could be the difference between the two?

There is no difference. You have made a mistake.

##
tt <- data.frame(read.csv(file="tt.csv", sep=""))  ## imports your data set
T.aov <- aov(PH~Community*Mowing*Water + Error(Block), data=tt)

summary(T.aov)  ## This matches your output

                       Df Sum Sq Mean Sq F value    Pr(>F)    
Mowing                  1 40.007  40.007 21.1747 0.0002215 ***
Water                   1 23.120  23.120 12.2370 0.0025673 ** 
Community:Mowing        1 21.060  21.060 11.1467 0.0036554 ** 
Community:Water         1  6.901   6.901  3.6524 0.0720478 .  
Mowing:Water            1  1.611   1.611  0.8527 0.3680090    
Community:Mowing:Water  1  0.858   0.858  0.4542 0.5089331    
Residuals              18 34.008   1.889                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

##
length(residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within))
[1] 24
length(proj(T.aov)[,"Residuals"])
[1] 24

## Details
residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within)

           9           10           11           12           13          
14           15           16 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118 
          17           18           19           20           21          
22           23           24 
 1.241311954  1.498811954 -0.616188046  0.483811954 -0.477827381
-1.660327381  2.684672619  1.924672619 
          25           26           27           28           29          
30           31           32 
-0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358
-0.751015358  0.836484642  1.181484642 

proj(T.aov)[,"Residuals"]
           9           10           11           12           13          
14           15           16           17           18           19          
20           21 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954 -0.616188046 
0.483811954 -0.477827381 
          22           23           24           25           26          
27           28           29           30           31           32 
-1.660327381  2.684672619  1.924672619 -0.465198010 -1.735198010 
1.502301990  0.839801990 -0.428515358 -0.751015358  0.836484642  1.181484642 

Regards, Mark.


Jean-Paul Maalouf wrote:
> 
> Thanks Mark and Richard for your propositions on how to extract residuals.
> 
> However, I've tried both solutions on my model, and I got different  
> residuals :
> If we consider the within residuals :
> 
> Mark's solution gave me a vector of 24 residuals :
> 
>> as.vector(residuals(aov(PH~Community*Mowing*Water +
>> Error(Block))$Within))
>   [1] -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882  
>   0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954  
> -0.616188046
> [12]  0.483811954 -0.477827381 -1.660327381  2.684672619  1.924672619  
> -0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358  
> -0.751015358
> [23]  0.836484642  1.181484642
> 
> and Richard's solution gave me a vector 32 residuals :
> 
> do.call(data.frame,proj(aov(PH~Community*Water*Mowing +
> Error(Block))))->proj
> proj$Within.Residuals
>   [1] -0.216875  1.745625 -1.174375 -0.354375  0.800625  0.460625  
> -0.799375 -0.461875 -0.404375 -0.274375  0.695625 -0.016875 -0.541875   
> 0.695625 -1.141875
> [16]  0.988125  0.589375  0.846875 -1.268125 -0.168125 -1.095625  
> -2.278125  2.066875  1.306875 -0.500625 -1.770625  1.466875  0.804375  
> -0.638125 -0.960625
> [31]  0.626875  0.971875
> 
> What could be the difference between the two?
> 
> Regards
> 
> JPM
> 
> 
> Quoting "Richard M. Heiberger" <rmh at temple.edu>:
> 
>> Jean-Paul Maalouf wrote:
>>> Do you have any idea on how can I verify preliminary assumptions in  
>>>  this model (normality of the residuals and variance homogeneity),   
>>> since R is not able to extract residuals?
>>
>> Of course, R extracts residuals.  Use the proj() function.  See ?proj
>> for the example
>> to get the projection of an aovlist object onto the terms of a linear
>> model
>>
>> proj(npk.aovE)
>>
>>
>> To get the projections into a simple data.frame, use
>> tmpE <- proj(npk.aovE)
>> tmpE.df <- do.call("data.frame", tmpE)
>> tmpE.df
>>
>> Mark Difford's solution effectively retrieved columns 3 and 10 from
>> tmpE.df
>>
>> Rich
> 
> 
> 
> -- 
> Jean-Paul Maalouf
> UMR 1202 BIOGECO
> Inra - Université Bordeaux 1
> Bâtiment B8, Avenue des Facultés
> 33405 Talence, France
> Tel : 05 40008772
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: http://www.nabble.com/Split-plot-analysis-problems-tp24589402p24632396.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list