Load libraries and data

library(here)
library(tidyverse)
library(glmmTMB)
library(ggeffects)
library(MuMIn)

dino <- read.csv(here("Week_08", "Data", "psittacosaurus.csv"))

Take a look at the data

ggplot(dino, aes(x = Age, y = Mass)) + 
  geom_point() + 
  labs(title = "Estimatation of dinosaur growth rates") +
  scale_x_continuous(limits = c(0, 11), 
                     breaks = seq(from = 0, to = 11, by = 1)) +
  theme_light()

Three models: Linear Model

# Mass = a + b*Age
# nonlinear least squares (NLS). In terms of the model structure, we fit a nonlinear relationship between the response and the predictor(s), and assume the data is normally distributed around the expected value of the response
linearmod = nls(Mass ~ a + b*Age, #formula as the predictor
                data = dino,
                start = list(a = 0.01, #intercept at age = 0 (extremely small)
                             b = 4/2)) #slope at age = 6 to 8 rise by 4 run by 2

summary(linearmod)

## 
## Formula: Mass ~ a + b * Age
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a -1.86634    0.36948  -5.051 2.82e-06 ***
## b  2.07211    0.09394  22.057  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.432 on 78 degrees of freedom
## 
## Number of iterations to convergence: 1 
## Achieved convergence tolerance: 8.333e-07

Exponential model

# Mass = a*exp(r*Age)
# nonlinear least squares (NLS). In terms of the model structure, we fit a nonlinear relationship between the response and the predictor(s), and assume the data is normally distributed around the expected value of the response
expmod = nls(Mass ~ a*exp(r*Age), #formula as the predictor
                data = dino,
                start = list(a = 0.01, #initial mass
                             r = log(29)/10)) #exponential growth rate. at age 10, weight = 29. 29 = exp(10r)
summary(expmod)

## 
## Formula: Mass ~ a * exp(r * Age)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  0.97643    0.09734   10.03  1.1e-15 ***
## r  0.32766    0.01088   30.11  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.519 on 78 degrees of freedom
## 
## Number of iterations to convergence: 10 
## Achieved convergence tolerance: 2.761e-06

confint(expmod)

## Waiting for profiling to be done...

##        2.5%     97.5%
## a 0.8108754 1.1592835
## r 0.3088246 0.3475626

The estimated growth rate for the exponential model is r = 0.32766. The confidence interval is between 0.31-0.35.

Logistic model

# Mass M / (1 + exp(-r(Age - Age0)))
logmod = nls(Mass ~ M / (1 + exp(-r * (Age - Age0))), #formula as the predictor
                data = dino,
                start = list(M = 30, #max mass 
                             r = 3, #slope from age 7 to 8
                             Age0 = 8)) #inflection point, age of fastest growth
summary(logmod)

## 
## Formula: Mass ~ M/(1 + exp(-r * (Age - Age0)))
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## M    42.99240    3.67708   11.69   <2e-16 ***
## r     0.51254    0.02854   17.96   <2e-16 ***
## Age0  9.12515    0.35587   25.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.04 on 77 degrees of freedom
## 
## Number of iterations to convergence: 10 
## Achieved convergence tolerance: 7.101e-06

confint(logmod)

## Waiting for profiling to be done...

##            2.5%      97.5%
## M    37.2608567 51.4490047
## r     0.4639475  0.5668192
## Age0  8.5454204  9.8756227

The estimated growth rate for the logistic model is r = 0.5124. The confidence interval is between 0.46-0.57.

Compare the models

plinearmod = predict(linearmod)
pexpmod = predict(expmod)
plogmod = predict(logmod)

plot(dino$Age, dino$Mass, 
     xlab = "Age", ylab = "Mass") 
lines(dino$Age, plinearmod, col = "blue") 
lines(dino$Age, pexpmod, col = "red") 
lines(dino$Age, plogmod, col = "green")

AICc

aicclinearmod = AICc(linearmod)
aiccexpmod = AICc(expmod)
aicclogmod = AICc(logmod)

cat("Linear AICc: ", aicclinearmod, "\n")

## Linear AICc:  373.5375

cat("Exponential AICc: ", aiccexpmod, "\n")

## Exponential AICc:  298.1695

cat("Logistic AICc: ", aicclogmod, "\n")

## Logistic AICc:  238.7308

Smaller AIC means that there is a better approxamation to the truth - given that the logistic model has the smallest AIC, we can assume this is the best model.

Is there evidence that this dinosaur has a maximum size? If so, what is the estimate for that size, and what is the confidence interval around that estimate? How does the estimated maximum size compare to the largest size in the data? How much stock do you put in the Mmax estimate, given the data we have? If this estimate of Mmax is true, about how big does this dinosaur get, relative to a human?

The logistic model predicts a maximum size for the dinosaur, which estimates a maximum of 42.99240 kg with a CI between 37.2608567 - 51.4490047. The estimate for this maximum size is larger than the actual data which is around 30 kg for the maximum size. This may be due to extrapolation. The average weight of a human is 62 kg and I’m surprised the dinosaur is smaller.