Background
Researchers are studying the impact of protein sources and protein
levels in the diet on the weight of rats. They feed the rats with diets
of beef, cereal and pork and use a low and high protein level for each
diet type. The researchers can include 60 rats in the experiment. Prior
to the experiment, the rats were divided in 10 homogeneous groups of 6
rats based on characteristics such as initial weight, appetite, etc.
Within each group a rat is randomly assigned to a diet. The rats are
fed during a month and the weight gain in grams is recorded for each
rat.
The researchers want to assess the effect of the type of diet and the
protein level on the weight of the rats.
In contrast to the previous exercise, we perform the analysis
for all three diets in the dataset.
Experimental
design
There are three explanatory variables in the experiment: the
factor diet type with two levels (beef and cereal), factor protein level
with levels (low and high) and a group blocking factor with 10
levels.
There are 6 treatments: beef-high, cereal-high, pork-high,
beef-low, cereal-low, pork-low protein.
The rats are the experimental units (the unit to which a
treatment is applied): in this design, there is a randomisation
restriction: Within a block, a rat is randomly assigned to a
diet.
The rats are the observational units (the unit on which the
response is measured): The weight is weighted for each rat.
The weight gain is the response variable.
The experiment is a randomized complete block (RCB)
design
Load libraries
Data import
diet <- read.table("https://raw.githubusercontent.com/statOmics/PSLS21/data/dietRats.txt",
header=TRUE)
head(diet)
Tidy data
diet <- diet %>%
mutate(block = as.factor(block),
protSource = as.factor(protSource),
protLevel = as.factor(protLevel)) %>%
mutate(protLevel = fct_relevel(protLevel, "l"))
Data exploration
Boxplot of the weight gain against protein source, protein level
with coloring according to block
Lineplot of the weight gain against protein source, protein level
with coloring and grouping according to block
Interpret the plots
Multivariate linear
regression analysis
Model
specification
Based on the data exploration, propose a sensible regression model to
analyse the data.
Assumptions
Check the assumptions of the linear regression model
Hypothesis
testing
Use the summary
function to get an initial test for the
parameters in the model.
Interpretation of the
regression parameters
To facilitate the interpretation of the different parameters in our
regression model, we can make use of the VisualizeDesign
function of the ExploreModelMatrix
R package. The first
argument of this function is the name of our target dataset, the second
argument is a model formula, which in this case is specified as a
~
followed by the explanatory variables in our model.
library(ExploreModelMatrix)
ExploreModelMatrix::VisualizeDesign(..., ~ ...)$plotlist
After seeing this, again think about the meaning of the parameters in
our model.
Testing the overall
(omnibus) effect of diet
By comparing a model containing diet effects to a model that does not
have diet effects, using anova.
Assessing the
interaction effect between protein source and protein level
Assessing specific
contrasts
Imagine that we are interested in assessing if there is an effect
of
- protein source in the low protein diets
- cereal versus beef
- pork verus beef
- pork versus cereal
- protein source in high protein diets
- cereal versus beef
- pork verus beef
- pork versus cereal
- protein level (high versus low) for
- beef diets
- cereal diets
- pork diets
- If the effect of the protein level differs between
- beef and cereal
- beef and pork
- cereal and pork diets
These effects of interest are so-called contrasts,
i.e. linear combinations of the parameters.
Step 1: translate these research questions into parameters of the
model (or combinations of multiple parameters).
Step 2: Assess the significance of the contrasts using the
multcomp
package. The contrasts are given as input in the
form of symbolic descriptions to the linfct
argument of the
glht
function.
Conclusion
Formulate a conclusion for the different research hypotheses.
LS0tCnRpdGxlOiAiRXhlcmNpc2UgOC42OiBCbG9ja2luZyBvbiB0aGUgcmF0IGRpZXQgZGF0YXNldCAoMikiICAgCmF1dGhvcjogIkxpZXZlbiBDbGVtZW50IGFuZCBKZXJvZW4gR2lsaXMiCmRhdGU6ICJzdGF0T21pY3MsIEdoZW50IFVuaXZlcnNpdHkgKGh0dHBzOi8vc3RhdG9taWNzLmdpdGh1Yi5pbykiICAKb3V0cHV0OgogICAgaHRtbF9kb2N1bWVudDoKICAgICAgY29kZV9kb3dubG9hZDogdHJ1ZSAgICAKICAgICAgdGhlbWU6IGNvc21vCiAgICAgIHRvYzogdHJ1ZQogICAgICB0b2NfZmxvYXQ6IHRydWUKICAgICAgaGlnaGxpZ2h0OiB0YW5nbwogICAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKLS0tCgojIEJhY2tncm91bmQgCgpSZXNlYXJjaGVycyBhcmUgc3R1ZHlpbmcgdGhlIGltcGFjdCBvZiBwcm90ZWluIHNvdXJjZXMgYW5kIHByb3RlaW4gbGV2ZWxzIGluIAp0aGUgZGlldCBvbiB0aGUgd2VpZ2h0IG9mIHJhdHMuIFRoZXkgZmVlZCB0aGUgcmF0cyB3aXRoIGRpZXRzIG9mIGJlZWYsIGNlcmVhbCAKYW5kIHBvcmsgYW5kIHVzZSBhIGxvdyBhbmQgaGlnaCBwcm90ZWluIGxldmVsIGZvciBlYWNoIGRpZXQgdHlwZS4gClRoZSByZXNlYXJjaGVycyBjYW4gaW5jbHVkZSA2MCByYXRzIGluIHRoZSBleHBlcmltZW50LiBQcmlvciB0byB0aGUgZXhwZXJpbWVudCwKdGhlIHJhdHMgd2VyZSBkaXZpZGVkIGluIDEwIGhvbW9nZW5lb3VzIGdyb3VwcyBvZiA2IHJhdHMgYmFzZWQgb24gCmNoYXJhY3RlcmlzdGljcyBzdWNoIGFzIGluaXRpYWwgd2VpZ2h0LCBhcHBldGl0ZSwgZXRjLgoKV2l0aGluIGVhY2ggZ3JvdXAgYSByYXQgaXMgcmFuZG9tbHkgYXNzaWduZWQgdG8gYSBkaWV0LiBUaGUgcmF0cyBhcmUgZmVkIGR1cmluZyAKYSBtb250aCBhbmQgdGhlIHdlaWdodCBnYWluIGluIGdyYW1zIGlzIHJlY29yZGVkIGZvciBlYWNoIHJhdC4gCgpUaGUgcmVzZWFyY2hlcnMgd2FudCB0byBhc3Nlc3MgdGhlIGVmZmVjdCBvZiB0aGUgdHlwZSBvZiBkaWV0IGFuZCB0aGUgcHJvdGVpbiAKbGV2ZWwgb24gdGhlIHdlaWdodCBvZiB0aGUgcmF0cy4gCgoqKkluIGNvbnRyYXN0IHRvIHRoZSBwcmV2aW91cyBleGVyY2lzZSwgd2UgcGVyZm9ybSB0aGUgYW5hbHlzaXMgZm9yIGFsbCB0aHJlZSoqCioqZGlldHMgaW4gdGhlIGRhdGFzZXQuKioKCiMgRXhwZXJpbWVudGFsIGRlc2lnbiAKCi0gVGhlcmUgYXJlIHRocmVlIGV4cGxhbmF0b3J5IHZhcmlhYmxlcyBpbiB0aGUgZXhwZXJpbWVudDogdGhlIGZhY3RvciBkaWV0IHR5cGUKd2l0aCB0d28gbGV2ZWxzIChiZWVmIGFuZCBjZXJlYWwpLCBmYWN0b3IgcHJvdGVpbiBsZXZlbCB3aXRoIGxldmVscyAKKGxvdyBhbmQgaGlnaCkgYW5kIGEgZ3JvdXAgYmxvY2tpbmcgZmFjdG9yIHdpdGggMTAgbGV2ZWxzLgoKLSBUaGVyZSBhcmUgNiB0cmVhdG1lbnRzOiBiZWVmLWhpZ2gsIGNlcmVhbC1oaWdoLCBwb3JrLWhpZ2gsIGJlZWYtbG93LCAKY2VyZWFsLWxvdywgcG9yay1sb3cgcHJvdGVpbi4gCgotIFRoZSByYXRzIGFyZSB0aGUgZXhwZXJpbWVudGFsIHVuaXRzICh0aGUgdW5pdCB0byB3aGljaCBhIHRyZWF0bWVudCBpcyBhcHBsaWVkKTogaW4gdGhpcyBkZXNpZ24sIHRoZXJlIGlzIGEgcmFuZG9taXNhdGlvbiByZXN0cmljdGlvbjogV2l0aGluIGEgYmxvY2ssIGEgcmF0IGlzIHJhbmRvbWx5IGFzc2lnbmVkIHRvIGEgZGlldC4KCi0gVGhlIHJhdHMgYXJlIHRoZSBvYnNlcnZhdGlvbmFsIHVuaXRzICh0aGUgdW5pdCBvbiB3aGljaCB0aGUgcmVzcG9uc2UgaXMgbWVhc3VyZWQpOiBUaGUgd2VpZ2h0IGlzIHdlaWdodGVkIGZvciBlYWNoIHJhdC4KCi0gVGhlIHdlaWdodCBnYWluIGlzIHRoZSByZXNwb25zZSB2YXJpYWJsZS4gCgotIFRoZSBleHBlcmltZW50IGlzIGEgcmFuZG9taXplZCBjb21wbGV0ZSBibG9jayAoUkNCKSBkZXNpZ24KCkxvYWQgbGlicmFyaWVzCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKIyBEYXRhIGltcG9ydAoKYGBge3J9CmRpZXQgPC0gcmVhZC50YWJsZSgiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL3N0YXRPbWljcy9QU0xTMjEvZGF0YS9kaWV0UmF0cy50eHQiLAogICAgICAgICAgICAgICAgICAgaGVhZGVyPVRSVUUpCmhlYWQoZGlldCkKYGBgCgojIFRpZHkgZGF0YQoKYGBge3J9CmRpZXQgPC0gZGlldCAlPiUKICBtdXRhdGUoYmxvY2sgPSBhcy5mYWN0b3IoYmxvY2spLAogICAgICAgICBwcm90U291cmNlID0gYXMuZmFjdG9yKHByb3RTb3VyY2UpLAogICAgICAgICBwcm90TGV2ZWwgPSBhcy5mYWN0b3IocHJvdExldmVsKSkgJT4lCiAgbXV0YXRlKHByb3RMZXZlbCA9IGZjdF9yZWxldmVsKHByb3RMZXZlbCwgImwiKSkKYGBgCgojIERhdGEgZXhwbG9yYXRpb24KCi0gQm94cGxvdCBvZiB0aGUgd2VpZ2h0IGdhaW4gYWdhaW5zdCBwcm90ZWluIHNvdXJjZSwgcHJvdGVpbiBsZXZlbCB3aXRoIGNvbG9yaW5nIGFjY29yZGluZyB0byBibG9jawoKYGBge3J9CmBgYAoKLSBMaW5lcGxvdCBvZiB0aGUgd2VpZ2h0IGdhaW4gYWdhaW5zdCBwcm90ZWluIHNvdXJjZSwgcHJvdGVpbiBsZXZlbCB3aXRoIGNvbG9yaW5nIGFuZCBncm91cGluZyBhY2NvcmRpbmcgdG8gYmxvY2sKCmBgYHtyfQpgYGAKCi0gSW50ZXJwcmV0IHRoZSBwbG90cwoKIyBNdWx0aXZhcmlhdGUgbGluZWFyIHJlZ3Jlc3Npb24gYW5hbHlzaXMKCiMjIE1vZGVsIHNwZWNpZmljYXRpb24KCkJhc2VkICBvbiB0aGUgZGF0YSBleHBsb3JhdGlvbiwgcHJvcG9zZSBhIHNlbnNpYmxlIHJlZ3Jlc3Npb24gbW9kZWwgdG8gYW5hbHlzZQp0aGUgZGF0YS4KCiMjIEFzc3VtcHRpb25zIAoKQ2hlY2sgdGhlIGFzc3VtcHRpb25zIG9mIHRoZSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbAoKIyMgSHlwb3RoZXNpcyB0ZXN0aW5nCgpVc2UgdGhlIGBzdW1tYXJ5YCBmdW5jdGlvbiB0byBnZXQgYW4gaW5pdGlhbCB0ZXN0IGZvciB0aGUgcGFyYW1ldGVycyBpbiB0aGUgCm1vZGVsLgoKIyMgSW50ZXJwcmV0YXRpb24gb2YgdGhlIHJlZ3Jlc3Npb24gcGFyYW1ldGVycwoKVG8gZmFjaWxpdGF0ZSB0aGUgaW50ZXJwcmV0YXRpb24gb2YgdGhlIGRpZmZlcmVudCBwYXJhbWV0ZXJzIGluIG91ciByZWdyZXNzaW9uCm1vZGVsLCB3ZSBjYW4gbWFrZSB1c2Ugb2YgdGhlIGBWaXN1YWxpemVEZXNpZ25gIGZ1bmN0aW9uIG9mIHRoZSAKYEV4cGxvcmVNb2RlbE1hdHJpeGAgUiBwYWNrYWdlLiBUaGUgZmlyc3QgYXJndW1lbnQgb2YgdGhpcyBmdW5jdGlvbiBpcyB0aGUgbmFtZQpvZiBvdXIgdGFyZ2V0IGRhdGFzZXQsIHRoZSBzZWNvbmQgYXJndW1lbnQgaXMgYSBtb2RlbCBmb3JtdWxhLCB3aGljaCBpbiB0aGlzCmNhc2UgaXMgc3BlY2lmaWVkIGFzIGEgYH5gIGZvbGxvd2VkIGJ5IHRoZSBleHBsYW5hdG9yeSB2YXJpYWJsZXMgaW4gb3VyIG1vZGVsLgoKYGBge3IsIGV2YWwgPSBGQUxTRX0KbGlicmFyeShFeHBsb3JlTW9kZWxNYXRyaXgpCkV4cGxvcmVNb2RlbE1hdHJpeDo6VmlzdWFsaXplRGVzaWduKC4uLiwgfiAuLi4pJHBsb3RsaXN0CmBgYAoKQWZ0ZXIgc2VlaW5nIHRoaXMsIGFnYWluIHRoaW5rIGFib3V0IHRoZSBtZWFuaW5nIG9mIHRoZSBwYXJhbWV0ZXJzIGluIG91ciBtb2RlbC4KCiMjIFRlc3RpbmcgdGhlIG92ZXJhbGwgKG9tbmlidXMpIGVmZmVjdCBvZiBkaWV0CgpCeSBjb21wYXJpbmcgYSBtb2RlbCBjb250YWluaW5nIGRpZXQgZWZmZWN0cyB0byBhIG1vZGVsIHRoYXQgZG9lcyBub3QgaGF2ZQpkaWV0IGVmZmVjdHMsIHVzaW5nIGFub3ZhLgoKYGBge3J9CmBgYAoKIyMgQXNzZXNzaW5nIHRoZSBpbnRlcmFjdGlvbiBlZmZlY3QgYmV0d2VlbiBwcm90ZWluIHNvdXJjZSBhbmQgcHJvdGVpbiBsZXZlbAoKCmBgYHtyfQpgYGAKCgojIyBBc3Nlc3Npbmcgc3BlY2lmaWMgY29udHJhc3RzCgpJbWFnaW5lIHRoYXQgd2UgYXJlIGludGVyZXN0ZWQgaW4gYXNzZXNzaW5nIGlmIHRoZXJlIGlzIGFuIGVmZmVjdCBvZgogCjEuIHByb3RlaW4gc291cmNlIGluIHRoZSBsb3cgcHJvdGVpbiBkaWV0cyAKCS0gY2VyZWFsIHZlcnN1cyBiZWVmCgktIHBvcmsgdmVydXMgYmVlZgoJLSBwb3JrIHZlcnN1cyBjZXJlYWwKCgoyLiBwcm90ZWluIHNvdXJjZSBpbiBoaWdoIHByb3RlaW4gZGlldHMgCgktIGNlcmVhbCB2ZXJzdXMgYmVlZgoJLSBwb3JrIHZlcnVzIGJlZWYKCS0gcG9yayB2ZXJzdXMgY2VyZWFsCgoKMy4gcHJvdGVpbiBsZXZlbCAoaGlnaCB2ZXJzdXMgbG93KSBmb3IKCS0gYmVlZiBkaWV0cwoJLSBjZXJlYWwgZGlldHMKCS0gcG9yayBkaWV0cwoKCjQuIElmIHRoZSBlZmZlY3Qgb2YgdGhlIHByb3RlaW4gbGV2ZWwgZGlmZmVycyBiZXR3ZWVuCgktIGJlZWYgYW5kIGNlcmVhbAoJLSBiZWVmIGFuZCBwb3JrCgktIGNlcmVhbCBhbmQgcG9yayBkaWV0cwoKVGhlc2UgZWZmZWN0cyBvZiBpbnRlcmVzdCBhcmUgc28tY2FsbGVkIAoqKmNvbnRyYXN0cywgaS5lLiBsaW5lYXIgY29tYmluYXRpb25zIG9mIHRoZSBwYXJhbWV0ZXJzKiouCgpTdGVwIDE6IHRyYW5zbGF0ZSB0aGVzZSByZXNlYXJjaCBxdWVzdGlvbnMgaW50byBwYXJhbWV0ZXJzIG9mIHRoZSBtb2RlbCAob3IKY29tYmluYXRpb25zIG9mIG11bHRpcGxlIHBhcmFtZXRlcnMpLgoKU3RlcCAyOiBBc3Nlc3MgdGhlIHNpZ25pZmljYW5jZSBvZiB0aGUgY29udHJhc3RzIHVzaW5nIHRoZSBgbXVsdGNvbXBgIHBhY2thZ2UuIApUaGUgY29udHJhc3RzIGFyZSBnaXZlbiBhcyBpbnB1dCBpbiB0aGUgZm9ybSBvZiBzeW1ib2xpYwpkZXNjcmlwdGlvbnMgdG8gdGhlIGBsaW5mY3RgIGFyZ3VtZW50IG9mIHRoZSBgZ2xodGAgZnVuY3Rpb24uCgojIENvbmNsdXNpb24KCkZvcm11bGF0ZSBhIGNvbmNsdXNpb24gZm9yIHRoZSBkaWZmZXJlbnQgcmVzZWFyY2ggaHlwb3RoZXNlcy4KCg==