R is a free software environment for statistical computing and graphics.
http://www.r-project.org/
Home
| Register
| Getting Started
| Schedule
| References
| FAQ
| Participants
| Discussion | Tom's site
Join our e-mailing list for other courses
For more information, please contact Paul Geissler (Paul_Geissler@usgs.gov).
Michael Crawley, 2007,The R Book, Chapter 11
See the text for descriptions. You can copy and paste the statements below into thr R Commander script window and execute them.
Anything after # on a line is a comment. I have added annotations as comments and shown the output as comments following each command.
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/yields.txt",header=T); attach(d); d
# sand clay loam
# 1 6 17 13
# 2 10 15 16
# . . .
d1=data.frame(stack(d)); colnames(d1)=c("y","soil"); attach(d1); d1
# y soil
# 1 6 sand
# 2 10 sand
# . . .
mean=by(y,soil,mean); var=by(y,soil,var); cbind(mean,var)
# mean var
# clay 11.5 15.388889
# loam 14.3 7.122222
# sand 9.9 12.544444
fligner.test(y~soil)
# Fligner-Killeen test of homogeneity of variances
# data: y by soil
# Fligner-Killeen:med chi-squared = 0.3651, df = 2, p-value = 0.8332 # no problems
plot(y~soil)
m=lm(y~soil); anova(m)
Analysis of Variance Table
# Response: y
# Df Sum Sq Mean Sq F value Pr(>F)
# soil 2 99.200 49.600 4.2447 0.02495 *
# Residuals 27 315.500 11.685
# ---
# Signif. codes: 0 |***| 0.001 |**| 0.01 |*| 0.05 |.| 0.1 | | 1
par(mfrow=c(2,2));plot(m)
a=aov(y~soil); model.tables(a,se=T)
# Tables of effects
# soil
# clay loam sand
# -0.4 2.4 -2.0 # difference from overall mean
# Standard errors of effects
# soil
# 1.081 # not appropriate for comparing effects
# replic. 10
#### another example
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/competition.txt",header=T); attach(d); d
# biomass clipping
# 1 551 n25
# 2 457 n25
# . . .
m=lm(biomass~clipping); anova(m)
# Analysis of Variance Table
# Response: biomass
# Df Sum Sq Mean Sq F value Pr(>F)
# clipping 4 85356 21339 4.3015 0.008752 **
# Residuals 25 124020 4961
par(mfrow=c(1,2)); plot(clipping,biomass); bargraph.CI(clipping,biomass)
# Note that the error bars show standard error NOT confidence intervals.
#### factorial experiments ########################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/growth.txt",header=T); attach(d); d
# supplement diet gain
# 1 supergain wheat 17.37125
# . . .
means=tapply(gain,list(diet,supplement),mean);means
# agrimore control supergain supersupp
# barley 26.34848 23.29665 22.46612 25.57530
# oats 23.29838 20.49366 19.66300 21.86023
# wheat 19.63907 17.40552 17.01243 19.66834par(mfrow=c(1,2))
barplot(means,beside=T,legend.text=levels(diet),ylim=c(0,35))
bargraph.CI(supplement,gain,group=diet,legend=T,ylim=c(0,35))
m=lm(gain~diet*supplement); anova(m)
# Analysis of Variance Table
# Response: gain
# Df Sum Sq Mean Sq F value Pr(>F)
# diet 2 287.171 143.586 83.5201 2.999e-14 ***
# supplement 3 91.881 30.627 17.8150 2.952e-07 ***
# diet:supplement 6 3.406 0.568 0.3302 0.9166
# Residuals 36 61.890 1.719
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m) # showes individual degrees of freedom coefficients
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 26.3485 0.6556 40.191 < 2e-16 ***
# diet[T.oats] -3.0501 0.9271 -3.290 0.002248 **
# diet[T.wheat] -6.7094 0.9271 -7.237 1.61e-08 ***
# supplement[T.control] -3.0518 0.9271 -3.292 0.002237 **
# supplement[T.supergain] -3.8824 0.9271 -4.187 0.000174 ***
# supplement[T.supersupp] -0.7732 0.9271 -0.834 0.409816
# diet[T.oats]:supplement[T.control] 0.2471 1.3112 0.188 0.851571
# diet[T.wheat]:supplement[T.control] 0.8183 1.3112 0.624 0.536512
# diet[T.oats]:supplement[T.supergain] 0.2470 1.3112 0.188 0.851652
# diet[T.wheat]:supplement[T.supergain] 1.2557 1.3112 0.958 0.344601
# diet[T.oats]:supplement[T.supersupp] -0.6650 1.3112 -0.507 0.615135
# diet[T.wheat]:supplement[T.supersupp] 0.8024 1.3112 0.612 0.544381
# Residual standard error: 1.311 on 36 degrees of freedom
# Multiple R-squared: 0.8607, Adjusted R-squared: 0.8182
# F-statistic: 20.22 on 11 and 36 DF, p-value: 3.295e-12
#### nested designs and split plots #####################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/splityield.txt",header=T); attach(d); d
# yield block irrigation density fertilizer
# 1 90 A control low N
# . . .
# Fields (blocks) were the largest units
# irrigation or control was applied to half of each field
# Each irrigation plot was split into 3 and a seedling density applied (low, medium, or high)
# Each density plot was split into 3 and a fertilizer was applied (N, P, or NP)
# Each level has a different error term because of restrictions on randomization.
m=aov(yield~irrigation*density*fertilizer+Error(block/irrigation/density)); summary(m)
# Error: block
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 3 194.444 64.815
# Error: block:irrigation
# Df Sum Sq Mean Sq F value Pr(>F)
# irrigation 1 8277.6 8277.6 17.590 0.02473 *
# Residuals 3 1411.8 470.6
#
# Error: block:irrigation:density
# Df Sum Sq Mean Sq F value Pr(>F)
# density 2 1758.36 879.18 3.7842 0.05318 .
# irrigation:density 2 2747.03 1373.51 5.9119 0.01633 *
# Residuals 12 2787.94 232.33
#
# Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
# fertilizer 2 1977.44 988.72 11.4493 0.0001418 ***
# irrigation:fertilizer 2 953.44 476.72 5.5204 0.0081078 **
# density:fertilizer 4 304.89 76.22 0.8826 0.4840526
# irrigation:density:fertilizer 4 234.72 58.68 0.6795 0.6106672
# Residuals 36 3108.83 86.36
#### another example
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/rats.txt",header=T); attach(d); d
# Glycogen Treatment Rat Liver
# 1 131 1 1 1
# . . .
# Three experimental treatments were applied to rats, with two rats per treatment.
# Right, center and left portions of each rat's liver was sampled.
# Two pieces of each was analyzed separately for glygogen using different assays.
Treatment=factor(Treatment); Rat=factor(Rat); Liver=factor(Liver)
# As above must provide error term for nested design.
m=aov(Glycogen~Treatment+Error(Treatment/Rat/Liver)); summary(m)
# Error: Treatment
# Df Sum Sq Mean Sq
# Treatment 2 1557.56 778.78
#
# Error: Treatment:Rat
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 3 797.67 265.89
#
# Error: Treatment:Rat:Liver
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 12 594.0 49.5
#
# Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 18 381.00 21.17
1-pf(778.78/265.89,2,3) # p for treatment = 0.1970992
#### You can obtain the same analysis more simply by using the means for each independent application of the treatment.
d1=aggregate(d,list(Treatment,Rat),mean); d1$Treatment=factor(d1$Treatment); d1
# Group.1 Group.2 Glycogen Treatment Rat Liver # Group.1=Treatment, Group.2=Rat, but will work with character values.
# 1 1 1 132.5000 1 1 2
# 2 2 1 149.6667 2 1 2
# 3 3 1 134.3333 3 1 2
# 4 1 2 148.5000 1 2 2
# 5 2 2 152.3333 2 2 2
# 6 3 2 136.0000 3 2 2
m1=aov(d1$Glycogen~d1$Treatment); summary(m1) # need to specify d1, as d is assumed
# Df Sum Sq Mean Sq F value Pr(>F)
# d1$Treatment 2 259.593 129.796 2.929 0.1971 # same F and p value as above
# Residuals 3 132.944 44.315
#### WRONG analysis (pseudo replication) by failing to note nested design.
m2=aov(d$Glycogen~d$Treatment); summary(m2) # WRONG pseudoreplication
# Df Sum Sq Mean Sq F value Pr(>F)
# d$Treatment 2 1557.56 778.78 14.498 3.031e-05 ***
# Residuals 33 1772.67 53.72 # Only 6 rats received independent treatments: 1 df mean, 2 df treatments, 3 df error
#### effect sizes ##################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/Daphnia.txt",header=T); attach(d); d
# Growth.rate Water Detergent Daphnia
# 1 2.919086 Tyne BrandA Clone1
# . . .
m=aov(Growth.rate~Water*Detergent*Daphnia); anova(m)
# Analysis of Variance Table
#Response: Growth.rate
# Df Sum Sq Mean Sq F value Pr(>F)
# Water 1 1.985 1.985 2.8504 0.0978380 .
# Detergent 3 2.212 0.737 1.0586 0.3754783
# Daphnia 2 39.178 19.589 28.1283 8.228e-09 ***
# Water:Detergent 3 0.175 0.058 0.0837 0.9686075
# Water:Daphnia 2 13.732 6.866 9.8591 0.0002587 ***
# Detergent:Daphnia 6 20.601 3.433 4.9302 0.0005323 ***
# Water:Detergent:Daphnia 6 5.848 0.975 1.3995 0.2343235
# Residuals 48 33.428 0.696
model.tables(m,"means",se=T) # Works with aov models but not lm
# Tables of means
# Grand mean
# 3.851905
#
# Water
# Water
# Tyne Wear
# 3.686 4.018
#
# Detergent
# Detergent
# BrandA BrandB BrandC BrandD
# 3.885 4.010 3.955 3.558
#
# Daphnia
# Daphnia
# Clone1 Clone2 Clone3
# 2.840 4.577 4.139
#
# Water:Detergent # : and not * incicates an interaction
# Detergent
# Water BrandA BrandB BrandC BrandD
# Tyne 3.662 3.911 3.814 3.356
# Wear 4.108 4.109 4.095 3.760
#
# Water:Daphnia
# Daphnia
# Water Clone1 Clone2 Clone3
# Tyne 2.868 3.806 4.383
# Wear 2.812 5.348 3.894
#
# Detergent:Daphnia
# Daphnia
# Detergent Clone1 Clone2 Clone3
# BrandA 2.732 3.919 5.003
# BrandB 2.929 4.403 4.698
# BrandC 3.071 4.773 4.019
# BrandD 2.627 5.214 2.834
#
# Water:Detergent:Daphnia
# , , Daphnia = Clone1
#
# Detergent
# Water BrandA BrandB BrandC BrandD
# Tyne 2.811 2.776 3.288 2.597
# Wear 2.653 3.082 2.855 2.656
#
# , , Daphnia = Clone2
#
# Detergent
# Water BrandA BrandB BrandC BrandD
# Tyne 3.308 4.191 3.621 4.106
# Wear 4.530 4.615 5.925 6.322
#
# , , Daphnia = Clone3
#
# Detergent
# Water BrandA BrandB BrandC BrandD
# Tyne 4.867 4.766 4.535 3.366
# Wear 5.140 4.630 3.504 2.303
#
# Standard errors for differences of means
# Water Detergent Daphnia Water:Detergent Water:Daphnia Detergent:Daphnia Water:Detergent:Daphnia
# 0.1967 0.2782 0.2409 0.3934 0.3407 0.4818 0.6814
# replic. 36 18 24 9 12 6 3
replications(Growth.rate~Water*Detergent*Daphnia,d) # number of replications per treatment
# Water Detergent Daphnia Water:Detergent Water:Daphnia
# 36 18 24 9 12
# Detergent:Daphnia Water:Detergent:Daphnia
# 6 3
plot(t) # You can plot these confidence intervals, but with this number of comparisons, it is unreadable.
#### multiple comparisons #############################################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/Fungi.txt",header=T); attach(d); d
# Habitat Fugus.yield
# 1 Birch 20.838291
# . . .
m=aov(Fugus.yield~Habitat); anova(m)
# Analysis of Variance Table
# Response: Fugus.yield
# Df Sum Sq Mean Sq F value Pr(>F)
# Habitat 15 7527.4 501.8 72.141 < 2.2e-16 ***
# Residuals 144 1001.7 7.0
t=TukeyHSD(m); t # must se aov not lm
# Tukey multiple comparisons of means
# 95% family-wise confidence level
# Fit: aov(formula = Fugus.yield ~ Habitat)
# $Habitat
# diff lwr upr p adj
# Ash-Alder 3.53292777 -0.5808096 7.6466651 0.1844088
# Aspen-Alder 12.78574402 8.6720067 16.8994814 0.0000000
# Beech-Alder 12.32365349 8.2099161 16.4373908 0.0000000
# Birch-Alder 14.11348150 9.9997441 18.2272189 0.0000000
# Cherry-Alder 10.29508769 6.1813503 14.4088250 0.0000000
# Chestnut-Alder 12.24107899 8.1273416 16.3548163 0.0000000
# Holmoak-Alder -1.44360558 -5.5573429 2.6701318 0.9975654
# Hornbeam-Alder 10.60271044 6.4889731 14.7164478 0.0000000
# Lime-Alder 19.19458205 15.0808447 23.3083194 0.0000000
# Oak-Alder 20.29457340 16.1808360 24.4083108 0.0000000
# Pine-Alder 14.34084715 10.2271098 18.4545845 0.0000000
# Rowan-Alder 6.29495226 2.1812149 10.4086896 0.0000410
# Spruce-Alder -2.15119456 -6.2649319 1.9625428 0.9036592
# Sycamore-Alder 2.80900108 -1.3047363 6.9227384 0.5644643
# Willow-Alder 2.77635167 -1.3373857 6.8900890 0.5848838
# Aspen-Ash 9.25281625 5.1390789 13.3665536 0.0000000
# Beech-Ash 8.79072572 4.6769884 12.9044631 0.0000000
# Birch-Ash 10.58055373 6.4668164 14.6942911 0.0000000
# Cherry-Ash 6.76215993 2.6484226 10.8758973 0.0000065
# Chestnut-Ash 8.70815122 4.5944139 12.8218886 0.0000000
# Holmoak-Ash -4.97653335 -9.0902707 -0.8627960 0.0042690
# Hornbeam-Ash 7.06978268 2.9560453 11.1835200 0.0000018
# Lime-Ash 15.66165428 11.5479169 19.7753916 0.0000000
# Oak-Ash 16.76164563 12.6479083 20.8753830 0.0000000
# Pine-Ash 10.80791938 6.6941820 14.9216567 0.0000000
# Rowan-Ash 2.76202449 -1.3517129 6.8757618 0.5938269
# Spruce-Ash -5.68412232 -9.7978597 -1.5703850 0.0003946
# Sycamore-Ash -0.72392669 -4.8376640 3.3898107 0.9999996
# Willow-Ash -0.75657610 -4.8703135 3.3571613 0.9999993
# Beech-Aspen -0.46209053 -4.5758279 3.6516468 1.0000000
# Birch-Aspen 1.32773748 -2.7859999 5.4414748 0.9990469
# Cherry-Aspen -2.49065633 -6.6043937 1.6230810 0.7546149
# Chestnut-Aspen -0.54466504 -4.6584024 3.5690723 1.0000000
# Holmoak-Aspen -14.22934960 -18.3430870 -10.1156123 0.0000000
# Hornbeam-Aspen -2.18303358 -6.2967709 1.9307038 0.8929860
# Lime-Aspen 6.40883803 2.2951007 10.5225754 0.0000264
# Oak-Aspen 7.50882938 3.3950920 11.6225667 0.0000003
# Pine-Aspen 1.55510312 -2.5586342 5.6688405 0.9946232
# Rowan-Aspen -6.49079177 -10.6045291 -2.3770544 0.0000191
# Spruce-Aspen -14.93693858 -19.0506759 -10.8232012 0.0000000
# Sycamore-Aspen -9.97674295 -14.0904803 -5.8630056 0.0000000
# Willow-Aspen -10.00939235 -14.1231297 -5.8956550 0.0000000
# Birch-Beech 1.78982801 -2.3239093 5.9035654 0.9787523
# Cherry-Beech -2.02856580 -6.1423031 2.0851716 0.9381871
# Chestnut-Beech -0.08257450 -4.1963119 4.0311629 1.0000000
# Holmoak-Beech -13.76725907 -17.8809964 -9.6535217 0.0000000
# Hornbeam-Beech -1.72094305 -5.8346804 2.3927943 0.9852674
# Lime-Beech 6.87092856 2.7571912 10.9846659 0.0000042
# Oak-Beech 7.97091991 3.8571826 12.0846573 0.0000000
# Pine-Beech 2.01719366 -2.0965437 6.1309310 0.9408763
# Rowan-Beech -6.02870123 -10.1424386 -1.9149639 0.0001123
# Spruce-Beech -14.47484805 -18.5885854 -10.3611107 0.0000000
# Sycamore-Beech -9.51465241 -13.6283898 -5.4009151 0.0000000
# Willow-Beech -9.54730182 -13.6610392 -5.4335645 0.0000000
# Cherry-Birch -3.81839381 -7.9321312 0.2953435 0.1011027
# Chestnut-Birch -1.87240252 -5.9861399 2.2413348 0.9682125
# Holmoak-Birch -15.55708708 -19.6708244 -11.4433497 0.0000000
# Hornbeam-Birch -3.51077106 -7.6245084 0.6029663 0.1925375
# Lime-Birch 5.08110055 0.9673632 9.1948379 0.0030585
# Oak-Birch 6.18109190 2.0673545 10.2948293 0.0000633
# Pine-Birch 0.22736564 -3.8863717 4.3411030 1.0000000
# Rowan-Birch -7.81852924 -11.9322666 -3.7047919 0.0000001
# Spruce-Birch -16.26467606 -20.3784134 -12.1509387 0.0000000
# Sycamore-Birch -11.30448042 -15.4182178 -7.1907431 0.0000000
# Willow-Birch -11.33712983 -15.4508672 -7.2233925 0.0000000
# Chestnut-Cherry 1.94599129 -2.1677461 6.0597286 0.9558589
# Holmoak-Cherry -11.73869328 -15.8524306 -7.6249559 0.0000000
# Hornbeam-Cherry 0.30762275 -3.8061146 4.4213601 1.0000000
# Lime-Cherry 8.89949435 4.7857570 13.0132317 0.0000000
# Oak-Cherry 9.99948571 5.8857484 14.1132231 0.0000000
# Pine-Cherry 4.04575945 -0.0679779 8.1594968 0.0592074
# Rowan-Cherry -4.00013544 -8.1138728 0.1136019 0.0661686
# Spruce-Cherry -12.44628225 -16.5600196 -8.3325449 0.0000000
# Sycamore-Cherry -7.48608662 -11.5998240 -3.3723493 0.0000003
# Willow-Cherry -7.51873603 -11.6324734 -3.4049987 0.0000003
# Holmoak-Chestnut -13.68468457 -17.7984219 -9.5709472 0.0000000
# Hornbeam-Chestnut -1.63836854 -5.7521059 2.4753688 0.9908681
# Lime-Chestnut 6.95350306 2.8397657 11.0672404 0.0000030
# Oak-Chestnut 8.05349441 3.9397571 12.1672318 0.0000000
# Pine-Chestnut 2.09976816 -2.0139692 6.2135055 0.9193954
# Rowan-Chestnut -5.94612673 -10.0598641 -1.8323894 0.0001526
# Spruce-Chestnut -14.39227354 -18.5060109 -10.2785362 0.0000000
# Sycamore-Chestnut -9.43207791 -13.5458153 -5.3183406 0.0000000
# Willow-Chestnut -9.46472732 -13.5784647 -5.3509900 0.0000000
# Hornbeam-Holmoak 12.04631603 7.9325787 16.1600534 0.0000000
# Lime-Holmoak 20.63818763 16.5244503 24.7519250 0.0000000
# Oak-Holmoak 21.73817898 17.6244416 25.8519163 0.0000000
# Pine-Holmoak 15.78445273 11.6707154 19.8981901 0.0000000
# Rowan-Holmoak 7.73855784 3.6248205 11.8522952 0.0000001
# Spruce-Holmoak -0.70758898 -4.8213263 3.4061484 0.9999997
# Sycamore-Holmoak 4.25260666 0.1388693 8.3663440 0.0349691
# Willow-Holmoak 4.21995725 0.1062199 8.3336946 0.0380910
# Lime-Hornbeam 8.59187160 4.4781343 12.7056090 0.0000000
# Oak-Hornbeam 9.69186296 5.5781256 13.8056003 0.0000000
# Pine-Hornbeam 3.73813670 -0.3756007 7.8518741 0.1207078
# Rowan-Hornbeam -4.30775819 -8.4214955 -0.1940208 0.0302058
# Spruce-Hornbeam -12.75390500 -16.8676424 -8.6401676 0.0000000
# Sycamore-Hornbeam -7.79370937 -11.9074467 -3.6799720 0.0000001
# Willow-Hornbeam -7.82635878 -11.9400961 -3.7126214 0.0000001
# Oak-Lime 1.09999135 -3.0137460 5.2137287 0.9999007
# Pine-Lime -4.85373490 -8.9674723 -0.7399975 0.0062590
# Rowan-Lime -12.89962979 -17.0133671 -8.7858924 0.0000000
# Spruce-Lime -21.34577661 -25.4595140 -17.2320393 0.0000000
# Sycamore-Lime -16.38558097 -20.4993183 -12.2718436 0.0000000
# Willow-Lime -16.41823038 -20.5319677 -12.3044930 0.0000000
# Pine-Oak -5.95372625 -10.0674636 -1.8399889 0.0001484
# Rowan-Oak -13.99962114 -18.1133585 -9.8858838 0.0000000
# Spruce-Oak -22.44576796 -26.5595053 -18.3320306 0.0000000
# Sycamore-Oak -17.48557232 -21.5993097 -13.3718350 0.0000000
# Willow-Oak -17.51822173 -21.6319591 -13.4044844 0.0000000
# Rowan-Pine -8.04589489 -12.1596322 -3.9321575 0.0000000
# Spruce-Pine -16.49204170 -20.6057791 -12.3783043 0.0000000
# Sycamore-Pine -11.53184607 -15.6455834 -7.4181087 0.0000000
# Willow-Pine -11.56449548 -15.6782328 -7.4507581 0.0000000
# Spruce-Rowan -8.44614681 -12.5598842 -4.3324095 0.0000000
# Sycamore-Rowan -3.48595118 -7.5996885 0.6277862 0.2019434
# Willow-Rowan -3.51860059 -7.6323379 0.5951368 0.1896363
# Sycamore-Spruce 4.96019563 0.8464583 9.0739330 0.0044944
# Willow-Spruce 4.92754623 0.8138089 9.0412836 0.0049788
# Willow-Sycamore -0.03264941 -4.1463868 4.0810879 1.0000000
pairwise.t.test(Fugus.yield,Habitat,p.adjust.method="none")
# Pairwise comparisons using t tests with pooled SD
# p adjustment methods include: bonferroni, holm, hochberg, hommel, Benjamini & Hochberg ("BH"), Benjamini & Yekutieli ("BY"), and none.
# data: Fugus.yield and Habitat
# Alder Ash Aspen Beech Birch Cherry Chestnut Holmoak Hornbeam Lime Oak Pine Rowan
# Ash 0.00323 - - - - - - - - - - - -
# Aspen < 2e-16 8.9e-13 - - - - - - - - - - -
# Beech < 2e-16 7.8e-12 0.69581 - - - - - - - - - -
# Birch < 2e-16 1.4e-15 0.26218 0.13135 - - - - - - - - -
# Cherry 5.9e-15 5.6e-08 0.03645 0.08761 0.00150 - - - - - - - -
# Chestnut < 2e-16 1.1e-11 0.64494 0.94428 0.11461 0.10116 - - - - - - -
# Holmoak 0.22299 4.3e-05 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 - - - - - -
# Hornbeam 1.3e-15 1.6e-08 0.06625 0.14673 0.00342 0.79462 0.16697 < 2e-16 - - - - -
# Lime < 2e-16 < 2e-16 2.3e-07 3.6e-08 3.0e-05 4.7e-12 2.5e-08 < 2e-16 2.0e-11 - - - -
# Oak < 2e-16 < 2e-16 2.4e-09 3.2e-10 5.6e-07 2.5e-14 2.2e-10 < 2e-16 1.1e-13 0.35260 - - -
# Pine < 2e-16 4.7e-16 0.18945 0.08938 0.84742 0.00079 0.07715 < 2e-16 0.00187 6.5e-05 1.3e-06 - -
# Rowan 3.6e-07 0.02057 1.7e-07 1.0e-06 6.3e-10 0.00090 1.4e-06 9.0e-10 0.00036 < 2e-16 < 2e-16 2.3e-10 -
# Spruce 0.07026 3.6e-06 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.54952 < 2e-16 < 2e-16 < 2e-16 < 2e-16 3.8e-11
# Sycamore 0.01855 0.54035 2.8e-14 2.6e-13 < 2e-16 2.7e-09 3.8e-13 0.00043 7.1e-10 < 2e-16 < 2e-16 < 2e-16 0.00365
# Willow 0.01993 0.52226 2.4e-14 2.2e-13 < 2e-16 2.3e-09 3.3e-13 0.00047 6.1e-10 < 2e-16 < 2e-16 < 2e-16 0.00335
# Spruce Sycamore
# Ash - -
# Aspen - -
# Beech - -
# Birch - -
# Cherry - -
# Chestnut - -
# Holmoak - -
# Hornbeam - -
# Lime - -
# Oak - -
# Pine - -
# Rowan - -
# Spruce - -
# Sycamore 4.6e-05 -
# Willow 5.1e-05 0.97796
# P value adjustment method: none
#### contrasts - pages 368-386 #############################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/competition.txt",header=T); attach(d); d
# biomass clipping
# 1 551 n25
# . . .
clipping=factor(clipping)
m=lm(biomass~clipping); anova(m)
# Df Sum Sq Mean Sq F value Pr(>F)
# clipping 4 85356 21339 4.3015 0.008752 **
# Residuals 25 124020 4961
summary(m)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 465.17 28.75 16.177 9.4e-15 ***
# clipping[T.n25] 88.17 40.66 2.168 0.03987 *
# clipping[T.n50] 104.17 40.66 2.562 0.01683 *
# clipping[T.r10] 145.50 40.66 3.578 0.00145 **
# clipping[T.r5] 145.33 40.66 3.574 0.00147 **
# Residual standard error: 70.43 on 25 degrees of freedom
# Multiple R-squared: 0.4077, Adjusted R-squared: 0.3129
# F-statistic: 4.302 on 4 and 25 DF, p-value: 0.008752
contrasts(clipping)
# [T.n25] [T.n50] [T.r10] [T.r5] # each is compared to the mean
# control 0 0 0 0
# n25 1 0 0 0
# n50 0 1 0 0
# r10 0 0 1 0
# r5 0 0 0 1
contrasts(clipping)=cbind(
c(4,-1,-1,-1,-1),
c(0, 1, 1,-1,-1),
c(0, 1,-1, 0, 0),
c(0, 0, 0, 1,-1));
contrasts(clipping) # orthoganal (independent) becasue each sums to 0 and cross products sum to 0.
# [,1] [,2] [,3] [,4] # [,1] is control vs others
# control 4 0 0 0 # [,2] is n vs r
# n25 -1 1 1 0 # [,3] is n25 vs n50
# n50 -1 1 -1 0 # [,4] is r5 vs r10
# r10 -1 -1 0 1 # Contrasts are orthogonal (uncorrelated) because each column sums to 0
# r5 -1 -1 0 -1 # and cross products all sum to 0.
m1=lm(biomass~clipping); summary(m1)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 561.80000 12.85926 43.688 < 2e-16 ***
# clipping1 -24.15833 6.42963 -3.757 0.000921 *** # control vs others
# clipping2 -24.62500 14.37708 -1.713 0.099128 . # n vs r
# clipping3 -8.00000 20.33227 -0.393 0.697313 # n25 vs n50
# clipping4 0.08333 20.33227 0.004 0.996762 # r5 vs r10
# Residual standard error: 70.43 on 25 degrees of freedom
# Multiple R-squared: 0.4077, Adjusted R-squared: 0.3129
# F-statistic: 4.302 on 4 and 25 DF, p-value: 0.008752
#### More contrasts from John Fox, 2002, An R and S-Plus Companion to Applied Regression pages 127-130. ####
data(Prestige); d=Prestige; attach(d); d;
# education income women prestige census type
# GOV.ADMINISTRATORS 13.11 12351 11.16 68.8 1113 prof
# . . .
contrasts(type)
# [T.prof] [T.wc] # default is treatment or indicator coding, each relative to the first, which is used as a baseline
# bc 0 0
# prof 1 0
# wc 0 1
contrasts(type)=contr.treatment(levels(type), base=2); contrasts(type)
# bc wc # Change base for comparisons.
# bc 1 0
# prof 0 0
# wc 0 1
contrasts(type)="contr.helmert"; contrasts(type)
# [,1] [,2]
# bc -1 -1
# prof 1 -1
# wc 0 2
contrasts(type)="contr.sum"; contrasts(type)
# [,1] [,2]
# bc 1 0
# prof 0 1
# wc -1 -1
contrasts(type)='contr.SAS'; contrasts(type)
# bc wc
# bc 1 0
# wc 0 1
# prof 0 0
#### set levels for ordered factor
type=ordered(type, levels=c("bc","wc","prof")); type # Levels: bc < wc < prof
contrasts(type)="contr.poly"; round(contrasts(type),3)
# .L .Q
# [1,] -0.707 0.408
# [2,] 0.000 -0.816
# [3,] 0.707 0.408
#### Multivariate Analysis of Variance ################################################
d=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/manova.txt",header=T); attach(d); d
# tear gloss opacity rate additive
# 1 6.5 9.5 4.4 Low Low
# . . .
y=cbind(tear,gloss,opacity); m=manova(y~rate*additive); summary(m)
# Df Pillai approx F num Df den Df Pr(>F)
# rate 1 0.6181 7.5543 3 14 0.003034 **
# additive 1 0.4770 4.2556 3 14 0.024745 *
# rate:additive 1 0.2229 1.3385 3 14 0.301782
# Residuals 16
summary.aov(m) # univariate anova
# Response tear :
# Df Sum Sq Mean Sq F value Pr(>F)
# rate 1 1.74050 1.74050 15.7868 0.001092 **
# additive 1 0.76050 0.76050 6.8980 0.018330 *
# rate:additive 1 0.00050 0.00050 0.0045 0.947143
# Residuals 16 1.76400 0.11025
#
# Response gloss :
# Df Sum Sq Mean Sq F value Pr(>F)
# rate 1 1.30050 1.30050 7.9178 0.01248 *
# additive 1 0.61250 0.61250 3.7291 0.07139 .
# rate:additive 1 0.54450 0.54450 3.3151 0.08740 .
#
# Response opacity :
# Df Sum Sq Mean Sq F value Pr(>F)
# rate 1 0.421 0.421 0.1036 0.7517
# additive 1 4.901 4.901 1.2077 0.2881
# rate:additive 1 3.961 3.961 0.9760 0.3379
# Residuals 16 64.924 4.058
U.S. Department of the Interior |
U.S. Geological Survey
URL: http://www.fort.usgs.gov/BRDScience/learnR10.aspx
Page Contact Information: Paul_Geissler@usgs.gov