Introduciton

This is a first-cut analysis of a greenhouse experiment conducted in the summer of 2014 in Davis, CA.

Data preparation

We’ll start by normalizing the data by species. Then, we’ll plot all the data together and do statistical analyes. In the normalization procedure, each data point is scaled by species so that the distribution of values for each species has a mean of 0 and a standard deviation of 1. After scaling, the data can be combined to look for overall treatment effects. We’ll also take care of a bunch of house cleaning along the way.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
library(reshape2)
setwd('C:/Users/AnacB1/Dropbox/12_psf_dens_experiment/')
source('04_Analysis/multiplot.R')
source('04_Analysis/theme_ba.R')
source('04_Analysis/theme_rot.R')
source('04_Analysis/summarySE.R')

d <- read.csv('02_Data/09_complete_cleaned_final.csv')
trt <- read.csv('01_ExpDesignMethods/01_exp_design/tags_final_n600.csv')
d <- d[!is.na(d$fc_repro_wt) & !is.na(d$fc_st_wt), ]
d$fc_wt <- with(d, fc_st_wt + fc_rt_wt + fc_repro_wt)
d$cp_wt <- with(d, cp_sh_wt + cp_rt_wt)
d$shoot_wt <- with(d, fc_st_wt + fc_repro_wt)
d$fc_wt_log <- log(d$fc_wt)
d$rs <- with(d, fc_rt_wt / shoot_wt)
d2 <- subset(d, pot != 172)
d3 <- merge(d2, trt)
d4 <- transform(d3, fc_wt.norm = ave(fc_wt, sp, FUN = scale))
d5 <- transform(d4, rs.norm = ave(rs, sp, FUN = scale))
d5$extract <- factor(d5$extract)
levels(d5$extract) <- c('large', 'small', 'sterile')
d5$comp <- factor(d5$comp)
levels(d5$comp) <- c('Alone', 'Lasthenia', 'Poa', 'Conspecific')
d5$sp <- factor(d5$sp)
levels(d5$sp) <- c('Daucus', 'Madia', 'Plantago', 'Potentilla', 'Trifolium')
d5$density <- factor(d5$density, levels = c('L', 'H'))
d5$comp <- factor(d5$comp, levels = c('Alone', 'Conspecific', 'Lasthenia', 'Poa'))
d5$extract <- factor(d5$extract, levels = c('sterile', 'small', 'large'))
d6 <- summarySE(d5, measurevar = 'fc_wt', groupvars = c('sp', 'density', 'extract', 'comp'))
d7 <- summarySE(d5, measurevar = 'fc_wt.norm', groupvars = c('density', 'extract', 'comp'))
d8 <- summarySE(d5, measurevar = 'rs', groupvars = c('sp', 'density', 'extract', 'comp'))
d9 <- summarySE(d5, measurevar = 'rs.norm', groupvars = c('density', 'extract', 'comp'))

A reminder on sample design

Let’s start with a reminder about the sample design.

# sample design
# 4 comp levels * 2 field densities * 3 extracts = 24 treatment combinations per species

# The number of reps varied by species due to seed availability and germination success
  # dp = 6; ms = 4; pa = 5; pe = 6; tb = 4


# so, our number of pots should be 
# 24 * 6 + 24 * 4 + 24 * 5 + 24 * 6 + 24 * 4 = 600

table(d5$comp, d5$extract, d5$density, d5$sp)
## , ,  = L,  = Daucus
## 
##              
##               sterile small large
##   Alone             6     6     6
##   Conspecific       6     6     6
##   Lasthenia         6     6     6
##   Poa               6     6     6
## 
## , ,  = H,  = Daucus
## 
##              
##               sterile small large
##   Alone             6     6     6
##   Conspecific       6     6     6
##   Lasthenia         6     6     6
##   Poa               6     6     6
## 
## , ,  = L,  = Madia
## 
##              
##               sterile small large
##   Alone             4     4     4
##   Conspecific       4     4     4
##   Lasthenia         4     4     4
##   Poa               4     4     4
## 
## , ,  = H,  = Madia
## 
##              
##               sterile small large
##   Alone             4     4     4
##   Conspecific       4     4     4
##   Lasthenia         4     4     4
##   Poa               4     4     4
## 
## , ,  = L,  = Plantago
## 
##              
##               sterile small large
##   Alone             6     6     6
##   Conspecific       6     6     6
##   Lasthenia         6     6     6
##   Poa               5     6     6
## 
## , ,  = H,  = Plantago
## 
##              
##               sterile small large
##   Alone             6     6     6
##   Conspecific       6     6     6
##   Lasthenia         6     6     6
##   Poa               6     6     6
## 
## , ,  = L,  = Potentilla
## 
##              
##               sterile small large
##   Alone             5     5     5
##   Conspecific       5     5     5
##   Lasthenia         5     5     5
##   Poa               5     5     5
## 
## , ,  = H,  = Potentilla
## 
##              
##               sterile small large
##   Alone             5     5     5
##   Conspecific       5     5     5
##   Lasthenia         5     5     5
##   Poa               5     5     5
## 
## , ,  = L,  = Trifolium
## 
##              
##               sterile small large
##   Alone             4     4     4
##   Conspecific       4     4     4
##   Lasthenia         3     4     4
##   Poa               3     4     4
## 
## , ,  = H,  = Trifolium
## 
##              
##               sterile small large
##   Alone             4     4     4
##   Conspecific       4     4     4
##   Lasthenia         4     3     4
##   Poa               4     4     4

In total, we should have 600 pots, but we ended up with data for 596 pots. As I recall, 4 of them got dumped over/spilled. Looks like the missing pots include 1 for plantago and three for Trifolium. But, in general, we have pretty complete dataset.

Plot of performance

OK, now that we are done with that, we can plot scaled biomass by treatment.

dodge <- position_dodge(width = .5) 


p1 <- ggplot(data = d7, aes(x = comp, y = fc_wt.norm, group = extract, 
  fill = density, shape =extract)) +
    geom_errorbar(aes(ymin = fc_wt.norm - se, ymax = fc_wt.norm + se), 
        color = 'black', width=.1, position = dodge) +
    geom_point(size = 3.5, position = dodge) +
    theme_ba() +
    scale_color_discrete(name="Field\nDensity") +
    scale_shape_discrete(name="Soil\nExtract") +
    xlab('Competition Treatment') +
    ylab('Biomass (scaled by focal species)') +
    ggtitle('Plant biomass by treatment') +
    scale_shape_manual(values = c(21, 22, 24)) +
    guides(fill = guide_legend(override.aes = list(shape = 21)))

p1

Looking at the results, it’s clear that performance goes down in the presence of competition, and performance is higher when growing in soils sampled from below patches with a high density of conspecifics in the field. The effect of extract looks to be non-existent. Let’s fit some models and see what is going on.

Models of performance

m1 <- aov(fc_wt.norm ~ comp + extract + density, data = d5)
summary(m1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## comp          3  132.3   44.11  58.874  < 2e-16 ***
## extract       2    0.9    0.45   0.605    0.546    
## density       1   16.4   16.41  21.897 3.57e-06 ***
## Residuals   589  441.3    0.75                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s look at the degrees of freedom: 589 in the residuals, 1 for density, 2 for extract, and 3 for competition. So 589 + 1 + 2 + 3 = 595, which is our n minus 1. So, that’s looking good. Does it seem right that we should be treating every pot independently? I think so.

As expected from looking at the plot, there is a significant effect of competition and density, but not extract. Let’s check for interactions.

full <- aov(fc_wt.norm ~ comp * extract * density, data = d5)
sub1 <- aov(fc_wt.norm ~ comp + extract * density, data = d5)
sub2 <- aov(fc_wt.norm ~ comp * extract + density, data = d5)
sub3 <- aov(fc_wt.norm ~ comp * density + extract, data = d5)

AIC(m1, full, sub1, sub2, sub3)
##      df      AIC
## m1    8 1528.324
## full 25 1558.722
## sub1 10 1531.930
## sub2 14 1539.308
## sub3 11 1533.548

Best model is clearly the model without interactions.

Did the extract treatment not work?

So, did our extract treatment just not work? We had good proof of concept with Trifolium as this photo suggests:

TB roots:

We also counted nodules for Trifolium. Let’s look at that data.

Trifolium nodulation

setwd('C:/Users/AnacB1/Dropbox/12_psf_dens_experiment/')
d <- read.csv('C:/Users/AnacB1/Dropbox/12_psf_dens_experiment/02_Data/10_trif_nodules.csv')
names(d)[1] <- 'pot'
trt <- read.csv('01_ExpDesignMethods/01_exp_design/tags_final_n600.csv')
d2 <- merge(d, trt)
d3 <- na.omit(d2)
d5 = d3
d5$extract <- factor(d5$extract)
levels(d5$extract) <- c('large', 'small', 'sterile')
d5$comp <- factor(d5$comp)
levels(d5$comp) <- c('Alone', 'Lasthenia', 'Poa', 'Conspecific')
d5$density <- factor(d5$density, levels = c('L', 'H'))
d5$comp <- factor(d5$comp, levels = c('Alone', 'Conspecific', 'Lasthenia', 'Poa'))
d5$extract <- factor(d5$extract, levels = c('sterile', 'small', 'large'))
d7 <- summarySE(d5, measurevar = 'total_nod', groupvars = c('density', 'extract', 'comp'))

p1 <- ggplot(data = d7, aes(x = comp, y = total_nod, group = extract, fill = density, shape =extract)) +
    geom_errorbar(aes(ymin = total_nod - se, ymax = total_nod + se), 
        color = 'black', width=.1, position = dodge) +
    geom_point(size = 3.5, position = dodge) +
    theme_ba() +
    scale_color_discrete(name="Field\nDensity") +
    scale_shape_discrete(name="Soil\nExtract") +
    xlab('Competition Treatment') +
    ylab('Number of nodules') +
    ggtitle('Trifolium nodulation by treatment') +
    scale_shape_manual(values = c(21, 22, 24)) +
    guides(fill = guide_legend(override.aes = list(shape = 21)))
## Scale for 'shape' is already present. Adding another scale for 'shape',
## which will replace the existing scale.
p1

Seems clear that nodulation is responding to our treatments.

m1 <- aov(total_nod ~ comp + extract + density, data = d7)
summary(m1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## comp         3  55085   18362   5.690  0.00692 ** 
## extract      2 116983   58492  18.126 6.09e-05 ***
## density      1  11340   11340   3.514  0.07813 .  
## Residuals   17  54857    3227                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
full <- aov(total_nod ~ comp * extract * density, data = d7)
sub1 <- aov(total_nod ~ comp + extract * density, data = d7)
sub2 <- aov(total_nod ~ comp * extract + density, data = d7)
sub3 <- aov(total_nod ~ comp * density + extract, data = d7)
AIC(m1, full, sub1, sub2, sub3)
##      df      AIC
## m1    8 269.7356
## full 25     -Inf
## sub1 10 271.0811
## sub2 14 257.1821
## sub3 11 271.5945
summary(sub2)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## comp          3  55085   18362  10.242  0.00163 ** 
## extract       2 116983   58492  32.626 2.37e-05 ***
## density       1  11340   11340   6.325  0.02873 *  
## comp:extract  6  35136    5856   3.266  0.04268 *  
## Residuals    11  19721    1793                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now we are talking! A signficant competition by extract interation. Looking back at the plot it looks like the extract effect is greatest in the absence of competition! Cool :)

There’s also a significant density effect, where nodulation is higher in low density soils…interesting.

Let’s re-do the biomass vs treatments model for each species seperately.

setwd('C:/Users/AnacB1/Dropbox/12_psf_dens_experiment/')
d <- read.csv('02_Data/09_complete_cleaned_final.csv')
trt <- read.csv('01_ExpDesignMethods/01_exp_design/tags_final_n600.csv')
d <- d[!is.na(d$fc_repro_wt) & !is.na(d$fc_st_wt), ]
d$fc_wt <- with(d, fc_st_wt + fc_rt_wt + fc_repro_wt)
d$cp_wt <- with(d, cp_sh_wt + cp_rt_wt)
d$shoot_wt <- with(d, fc_st_wt + fc_repro_wt)
d$fc_wt_log <- log(d$fc_wt)
d$rs <- with(d, fc_rt_wt / shoot_wt)
d2 <- subset(d, pot != 172)
d3 <- merge(d2, trt)
d4 <- transform(d3, fc_wt.norm = ave(fc_wt, sp, FUN = scale))
d5 <- transform(d4, rs.norm = ave(rs, sp, FUN = scale))
d5$extract <- factor(d5$extract)
levels(d5$extract) <- c('large', 'small', 'sterile')
d5$comp <- factor(d5$comp)
levels(d5$comp) <- c('Alone', 'Lasthenia', 'Poa', 'Conspecific')
d5$sp <- factor(d5$sp)
levels(d5$sp) <- c('Daucus', 'Madia', 'Plantago', 'Potentilla', 'Trifolium')
d5$density <- factor(d5$density, levels = c('L', 'H'))
d5$comp <- factor(d5$comp, levels = c('Alone', 'Conspecific', 'Lasthenia', 'Poa'))
d5$extract <- factor(d5$extract, levels = c('sterile', 'small', 'large'))


d5L <- split(d5, d5$sp)

mods <- lapply(d5L, function(x){
  #x <- d5L[[1]]
  m1 <- aov(fc_wt ~ comp + extract + density, data = x)
  return(summary(m1)) 
})

mods
## $Daucus
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## comp          3  13.95   4.651  19.562 1.29e-10 ***
## extract       2   0.16   0.080   0.335    0.716    
## density       1  16.66  16.658  70.057 5.96e-14 ***
## Residuals   137  32.58   0.238                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Madia
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## comp         3 12.224   4.075  12.267 8.64e-07 ***
## extract      2  0.174   0.087   0.262 0.769898    
## density      1  4.360   4.360  13.126 0.000484 ***
## Residuals   89 29.563   0.332                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Plantago
##              Df Sum Sq Mean Sq F value Pr(>F)    
## comp          3  59.61  19.871  47.816 <2e-16 ***
## extract       2   0.50   0.250   0.601  0.550    
## density       1   0.10   0.096   0.230  0.632    
## Residuals   136  56.52   0.416                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Potentilla
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## comp          3  24.73   8.244   8.328 4.73e-05 ***
## extract       2   0.48   0.242   0.245    0.783    
## density       1   0.22   0.223   0.226    0.636    
## Residuals   113 111.86   0.990                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Trifolium
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## comp         3  5.993  1.9976  14.464 1.04e-07 ***
## extract      2  1.480  0.7398   5.357  0.00642 ** 
## density      1  0.327  0.3272   2.369  0.12741    
## Residuals   86 11.877  0.1381                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks like extract has a signficant effect on biomass for Triolium only..

Next steps

Root to shoot ratio….