1 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 this exercise we will perform the data exploration using all diets, but, to keep the data analysis simple we will only assess the beef and cereal diets.

2 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 and observational units.

  • The weight gain is the response variable.

  • The experiment is a randomized complete block (RCB) design

Load libraries

library(tidyverse)

3 Data import

diet <- read.table("https://raw.githubusercontent.com/statOmics/PSLS21/data/dietRats.txt",
                   header=TRUE)
head(diet)

4 Tidy data

diet <- diet %>%
  mutate(block = as.factor(block),
         protSource = as.factor(protSource),
         protLevel = as.factor(protLevel)) %>%
  mutate(protLevel = fct_relevel(protLevel, "l"))

5 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

6 Data filtering

Filter the data to only use the beef and cereal diet.

7 Multivariate linear regression analysis

7.1 Model specification

Based on the data exploration, propose a sensible regression model to analyse the data.

7.2 Assumptions

Check the assumptions of the linear regression model

7.3 Hypothesis testing

Use the summary function to get an initial test for the parameters in the model.

7.4 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.

7.5 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.

7.6 Assessing the interaction effect between protein source and protein level

7.7 Assessing specific contrasts

Imagine that we are interested in assessing if there is an effect of

  1. protein source in the low protein diets

  2. protein source in high protein diets

  3. protein level for beef diets,

  4. protein level for cereal diets,

  5. if the effect of the protein level differs between beef and cereal

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.

8 Conclusion

Formulate a conclusion for the different research hypotheses.

LS0tCnRpdGxlOiAiRXhlcmNpc2UgOC41OiBCbG9ja2luZyBvbiB0aGUgcmF0IGRpZXQgZGF0YXNldCIgICAKYXV0aG9yOiAiTGlldmVuIENsZW1lbnQgYW5kIEplcm9lbiBHaWxpcyIKZGF0ZTogInN0YXRPbWljcywgR2hlbnQgVW5pdmVyc2l0eSAoaHR0cHM6Ly9zdGF0b21pY3MuZ2l0aHViLmlvKSIgIApvdXRwdXQ6CiAgICBodG1sX2RvY3VtZW50OgogICAgICBjb2RlX2Rvd25sb2FkOiB0cnVlICAgIAogICAgICB0aGVtZTogY29zbW8KICAgICAgdG9jOiB0cnVlCiAgICAgIHRvY19mbG9hdDogdHJ1ZQogICAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQotLS0KCiMgQmFja2dyb3VuZCAKClJlc2VhcmNoZXJzIGFyZSBzdHVkeWluZyB0aGUgaW1wYWN0IG9mIHByb3RlaW4gc291cmNlcyBhbmQgcHJvdGVpbiBsZXZlbHMgaW4gCnRoZSBkaWV0IG9uIHRoZSB3ZWlnaHQgb2YgcmF0cy4gVGhleSBmZWVkIHRoZSByYXRzIHdpdGggZGlldHMgb2YgYmVlZiwgY2VyZWFsIAphbmQgcG9yayBhbmQgdXNlIGEgbG93IGFuZCBoaWdoIHByb3RlaW4gbGV2ZWwgZm9yIGVhY2ggZGlldCB0eXBlLiAKVGhlIHJlc2VhcmNoZXJzIGNhbiBpbmNsdWRlIDYwIHJhdHMgaW4gdGhlIGV4cGVyaW1lbnQuIFByaW9yIHRvIHRoZSBleHBlcmltZW50LAp0aGUgcmF0cyB3ZXJlIGRpdmlkZWQgaW4gMTAgaG9tb2dlbmVvdXMgZ3JvdXBzIG9mIDYgcmF0cyBiYXNlZCBvbiAKY2hhcmFjdGVyaXN0aWNzIHN1Y2ggYXMgaW5pdGlhbCB3ZWlnaHQsIGFwcGV0aXRlLCBldGMuCgpXaXRoaW4gZWFjaCBncm91cCBhIHJhdCBpcyByYW5kb21seSBhc3NpZ25lZCB0byBhIGRpZXQuIFRoZSByYXRzIGFyZSBmZWQgZHVyaW5nIAphIG1vbnRoIGFuZCB0aGUgd2VpZ2h0IGdhaW4gaW4gZ3JhbXMgaXMgcmVjb3JkZWQgZm9yIGVhY2ggcmF0LiAKClRoZSByZXNlYXJjaGVycyB3YW50IHRvIGFzc2VzcyB0aGUgZWZmZWN0IG9mIHRoZSB0eXBlIG9mIGRpZXQgYW5kIHRoZSBwcm90ZWluIApsZXZlbCBvbiB0aGUgd2VpZ2h0IG9mIHRoZSByYXRzLiAKCkluIHRoaXMgZXhlcmNpc2Ugd2Ugd2lsbCBwZXJmb3JtIHRoZSBkYXRhIGV4cGxvcmF0aW9uIHVzaW5nIGFsbCBkaWV0cywgYnV0LCAKdG8ga2VlcCB0aGUgZGF0YSBhbmFseXNpcyBzaW1wbGUgd2Ugd2lsbCBvbmx5IGFzc2VzcyB0aGUgYmVlZiBhbmQgY2VyZWFsIGRpZXRzLgoKIyBFeHBlcmltZW50YWwgZGVzaWduIAoKLSBUaGVyZSBhcmUgdGhyZWUgZXhwbGFuYXRvcnkgdmFyaWFibGVzIGluIHRoZSBleHBlcmltZW50OiB0aGUgZmFjdG9yIGRpZXQgdHlwZQp3aXRoIHR3byBsZXZlbHMgKGJlZWYgYW5kIGNlcmVhbCksIGZhY3RvciBwcm90ZWluIGxldmVsIHdpdGggbGV2ZWxzIAoobG93IGFuZCBoaWdoKSBhbmQgYSBncm91cCBibG9ja2luZyBmYWN0b3Igd2l0aCAxMCBsZXZlbHMuCgotIFRoZXJlIGFyZSA2IHRyZWF0bWVudHM6IGJlZWYtaGlnaCwgY2VyZWFsLWhpZ2gsIHBvcmstaGlnaCwgYmVlZi1sb3csIApjZXJlYWwtbG93LCBwb3JrLWxvdyBwcm90ZWluLiAKCi0gVGhlIHJhdHMgYXJlIHRoZSBleHBlcmltZW50YWwgYW5kIG9ic2VydmF0aW9uYWwgdW5pdHMuIAoKLSBUaGUgd2VpZ2h0IGdhaW4gaXMgdGhlIHJlc3BvbnNlIHZhcmlhYmxlLiAKCi0gVGhlIGV4cGVyaW1lbnQgaXMgYSByYW5kb21pemVkIGNvbXBsZXRlIGJsb2NrIChSQ0IpIGRlc2lnbgoKTG9hZCBsaWJyYXJpZXMKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgojIERhdGEgaW1wb3J0CgpgYGB7cn0KZGlldCA8LSByZWFkLnRhYmxlKCJodHRwczovL3Jhdy5naXRodWJ1c2VyY29udGVudC5jb20vc3RhdE9taWNzL1BTTFMyMS9kYXRhL2RpZXRSYXRzLnR4dCIsCiAgICAgICAgICAgICAgICAgICBoZWFkZXI9VFJVRSkKaGVhZChkaWV0KQpgYGAKCiMgVGlkeSBkYXRhCgpgYGB7cn0KZGlldCA8LSBkaWV0ICU+JQogIG11dGF0ZShibG9jayA9IGFzLmZhY3RvcihibG9jayksCiAgICAgICAgIHByb3RTb3VyY2UgPSBhcy5mYWN0b3IocHJvdFNvdXJjZSksCiAgICAgICAgIHByb3RMZXZlbCA9IGFzLmZhY3Rvcihwcm90TGV2ZWwpKSAlPiUKICBtdXRhdGUocHJvdExldmVsID0gZmN0X3JlbGV2ZWwocHJvdExldmVsLCAibCIpKQpgYGAKCiMgRGF0YSBleHBsb3JhdGlvbgoKLSBCb3hwbG90IG9mIHRoZSB3ZWlnaHQgZ2FpbiBhZ2FpbnN0IHByb3RlaW4gc291cmNlLCBwcm90ZWluIGxldmVsIHdpdGggY29sb3JpbmcgYWNjb3JkaW5nIHRvIGJsb2NrCgpgYGB7cn0KCmBgYAoKLSBMaW5lcGxvdCBvZiB0aGUgd2VpZ2h0IGdhaW4gYWdhaW5zdCBwcm90ZWluIHNvdXJjZSwgcHJvdGVpbiBsZXZlbCB3aXRoIGNvbG9yaW5nIGFuZCBncm91cGluZyBhY2NvcmRpbmcgdG8gYmxvY2sKCmBgYHtyfQoKYGBgCgotIEludGVycHJldCB0aGUgcGxvdHMKCiMgRGF0YSBmaWx0ZXJpbmcKCkZpbHRlciB0aGUgZGF0YSB0byBvbmx5IHVzZSB0aGUgYmVlZiBhbmQgY2VyZWFsIGRpZXQuCgpgYGB7cn0KCmBgYAoKIyBNdWx0aXZhcmlhdGUgbGluZWFyIHJlZ3Jlc3Npb24gYW5hbHlzaXMKCiMjIE1vZGVsIHNwZWNpZmljYXRpb24KCkJhc2VkICBvbiB0aGUgZGF0YSBleHBsb3JhdGlvbiwgcHJvcG9zZSBhIHNlbnNpYmxlIHJlZ3Jlc3Npb24gbW9kZWwgdG8gYW5hbHlzZQp0aGUgZGF0YS4KCiMjIEFzc3VtcHRpb25zIAoKQ2hlY2sgdGhlIGFzc3VtcHRpb25zIG9mIHRoZSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbAoKIyMgSHlwb3RoZXNpcyB0ZXN0aW5nCgpVc2UgdGhlIGBzdW1tYXJ5YCBmdW5jdGlvbiB0byBnZXQgYW4gaW5pdGlhbCB0ZXN0IGZvciB0aGUgcGFyYW1ldGVycyBpbiB0aGUgCm1vZGVsLgoKIyMgSW50ZXJwcmV0YXRpb24gb2YgdGhlIHJlZ3Jlc3Npb24gcGFyYW1ldGVycwoKVG8gZmFjaWxpdGF0ZSB0aGUgaW50ZXJwcmV0YXRpb24gb2YgdGhlIGRpZmZlcmVudCBwYXJhbWV0ZXJzIGluIG91ciByZWdyZXNzaW9uCm1vZGVsLCB3ZSBjYW4gbWFrZSB1c2Ugb2YgdGhlIGBWaXN1YWxpemVEZXNpZ25gIGZ1bmN0aW9uIG9mIHRoZSAKYEV4cGxvcmVNb2RlbE1hdHJpeGAgUiBwYWNrYWdlLiBUaGUgZmlyc3QgYXJndW1lbnQgb2YgdGhpcyBmdW5jdGlvbiBpcyB0aGUgbmFtZQpvZiBvdXIgdGFyZ2V0IGRhdGFzZXQsIHRoZSBzZWNvbmQgYXJndW1lbnQgaXMgYSBtb2RlbCBmb3JtdWxhLCB3aGljaCBpbiB0aGlzCmNhc2UgaXMgc3BlY2lmaWVkIGFzIGEgYH5gIGZvbGxvd2VkIGJ5IHRoZSBleHBsYW5hdG9yeSB2YXJpYWJsZXMgaW4gb3VyIG1vZGVsLgoKYGBge3IsIGV2YWwgPSBGQUxTRX0KbGlicmFyeShFeHBsb3JlTW9kZWxNYXRyaXgpCkV4cGxvcmVNb2RlbE1hdHJpeDo6VmlzdWFsaXplRGVzaWduKC4uLiwgfiAuLi4pJHBsb3RsaXN0CmBgYAoKQWZ0ZXIgc2VlaW5nIHRoaXMsIGFnYWluIHRoaW5rIGFib3V0IHRoZSBtZWFuaW5nIG9mIHRoZSBwYXJhbWV0ZXJzIGluIG91ciBtb2RlbC4KCiMjIFRlc3RpbmcgdGhlIG92ZXJhbGwgKG9tbmlidXMpIGVmZmVjdCBvZiBkaWV0CgpCeSBjb21wYXJpbmcgYSBtb2RlbCBjb250YWluaW5nIGRpZXQgZWZmZWN0cyB0byBhIG1vZGVsIHRoYXQgZG9lcyBub3QgaGF2ZQpkaWV0IGVmZmVjdHMsIHVzaW5nIGFub3ZhLgoKYGBge3J9CgpgYGAKCiMjIEFzc2Vzc2luZyB0aGUgaW50ZXJhY3Rpb24gZWZmZWN0IGJldHdlZW4gcHJvdGVpbiBzb3VyY2UgYW5kIHByb3RlaW4gbGV2ZWwKCgpgYGB7cn0KCmBgYAoKCiMjIEFzc2Vzc2luZyBzcGVjaWZpYyBjb250cmFzdHMKCkltYWdpbmUgdGhhdCB3ZSBhcmUgaW50ZXJlc3RlZCBpbiBhc3Nlc3NpbmcgaWYgdGhlcmUgaXMgYW4gZWZmZWN0IG9mCiAKMS4gcHJvdGVpbiBzb3VyY2UgaW4gdGhlIGxvdyBwcm90ZWluIGRpZXRzCgoyLiBwcm90ZWluIHNvdXJjZSBpbiBoaWdoIHByb3RlaW4gZGlldHMgCgozLiBwcm90ZWluIGxldmVsIGZvciBiZWVmIGRpZXRzLCAKCjQuICBwcm90ZWluIGxldmVsIGZvciBjZXJlYWwgZGlldHMsIAoKNS4gaWYgdGhlIGVmZmVjdCBvZiB0aGUgcHJvdGVpbiBsZXZlbCBkaWZmZXJzIGJldHdlZW4gYmVlZiBhbmQgY2VyZWFsIAoKU3RlcCAxOiB0cmFuc2xhdGUgdGhlc2UgcmVzZWFyY2ggcXVlc3Rpb25zIGludG8gcGFyYW1ldGVycyBvZiB0aGUgbW9kZWwgKG9yCmNvbWJpbmF0aW9ucyBvZiBtdWx0aXBsZSBwYXJhbWV0ZXJzKS4KClN0ZXAgMjogQXNzZXNzIHRoZSBzaWduaWZpY2FuY2Ugb2YgdGhlIGNvbnRyYXN0cyB1c2luZyB0aGUgYG11bHRjb21wYCBwYWNrYWdlLiAKVGhlIGNvbnRyYXN0cyBhcmUgZ2l2ZW4gYXMgaW5wdXQgaW4gdGhlIGZvcm0gb2Ygc3ltYm9saWMKZGVzY3JpcHRpb25zIHRvIHRoZSBgbGluZmN0YCBhcmd1bWVudCBvZiB0aGUgYGdsaHRgIGZ1bmN0aW9uLgoKIyBDb25jbHVzaW9uCgpGb3JtdWxhdGUgYSBjb25jbHVzaW9uIGZvciB0aGUgZGlmZmVyZW50IHJlc2VhcmNoIGh5cG90aGVzZXMuCgoK