OCN 683 - Homework 9
Sarah Pennington
We would like to know how fast dinosaurs grew (obviously). Some clever folks have extracted size-at-age data from fossils using growth rings in the bones. For the species Psittacosaurus lujiatunensis: .
Fit the following three curves: Linear, exponential, logistic
Fit these three curves to the dataset ‘psittacosaurus.csv’, using nls() as described in lecture. You’ll need to supply start values for the parameters, which may be tricky. If you get error messages when trying to fit the model, use trial-and-error, or try plotting what the curve looks like with different parameter values to come up with reasonable guesses. For each model report the coefficient estimates and confidence intervals, and plot the fitted curves on top of the raw data.
dino <- read_csv("psittacosaurus.csv")## Rows: 80 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): Age, Mass
##
## ℹ 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.ggplot(dino, aes(x = Age, y = Mass)) + geom_point() +theme_clean() + scale_x_continuous(limits = c(0, 11), breaks = seq(0, 11, by = 1))
#Fit three curves
#a=initial slop, h= asymptote from lecture example
###Linear model
linear <- nls(Mass ~ a + (b*Age), data= dino, start = list(a = 0.01, b = 10/7))
plot(linear)
summary(linear)##
## 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: 3.703e-07#coef: a = -1.86634, b = 2.07211
confint(linear)## Waiting for profiling to be done...## 2.5% 97.5%
## a -2.601915 -1.130762
## b 1.885089 2.259139#confidence intervals
# 2.5% 97.5%
#a -2.601915 -1.130762
#b 1.885089 2.259139
predicted_linear <- predict(linear, newdata = dino)
###Exponential model
expo <- nls(Mass ~ exp(r*Age)*a, data = dino, start = list (a = 0.005, r = 0.5))
plot(expo)
summary(expo)##
## Formula: Mass ~ exp(r * Age) * a
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.97642 0.09734 10.03 1.1e-15 ***
## r 0.32767 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: 12
## Achieved convergence tolerance: 7.481e-06#coef: a = 0.97642, r = 0.32767
confint(expo)## Waiting for profiling to be done...## 2.5% 97.5%
## a 0.8108754 1.1592835
## r 0.3088246 0.3475626# 2.5% 97.5%
#a 0.8108754 1.1592835
#r 0.3088246 0.3475626
predicted_expo <- predict(expo, newdata = dino)
###Logistic model
logistic <- nls(Mass ~ Mmax / (1+ exp(-r *(Age - Age0))), data = dino, start = list (r = 1, Age0 = 8, Mmax = 30))
plot(logistic)
summary(logistic)##
## Formula: Mass ~ Mmax/(1 + exp(-r * (Age - Age0)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## r 0.51254 0.02854 17.96 <2e-16 ***
## Age0 9.12513 0.35587 25.64 <2e-16 ***
## Mmax 42.99217 3.67703 11.69 <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: 9
## Achieved convergence tolerance: 2.767e-06#coef: r = 0.51254, Age0 = 9.12513, Mmax = 42.99217
#confint(logistic)
# 2.5% 97.5%
#r 0.4639475 0.5668192
#Age0 8.5454204 9.8756227
#Mmax 37.2608575 51.4490031
predicted_log <- predict(logistic, newdata = dino)
###Plot all 3 models
plot(dino$Age, dino$Mass,
xlab = "Age", ylab = "Mass") +
lines(dino$Age, predicted_linear, col = "blue", lwd = 1) + lines(dino$Age, predicted_expo, col = "red", lwd = 1) + lines(dino$Age, predicted_log, col = "green", lwd = 1)
## integer(0)Compare the three models using AICc. Which model is the best? What are the ΔAICc values and the Akaike weights for the three models? How do you interpret these results in terms of the relative support for the three models?
#Compare the three models using AICc. Which model is the best?
library(MuMIn)##
## Attaching package: 'MuMIn'## The following object is masked from 'package:arm':
##
## coefplotAICc(linear) #373.5375## [1] 373.5375AICc(expo) #298.1695## [1] 298.1695AICc(logistic) #238.7308## [1] 238.7308 ###The logistic model is the best.
# What are the ΔAICc values and the Akaike weights for the three models?
AIC(linear) #373.2217## [1] 373.2217AIC(expo) #297.8537## [1] 297.8537AIC(logistic) #238.1975## [1] 238.1975lin <- AICc(linear) - AIC(linear) #substantial
expo <- AICc(expo) - AIC(expo) #substantial
log <- AICc(logistic) - AIC(logistic) #substantial
#use delta-AIC to calculate the model likelihoods
model.like.1 = exp(-0.5*(AICc(linear) - AICc(logistic)))
model.like.2 = exp(-0.5*(AICc(expo) - AICc(logistic)))## Error in UseMethod("logLik"): no applicable method for 'logLik' applied to an object of class "c('double', 'numeric')"model.like.3 = exp(-0.5*(AICc(logistic) - AICc(logistic)))
#sum the model likelihoods for standardization
summed.likes = sum(c(model.like.1, model.like.2, model.like.3))## Error: object 'model.like.2' not found#calculate the Akaike weights
weight1 = model.like.1/summed.likes## Error: object 'summed.likes' not foundweight2 = model.like.2/summed.likes## Error: object 'model.like.2' not foundweight3 = model.like.3/summed.likes## Error: object 'summed.likes' not found#How do you interpret these results in terms of the relative support for the #three models?
### Based on these results, the logistic model has the most support as the one that fits the data the best. What is the estimated exponential growth rate (r) for the exponential and logistic models? What is the confidence interval on this parameter for the two models? For exponential growth, the doubling time is log(2)/r. How long does it take the dinosaur to double in size, based on the two models?
###Estimated growth rate (r) is r = 0.51254 for logistic and r = 0.32767 for exponential.
### For the logistic model the confidence interval for r is 0.46- 0.57. For the exponential model the confidence interval is 0.31 - 0.35.
dt <- log(2)/0.32767 #2.11
### According to the exponential model, it takes approx 2 years for the dinosaur to double in size. #####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 does predict a maximum size for the dinosaur.
###The logistic models estimate for the maximum size of the dinosaur is 42.99kg with a 95% confidence interval of 37.2608575- 51.4490031kg.
###The largest size in the data is 29.67kg.
###The estimate for maximum size is larger than the largest sized dinosaur observed in the data, so the maximum size is an extrapolation that is hard to take seriously without other evidence that the dinosaurs can grow that large.
###Quick google search says the average kg weight of a human is 62kg, so this dinosaud at its largest observed size is about half the mass of the average person.