Background
Researchers want to assess the effect of three different diets on the
weight gain of fish. They have set up an experiment with 12 different
tanks of fish. Each tank contains the same number of fish. The weight of
6 fish in each tank was measured at the beginning and the end of the
experiment. The researchers recorded the weight gain.
Analysis
Pseudoreplication
The diets are randomized to the tanks, hence, the tanks are the
experimental units.
We measure the weight of fish, they are the observational
unit.
We measure multiple fish for each tank, hence the fish of the
same tank are pseudoreplicates. Indeed, fish from the
same tank are exposed to more similar conditions and their measurements
will be more similar and are thus not independent.
If we do not account for this pseudoreplication
our (regression) analysis and consider all fish as independent repeats,
the standard error on the estimate for the diet effect will be
underestimated, resulting in overly liberal inference.
To account for pseudoreplication in the data, we can average over the
pseudoreplicates because we assess the same number fish in each tank.
Hence, the tank averages will have an equal precision.
Note, that prior to averaging, the experimental unit was the tank,
while the observational unit was the fish inside the tank. After
averaging, the tank will be both the experimental and the observational
unit of the experiment.
Averaging can be achieved with the group_by
and
summarise
functions of the dplyr
R
package.
fish <- fish %>%
group_by(Tank, Diet) %>%
summarise(aveWtGain = mean(WtGain))
## `summarise()` has grouped output by 'Tank'. You can override using the
## `.groups` argument.
| | | | |
---|
1 | 1 | 2.833333 | | |
2 | 1 | 3.250000 | | |
3 | 1 | 5.700000 | | |
4 | 1 | 5.133333 | | |
5 | 2 | 29.600000 | | |
6 | 2 | 20.000000 | | |
7 | 2 | 18.550000 | | |
8 | 2 | 30.100000 | | |
9 | 3 | 45.966667 | | |
10 | 3 | 47.116667 | | |
Now, the weight gain values are averaged over de six fish in each of
the tanks.
Model
assumptions
List assumptions:
- The observations are independent
- Linearity between the response and predictor variable
- The residues of the model must be normally distributed
- Homoscedasticity of the data
After averaging the weight gain values over de six fish in each of
the tanks, we dealt with the pseudo-replication within tank and we
expect the tanks to be independent repeats.
To assess the remaining assumptions, we first fit a linear regression
model.
lmDiet <- lm(aveWtGain ~ Diet, fish)
par(mfrow=c(2,2))
plot(lmDiet)

We see some undesirable patterns in the diagnostic plots of the
linear model.
While the smoother for assessing linearity in the data is flat and
centered around zero, there is a much larger spread around the smoother
for the observations (tanks) of the second diet type (tanks 5-8). These
tanks are also flagged in the QQ-plot. Finally, the same tanks seam to
have a larger variability in Diet 2 and/or additional tank/batch effects
seem to occur. The tanks seem to cluster per two tanks and do not seem
to be independent. We therefore have to be careful with the
interpretation of the results.
Overall effect of
diet
Test (ANOVA)
lmDiet <- lm(aveWtGain ~ Diet, fish)
library(car)
anova_diet <- Anova(lmDiet, type=3)
anova_diet
| | | | |
---|
(Intercept) | 71.5434 | 1 | 4.914193 | 5.384469e-02 |
Diet | 3879.4074 | 2 | 133.234894 | 2.059484e-07 |
Residuals | 131.0267 | 9 | NA | NA |
There is an extremely significant effect of the diet on the average
weight gain of fish (p << 0.001). Note, that we have to be
cautious with the interpretation due to the violation of the model
assumptions
Post-hoc
analysis
library(multcomp)
multComp <- glht(model = lmDiet,
linfct = mcp(Diet = "Tukey"))
sumDiet <- summary(multComp)
sumDiet
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = aveWtGain ~ Diet, data = fish)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 20.333 2.698 7.536 <1e-04 ***
## 3 - 1 == 0 44.000 2.698 16.308 <1e-04 ***
## 3 - 2 == 0 23.667 2.698 8.772 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confDiet <- confint(multComp)
confDiet$confint
## Estimate lwr upr
## 2 - 1 20.33333 12.80085 27.86581
## 3 - 1 44.00000 36.46752 51.53248
## 3 - 2 23.66667 16.13419 31.19915
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.791865
Conclusion
There is an extremely significant effect of the diet on the average
weight gain of fish (p << 0.001). The average weight gain of fish
is extremely significantly different between all three diets (all p
<< 0.001).
The weight gain of fish in diet 3 is on average 44g and 23.7g higher
than that for fish that were fed with diet 1 and 2, respectively (95% CI
[36.5, 51.5] and [16.1, 31.2], respectively). The weight gain of fish in
diet 2 is on average 20.3g higher than that for fish that were fed with
diet 1 (95% CI [12.8, 27.9]).
Note, however, that the assumptions of test are
violated, i.e. larger variability in Diet 2 and/or additional
tank/batch effects in Diet 2. Indeed, the tanks seem to cluster per two
tanks, which may indicate additional batch effects. We therefore
have to be very careful with the interpretation of the
results.
LS0tCnRpdGxlOiAiRXhlcmNpc2UgOC54OiBUaGUgZmlzaCB0YW5rIGRhdGFzZXQgLSBzb2x1dGlvbiIgICAKYXV0aG9yOiAiTGlldmVuIENsZW1lbnQgYW5kIEplcm9lbiBHaWxpcyIKZGF0ZTogInN0YXRPbWljcywgR2hlbnQgVW5pdmVyc2l0eSAoaHR0cHM6Ly9zdGF0b21pY3MuZ2l0aHViLmlvKSIgIApvdXRwdXQ6CiAgICBodG1sX2RvY3VtZW50OgogICAgICBjb2RlX2Rvd25sb2FkOiB0cnVlICAgIAogICAgICB0aGVtZTogY29zbW8KICAgICAgdG9jOiB0cnVlCiAgICAgIHRvY19mbG9hdDogdHJ1ZQogICAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQogICAgcGRmX2RvY3VtZW50OgogICAgICB0b2M6IHRydWUKICAgICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlCiAgICAgIGxhdGV4X2VuZ2luZTogeGVsYXRleAoKLS0tCgojIEJhY2tncm91bmQgCgpSZXNlYXJjaGVycyB3YW50IHRvIGFzc2VzcyB0aGUgZWZmZWN0IG9mIHRocmVlIGRpZmZlcmVudCBkaWV0cyBvbiB0aGUgd2VpZ2h0IApnYWluIG9mIGZpc2guIFRoZXkgaGF2ZSBzZXQgdXAgYW4gZXhwZXJpbWVudCB3aXRoIDEyIGRpZmZlcmVudCB0YW5rcyBvZiBmaXNoLiAKRWFjaCB0YW5rIGNvbnRhaW5zIHRoZSBzYW1lIG51bWJlciBvZiBmaXNoLiBUaGUgd2VpZ2h0IG9mIDYgZmlzaCBpbiBlYWNoIHRhbmsgCndhcyBtZWFzdXJlZCBhdCB0aGUgYmVnaW5uaW5nIGFuZCB0aGUgZW5kIG9mIHRoZSBleHBlcmltZW50LiBUaGUgcmVzZWFyY2hlcnMgCnJlY29yZGVkIHRoZSB3ZWlnaHQgZ2Fpbi4KCiMgRXhwZXJpbWVudGFsIGRlc2lnbgoKLSBUaGUgZXhwbGFuYXRvcnkgdmFyaWFibGUgaW4gdGhlIGV4cGVyaW1lbnQgaXMgdGhlIGZhY3RvciB0eXBlIG9mIGRpZXQKCi0gVGhlcmUgYXJlIDMgZGlldCB0eXBlczogZGlldCAxLCBkaWV0IDIsIGFuZCBkaWV0IDMuCgotIFRoZSBleHBlcmltZW50YWwgdW5pdCBpcyB0aGUgdGFuazogdGhlIGRpZXQgdHJlYXRtZW50cyBhcmUgYWRvcHRlZCBvbiB0aGUgCnRhbmtzLCBpLmUuLCB0aGUgdGFua3MgYXJlIHJhbmRvbWl6ZWQgdG8gdGhlIGRpZXQgdHJlYXRtZW50cy4KCi0gVGhlIHdlaWdodCBnYWluIGlzIHRoZSByZXNwb25zZSBhbmQgaXQgaXMgbWVhc3VyZWQgb24gZWFjaCBmaXNoIHdoaWNoIGFyZSAKdGhlIG9ic2VydmF0aW9uYWwgdW5pdHMuCgpJbXBvcnQgbGlicmFyaWVzCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKIyBEYXRhIGltcG9ydAoKYGBge3J9CmZpc2ggPC0gcmVhZC50YWJsZSgiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL3N0YXRPbWljcy9QU0xTMjEvZGF0YS9maXNoVGFuay50eHQiLGhlYWRlcj1UUlVFKQpoZWFkKGZpc2gpCmBgYAoKIyBUaWR5IGRhdGEKCmBgYHtyfQpmaXNoIDwtIGZpc2ggJT4lCiAgICBtdXRhdGUoVGFuayA9IGFzLmZhY3RvcihUYW5rKSwgRGlldCA9IGFzLmZhY3RvcihEaWV0KSkKYGBgCgojIERhdGEgZXhwbG9yYXRpb24KCmBgYHtyfQpmaXNoICU+JQogIGdncGxvdChhZXMoeCA9IFRhbmssIHkgPSBXdEdhaW4sIGZpbGw9RGlldCkpICsKICAgIGdlb21fYm94cGxvdChvdXRsaWVyLnNoYXBlPU5BKSArIAogICAgZ2VvbV9qaXR0ZXIod2lkdGggPSAwLjIpICsKICAgIGdndGl0bGUoIldlaWdodCBnYWluIHBlciB0YW5rIikgKwogICAgeWxhYigiV2VpZ3RoIGdhaW4gKGcpIikgKyAKICAgIHN0YXRfc3VtbWFyeShmdW4gPSBtZWFuLCBnZW9tPSJwb2ludCIsIHNoYXBlPTUsIHNpemU9MywgY29sb3I9ImJsYWNrIiwgZmlsbD0iYmxhY2siKSArCiAgICBzY2FsZV9maWxsX2JyZXdlcihwYWxldHRlPSJSZEd5IikgKwogICAgdGhlbWVfYncoKQpgYGAKCiMgQW5hbHlzaXMKCiMjIFBzZXVkb3JlcGxpY2F0aW9uCgotIFRoZSBkaWV0cyBhcmUgcmFuZG9taXplZCB0byB0aGUgdGFua3MsIGhlbmNlLCB0aGUgdGFua3MgYXJlIHRoZSBleHBlcmltZW50YWwgdW5pdHMuIAotIFdlIG1lYXN1cmUgdGhlIHdlaWdodCBvZiBmaXNoLCB0aGV5IGFyZSB0aGUgb2JzZXJ2YXRpb25hbCB1bml0LiAKLSBXZSBtZWFzdXJlIG11bHRpcGxlIGZpc2ggZm9yIGVhY2ggdGFuaywgaGVuY2UgIHRoZSBmaXNoIG9mIHRoZSBzYW1lIHRhbmsgYXJlICoqcHNldWRvcmVwbGljYXRlcyoqLiBJbmRlZWQsIGZpc2ggZnJvbSB0aGUgc2FtZSB0YW5rIGFyZSBleHBvc2VkIHRvIG1vcmUgc2ltaWxhciBjb25kaXRpb25zIGFuZCB0aGVpciBtZWFzdXJlbWVudHMgd2lsbCBiZSBtb3JlIHNpbWlsYXIgYW5kIGFyZSB0aHVzIG5vdCBpbmRlcGVuZGVudC4gCgotIElmIHdlIGRvIG5vdCBhY2NvdW50IGZvciB0aGlzIAoqKnBzZXVkb3JlcGxpY2F0aW9uKiogb3VyIChyZWdyZXNzaW9uKSBhbmFseXNpcyBhbmQgY29uc2lkZXIgYWxsIGZpc2ggYXMgaW5kZXBlbmRlbnQgcmVwZWF0cywgdGhlIHN0YW5kYXJkIGVycm9yIG9uIHRoZSBlc3RpbWF0ZSBmb3IgdGhlIGRpZXQgZWZmZWN0IHdpbGwgYmUgdW5kZXJlc3RpbWF0ZWQsIHJlc3VsdGluZyBpbiBvdmVybHkgbGliZXJhbCBpbmZlcmVuY2UuCgpUbyBhY2NvdW50IGZvciBwc2V1ZG9yZXBsaWNhdGlvbiBpbiB0aGUgZGF0YSwgd2UgY2FuIGF2ZXJhZ2Ugb3ZlciB0aGUgcHNldWRvcmVwbGljYXRlcyBiZWNhdXNlIHdlIGFzc2VzcyB0aGUgc2FtZSBudW1iZXIgZmlzaCBpbiBlYWNoIHRhbmsuIApIZW5jZSwgdGhlIHRhbmsgYXZlcmFnZXMgd2lsbCBoYXZlIGFuIGVxdWFsIHByZWNpc2lvbi4gCgpOb3RlLCB0aGF0IHByaW9yIHRvIGF2ZXJhZ2luZywgdGhlIGV4cGVyaW1lbnRhbCB1bml0IHdhcyB0aGUgdGFuaywgd2hpbGUgdGhlCm9ic2VydmF0aW9uYWwgdW5pdCB3YXMgdGhlIGZpc2ggaW5zaWRlIHRoZSB0YW5rLiBBZnRlciBhdmVyYWdpbmcsIHRoZSB0YW5rCndpbGwgYmUgYm90aCB0aGUgZXhwZXJpbWVudGFsIGFuZCB0aGUgb2JzZXJ2YXRpb25hbCB1bml0IG9mIHRoZSBleHBlcmltZW50LgogCkF2ZXJhZ2luZyBjYW4gYmUgYWNoaWV2ZWQgd2l0aCB0aGUgYGdyb3VwX2J5YCBhbmQgYHN1bW1hcmlzZWAgZnVuY3Rpb25zIG9mIAp0aGUgYGRwbHlyYCBSIHBhY2thZ2UuCgpgYGB7cn0KZmlzaCA8LSBmaXNoICU+JSAKICAgIGdyb3VwX2J5KFRhbmssIERpZXQpICU+JQogICAgc3VtbWFyaXNlKGF2ZVd0R2FpbiA9IG1lYW4oV3RHYWluKSkKZmlzaApgYGAKCk5vdywgdGhlIHdlaWdodCBnYWluIHZhbHVlcyBhcmUgYXZlcmFnZWQgb3ZlciBkZSBzaXggZmlzaCBpbiBlYWNoIG9mIHRoZSB0YW5rcy4KCiMjIE1vZGVsIGFzc3VtcHRpb25zCgpMaXN0IGFzc3VtcHRpb25zOgoKMS4gVGhlIG9ic2VydmF0aW9ucyBhcmUgaW5kZXBlbmRlbnQKMi4gTGluZWFyaXR5IGJldHdlZW4gdGhlIHJlc3BvbnNlIGFuZCBwcmVkaWN0b3IgdmFyaWFibGUKMy4gVGhlIHJlc2lkdWVzIG9mIHRoZSBtb2RlbCBtdXN0IGJlIG5vcm1hbGx5IGRpc3RyaWJ1dGVkCjQuIEhvbW9zY2VkYXN0aWNpdHkgb2YgdGhlIGRhdGEKCkFmdGVyIGF2ZXJhZ2luZyB0aGUgd2VpZ2h0IGdhaW4gdmFsdWVzIG92ZXIgZGUgc2l4IGZpc2ggaW4gZWFjaCBvZiB0aGUgdGFua3MsCndlIGRlYWx0IHdpdGggdGhlIHBzZXVkby1yZXBsaWNhdGlvbiB3aXRoaW4gdGFuayBhbmQgd2UgZXhwZWN0IHRoZSB0YW5rcyB0byBiZSBpbmRlcGVuZGVudCByZXBlYXRzLiAKClRvIGFzc2VzcyB0aGUgcmVtYWluaW5nIGFzc3VtcHRpb25zLCB3ZSBmaXJzdCBmaXQgYSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbC4KCmBgYHtyfQpsbURpZXQgPC0gbG0oYXZlV3RHYWluIH4gRGlldCwgZmlzaCkKcGFyKG1mcm93PWMoMiwyKSkKcGxvdChsbURpZXQpCmBgYAoKV2Ugc2VlIHNvbWUgdW5kZXNpcmFibGUgcGF0dGVybnMgaW4gdGhlIGRpYWdub3N0aWMgcGxvdHMgb2YgdGhlIGxpbmVhciBtb2RlbC4KCldoaWxlIHRoZSBzbW9vdGhlciBmb3IgYXNzZXNzaW5nIGxpbmVhcml0eSBpbiB0aGUgZGF0YSBpcyBmbGF0IGFuZCBjZW50ZXJlZAphcm91bmQgemVybywgdGhlcmUgaXMgYSBtdWNoIGxhcmdlciBzcHJlYWQgYXJvdW5kIHRoZSBzbW9vdGhlciBmb3IgdGhlIApvYnNlcnZhdGlvbnMgKHRhbmtzKSBvZiB0aGUgc2Vjb25kIGRpZXQgdHlwZSAodGFua3MgNS04KS4gVGhlc2UgdGFua3MgYXJlCmFsc28gZmxhZ2dlZCBpbiB0aGUgUVEtcGxvdC4gRmluYWxseSwgdGhlIHNhbWUgdGFua3Mgc2VhbSB0byBoYXZlIGEgbGFyZ2VyIHZhcmlhYmlsaXR5IGluIERpZXQgMiBhbmQvb3IgYWRkaXRpb25hbCB0YW5rL2JhdGNoIGVmZmVjdHMgc2VlbSB0byBvY2N1ci4gVGhlIHRhbmtzIHNlZW0gdG8gY2x1c3RlciBwZXIgdHdvIHRhbmtzIGFuZCBkbyBub3Qgc2VlbSB0byBiZSBpbmRlcGVuZGVudC4gV2UgdGhlcmVmb3JlIGhhdmUgdG8gYmUgY2FyZWZ1bCB3aXRoIHRoZSBpbnRlcnByZXRhdGlvbiBvZiB0aGUgcmVzdWx0cy4KCgojIyBPdmVyYWxsIGVmZmVjdCBvZiBkaWV0CgojIyMgVGVzdCAoQU5PVkEpCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbG1EaWV0IDwtIGxtKGF2ZVd0R2FpbiB+IERpZXQsIGZpc2gpCmxpYnJhcnkoY2FyKQphbm92YV9kaWV0IDwtIEFub3ZhKGxtRGlldCwgdHlwZT0zKQphbm92YV9kaWV0CmBgYAoKVGhlcmUgaXMgYW4gZXh0cmVtZWx5IHNpZ25pZmljYW50IGVmZmVjdCBvZiB0aGUgZGlldCBvbiB0aGUgYXZlcmFnZSB3ZWlnaHQgZ2FpbiBvZiBmaXNoIChwIDw8IDAuMDAxKS4gCipOb3RlLCB0aGF0IHdlIGhhdmUgdG8gYmUgY2F1dGlvdXMgd2l0aCB0aGUgaW50ZXJwcmV0YXRpb24gZHVlIHRvIHRoZSB2aW9sYXRpb24gb2YgdGhlIG1vZGVsIGFzc3VtcHRpb25zKgoKIyMgUG9zdC1ob2MgYW5hbHlzaXMKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KG11bHRjb21wKQptdWx0Q29tcCA8LSBnbGh0KG1vZGVsID0gbG1EaWV0LCAKICAgICAgICAgICAgICAgICBsaW5mY3QgPSBtY3AoRGlldCA9ICJUdWtleSIpKQpzdW1EaWV0IDwtIHN1bW1hcnkobXVsdENvbXApCnN1bURpZXQKYGBgCgpgYGB7cn0KY29uZkRpZXQgPC0gY29uZmludChtdWx0Q29tcCkKY29uZkRpZXQkY29uZmludApgYGAKCiMjIENvbmNsdXNpb24KClRoZXJlIGlzIGFuIGV4dHJlbWVseSBzaWduaWZpY2FudCBlZmZlY3Qgb2YgdGhlIGRpZXQgb24gdGhlIGF2ZXJhZ2Ugd2VpZ2h0IGdhaW4gb2YgZmlzaCAocCA8PCAwLjAwMSkuIApUaGUgYXZlcmFnZSB3ZWlnaHQgZ2FpbiBvZiBmaXNoIGlzIGV4dHJlbWVseSBzaWduaWZpY2FudGx5IGRpZmZlcmVudCBiZXR3ZWVuIGFsbCB0aHJlZSBkaWV0cyAoYWxsIHAgPDwgMC4wMDEpLiAKClRoZSB3ZWlnaHQgZ2FpbiBvZiBmaXNoIGluIGRpZXQgMyBpcyBvbiBhdmVyYWdlIGByIHJvdW5kKHN1bURpZXQkdGVzdCRjb2VmZmljaWVudHNbMl0sZGlnaXRzPTEpYGcgYW5kIGByIHJvdW5kKHN1bURpZXQkdGVzdCRjb2VmZmljaWVudHNbM10sZGlnaXRzPTEpYGcgaGlnaGVyIHRoYW4gdGhhdCBmb3IgZmlzaCB0aGF0IHdlcmUgZmVkIHdpdGggZGlldCAxIGFuZCAyLCByZXNwZWN0aXZlbHkgKDk1JSBDSQpbYHIgcm91bmQoY29uZkRpZXQkY29uZmludFsyLC0xXSwxKWBdIGFuZCBbYHIgcm91bmQoY29uZkRpZXQkY29uZmludFszLC0xXSwxKWBdLCByZXNwZWN0aXZlbHkpLgpUaGUgd2VpZ2h0IGdhaW4gb2YgZmlzaCBpbiBkaWV0IDIgaXMgb24gYXZlcmFnZSBgciByb3VuZChzdW1EaWV0JHRlc3QkY29lZmZpY2llbnRzWzFdLGRpZ2l0cz0xKWBnIGhpZ2hlciB0aGFuIHRoYXQgZm9yIGZpc2ggdGhhdCB3ZXJlIGZlZCB3aXRoIGRpZXQgMSAoOTUlIENJCltgciByb3VuZChjb25mRGlldCRjb25maW50WzEsLTFdLDEpYF0pLgoKKipOb3RlLCBob3dldmVyLCB0aGF0IHRoZSBhc3N1bXB0aW9ucyBvZiB0ZXN0IGFyZSB2aW9sYXRlZCwqKiBpLmUuIGxhcmdlciAKdmFyaWFiaWxpdHkgaW4gRGlldCAyIGFuZC9vciBhZGRpdGlvbmFsIHRhbmsvYmF0Y2ggZWZmZWN0cyBpbiBEaWV0IDIuIEluZGVlZCwgdGhlIHRhbmtzIApzZWVtIHRvIGNsdXN0ZXIgcGVyIHR3byB0YW5rcywgd2hpY2ggbWF5IGluZGljYXRlIGFkZGl0aW9uYWwgYmF0Y2ggZWZmZWN0cy4gCioqV2UgdGhlcmVmb3JlIGhhdmUgdG8gYmUgdmVyeSBjYXJlZnVsIHdpdGggdGhlIGludGVycHJldGF0aW9uIG9mIHRoZSByZXN1bHRzLioqCg==