Advanced Statistics with R: Mixed Model – Exam

The exam consists in three datasets to analyze. Try to follow the subsequent indications to analyze the datasets, and then send back a .docx (or similar) file including:

?(all) the lines of the script you run in R;

?the output for each line;

?the information requested for each dataset, indicated as Q1, Q2… within each dataset

N.B.: remember to load the libraries lme4 and car (and other libraries if necessary), and remember to convert in factors the effects that are written as number but they should be treated as categorical

N.B.: if you prefer, you can directly modify this file for answering

1.Pines

The dataset “pines.txt” consists of the tree volume (Volume: ft3/acre) of pine trees included within each one of the 4 plots (Block) considered in the loblolly pine midrotation fertilization experiment of Gumpertz and Brownie (1993; Can J For Res 23: 625-639). The data, reported in the following Table 1, included the tree volume after two years of experiment (years 4 and 6 not considered). Each block received two treatments, Nitrogen (N, 4 levels) and Phosphorous (P, 3 levels), as reported in Table 1.  

After imported the dataset in R, let’s write a mixed model (M1) considering the Volume as response variable and the main effects (= no interactions) of N, P and Block; (Q1): which effects are fixed and which is random?

After run the summary of the model, let’s plot the residuals for checking the normality of distribution using a qqnorm or a qqPlot statistics.

Let run an ANOVA on fixed effects; (Q2): what effect(s) is/are significant?

Let’s write a second model, M2, also considering the interaction between N and P; run the summary, the ANOVA on fixed effects, and make a model comparison using the AIC statistics; (Q3): which is the model that fits better the data?

Let’s write a third model, M3, by including all the previous effects of M2 as fixed, run the summary, the ANOVA, and compare the ANOVA results with the ones of M2; (Q4): do you find similar significances for fixed effects in models M2 and M3 (for the effects that were accounted in both); what model would you like to choose for analyzing your data (you can also look at model fitting)?

2.Guppies

The dataset “guppies.txt” includes the experimental measures about the sexual behavior of 224 male guppies (Poecilia reticulata). The behaviors were measured for 10 min for each male within experimental tanks, in which each male was put in front of a receptive female. Behaviors included the “Display”, a visual exhibition of the male in front of the female, and a sneaky attempt to mate called “GT”, both included in the dataset and measured in sec. The dataset also included the dams and the sires of each male guppy, that were mated following the scheme in the picture: each sire was mated with two dams, and each dam had a number of offspring (the males, called “ID” in the dataset). This data structure allows to find out the heritability of the sexual behavior (the quote of the behavior that is transmitted to progeny).

After imported the dataset in R, let’s write a mixed model (M1) for the “Display” behavior including the sire and the dam as random effects. Take care about the data structure: (Q1) which is the upper level effect? And which is the nested effect? (N.B. in this dataset both the upper level and the nested effect should be indicated as random effect).

Then look at the summary and (Q2) calculate the heritability of Display as ratio of the sire (intercept) variance on the sum of all the variances for random effects (sire, dam and residual). N.B.: since there is a lot of residual variance, residuals of the model are partly without the range of the normal distribution and are correlated with observation.

Moving from M1, let’s write a model M2 without including the dam effect, and compare the two models using AIC statistics and the likelihood ratio test: (Q3) which is the best fitted model? Does the introduction of the dam effect significantly improve the model fitting?

Let’s run the same analyses but for “GT” behavior and compare the two models including the dam effect or not (M3 vs. M4); (Q4): does the introduction of the dam effect significantly improve the model fitting?

3.Growth of tpr foals

The dataset “growth_tpr.txt” includes the periodical measures of the body weight (weight) of 43 foals (ID, matr or name columns), 13 females and 30 males (sex effect) of the CAITPR breed. After arrived at experimental station at an average age of 14.1 days, the foals were fed two diets (diet: HP: high protein, 21 foals; LP: low protein, 22 foals) and measured 7 times, every 4-5 weeks (date_num). The experiment is aimed to look if a reduction of protein in diet (LP) is able to provide the same growth of the normal diet (HP), in both sexes. The following table summarizes the data structure:

Let’s write a mixed model (M1) including the body weight of the foals as response variable (y) and the fixed effects of diet, sex, date_num and the random effect of ID.

Run the summary and the ANOVA on fixed effects; (Q1): what effects are significant? Does the protein reduction affect the growth of the foals?

Let’s write a further mixed model (M2) by including the date_num also as random slope for each individual (ID), run the summary and the ANOVA on fixed effects; (Q2): do you find some differences?

Let’s compare M1 and M2 using model fitting statistics (AIC, likelihood ratio test); (Q3): what is the model that better fits the data?