Improving the models: What’s next?

What’s next? Any aspects of the experimental design missing from the model?

We can construct models that includes any missing factors.

We can also make an individual-level random effect to use for overdispersed poisson:

d1$obs <- 1:length(d1$y) ## makes a unique number for each row in dataset

We can also create other models that include different random effects that account for blocks:

r0 <- glmmTMB(y ~ spray * lead, data=d1, family="gaussian")              ### gaussian glm just to see how bad it really is
r1 <- glmmTMB(y ~ spray * lead, data=d1, family="poisson")               ### poisson glm
r2 <- glmmTMB(y ~ spray * lead, data=d1, family="nbinom2")               ### nb glm
r3 <- glmmTMB(y ~ spray * lead + (1|obs), data=d1, family="poisson")     ### overdispersed poisson glm
r4 <- glmmTMB(y ~ spray * lead + (1|block), data=d1, family="nbinom2")   ### nb w/ block
r5 <- glmmTMB(y ~ spray * lead + (1|obs) + (1|block), data=d1, family="poisson")   ### OD poisson w/ block

Six models above differ only distribution and random effects. Fixed effects are the same. We can now use AIC for selecting the most appropriate distribution and random effects:

#### for model selection, use AIC or likihood ratio test
model.sel(r0,r1,r2,r3,r4,r5) ## from MuMIn

## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(led) cnd(spr) cnd(led:spr)             family
## r4     0.2835          +        +        +            +       nbinom2(log)
## r5     0.1399          +        +        +            +       poisson(log)
## r3     0.1042          +        +        +            +       poisson(log)
## r2     0.3365          +        +        +            +       nbinom2(log)
## r1     0.3365          +        +        +            +       poisson(log)
## r0     1.4000          +        +        +            + gaussian(identity)
##       random df    logLik   AICc  delta weight
## r4      c(b)  6 -1497.706 3007.5   0.00  0.903
## r5 c(o)+c(b)  6 -1499.931 3011.9   4.45  0.097
## r3      c(o)  5 -1518.031 3046.1  38.63  0.000
## r2            5 -1521.484 3053.0  45.54  0.000
## r1            4 -1558.739 3125.5 118.03  0.000
## r0            5 -1929.308 3868.7 861.18  0.000
## Models ranked by AICc(x) 
## Random terms: 
## c(b) = 'cond(1 | block)'
## c(o) = 'cond(1 | obs)'

AICtab(r0,r1,r2,r3,r4,r5,base=T,logLik=T,weights=T)    ## from bbmle

##    logLik  AIC     dLogLik dAIC    df weight
## r4 -1497.7  3007.4   431.6     0.0 6  0.903 
## r5 -1499.9  3011.9   429.4     4.5 6  0.097 
## r3 -1518.0  3046.1   411.3    38.7 5  <0.001
## r2 -1521.5  3053.0   407.8    45.6 5  <0.001
## r1 -1558.7  3125.5   370.6   118.1 4  <0.001
## r0 -1929.3  3868.6     0.0   861.2 5  <0.001

Now for further model interpretation you can examine residuals, fixed effects using Anova(), emmeans, etc.

Binomial

In this section we look at the binomial distribution. Let’s load in the data we will use:

data("wheatley.carrot")
?wheatley.carrot

Let’s clean it up a bit:

dat1 <- wheatley.carrot %>% filter(insecticide!='nil')
dat1$propdmg <- dat1$damaged/dat1$total
head(dat1) 

##   treatment insecticide depth rep damaged total    propdmg
## 1       T01    diazinon   1.0  R1     120   187 0.64171123
## 2       T02    diazinon   2.5  R1      60   184 0.32608696
## 3       T03    diazinon   5.0  R1      35   179 0.19553073
## 4       T04    diazinon  10.0  R1       5   178 0.02808989
## 5       T05    diazinon  25.0  R1      66   182 0.36263736
## 6       T06  disulfoton   1.0  R1      97   187 0.51871658

hist(dat1$propdmg)

Now let’s plot out the data. Note that for plotting we can use the “beta_family” although for analysis we use “binomial”.

ggplot(dat1, aes(x=log(depth),
                 y=(damaged/total), color=insecticide)) +
  geom_point(size=3) + 
  theme_bw(base_size = 16)

ggplot(dat1, aes(x=log(depth), y=(damaged/total), color=insecticide)) +
  geom_point(size=3) +
  geom_smooth(method='glm',
              method.args=list(family="beta_family"),
              formula = y~x) +
  theme_bw(base_size = 16)

ggplot(dat1, aes(x=log(depth), y=(damaged/total),
                 color=insecticide)) + geom_point() +
  geom_smooth(method='glm', 
              method.args=list(family="beta_family"), 
              formula = y~x+I(x^2))

Now let’s build two different binomical GLMMs with different set ups:

mod1 <- glmmTMB(cbind(damaged,total-damaged) ~ insecticide * depth + (1|rep), data=dat1, family='binomial')

mod2 <- glmmTMB(cbind(damaged,total-damaged) ~ insecticide * depth + I(depth^2) + (1|rep), data=dat1, family='binomial')

Note above that mod2 is applying a non-linear modification where we include depth and depth squared. The “^” is another way of saying raise depth to the 2nd power. The I() isolate this and tells the R interpreter that the model is including a second order fixed effect. For more info check ?formula.

We can plot residuals:

plot(simulateResiduals(mod1))

plot(simulateResiduals(mod2))

Which of the two is better?

We can also run ANOVAs on the models:

Anova(mod1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: cbind(damaged, total - damaged)
##                    Chisq Df Pr(>Chisq)    
## insecticide       91.449  1  < 2.2e-16 ***
## depth              2.543  1  0.1107873    
## insecticide:depth 14.525  1  0.0001383 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(mod2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: cbind(damaged, total - damaged)
##                    Chisq Df Pr(>Chisq)    
## insecticide       111.62  1  < 2.2e-16 ***
## depth             563.71  1  < 2.2e-16 ***
## I(depth^2)        594.17  1  < 2.2e-16 ***
## insecticide:depth  15.85  1  6.858e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##Beta

In this section we show an example of applying the beta family using true proportions (or percentages). First let’s load data:

b1 <-read_csv("LynnETAL_2022_Ecography_sodCaseStudy.csv") %>%
  filter(species=='SALA')
head(b1)

## # A tibble: 6 × 18
##   study        family speci…¹ species meas_…² site  latit…³ longi…⁴ o_lat o_long
##   <chr>        <chr>  <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl> <dbl>  <dbl>
## 1 RivkinETAL_… Alism… Sagitt… SALA    2015    ACT-…    44.5   -77.3    NA     NA
## 2 RivkinETAL_… Alism… Sagitt… SALA    2015    BPV-…    45.7   -77.6    NA     NA
## 3 RivkinETAL_… Alism… Sagitt… SALA    2015    BSM-…    44.8   -76.1    NA     NA
## 4 RivkinETAL_… Alism… Sagitt… SALA    2015    BTC-…    44.2   -76.8    NA     NA
## 5 RivkinETAL_… Alism… Sagitt… SALA    2015    CBT-…    47.4   -79.7    NA     NA
## 6 RivkinETAL_… Alism… Sagitt… SALA    2015    CFR-…    44.7   -77.1    NA     NA
## # … with 8 more variables: mean_herb <dbl>, na_wdcpc <dbl>, mat <dbl>,
## #   iso <dbl>, dtr <dbl>, precip <dbl>, aet <dbl>, notes <chr>, and abbreviated
## #   variable names ¹​species_name, ²​meas_year, ³​latitude, ⁴​longitude

Let’s examine a histogram:

hist(b1$mean_herb) ## examine histogram. Data are entered as percentage (0-100%)

Let’s plot out the data. Note: data was collected at 42 unique sites (one data point per site).

ggplot(b1 , aes(x=latitude, y=mean_herb)) +
  geom_point()

Now we need to convert to proportion and then do a small transformation to remove 0 and 1:

b1$mean_herb1 <- (b1$mean_herb/100)
b1$mean_herb1 <- (b1$mean_herb1*(length(b1$mean_herb1)-1)+.5)/length(b1$mean_herb1)

Let’s plot this:

ggplot(b1, aes(x=latitude, y=mean_herb1)) + geom_point() +
  geom_smooth(method='glm', 
              method.args=list(family="beta_family"),
              formula = y~x)

Now let’s build a beta GLMM:

bm1 <- glmmTMB(mean_herb1 ~ latitude, data=b1, family='beta_family')
Anova(bm1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: mean_herb1
##           Chisq Df Pr(>Chisq)   
## latitude 9.1383  1   0.002503 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check the summary:

summary(bm1)

##  Family: beta  ( logit )
## Formula:          mean_herb1 ~ latitude
## Data: b1
## 
##      AIC      BIC   logLik deviance df.resid 
##   -126.9   -121.7     66.5   -132.9       39 
## 
## 
## Dispersion parameter for beta family (): 30.5 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  5.69631    2.58221   2.206   0.0274 * 
## latitude    -0.17525    0.05797  -3.023   0.0025 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot and check residuals:

plot(simulateResiduals(bm1)) ## QQ plot is good. Resid vs pred is wiggle, probably ok.

Problem set

1. Problem set with real data (Rivkin et al. 2018 Am J Bot)

Herbivory data was collected on 43 populations of an aquatic plant across a latitudinal gradient in Canada. At each population, many plants (~5-15) were examined for herbivory damage. Some additional covariates were recorded, such as Competition around the plant (1-3 from less to more) and plant height (cm).

Read in data:

d1 <-read.csv("ajb21098-sup-0002-appendixs2.csv")
head(d1)

##   Population Latitude Longitude Competition Individual PlantHeight LeafDamage
## 1     ACT-ON 44.54829 -77.32384           2          4        41.0 0.02000000
## 2     ACT-ON 44.54829 -77.32384           2          6        43.0 0.03000000
## 3     ACT-ON 44.54829 -77.32384           2         18        17.5 0.04000000
## 4     ACT-ON 44.54829 -77.32384           2          2        37.0 0.06666667
## 5     ACT-ON 44.54829 -77.32384           2          3        48.0 0.07000000
## 6     ACT-ON 44.54829 -77.32384           2         12        39.5 0.08000000
##   LarvalAbundance
## 1               0
## 2               0
## 3               0
## 4               1
## 5               0
## 6               0

hist(d1$LeafDamage)

Let’s remove 0’s and 1’s (see Smithson & Verkuilen 2006 or Douma & Weedon 2018):

d1$LeafDamage <- (d1$LeafDamage*(length(d1$LeafDamage)-1)+.5)/length(d1$LeafDamage)

Now plot data:

ggplot(d1 , aes(x=Latitude, y=LeafDamage)) +
  geom_point() +
  geom_smooth(method='glm', 
              method.args=list(family="beta_family"),
              formula = y~x)

Question: Does herbivory increase towards the equator? How do residuals look? Any way to improve them?

2. Problem set with real data (Rivkin et al. 2018 Am J Bot), Rats dataset with only females. Half treated with a drug and the other half were untreated. Then they checked for tumors.Does the drug reduce probability of developing tumor? After 50 days? After 100 days?

Load data:

library(survival)
rats <- rats %>% filter(sex=="f")
?rats

## Help on topic 'rats' was found in the following packages:
## 
##   Package               Library
##   aods3                 /Users/leoohyama/Library/R/x86_64/4.1/library
##   survival              /Library/Frameworks/R.framework/Versions/4.1/Resources/library
## 
## 
## Using the first match ...

head(rats)

##   litter rx time status sex
## 1      1  1  101      0   f
## 2      1  0   49      1   f
## 3      1  0  104      0   f
## 4      3  1  104      0   f
## 5      3  0  102      0   f
## 6      3  0  104      0   f

Plot data:

ggplot(rats , aes(x=time, y=status)) + 
  geom_point() +
  geom_smooth(method='glm', 
              method.args=list(family="binomial"), 
              formula = y~x) + 
  facet_wrap(~rx)