The aim of these exercises is to …
For each of the examples below, you should follow the sequence we’ve used previously:
What is the biological question?
Is the predictor continuous or categorical?
Write out the linear model corresponding to this question.
What distribution do you expect the response variable to follow?
What are the assumptions behind the statistical model you’ll fit?
Fit the model
How will you assess whether the model fits well?
Can you detect an effect of the predictor?
How do you measure the effect?
What do you conclude (including any cautions)
Li et al. (2015) investigated ways of revegetating damaged freshwater enviroments with vegetation. They used a suite of 4 plant species, and focused on the combined effects of propagule pressure (i.e., the initial seeding event) and water depth on the aquatic vegetation that developed. They reported data for each plant species separately and considered overall plant “community”, putting all species together.
The four plant species can propagate clonally, by shoot regeneration, and they used shoot fragments to “seed” experimental pots. Each pot had 1, 2, or 4 fragments of each of the four species (i.e. there were three propagule pressures). The pots, filled with natural lake sediment, were placed in larger plastic tanks, which were filled to a depth of 30 or 70 cm. There were 5 tanks for each depth. The tanks were outside in full sun, and plants were allowed to grow for 7 weeks, after which time plants were harvested and measured.
For each species in a pot, five shoots were selected and their length measured, the number of nodes (indicating clonal growth) counted, and dry biomass calculated. The total biomass for that species was also measured and used to estimate total nodes and total shoot length. Values for each species were then used to calulate total biomass, total shoot length, and total number of nodes for that pot.
Their data can be found in Table S1 of the paper, but we’ve extracted the data and renamed some of the columns to make things clearer. Our file is here.
df <-read.csv("../data/li.csv")
set_sum_contrasts()## setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))head(df, 10)## depth tank pressure biomass nodes slength hvbiomass
## 1 30 1 low 3.204 1161.9653 14.095498 2.606
## 2 30 1 medium 3.523 1456.9005 16.637120 2.496
## 3 30 1 high 7.077 1679.0204 19.736082 5.690
## 4 30 2 low 3.969 1257.9132 15.239968 2.789
## 5 30 2 medium 4.470 1402.1528 20.200615 3.609
## 6 30 2 high 8.041 2172.3806 28.823652 5.935
## 7 30 3 low 0.530 371.1823 4.193746 0.339
## 8 30 3 medium 1.473 756.2558 8.357996 1.019
## 9 30 3 high 1.719 1191.0750 15.192225 1.532
## 10 30 4 low 1.879 814.2920 10.042774 1.521df$depth<-as.factor(df$depth)
df$tank<-as.factor(df$tank)Remember to go through the steps at the start of these exercises to make sure you have identified the design and assumptions clearly
Seitz, Lipcius, and Hines (2016) studied the estuarine clam Macoma balthica and set up a field experiment to examine the effects of two different sediment types (shallow-mud with high food availability and muddy-sand with low food availability) and predation by fish and crabs on the growth rate of clams in upper Chesapeake Bay. They used a split-plot design where there were eight shallow-mud and four muddy-sand sites, with more shallow-mud sites because of greater variability; sediment type was the between-plots fixed factor with random sites (“plots”) nested in each type. Within each site, there was one sub-plot (0.25m2) that excluded predators with a cage and a second sub-plot that had no cage; predation was the within-plots fixed factor. Ten individually marked and measured clams were transplanted into each sub-plot and retrieved 20-22 days later to record growth.
The data are available from dryad. The data for this example are in GrowSplitPlot.xlsx, but it needs quite a bit of work for us to work with it in R. Our tidied version is here
Load the data and have a quick look. The variable names should be clear
df <- read.csv("../data/seitz.csv")
df$Plot <-as.factor(df$Plot)
head(df,10)## Habitat Plot Growth Predation
## 1 Mud 1 0.7 Exclusion
## 2 Mud 1 0.7 Exclusion
## 3 Mud 1 0.2 Exclusion
## 4 Mud 1 1.1 Exclusion
## 5 Mud 1 0.9 Exclusion
## 6 Mud 1 0.3 Exclusion
## 7 Mud 1 0.0 Exclusion
## 8 Mud 2 0.1 Exclusion
## 9 Mud 2 0.6 Exclusion
## 10 Mud 2 0.1 ExclusionExamine this effect using a linear mixed model approach and REML
If we follow the example of Box 11.3, there is some variation in opinion as to whether we’d try a slightly simpler model, omitting the Predation x Plot interaction, which doesn’t account for much variation.
What do you think of this option, and would it make a difference?
Simpler is better when explaining your analysis to an audience!
Dalziel et al. (2015) were interested in disentangling the processes of acclimation and adaptation to changing environments, and examined levels of phenotypic plasticity in a small fish, the lake whitefish. This fish has developed two ecotypes - phenotypes that feed in different ways. Dwarf ecotypes are actively swimming and feeding, but slow-growing, while the ancestral normal phenotype is more sedentary. It’s thought that these ecotypes evolved separately under past climates, but then came together and can co-occur. There are anatomical and physiological differences between the ecotypes that are linked to swimming ability, and Dalziel and her colleagues were interested in how these ecotypes responed to a changed environment with faster-flowing water, particularly any role of phenotypic plasticity.
Phenotypic plasticity is best assessed by placing different phenotypes into the same set of varied environments and measuring their performance. They did this for whitefish by using roughly size-matched, laboratory-fertilized and reared fish from each ecotype, which were then placed in two environments, still water and fast-flowing water. Fish were kept in these environments for approximately 3 months, then sampled.
Flow environments were produced in 1 m high circular tanks, .6 m in diameter, with a .16 m diameter cylinder in the centre of the tank, creating a circular swimming “race track”. Within this track, water flowed at 7.6 cm/s for 6 h/day or 0.5 cm/s. Each tank contained 8 fish of each ecotype.
At the end of the experiment, fish were sacrficed, and the researchers took hematocrit samples, measured heart morphology and took samples of red and white muscle. Red mussel samples were used for enzyme assays and histology, and white muscle to obtain mitochondria, whose performance was assessed.
They provided data and R code through dryad. We’ll focus on their data focusing on mitochondrial function in white muscle. These data are in their file Fig5_MitoRespiration.csv. If you want to explore some of their other analyses, you can use their R script, which allows easy selection of data files.
plasticity<-read.csv("../data/Fig5_MitoRespiration.csv",header=TRUE,stringsAsFactors = FALSE, na.strings = c("NA",""))
head(plasticity, 10)## indv ecotype treatment group tank mass length mitoprotein.ml
## 1 WP1 dwarf swim DS 1 3.2197 7.580 6.25
## 2 WP2 dwarf control DC 2 4.1880 7.959 5.26
## 3 WP3 normal swim NS 3 3.2915 6.501 4.49
## 4 WP4 normal control NC 4 4.6945 7.840 3.40
## 5 WP5 normal control NC 5 4.2127 7.582 3.49
## 6 WP6 normal swim NS 6 3.7213 7.133 4.16
## 7 WP7 normal control NC 7 3.9910 7.766 5.28
## 8 WP8 dwarf swim DS 8 5.8855 9.091 4.84
## 9 WP9 normal swim NS 8 4.5422 7.742 2.92
## 10 WP10 normal control NC 7 4.8839 NA 4.22
## mitoprotein.gWM state2.ml state2.mg state3.ml state3.mg state4.ml state4.mg
## 1 0.622 74.0 11.8 404.2 64.7 137.8 22.0
## 2 0.526 42.3 8.1 344.9 65.6 128.7 24.5
## 3 0.458 24.3 5.4 261.3 58.2 122.5 27.3
## 4 0.260 33.5 9.9 331.9 97.6 90.0 26.5
## 5 0.347 28.3 8.1 266.1 76.2 63.7 18.3
## 6 0.435 NA NA NA NA NA NA
## 7 0.541 27.7 5.2 249.2 47.2 77.4 14.7
## 8 0.484 53.4 11.0 548.8 113.4 69.0 14.2
## 9 0.299 34.9 12.0 248.3 85.0 61.8 21.2
## 10 0.427 39.3 9.3 306.7 72.7 59.1 14.0
## RCR ACR flux1.4.ml flux1.4.mg flux2.4.ml flux2.4.mg flux.cox.ml flux.cox.mg
## 1 2.9 5.5 350.5 56.1 89.8 14.4 955.9 152.9
## 2 2.7 8.1 404.8 77.0 121.2 23.0 388.4 73.8
## 3 2.1 10.8 286.4 63.8 98.0 21.8 264.0 58.8
## 4 3.7 9.9 405.9 119.4 96.0 28.2 750.7 220.8
## 5 4.2 9.4 310.8 89.1 97.3 27.9 608.3 174.3
## 6 NA NA NA NA NA NA NA NA
## 7 3.2 9.0 388.8 73.6 47.9 9.1 671.8 127.2
## 8 8.0 10.3 675.8 139.6 171.1 35.4 1143.5 236.3
## 9 4.0 7.1 281.3 96.3 59.7 20.5 687.1 235.3
## 10 5.2 7.8 374.8 88.8 94.9 22.5 748.6 177.4
## notes1
## 1 mitoprotein.ml=mitochondrial protein per mL mito resuspension (from BSA assay)
## 2 mitoprotein.gWM1= mitochondrial protein per gram tissue based upon calculations from uL buffer added
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## 7 <NA>
## 8 <NA>
## 9 <NA>
## 10 <NA>Look at data measuring fluxes of four oxidative complexes, I-IV. These are the variables flux1.4.mg, flux2.4mg, and flux.cox.mg
flux1.4.mg is total flux across the four complexes, flux2.4.mg is flux across complexes II-IV, and flux.cox.mg is flux across complex IV
#Make sure categorical predictors are treated the right way
plasticity$treatment <-as.factor(plasticity$treatment)
plasticity$ecotype <-as.factor(plasticity$ecotype)
plasticity$tank <-as.factor(plasticity$tank)Morrissette-Boileau et al. (2018) examined factors affecting a shrub (dwarf birch) that is potentially expanding its range northwards with warming. They were interested in modifying factors - enhanced nutrients, warming, and browsing by caribou.
Caribou were not manipulated directly, but were excluded by fencing them out, and their foraging simulated by manual removal of 0, 25, or 75% of available shoots once each spring.
Nutrients were manipulated in 4 x 24 m plots, half of which received a nutrient supplement at bud-burst each year.
Each 4 x 24 area was divided into 4 x 4 m sub-plots. Each of these sub-plots had a 1 m x 1m area that was the subject of a combination of warming and caribou treatments.
Global warming was simulated using open-top hexagonal chambers that trap solar energy. They were used to establish two levels, ambient and warmed. Warming was applied to the target 1 x 1 m area thoughout the growing season.
The experiment ran for 5 years, after which time the authors recorded the total leafy biomass from each birch stem in the 1 x 1 target area. This time was the start of the 2014 growing season.
Their data are available from dryad. The relevant files are inside_productivity.csv and a very clear README_for_inside_productivity.txt file. We won’t explain the variables, and just point you to the readme. Figure 1 of the paper also has a nice diagram of the experiment.
Start by taking a quick look at the data file. The variables you may be working with for now are block, fertilization, temp, browsing, and biomass. Note that this file has a semicolon rather than comma as the delimiter, so make sure you take account of that when importing the data
df <- read_delim("../data/inside_productivity.csv",
delim = ";", escape_double = FALSE, trim_ws = TRUE)## Rows: 60 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (4): Block, Fertilization, Temp, Parcelle_id
## dbl (4): Biomass, Browsing, Biomass_beg2009, Biomass_end2009
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.set_sum_contrasts()## setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))head(df, 10)## # A tibble: 10 × 8
## Biomass Block Fertilization Temp Browsing Parcelle_id Biomass_beg2009
## <dbl> <chr> <chr> <chr> <dbl> <chr> <dbl>
## 1 535. A N H 0 ANH0 41
## 2 471. A N H 1 ANH1 34
## 3 244. A N H 3 ANH3 20
## 4 1172. A N L 0 ANL0 25
## 5 776. A N L 1 ANL1 27
## 6 928. A N L 3 ANL3 32
## 7 1274. A S H 0 ASH0 30
## 8 852. A S H 1 ASH1 28
## 9 1219. A S H 3 ASH3 22
## 10 1139. A S L 0 ASL0 17
## # ℹ 1 more variable: Biomass_end2009 <dbl>df$Browsing <- as.factor(df$Browsing)
df$Block <- as.factor(df$Block)
df$Temp <- as.factor(df$Temp)
df$Fertilization <-as.factor(df$Fertilization)Make sure you think carefully about scales of experimental units, which factors are between- and which are within-plots, and which are fixed vs random.
Make sure you include initial checking of assumptions
In the paper, the authors used initial biomass (Biomass_end2009) as a covariate, to account for initial variation in the number of stems occurring in sub-plots. Does inclusion of initial biomass improve the model fit?
Examine the residuals, and see if you think the response variable needs transformation, and if you do a transform, does the model fit any better?
In the discussion of the paper, there is discussion over the success of the method used to produce Warming. The technique has been used before, but data loggers in place during the experiment suggested no consistent warming was produced.
Would a simpler model omitting Warming lead you to any different conclusions?
The paper includes other response variables. One you could get practice with is gr_experiment.csv. This file has data for 3-6 shoots from each sub-plot, which underwent dendrochronological analysis. The shoots were sectioned and growth increments for each year between 1966-2013 recorded. These increments for each shoot are repeated measures, with the individual shoots as the “subjects”.
Think about how you might analyse this data set. If you want, you could try analyzing the full data set. An alternative could be to look at some slightly simpler ways:
If you analyze a single ring or an average, you’ll then have multiple measurements for each sub-plots. Are these values independent experimental units, or should they be averaged?
We’ve just considered the 4 x 24 m plots in which nutrients were in which nutrients were manipulated, but in the actual experiment, these plots were themselves grouped within 5 caribou fences, each with 12 plots.
You could challenge yourself by thinking about how you’d incorporate this aspect into your model & analysis, including the complications that might arise from adding a second random factor.