16 Complex models

16.1 Designing a Model

We are introduced to the fruitfly dataset Partridge and Farquhar (1981)⁷. From our understanding of sexual selection and reproductive biology in fruit flies, we know there is a well established 'cost' to reproduction in terms of reduced longevity for female fruitflies. The data from this experiment is designed to test whether increased sexual activity affects the lifespan of male fruitflies.

The flies used were an outbred stock, sexual activity was manipulated by supplying males with either new virgin females each day, previously mated females ( Inseminated, so remating rates are lower), or provide no females at all (Control). All groups were otherwise treated identically.

• type: type of female companion (virgin, inseminated, control(partners = 0))

• longevity: lifespan in days

• thorax: length of thorax in micrometres (a proxy for body size)

• sleep: percentage of the day spent sleeping

16.2 Hypothesis

Before you start any formal analysis you should think clearly about the sensible parameters to test. In this example, we are most interested in the effect of sexual activity on longevity. But it is possible that other factors may also affect longevity and we should include these in our model as well, and we should think hard about what terms might reasonably be expected to interact with sexual activity to affect longevity.

Task

Which terms and interactions do you think we should include in our model?

Solution

In this exercise I have just asked you to try and think logically about suitable predictors. For a more formal investigation you should support this with evidence where possible

• type - should definitely be included.

• thorax - the size of the flies could determine longevity. Carreira et al (2009)⁸

• sleep - sleep could easily help determine longevity. Thompson et al (2020)⁹

• type:sleep - the amount that sleep (rest) helps promote longevity could change depending on how much activity the fly engages in when awake. Chen et al (2017)¹⁰

Other interactions could be included but you should have a strong reason for them.

16.3 Checking the data

You should now import, clean and tidy your data. Making sure it is in tidy format, all variables have useful names, and there are no mistakes, missing data or typos.

Based on the variables you have decided to test you should start with some simple visualisations, to understand the distribution of your data, and investigate visually the relationships you wish to test.

This is a full two-by-two plot of the entire dataset, but you should try and follow this up with some specific plots.

GGally::ggpairs(fruitfly)

16.4 Activity 1: Think about your data

Task

Think carefully about the plots you should make to investigate the potential differences and relationships you wish to investigate - examples hidden behing this dropdown

Solution

In this first figure - we can investigate whether there is an obvious difference in the longevities of males across the three treatments

colours <- c("cyan", "darkorange", "purple")

fruitfly %>% 
  ggplot(aes(x = longevity, y = type, fill = type))+
  geom_density_ridges(alpha = 0.5)+
  scale_fill_manual(values = colours)+
  theme_minimal()+
  theme(legend.position = "none")

Figure 16.1: A density distribution of longevity across the three sexual activity treatments

With careful thought we can construct figures that help us investigate the distribution of multiple potential predictors at once. Here we can overlay body size onto the three treatments.

Q Does it like size affects longevity? Yes No

Q Does it look like size affects longevity differently between treatment groups? Yes No

Here it does look as though larger flies have a longer lifespan than smaller flies. But there appears to be little difference in the distribution of body sizes by treatment. This does not mean we can't test this in our model. But we should have a clear idea of which will have largest probable effect sizes

fruitfly %>% 
  ggplot(aes(x=type, y = longevity, fill = type))+
  geom_boxplot(alpha = 0.4)+
  geom_jitter(aes(size = thorax),
              alpha = 0.6,
              width = 0.4,
              shape = 21)+
  scale_fill_manual(values = colours)+
    guides(fill = "none",
         colour = "none")+
  theme_minimal()

Figure 16.2: A boxplot of longevity across three treatments of sexual activity. Individual points represent individual males, where the size of each point represents thorax length.

We are also interested in the potential effect of sleep on activity, we can construct a scatter plot of sleep against longevity, while including body size as a covariate.

In these faceted plots - Are the trendlines moving in the same direction? Yes No

Investigating these can help us determine is there much evidence of a potential effect of sleep, is there evidence for whether this might be an additive effect or one which interacts with the the three treatments?

fruitfly %>% 
ggplot(aes(x=sleep, y = longevity, fill = type))+
  geom_point(aes(size = thorax,
                 fill = type),
             shape = 21,
             alpha = 0.4)+
  geom_smooth(method = "lm", colour = "black",
              se = FALSE)+
  scale_fill_manual(values = colours)+
  guides(fill = "none",
         colour = "none")+
  facet_wrap(~ type)+
  theme_minimal()

Figure 3.3: A scatter plot of proportion of time spent sleeping against longevity with a linear model trendline. Points represent individual flies, and the size of each point represents thorax length (mm).

16.5 Designing a model

When you include an interaction term, the numbers produced from this are how much more or less the mean estimate is than if you just combined the main effects.

# a full model
flyls1 <- lm(longevity ~ type + thorax + sleep + type:sleep, data = fruitfly)

flyls1 %>% 
  broom::tidy()

term	estimate	std.error	statistic	p.value
(Intercept)	-57.5275383	11.3554560	-5.0660703	0.0000015
typeInseminated	7.9883828	5.3412012	1.4956154	0.1374236
typeVirgin	-10.9075381	5.4745755	-1.9923989	0.0486358
thorax	142.5090010	13.4115350	10.6258531	0.0000000
sleep	0.0904459	0.1885893	0.4795919	0.6324053
typeInseminated:sleep	-0.1965054	0.2082301	-0.9436937	0.3472544
typeVirgin:sleep	-0.1124276	0.2166543	-0.5189260	0.6047842

Because we have included an interaction effect the number of terms is quite long and takes more consideration to understand. We can see for the individual estimates that it does not appear that the interaction is having a strong effect (estimate) and this does not appear to be different from a null hypothesis of no interaction effect. But we we should use an F test to look at the overall effect to be sure.

Task

From the model summary table could you say what the mean longevity of a male with a 0.79mm thorax, that sleeps for 22% of the day and is paired with virgin females would be?

Solution

# intercept
coef(flyls1)[1] + 
  
# 1*coefficient for virgin treatment  
coef(flyls1)[3] + 
  
# 0.79 * coefficient for thorax size  
(coef(flyls1)[4]*0.79) + 
  
# 22 * coefficient for sleep  
(coef(flyls1)[5]*22) + 

# 22 * 1 * coefficient for interaction
(coef(flyls1)[7]*22*1)

## typeVirgin:sleep 
##        -2.473406

16.6 Model checking & collinearity

Before we start playing with the terms in our model, we should check to see if this is even a good way of fitting and measuring our data. We should check the assumptions of our model are being met.

performance::check_model(flyls1)

16.7 Activity 2: Model checking

Question - IS the assumption of homogeneity of variance met? Yes No

Solution

• Mostly - the reference line is fairly flat (there is a slight curve).

• It looks as though there might be some increasing heterogeneity with larger values, though very minor.

VERDICT, pretty much ok, should be fine for making inferences.

With a slight curvature this could indicate that you might get a better fit with a transformation, or perhaps that there is a missing variable that if included in the model would improve the residuals. In this instance I wouldn't be overly concerned. See here for a great explainer on intepreting residuals¹¹.

Question - ARE the residuals normally distributed? Yes No

Solution

Yes - the QQplot looks pretty good, a very minor indication of a right skew, but nothing to worry about.

Interpreting QQ plots

Question - IS their an issue with Collinearity? Yes No

Solution

This graph clearly shows there is collinearity. But this is not unusual when we include an interaction term, if we see evidence of collinearity in terms that are not part of an interaction then we should take another look¹².

What can you do about collinearity in main effects? 1) Nothing 2) Transform 3) Drop one of the terms.

The check_performance() function produces a visual summary of a Variance Inflation Factor produced from the vif() function. This is a measure of the standard error of each estimated coefficient. If this is very larger (greater than 5 or 10), this indicates the model has problems estimating the coefficient. This does not affect model predictions, but makes it more difficult to determine the estimate change from a predictor.

car::vif(flyls1)

##                 GVIF Df GVIF^(1/(2*Df))
## type       12.478906  2        1.879508
## thorax      1.052967  1        1.026142
## sleep       8.750764  1        2.958169
## type:sleep 38.749001  2        2.494969

16.8 Data transformations

The most common issues when trying to fit simple linear regression models is that our response variable is not normal which violates our modelling assumption. There are two things we can do in this case:

• Variable transformation e.g lm(sqrt(x) ~ y, data = data)

  • Can sometimes fix linearity

  • Can sometimes fix non-normality and heteroscedasticity (i.e non-constant variance)

• Generalized Linear Models (GLMs) to change the error structure (i.e the assumption that residuals need to be normal - see next week.)

16.8.1 BoxCox

The BoxCox gets its name from its two inventors, George Box and David Cox. Implemented by the MASS package, when applied to a linear model it sytematically applies transformations by raising the y variable to a power (lambda).

The R output for the MASS::boxcox() function plots a maximum likelihood curve (with a 95% confidence interval - drops down as dotted lines) for the best transformation for fitting the data to the model.

Table 16.1: Common Box-Cox Transformations

lambda value	transformation
0.0	log(Y)
0.5	sqrt(Y)
1.0	Y
2.0	Y^1

# run this, pick a transformation and retest the model fit
MASS::boxcox(flyls1)

Figure 16.3: standard curve fitted by maximum likelihood, dashed lines represent the 95% confidence interval range for picking the 'best' transformation for the dependent variable

Question - Does the fit of the model improve with a square root transformation? Yes No

Solution

Not really, despite the suggestion that a sqrt transformation would improve the model, residual fits are not really any better - so we might as well stick with the original scale.

16.9 Model selection

# use drop1 function to remove top-level terms
drop1(flyls1, test = "F")

	Df	Sum of Sq	RSS	AIC	F value	Pr(>F)
<none>	NA	NA	14994.43	612.3900	NA	NA
thorax	1	14347.4733	29341.90	694.3073	112.9087541	0.0000000
type:sleep	2	130.1431	15124.57	609.4702	0.5120865	0.6005695

Based on this ANOVA table, we do not appear to have a strong rationale for keeping the interaction term in the model (AIC or F-test). Therefore we can confidently remove the interaction, simplifying our model and making interpretation easier.

flyls2 <- lm(longevity ~ type + thorax + sleep, data = fruitfly)

drop1(flyls2, test = "F")

	Df	Sum of Sq	RSS	AIC	F value	Pr(>F)
<none>	NA	NA	15124.57	609.4702	NA	NA
type	2	7576.86233	22701.43	656.2337	30.057833	0.0000000
thorax	1	15282.82102	30407.39	694.7659	121.255596	0.0000000
sleep	1	86.27949	15210.85	608.1813	0.684551	0.4096663

Question - Should we drop sleep from this model? Yes No

Solution

There is good reason to remove non-significant interaction terms from a model, they complicate estimates and make interpretations more difficult. For main effects things are a little more ambiguous.

When the main aim is prediction, it makes sense to be cautious and retain non-significant terms, as extra terms make no difference to the R^2 of a model.

When the focus is on hypothesis testing, then removal of non-significant terms can help produce a 'true' model, but this is optional. Generally speaking it is often simpler to leave main effects in the model (you should have carefully considered the terms which were included in the first place).

In this example we can also see that AIC has not really changed - so the quality of the model is also not improved vby dropping this term.

16.10 Posthoc

Using the emmeans package is a very easy way to produce the estimate mean values (rather than mean differences) for different categories emmeans. If the term pairwise is included then it will also include post-hoc pairwise comparisons between all levels with a tukey test contrasts.

emmeans::emmeans(flyls2, specs = pairwise ~ type + thorax + sleep)

## $emmeans
##  type        thorax sleep emmean   SE  df lower.CL upper.CL
##  Control      0.821  23.5   61.3 2.26 120     56.8     65.8
##  Inseminated  0.821  23.5   64.9 1.59 120     61.8     68.1
##  Virgin       0.821  23.5   48.0 1.59 120     44.9     51.2
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                            estimate   SE  df t.ratio
##  Control 0.82096 23.464 - Inseminated 0.82096 23.464    -3.63 2.77 120 -1.309 
##  Control 0.82096 23.464 - Virgin 0.82096 23.464         13.25 2.76 120  4.796 
##  Inseminated 0.82096 23.464 - Virgin 0.82096 23.464     16.87 2.25 120  7.508 
##  p.value
##  0.3929 
##  <.0001 
##  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

For continuous variables (sleep and thorax) - emmeans has set these to the mean value within the dataset, so comparisons are constant between categories at the average value of all continuous variables.

16.11 Activity 3: Write-up

Task

Can you write a Results section?

Solution

I tested the hypothesis that sexual activity is costly for male Drosophila melanogaster fruitflies. Previous research indicated that sleep deprived males are less attractive to females, this would indicate that levels of sexual activity might be affected by sleep and impact the effect on longevity, as such this was included as an interaction term in the full model. Body size is also know to affect lifespan, as such this was included as a covariate in the mode.

There was a small interaction effect of decreased lifespan with increasing sleep in the treatment groups compared to control in our samples, but this was not significantly different from no effect (F_2,118 = 0.512, P = 0.6), and was therefore dropped from the full model (Table 15.1).

library(kableExtra)
flyls2 %>% broom::tidy(conf.int = T) %>% 
 select(-`std.error`) %>% 
mutate_if(is.numeric, round, 2) %>% 
kbl(col.names = c("Predictors",
                    "Estimates",
                    "Z-value",
                    "P",
                    "Lower 95% CI",
                    "Upper 95% CI"),
      caption = "Linear model coefficients", 
    booktabs = T) %>% 
   kable_styling(full_width = FALSE, font_size=16)

Table 16.2: Linear model coefficients

Predictors	Estimates	Z-value	P	Lower 95% CI	Upper 95% CI
(Intercept)	-56.05	-5.01	0.00	-78.18	-33.91
typeInseminated	3.63	1.31	0.19	-1.86	9.11
typeVirgin	-13.25	-4.80	0.00	-18.71	-7.78
thorax	144.43	11.01	0.00	118.46	170.40
sleep	-0.05	-0.83	0.41	-0.18	0.07

There was a significant overall effect of treatment on male longevity (Linear model: F_2,120 = 30.1, P < 0.001), with males paired to virgin females having the lowest mean longevity (48 days, [95%CI: 44.9 - 51.2]) (when holding body size and sleep constant), compared to control males (61.3 days [56.8 - 65.8]) and males paired with inseminated females (64.9 days [61.8 - 68.1 days]).

Post hoc analysis showed that these differences were statistically significant for males paired with control females compared to the inseminated (Tukey test: t₁₂₀ = 4.8, P < 0.001) and virgin groups (t₁₂₀ = 7.5, P < 0.001), but there was no overall evidence of a difference between inseminated and virgin groups (t₁₂₀ = -1.309 P < 0.3929) (Figure 19.4).

Comparing the treatment effects against other predictors of longevity such as body size and sleep, I found that sleep had a very small effect on longevity (mean change -0.05 days [-0.18 - 0.07]) which was not significantly different from no effect (Linear model: F_1,120 = 0.68, P = 0.41). Body size (taken from thorax length) was a significant predictor of longevity (F_1,120 = 121, P < 0.001), with each 0.1 mm increase in body size adding 14.4 days to the individual lifespan [11.8 - 17]. It appears as though body size has a stronger effect on longevity than treatment, indicating that while there is a measurable cost of sexual activity to males, it may be less severe than in females (not compared here), and less severe than other measurable predictors.

# optional extras for adding animal silhouette
library(rphylopic)
# search for the right image in the database
fruitfly_phylo <- name_search(text = "Drosophila melanogaster", options = "namebankID")[[1]] 

# name_images(uuid = fruitfly$uid[1])

# extract the ID of thr right image
fruitfly_2 <- name_images(uuid = fruitfly_phylo$uid[1])$same[[2]]$ui

# get the actual image
fruitfly_pic <- image_data(fruitfly_2, size = 256)[[1]]

# density plot of thorax length by treatment
marginal1 <- fruitfly %>% 
  ggplot()+
  geom_density(aes(x = thorax, fill = type),
               alpha = 0.5)+
  scale_fill_manual(values = colours)+
  theme_void()+
  theme(legend.position = "none")

# density plot of longevity by treatment
marginal2 <- fruitfly %>% 
  ggplot()+
  geom_density(aes(x = longevity, fill = type),
               alpha = 0.5)+
  scale_fill_manual(values = colours)+
  theme_void()+
  coord_flip()


# create a new dataset where all sleep values are set to a constant (mean) values - see the emmeans table

fruitfly_2 <- fruitfly %>% 
  mutate(sleep = mean(sleep))

# use the final model to produce model predictions set to the new constant sleep dataframe
model_plot <- broom::augment(flyls2, newdata = fruitfly_2, interval = "confidence") %>% 
  ggplot(aes(x=thorax, y = .fitted, colour = type, fill = type))+
  geom_ribbon(aes(ymin=.lower, ymax=.upper), alpha = 0.2)+
    geom_line(linetype = "dashed", show.legend = FALSE)+
  geom_point(data = fruitfly, aes(x = thorax, y = longevity),
             show.legend = FALSE)+
  scale_colour_manual(values = colours)+
  scale_fill_manual(values = colours)+
  labs(y = "Lifespan in days",
       x = "Thorax length (mm)",
       fill = "Type of female exposure")+
  guides(colour = "none")+
  theme_classic()+
  theme(legend.position = "none")+
  add_phylopic(fruitfly_pic, alpha = 0.8, x = 0.7, y = 80, ysize = 10)



layout <- "
AAA#
BBBC
BBBC
BBBC"


marginal1+model_plot+marginal2 +plot_layout(design = layout)

Figure 9.2: A scatter plot of longevity against body size across three treatments of differening male sexual activity. Fitted model slopes are from the reduced linear model (main effects only of thorax size, sleep and treatment group), with 95% confidence intervals, circles are individual data points. Marginal plots are density plot distributions for thorax length and longevity split by treatments.

16.12 Summary

In this chapter we have worked with our scientific knowledge to develop testable hypotheses and built statistical models to formally assess them. We now have a working pipeline for tackling complex datasets, developing insights and producing and explaining robust linear models.

16.12.1 Checklist

• Think carefully about the hypotheses to test, use your scientific knowledge and background reading to support this

• Import, clean and understand your dataset: use data visuals to investigate trends and determine if there is clear support for your hypotheses

• Fit a linear model, including interaction terms with caution

• Investigate the fit of your model, understand that parameters may never be perfect, but that classic patterns in residuals may indicate a poorly fitting model - sometimes this can be fixed with careful consideration of missing variables or through data transformation

• Test the removal of any interaction terms from a model, look at AIC and significance tests

• Make sure you understand the output of a model summary, sense check this against the graphs you have made

• The direction and size of any effects are the priority - produce estimates and uncertainties. Make sure the observations are clear.

• Write-up your significance test results, taking care to report not just significance (and all required parts of a significance test). Do you know what to report? Within a complex model - reporting t will indicate the slope of the line for that single term against the intercept, F is the overall effect of a predictor across all levels, post-hoc if you wish to compare across all levels.

• Well described tables and figures can enhance your results sections - take the time to make sure these are informative and attractive.

16.13 Supplementary code

sjPlot A really nice package that helps produce model summaries for you automatically

library(sjPlot)
tab_model(flyls2)

	longevity
Predictors	Estimates	CI	p
(Intercept)	-56.05	-78.18 – -33.91	<0.001
type [Inseminated]	3.63	-1.86 – 9.11	0.193
type [Virgin]	-13.25	-18.71 – -7.78	<0.001
thorax	144.43	118.46 – 170.40	<0.001
sleep	-0.05	-0.18 – 0.07	0.410
Observations	125
R2 / R2 adjusted	0.605 / 0.591

library(gtsummary)
tbl_regression(flyls2)

Characteristic	Beta	95% CI1	p-value
type
Control	—	—
Inseminated	3.6	-1.9, 9.1	0.2
Virgin	-13	-19, -7.8	<0.001
thorax	144	118, 170	<0.001
sleep	-0.05	-0.18, 0.07	0.4
1 CI = Confidence Interval