Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
week_06_exercise_v4.R
#######################################################
# BIOSCI220 WEEK 6 LAB EXERCISE
#######################################################
################################################
# Note for people using online R
################################################
#
# In previous labs, people who couldn't get R
# installed on a home machine could use the
# "R Snippets" website to run the lab code.
#
# However, I have discovered that R Snippets
# puts a 10-second limit on R runs, and our
# Maximum Likelihood searches in this lab
# will often take a little longer than this.
#
# Therefore, we have to use something else.
# The best option I have found is RStudio Cloud.
#
# I have run all of the lab code in RStudio Cloud,
# and it works just like regular RStudio.
#
# Also, it remembers your previous actions, so
# you no longer have to do the "copy-paste the
# whole script down to your question" business
# that R Snippets required. (However, this can
# still be a useful tactic for "resetting"
# everything, if you are getting errors due to
# having accidentally skipped some steps.)
#
# However, you do need to register for a free
# account. I have put instructions below.
#
################################################
# Instructions for using RStudio Cloud:
################################################
#
# 1. Go to https://rstudio.cloud/
#
# 2. Register for a free account, e.g. using a gmail address.
#
# 3. Wait awhile for your Workspace to load. Be patient.
#
# 4. Click the blue "New Project" button in the upper right.
#
# 5. Wait awhile for your new project to load. Be patient.
#
# 6. Click File -> New File -> R script
#
# 7. Paste your R script into this window. You should be able to
# run the script's code by either (a) copy-pasting into
# the Console, or (b) clicking the "Run" button.
#
# 8. Go to the Console. Type
#
# install.packages("deSolve")
#
# ...to get deSolve installed.
#
# 9. Type library(deSolve). If no errors are returned, then
# deSolve is working, and you can start the lab.
################################################
#######################################################
# Introduction
#######################################################
# In week 5, we got a lot of experience working with SIR models,
# and simulating possible futures, given different values of R0
# for different levels of lockdown, different timings of interventions,
# different numbers of starting infections, etc.
#
# However, many students came across one big difficulty -- how do we
# know what parameters to use? In some cases, we can justify parameters
# with background knowledge. For example, it appears, from medical
# observations, that COVID-19 patients are infectious for about 14 days.
# We don't necessarily need fancy inference methods for that.
#
# However, some parameters, like R0, are fairly uncertain -- in the
# media, you will see scientists quoting numbers between 2 and 6.
# Furthermore, the R0 will (hopefully!) change depending on how well
# a country's leadership and population enact the various control
# measures. We would like to *measure* this, if possible.
#
# Regarding the number of starting infections, we rarely know
# exactly how many infections were imported into a country at
# the start of an epidemic. We can guesstimate based on the number
# of early cases, but we might want to use a more formal method.
#
# Furthermore, we would like to have some measure of uncertainty
# about these parameters. Are we highly confident that R0 is 3.0,
# plus or minus 0.1? Or are we highly uncertain, and we can only
# put R0 between 2 and 6?
#
# Answers to these questions can be found, pretty easily, with
# MODEL-BASED INFERENCE.
#
# In this lab, we will perform some model-based inference on SIR
# models. This will give you a decent taste of what real-life
# model-fitting is like -- although I hasten to add that we are
# using quite simple models, and professional epidemiologists
# would use somewhat more complex models to help fix parts of
# our models that are unrealistic (for example, our models allow
# there to be 0.0639 of an infection, whereas in real life there
# is either 0 or 1 infections).
#
# In a quick appendix, I connect up the model-based, Maximum
# Likelihood inference we have studied, to Ordinary Least Squares
# (OLS) or linear regression, which you probably remember
# from introductory statistics. It turns
# out that a great deal of statistics is tied to likelihood,
# Maximum Likelihood, and model-based inference in one
# form or another, although often the connections are not
# always apparent in introductory courses.
#######################################################
#######################################################
# Likelihood of a case count
#######################################################
# In week 4, we introduced the concept of likelihood, where
# "likelihood" means the probability of the data, given a model.
# We can write this definition like this:
#
# L = P(data | model)
#
# The "|" means "given".
#
# Sometimes, we write the likelihood as a function of the model,
# and of the parameters. The list of parameters is often abbreviated
# as the Greek symbol "theta". We cannot use greek symbols in
# R scripts, which have to be plain-text, so we will just write
# "theta".
#
# L = P(data | model, theta)
#
# ...this just means that the likelihood (L) is the probability (P)
# of the data, given a model and particular parameter values stored
# in theta.
#
#
# As we saw in week 4, we can try to learn what parameter values
# maximize the probability of the data. This is the method of
# "Maximum Likelihood" (ML), and is one common method to attempt
# to estimate likely values of parameters. It is commonly used
# in hundreds of different scientific fields.
#
# Let's set up a ML analysis, in order to fit SIR model parameters
# to real data.
#######################################################
# NOTE: AS IN PREVIOUS WEEKS, IT IS CRUCIAL TO COPY-
# PASTE *ALL* OF THE CODE, WITHOUT SKIPPING STEPS.
#
# IT IS A BIT TEDIOUS, COMPARED TO E.G. LOADING FUNCTIONS
# FROM A WEBSITE, BUT AT LEAST IT KEEPS US ALL
# ON THE SAME PAGE, AND MAKES THE LAB PORTABLE
#
# In R on a computer, you can run the code in chunks.
#
# This should also be true on RStudio Cloud.
#
#######################################################
# We are using the deSolve library again; see Week 5 for install instructions
library(deSolve)
# Let's start with the New Zealand COVID-19 "active cases"
# data. I acquired this data from:
#
# Source: Hancock, Faith (2020). Covid-19 in NZ. NZ Herald, April 13, 2020.
# https://www.newsroom.co.nz/2020/03/27/1101216/covid-19-in-nz-the-numbers
#
# ...which gets its info from the Ministry of Health. These are
# basically the numbers that are announced by the government
# in a press conference at 1 pm each day.
#
# I have the data saved as a text file, and normally I would
# load it with:
################################################
## EXAMPLE CODE FOR LOADING A TEXT FILE IN R
################################################
## Set the working directory (wd) to wherever you have the file saved:
# getwd() # getwd = "get working directory", to check what
# working directory you are in
## You would type the full path for the directory containing your file here:
# wd = "/drives/GDrive/__classes/BIOSCI220/__week06/"
#
## Then, set R to make this your working directory
# setwd(wd) # setwd = "set working directory:
# getwd() # getwd = "get working directory", to check what
# working directory you are in
## Load the file into table "dat"
## fn = filename
# fn = "2020-04-13_NZ_daily_cases.txt"
# dat = read.table(file=fn, header=TRUE, stringsAsFactors=FALSE)
# dat
# dim(dat)
################################################
# WHY WE AREN'T TRYING TO LOAD FILES IN THIS WORKSHEET
#
# Getting everyone to save the file somewhere they can
# find it on their computer systems, then getting everyone to
# successfully reference that exact file path in their R
# script, is always a huge pain. This is why people say,
# "the hardest part about teaching R is getting students to
# successfully load their data."
#
# It would be even harder to do this while teaching remotely,
# so we'll "cheat", and just have the data pasted straight
# into the R script.
#
# In real life R analyses, you will have to put some effort
# into saving your data files into a directory, setting
# R's working directory with setwd(), and loading
# the data file with a function like read.table().
#
################################################
datafile_as_a_string = '
# Source: Hancock, Faith (2020). Covid-19 in NZ. NZ Herald, April 13, 2020.
# https://www.newsroom.co.nz/2020/03/27/1101216/covid-19-in-nz-the-numbers
Date Total_cases Active_cases
2020-02-28 1 1
2020-02-29 1 1
2020-03-01 1 1
2020-03-02 1 1
2020-03-03 1 1
2020-03-04 2 2
2020-03-05 3 3
2020-03-06 3 3
2020-03-07 5 5
2020-03-08 5 5
2020-03-09 5 5
2020-03-10 5 5
2020-03-11 5 5
2020-03-12 5 5
2020-03-13 5 5
2020-03-14 6 6
2020-03-15 8 8
2020-03-16 8 8
2020-03-17 12 12
2020-03-18 20 20
2020-03-19 28 28
2020-03-20 39 39
2020-03-21 52 52
2020-03-22 66 66
2020-03-23 102 102
2020-03-24 155 143
2020-03-25 205 183
2020-03-26 283 256
2020-03-27 368 331
2020-03-28 451 401
2020-03-29 514 457
2020-03-30 589 525
2020-03-31 647 572
2020-04-01 708 625
2020-04-02 797 704
2020-04-03 868 764
2020-04-04 950 822
2020-04-05 1039 882
2020-04-06 1106 929
2020-04-07 1160 918
2020-04-08 1210 927
2020-04-09 1239 921
2020-04-10 1283 908
2020-04-11 1312 886
2020-04-12 1330 855
2020-04-13 1349 798
2020-04-14 1366 729
2020-04-15 1386 649
2020-04-16 1401 622
2020-04-17 1409 582
'
# Here, we will load the data by reading in the text string we
# created above
dat = read.table(text=datafile_as_a_string, header=TRUE, stringsAsFactors=FALSE)
# Look at the full data:
# Look at the head/tail of the data
head(dat)
tail(dat)
# Looking at dat, please answer the questions below:
#
# QUESTION #1: How many rows are in the dataset loaded into 'dat'?
#
# QUESTION #2: How many columns are in the dataset loaded into 'dat'?
#
# QUESTION #3: What is the first day data are reported? Put your answer in
# YYYY-MM-DD format.
#
# QUESTION #4: What is the last day data are reported? Put your answer in
# YYYY-MM-DD format.
#
# Hint: You can also at the dimensions of dat, with dim().
# R *ALWAYS* refers to rows, then columns
dim(dat)
# Or, you can use "nrow" and "ncol"
nrow(dat)
ncol(dat)
#######################################################
# Let's plot this data
#######################################################
# Let's make x-values just be 1-46
xvals = 1:nrow(dat)
plot(x=xvals, y=dat$Total_cases, col="grey70", pch=19, xlab="Day", ylab="Number of cases")
points(x=xvals, y=dat$Active_cases, col="red", pch=19)
legend(x="topleft", legend=c("Total cases", "Active cases"), col=c("grey50", "red"), pch=c(19,19))
# QUESTION #5: Why are "Total cases" and "Active cases" the
# same at the start, but then diverge later in the plot?
# Distinguishing "classes" (different data types) in R
#
# It might be nice to make this plot have dates on the
# x-axis.
#
# What happens if we type this?
#
# plot(x=dat$Date, y=dat$Total_cases, col="grey70", pch=19, xlab="Day", ylab="Number of cases")
#
# (I am commenting this out, as it will "crash" peoples'
# R-Snippets runs. Uncomment it to try it.)
#
# The error message returned is:
#
# Error in plot.window(...) : need finite 'xlim' values
#
# This error occurs because the x-values are not numbers
# or dates, they are character strings representing
# dates, and character strings don't provide a way for
# R to determine the finite numbers determining the
# limits of the x-axis.
#
# Although the dates in dat$Date look like dates to us humans,
# R just sees them as character strings. In R, character strings
# have quote marks around them. You can see this by looking at a
# single value from the column:
dat$Date[1] # the first value in column "Date"
dat$Date[2] # the second value in column "Date"
dat$Date[3] # the third value in column "Date"
# We can do the plot with calendar dates. To do this,
# we have to convert the "Dates" column
# - from "character" format (meaning a string of text)
# - to "Date" format (meaning a number R will recognize as a calendar date)
# Convert the dates to R's "Date" format, then order d by date
#
# Functions like "as.character", and "as.Date", can convert between
# these formats (also known as data types or variable classes).
#
# You can check the format/class of a variable or column with the class() function.
# The initial class for column "Date"
class(dat$Date)
# QUESTION #6: What was the initial class of dat$Date? (put the exact text returned,
# i.e. the text inside the quotes, but not the quotes themselves)
dat$Date = as.Date(dat$Date, format="%Y-%m-%d")
# The new class for column "Date"
class(dat$Date)
# QUESTION #7: What was the new class of dat$Date? (put the exact text returned,
# i.e. the text inside the quotes, but not the quotes themselves)
# Now we can plot the data, with dates on the x-axis:
plot(x=dat$Date, y=dat$Total_cases, col="grey70", pch=19, xlab="Day", ylab="Number of cases")
points(x=dat$Date, y=dat$Active_cases, col="red", pch=19)
legend(x="topleft", legend=c("Total cases", "Active cases"), col=c("grey70", "red"), pch=c(19,19))
# Compare this plot to the "Total and active cases plot" available here:
# https://www.newsroom.co.nz/2020/03/27/1101216/covid-19-in-nz-the-numbers
#
# QUESTION #8: What is different about the NZ Herald plot, compared to
# your R plot?
#
# QUESTION #9: What would you have to do to make your R plot match the NZ Herald plot?
#
#######################################################
# Maximum likelihood inference of R0,
# pre- and post-lockdown
#######################################################
# You will recall the several steps required to do
# maximum likelihood inference:
#
# PART 1. Be able to calculate the likelihood of a single
# datum, given a model and parameter values.
#
# PART 2. Calculate the likelihood for all data points under
# a model+parameters.
#
# NOTE: For raw likelihoods, this means MULTIPLYING the likelihoods
# of the individual datapoints together. (In other words,
# getting the product of all of the individual likelihoods.)
# Remember the coin-flipping example: if the likelihood of
# a single heads is 0.7, then the likelihood of two Heads is
# 0.7 * 0.7, or 0.7^2.
#
# For log-likelihoods, we ADD UP the log-likelihoods of
# all of the individual data. This is true because adding
# logged likelihoods is equivalent to multiplying raw likelihoods.
#
# Proof of this:
#
# If Prob_heads is 0.7, then the raw likelihood of 2 Heads is:
0.7 * 0.7
# If Prob_heads is 0.7, then the log-likelihood of 2 Heads is:
log(0.7) + log(0.7)
# ...which is the same as...
log(0.7 * 0.7)
log(0.49)
# ...so, the log operation turns multiplication into addition.
#
# (The log operation also turns power operations into multiplication,
# so these are equivalent:
log(0.7^2)
2*log(0.7)
# ...remember, whenever you find yourself confused about some bit of
# math like this, just experiment with it in R to figure it out.)
#
# PART 3. Optimize the parameters in some way, in order
# to find the parameter values that maximize
# the likelihood.
#
# The "maximum likelihood parameter values" are the
# parameter values that maximize the probability of
# the observed data. They are the "best fit" parameters.
# The method of Maximum Likelihood is one very common,
# very simple method of trying to learn parameter
# values from the data. Learning parameters from the
# data is also called "estimation" or "inference".
# A vast amount of statistics is, deep down, based on
# Maximum Likelihood methods, even if they are
# (unfortunately) not usually emphasized in
# introductory classes.
#
#######################################################
# Here, we will use ML to infer R0 ("R-naught") for
# the NZ Covid epidemic.
#
# PLEASE NOTE that this is a practice exercise, using
# simplified models. "Real-life" epidemiology models
# would be somewhat more complex, and allow, for
# example, multiple different introductions of Covid
# cases on multiple different days.
#
# I think the models we are using are very instructive --
# I'm just saying that you shouldn't go around claiming
# you are a professional epidemiologist based on this
# exercise.
#######################################################
#######################################################
# Part 1. Be able to calculate the likelihood of a single
# datum, given a model and parameter values.
#######################################################
# Let's start with calculating the likelihood of a
# single data point.
#
# We will do this with the "likelihood_single_point()"
# function. Following Funk (2019), we will calculate
# the probability of a particular observed
# count with a Poisson distribution. The Poisson
# distribution is often used for count data,
# for several reasons:
#
# * It has only one parameter: lambda, which is equal
# to both the expected mean and the variance
# * It has a minimum of zero (so it is impossible
# to produce negative counts).
#
# So, all that likelihood_single_point() does is
# take a model projection for a particular day
# ("model.point"), set that as the lambda (expected
# mean) of the the dpois() function, and then get
# the probability of the observed case count
# ("data.point") for the day.
#
# The "log=TRUE" statement means that the
# log-likelihood is returned, instead of the
# regular or "raw" likelihood.
#
# Let's load this function, then use it to
# calculate a few likelihoods for single
# points.
# Function definition for: likelihood_single_point
# This function returns the likelihood or log-likelihood of a
# single day's case count (data.point), given a prediction
# for that day from an SIR model (model.point)
likelihood_single_point <- function(data.point, model.point, log=TRUE)
{
# Get the probability density of the count, assuming
# the observed count is a poisson process with an
# average value matching the true process.
#
# This allows for variation, e.g. if a case is reported
# days after it occurs.
likelihood = dpois(x=data.point, lambda=model.point, log=log)
# Return the likelihood value
return(likelihood)
} # END of likelihood_single_point
# Let's say that a particular set of model parameters
# predict that, on average, 10 active cases will
# be observed. What is the likelihood for an observation
# of exactly 10 cases?
# What is the likelihood for an observation
# of exactly 10 cases?
likelihood_single_point(data.point=10, model.point=10, log=FALSE)
# What is the log-likelihood for an observation
# of exactly 10 cases?
likelihood_single_point(data.point=10, model.point=10, log=TRUE)
# Now that you have seen how this code works, please answer the questions
# below.
# QUESTION #10: What is the log-likelihood for an observation of
# 8 active cases, if the model predicts 10 cases?
# QUESTION #11: What is the log-likelihood for an observation of
# 6 active cases, if the model predicts 10 cases?
# QUESTION #12: What is the log-likelihood for an observation of
# 4 active cases, if the model predicts 10 cases?
# QUESTION #13: Does the log-likelihood get higher or lower as
# the number of observed cases gets further from the model
# projection?
# NOTE: FOR QUESTIONS #14-15, MAKE SURE YOU READ AND UNDERSTAND
# THE NOTES ABOVE ON MULTIPLYING RAW LIKELIHOOD, VERSUS SUMMING
# LOG-LIKELIHOOD:
#
# "For raw likelihoods, this means MULTIPLYING the likelihoods
# of the individual datapoints together..."
# "For log-likelihoods, we ADD UP the log-likelihoods of
# all of the individual data..."
# QUESTION #14: What is the total log-likelihood of the 3
# observed counts (8, 6, and 4), assuming the model prediction
# was 10 for each?
# QUESTION #15: What is the total log-likelihood of the 3
# observed counts (8, 6, and 4), assuming the model prediction
# was 6 for each?
# QUESTION #16: Which model prediction yields a higher
# total log-likelihood for the 3 observed counts, the
# prediction of 6, or 10?
###############################################################
# Part 2. Calculate the likelihood for all data points under
# a model+parameters, and multiply all of these
# likelihoods together to get their product.
# (Or, more commonly, add together the logs of the
# likelihoods, to get their sum.)
###############################################################
# Rather than summing log-likelihoods by hand, we can just
# have another function. This way, we can take a series
# of observed counts, and a model prediction, and get the
# total log-likelihood easily.
#
# I have done this in likelihood_time_series(), which
# just runs the likelihood_single_point() function across
# each pair of data.points and model.points.
likelihood_time_series <- function(data.points, model.points, log=TRUE)
{
# Error check - return a stop() message if check fails
if (length(data.points) != length(model.points))
{ stop("ERROR in likelihood_time_series(): the lengths of data.points and model.points must match.") }
# Calculate the likelihoods of the data.points
likelihoods = mapply(FUN=likelihood_single_point, data.point=data.points, model.point=model.points, MoreArgs=list(log=log))
# If you calculated log-likelihoods, then sum them
# If you calculated raw likelihoods, then multiply them
if (log == TRUE)
{
total_likelihood = sum(likelihoods)
} else {
total_likelihood = prod(likelihoods) # prod=take the product
} # END if/else statement: if (log == TRUE)
# Return that value
return (total_likelihood)
} # END of likelihood_time_series
# We can use this function as another way to answer our question, above.
likelihood_time_series(data.points=c(8,6,4), model.points=c(10,10,10), log=TRUE)
# Now, in order to fit a series of observations to a
# model projection, we need the functions we used last
# week to make model projections. As before, load
# these functions by pasting them in ALL AT ONCE,
# i.e. from the function name, all the way through to the
# close of the function at the "}" character.
#######################################################
# Functions from Week 5
#######################################################
library(deSolve)
SIR_ode <- function(time, state, parameters)
{
# The above inputs are:
# time = a time point
# state = an R list of the 3 state values, i.e. the
# population sizes of S, I, and R
# parameters = an R list of the parameters R0 and D_inf
# Convert the input parameters "R0" and "D_inf"
# into the rates "beta" and "nu"
beta <- parameters[["R0"]] / parameters[["D_inf"]]
nu <- 1 / parameters[["D_inf"]]
# Extract the current values of the states S, I, and R
S <- state[["S"]]
I <- state[["I"]]
R <- state[["R"]]
# Calculate N, the total of S+I+R,
# i.e. the population size
N <- S + I + R
# Write out the system of
# Ordinary Differential Equations (ODEs)
dS <- -beta * I/N * S
dI <- beta * I/N * S - nu * I
dR <- nu * I
# Return the current rates of change of S, I, and R:
return(list(c(dS, dI, dR)))
} # End of function SIR_ode
simulate_SIR <- function(theta, init.state, times)
{
trajectory <- data.frame(ode(y = init.state,
times = times,
func = SIR_ode,
parms = theta,
method = "ode45"))
return(trajectory)
} # End of simulate_SIR()
simulate_SIR_changes <- function(thetas_states_df, times)
{
# Loop through the rows
trajectory = NULL
for (rn in 1:nrow(thetas_states_df))
{
if (rn == 1)
{
init.state = c(S=thetas_states_df$S[rn], I=thetas_states_df$I[rn], R=thetas_states_df$R[rn])
theta = list(R0=thetas_states_df$R0[rn], D_inf=thetas_states_df$D_inf[rn])
} else {
last_row = trajectory[nrow(trajectory),]
init.state = c(S=last_row$S+thetas_states_df$S[rn], I=last_row$I+thetas_states_df$I[rn], R=last_row$R+thetas_states_df$R[rn])
theta = list(R0=thetas_states_df$R0[rn], D_inf=thetas_states_df$D_inf[rn])
}
# Read in the parameters for the initial, versus later, time-bins
if (rn >= nrow(thetas_states_df))
{
start_time = thetas_states_df$time[rn]
end_time = max(times)
} else {
start_time = thetas_states_df$time[rn]
end_time = thetas_states_df$time[rn+1]
}
TF1 = times >= start_time
TF2 = times <= end_time
TF = (TF1 + TF2) == 2
tmp_times = times[TF]
tmp_trajectory = simulate_SIR(theta=theta, init.state=init.state, times=tmp_times)
if (rn == 1)
{
trajectory = rbind(trajectory, tmp_trajectory)
} else {
trajectory_minus_last_state = trajectory[-nrow(trajectory),]
trajectory = rbind(trajectory_minus_last_state, tmp_trajectory)
} # if (rn == 1)
} # END for (rn in 1:nrow(thetas_states_df))
return(trajectory)
} # End simulate_SIR_changes()
#######################################################
# Code to plot a model projection, observed counts,
# and the total log-likelihood
#######################################################
# Let's get our active case counts into a variable, "data.points"
data.points = dat$Active_cases
# ...and, let's make a model projection like we did last week:
times = 1:length(data.points)
time = c(1, 38) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case;
# Assume 24+14=38 days to see effect.
R0 = c(3.0, 0.3)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(1, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_df, times=1:length(data.points))
# Plot the results
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'Model projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
# Print the total log-likelihood in the plot title
# NOTE how we take the model projection of # infectious, trajectory$I, and
# store that in model.points:
model.points = trajectory$I
total_log_likelihood = likelihood_time_series(data.points=data.points, model.points=model.points, log=TRUE)
titletxt = paste0("Total log-likelihood = ", round(total_log_likelihood, 2))
title(titletxt)
#######################################################
# Now, that model projection (red curve) doesn't look very close to your
# observations (the circles), does it? Can we do better?
#
# Let's spend a few minutes changing the two "R0" values, and/or the
# first "I" value, and re-running the code above.
#
# For example, I've made another attempt:
#
#######################################################
# Code to plot a model projection, observed counts,
# and the total log-likelihood
#######################################################
# Let's get our active case counts into a variable, "data.points"
data.points = dat$Active_cases
# ...and, let's make a model projection like we did last week:
times = 1:length(data.points)
time = c(1, 38) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case;
# Assume 24+14=38 days to see effect.
R0 = c(3.5, 0.4)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(2, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_df, times=1:length(data.points))
# Plot the results
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'Model projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
# Print the total log-likelihood in the plot title
# NOTE how we take the model projection of # infectious, trajectory$I, and
# store that in model.points:
model.points = trajectory$I
total_log_likelihood = likelihood_time_series(data.points=data.points, model.points=model.points, log=TRUE)
titletxt = paste0("Total log-likelihood = ", round(total_log_likelihood, 2))
title(titletxt)
#######################################################
# QUESTION #17: Spend about 10 minutes re-running the
# code chunk above, changing the two "R0" values,
# and/or the first "I" value.
#
# Answer these questions in the answer box:
#
# What is the highest log-likelihood you found?
#
# What parameter values did you use to find this
# log-likelihood (lnL)?
#
# How hard was it to find parameter values that
# made a curve that matched the data?
#
# I will not mark down for finding the "wrong" answer, this
# is about getting an intuition about how
# changing parameters -> changes model projections ->
# changes fit between line & data -> changes the
# log-likelihood.
#
# Report your answer like this:
# The highest log-likelihood (lnL) I found was lnL = -12500.74
#
# What parameter values did you use to find this
# log-likelihood (lnL)?
#
# R0 = c(3.0, 0.3)
# I = c(1, 0)
#
# How hard was it to find parameter values that
# made a curve that roughly matched the data?
#
# I found it very easy/very difficult to find
# the parameter values because...
#
##########################################################
# Part 3. Optimize the parameters in some way, in order
# to find the parameter values that maximize
# the likelihood.
##########################################################
# Hopefully, you just saw that "optimization by hand"
# is pretty difficult, although possible if you
# had no other options.
#
# Fortunately, we have other options!
#
# The function below, "params_to_likelihood_v1()", will help us
# automate this process of finding good parameter values, by:
#
# (1) taking the free parameters you want to infer ("params")
# (2) loading them into simulate_SIR_changes() to get a trajectory (model.points)
# (3) calculating the likelihood of the observed data.points against the model.points
#
# It takes parameters, simulates the model projection,
# then calculates the total log-likelihood - this just
# combines everything we have done separately above!
#
# It takes as inputs:
#
# params - starting values for the parameters that you want to infer
# data.points - the observed active cases, as before
# thetas_states_df - the data.frame holding time,R0, D_inf, S, I, R
# delay_val - if you think there is a delay between an intervention,
# and seeing the effects, put it here. Default is 14.
# printvals - if TRUE, prints the parameter values you are trying.
#
params_to_likelihood_v1 <- function(params, data.points, thetas_states_df, delay_val=14, printvals=TRUE, makePlots=FALSE)
{
# The times are just the length of the data.points
times = 1:length(data.points)
# (1) input the free parameters
thetas_states_temp = thetas_states_df
TF = thetas_states_temp == "free"
thetas_states_temp[TF] = params
# Error check
if ((max(thetas_states_temp$time)+delay_val) > max(times))
{
stop("ERROR in params_to_likelihood_v1: there is a 'thetas_states_df$time'+delay_val greater than found in 'times'")
}
# Make sure everything in thetas_states_df is numeric, instead of character
for (i in 1:ncol(thetas_states_temp))
{
thetas_states_temp[,i] = as.numeric(thetas_states_temp[,i])
}
# Ad2 the delay (allows time for interventions to show up
# in the detections)
thetas_states_temp$time[2:length(thetas_states_temp$time)] = thetas_states_temp$time[2:length(thetas_states_temp$time)] + delay_val
# Calculate the trajectory given the parameters,
# input into model.points
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_temp, times)
model.points = trajectory$I
lnL = likelihood_time_series(data.points, model.points, log=TRUE)
if (is.finite(lnL)== FALSE)
{
lnL=-10e100
}
# Print, if desired
if (printvals == TRUE)
{
print(thetas_states_temp)
cat("log-likelihood = ", lnL, "\n", sep="")
}
# Make plots, if desired
if (makePlots == TRUE)
{
plot_trajectory_lines(trajectory, theta=NULL, thetas_states_df=thetas_states_df, title_prefix="Plot #5:")
}
return(lnL)
} # END params_to_likelihood_v1
# We can use params_to_likelihood_v1(), in conjunction with with
# the ML optimization function optim() (we used this in Week 4!),
# to try lots of parameter values, and gradually find the
# parameters that maximize the likelihood.
#
# Let's try it!
#
# The major difference from what we did above is that, after
# setting up thetas_states_df, we write "free" in certain
# cells of the table. These "free" parameters are then filled
# in with the values from "params" -- and the optim() function
# will change the values in params, in an attempt to find the
# best parameters.
#
#######################################################
# Fitting a 2-regime model
#######################################################
time = c(1, 24) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case
R0 = c(3.0, 0.3)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(1, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# Make some parameters into "free" parameters, to be inferred
thetas_states_df$R0 = c("free", "free")
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(2.5, 1.0, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, delay_val=14, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 14
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, 3-regime model\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# Congratulations! You have run a Maximum Likelihood inference
# on New Zealand's Covid-19 active case count data.
# Using the plot, and the variables & screen output above, please answer
# the questions below.
#
# QUESTION #18. What was the log-likelihood obtained for the
# Maximum Likelihood parameters? Provide the answer to 2
# decimal places.
#
# QUESTION #19. Was this lnL higher than the one you found
# when you searched the parameter space "by hand"?
#
# QUESTION #20. During the ML search, the thetas_states table,
# and the log-likelihood, were printed to screen at each
# step of the search. Which option best describes what this
# looks like?
#
# NOTE: FOR QUESTIONS 21-24, AND ANY OTHER QUESTIONS ASKING ABOUT
# THE "first", "second", etc. R0 parameters, see lecture slide
# #101. And this example:
#
# If your thetas_states_ML table looks like this:
#
# time R0 D_inf S I R
# 1 1 2.574564 14 5e+06 1.028189 0
# 2 38 1.007426 14 0e+00 0.000000 0
#
# ...then:
# the 1st R0 parameter is: 2.574564
# the 2nd R0 parameter is: 1.007426
# the starting number of Infectious is: 1.028189
#
#
# QUESTION #21. What was the first R0 parameter, as inferred by ML?
# (This parameter was operating from Days 1-36.)
#
# QUESTION #22. What was the second R0 parameter, as inferred by ML?
# (This parameter was operating from Day 37 to the end.)
#
# QUESTION #23. What was the starting number of Infected, as inferred by ML?
# (This parameter was operating from Day 37 to the end.)
#
# QUESTION #24. Is this "starting number of Infected" a plausible number in
# real life, or is this more of a "nuisance parameter" in our analysis?
#
# You may have noticed that even our ML-fitted model parameters
# did not produce an excellent match to our data points.
#
# For example, the projected curve first rises faster than
# our observed case counts, but then it gets lower than the
# counts, and then it gets higher than the counts again.
#
# This is probably the result of trying to fit a single
# R0 parameter, which will produce a curve increasing
# exponentially by a constant proportion, to a dataset
# where, possibly, multiple things are going on.
#
# For example, what if, before the Level 4 lockdown happened,
# New Zealanders were already reducing transmission of the virus,
# for example through:
#
# * Level 2 and Level 3 restrictions
# * Increased handwashing & use of sanitizer
# * Stopping handshaking etc.
#
# This suggests that we might end up with a better-fitting
# model, and better parameter estimates, if we fit a
# 3-regime model, instead of a 2-regime model.
#
# Let's try that!
#######################################################
# Fitting a 3-regime model
#######################################################
time = c(1, 16, 24) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case
R0 = c(3.0, 1.0, 0.3)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(1, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# Make some parameters into "free" parameters, to be inferred
thetas_states_df$R0 = c("free", "free", "free")
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(2.5, 1.0, 0.3, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 14
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, 3-regime model\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# QUESTION #25. Visually, in the 3-regime model, does it look like
# this curve fits the data better?
#
# QUESTION #26. In the 3-regime model, what was the log-likelihood obtained
# for the Maximum Likelihood parameters? Provide the answer to 2
# decimal places.
#
# QUESTION #27. In the 3-regime model, is the lnL higher than the one you found
# for the 2-regime model?
#
# QUESTION #28. What was the first R0 parameter, as inferred by ML?
# (This parameter was operating from Days 1-30.)
#
# QUESTION #29. What was the second R0 parameter, as inferred by ML?
# (This parameter was operating from Days 31-38.)
#
# QUESTION #30. What was the third R0 parameter, as inferred by ML?
# (This parameter was operating from Day 39-end.)
#
# QUESTION #31. What was the starting number of Infected, as inferred by ML?
# (note: as we saw before, this is a nuisance parameter, not intended
# as a real estimate in this simplified model)
#
# QUESTION #32. What do our inferred R0 values suggest about the
# coronavirus's natural spreading ability, and about the impact of
# sanitation and physical distancing measures?
#
#######################################################
# Statistical model comparison
#
# One huge advantage of likelihood-based inference is
# that the likelihoods can be used to statistically
# compare model fits.
#
# One very common method for comparing ML-fitted models
# is Akaike Information Criterion (AIC). AIC has been
# used in tens of thousands of research papers to
# compare multiple model fits. See the lecture slides/
# audio, and the Week 6 reading, for an introduction to
# AIC.
#
# The very-very short version of AIC is that AIC
# *estimates the relative distance from the truth*.
#
# We never know what the "true model" is in real life,
# since any true biological history that produced observations
# is preposterously complex. However, we can at least
# compare the *relative* fit of models to data, using AIC.
#
# AIC is a simple function of likelihood:
#
# AIC = -2 * (lnL - nparams)
#
# lnL = the maximized log-likelihood under some model
# nparams = the number of free parameters that you fit with ML
#
# Basically, the AIC is a form of penalized likelihood: we
# favor models that give higher lnLs, but then we penalize
# models base on the number of free parameters that we fit.
#
# This helps enforce parsimony in our explanations. If we
# put no penalty on extra free parameters, we could always
# fit models with higher and higher likelihood, until we
# ended up with a free parameter for every data point, at
# which point we fit all the data perfectly, but we have
# explained nothing at all, because the model just restates
# the dataset.
#
# Some details:
#
# The "-2" in the AIC formula needs explanation. The "-"
# just converts the log-likelihood (which is always negative,
# if the data are discrete) from negative to positive. This
# makes sense because AIC is a distance, and distances
# should be positive.
#
# The "2" part in AIC is actually an arbitrary convention -
# it would work the same if we used 1, 3, etc. The number "2"
# was chosen because another famous statistical test, the
# Likelihood Ratio Test, uses 2*(difference in log-likelihoods)
# as the test statistic, often called the "deviance".
#
# The rationale for 1 parameter = 1 unit of log-likelihood
# is complex, but it can be derived from a view of the world
# where we think that processes are controlled by a series of
# parameters of decreasing importance, such that there are
# a few major controlling variables, and an exponential
# decay in the importance of other controlling variables.
#
# Other rationales can produce other penalties, creating
# other model comparison criteria, such as BIC (Bayesian
# Information Criterion). Don't worry about these for now.
#
# AIC model weights
#
# AIC is just a number that doesn't mean much by itself.
# AIC values are used to *compare* multiple models. We
# do this with AIC weights, which assign each model a
# weight, where all of the weights sum to 1.
#
# For example, if you had these weights:
#
# Model 1 AIC weight = 0.05
# Model 2 AIC weight = 0.9
# Model 3 AIC weight = 0.05
#
# You would say that "Model 2 acquired 90% of the AIC
# weight, and was the best-fitting model in our
# model set. It fits about 9 times better than the
# other models combined, which is moderately good
# evidence that it is a better fit. Models 1 and 3
# have equal model weight, so they are approximately
# equivalent fits to the data."
#
##############################################
# Here is how to calculate AIC weights:
##############################################
#
# 1. Calculate AIC values for all the models
# in your model set.
#
# 2. Find the AIC-best model. This is the model
# with the *LOWEST* AIC (i.e., the smallest
# relative distance from the truth).
#
# 3. For all the models, calculate the AIC difference
# from the best model. (The AIC difference between
# the best model and itself is 0.0.). This
# number is called "deltaAIC" (in science, "delta"
# often means "change" or "difference", so "deltaAIC"
# just means "difference in AIC".
#
# 4. For each model, calculate exp(-0.5 * deltaAIC). These
# are the relative likelihoods.
#
# 5. Sum the relative likelihoods, and divide each
# individual relative likelihood by the sum of
# relative likelihoods. This is the AIC weight.
# The AIC weights will sum to 1.0.
#
# 6. If desired, convert these fractions (numbers
# between 0 and 1) to percentages, by multiplying
# by 100, or by clicking the "%" option in Excel.
#
###############################################
###############################################
# Advantages and disadvantages of AIC
###############################################
# The advantages of the AIC framework include:
#
# * Multiple models (more than 2) can be compared easily.
# Traditional statistical tests, which give you p-values,
# typically only compare two models at a time, namely
# a null model and an alternate.
#
# * There is no concept of "p-value cutoffs" or
# "statistical signficance" in AIC-based analyses.
# The proponents of AIC acknowledge that evidence is
# continuous. If AIC favours a model 1 over model 2
# 60%-40%, this is a tiny bit of evidence that the
# model 1 is closer to the truth. If AIC favours a
# model 1 99.999% to 0.001%, then this is strong
# evidence that model 1 closer to the truth.
#
# * AIC assumes that we never have access to the
# complete, true model that produced our data. It
# acknowledges that the "true" model behind most
# empirical data would be fantastically complex. So
# all we are really trying to do is find better
# simple approximations of a very complex reality.
# AIC is one easy way to measure the *relative*
# distance of models from the truth. We never know
# the "absolute" distance from the truth, because
# that would require knowing the true model, which
# we never do.
#
# * Assuming you have maximized lnL values for each
# model, the formula for AIC, and AIC weights, is very
# simple. The calculations can be done by hand, or
# in Excel or Google Sheets (or R). All you have to
# be able to do is count the number of free parameters
# in each model.
#
# * AIC can be used very broadly for comparing models.
# This includes not just models where the lnL is
# explicitly reported (like those above), but
# also in linear regression and other linear models.
# Here, model fit is more commonly reported in
# terms of "deviance" or "RSS" (residual sum
# of squares, the sum of the squared errors between
# the line of best fit, and the data). These are
# all proportional to -1 * log-likelihood, under
# the assumption that the errors are normally
# distributed. So, the "least squares" method can
# also be interpret as a Maximum Likelihood method!
#
# (Unfortunately, many introductory students only
# ever learn about least-squares, when Maximum
# Likelihood has much broader application to
# all sorts of models, not just cases where
# a normal distribution of errors can be assumed.)
#
# AIC has some limitations as well:
#
# * If you are literally doing a randomized controlled
# experiment, where you have two hypotheses (models)
# that exhaust the space of models, namely "no effect"
# (the null model) and "effect" (the alternative model),
# then the traditional p-value framework can make a
# lot of sense. This is especially the case where
# further decisions will be based on the result, so
# the effect of false-positives and false-negatives
# needs serious consideration.
#
# * If you *do* have the true model in your set, for
# example when your data come from a computer simulation
# where you know the true model, then a criterion
# like BIC, which has a stronger penalty against
# extra parameters.
#
# * If you have a dataset where the number of data is
# small (typically <30 data points), and/or the number
# of free parameters is large, it is more appropriate
# to use "AICc", which is sample-size corrected AIC.
# The formula for this is a little more complex, so
# we will ignore it here.
#
# * There are all kinds of intense philosophical debates
# in statistics and philosophy of science about
# models, inference, and the best way to measure
# support for models/hypotheses. AIC is just one
# choice amongst these. I would put it roughly in
# the "likelihoodist" school. Other major schools of
# thought are "frequentism" (which includes p-values
# and worries about the long-term error rate) and
# Bayesianism (which explicitly puts prior probability
# distributions on parameter values, and on models.)
#
# * It is easy forget that "all models are wrong", and
# get very attached to your best model. This is a
# challenge in all of statistics & data analysis,
# however. The famous statement from George Box is
# "All models are wrong, but some models are useful."
# Constant critical thinking is required, as I have
# tried to encourage!
#######################################################
#######################################################
# EXAMPLE USE OF AIC TO COMPARE MODELS
#######################################################
# Let's compare four models for our NZ Covid-19
# active cases dataset. Here, we run the 2-regime
# and 3-regime models we used before, but I also
# include a 1-regime model, and a 4-regime model.
#######################################################
# Model 1: A 1-regime model
#######################################################
time = c(1, 49) # Day 49 is the 2nd-to-last day, days 49-50 will have
# a fixed R0 of 2.5
R0 = c(3.0, 2.5)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(1, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# NOTE: ONLY 2 FREE PARAMETERS
thetas_states_df$R0[1] = "free"
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(2.5, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, delay_val=0, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 0
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, Model 1 (a 1-regime model)\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# Save this model's parameters and log-likelihood
thetas_states_ML_model1 = thetas_states_ML
total_lnL_Model1 = ML_results$value
#######################################################
# Model 2: A 2-regime model
#######################################################
time = c(1, 24) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case
R0 = c(3.0, 0.3)
D_inf = c(14.0, 14.0)
S = c(5000000, 0)
I = c(1, 0)
R = c(0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# Make some parameters into "free" parameters, to be inferred
thetas_states_df$R0 = c("free", "free")
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(3.0, 0.3, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, delay_val=14, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 14
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, Model 2 (a 2-regime model)\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# Save this model's parameters and log-likelihood
thetas_states_ML_model2 = thetas_states_ML
total_lnL_Model2 = ML_results$value
#######################################################
# Model 3: A 3-regime model
#######################################################
time = c(1, 16, 24) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case
R0 = c(3.0, 1.0, 0.3)
D_inf = c(14.0, 14.0, 14.0)
S = c(5000000, 0, 0)
I = c(1, 0, 0)
R = c(0, 0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# Make some parameters into "free" parameters, to be inferred
thetas_states_df$R0 = c("free", "free", "free")
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(3.0, 2.0, 0.3, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, delay_val=14, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 14
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, Model 3 (a 3-regime model)\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# Save this model's parameters and log-likelihood
thetas_states_ML_model3 = thetas_states_ML
total_lnL_Model3 = ML_results$value
#######################################################
# Model 4: A 4-regime model
#######################################################
time = c(1, 16, 20, 24) # NZ Level 4 lockdown went in on March 26, day 24 since 1st case
R0 = c(3.0, 1.0, 1.0, 0.3)
D_inf = c(14.0, 14.0, 14.0, 14.0)
S = c(5000000, 0, 0, 0)
I = c(1, 0, 0, 0)
R = c(0, 0, 0, 0)
thetas_states_table = cbind(time, R0, D_inf, S, I, R)
thetas_states_df = as.data.frame(thetas_states_table, stringsAsFactors=FALSE)
# Make some parameters into "free" parameters, to be inferred
thetas_states_df$R0 = c("free", "free", "free", "free")
thetas_states_df$I[1] = "free"
thetas_states_df
# Run the params_to_likelihood_v1() function once, to see your starting likelihood
data.points = dat$Active_cases
params = c(3.0, 2.0, 1.0, 0.3, 0.03) # starting parameter values
params_to_likelihood_v1(params, data.points, thetas_states_df, delay_val=14, printvals=TRUE, makePlots=FALSE)
# Running the Maximum Likelihood search with the optim() function.
# LOOK AT THE OUTPUT THAT PRINTS TO SCREEN!!
delay_val = 14
ML_results = optim(par=params, fn=params_to_likelihood_v1, data.points=data.points, thetas_states_df=thetas_states_df, delay_val=delay_val, printvals=TRUE, method="L-BFGS-B", lower=0.0, control=list(fnscale=-1))
# Graphing the ML model
# Take the learned parameters from "ML_results", put them
# into a theta_states data.frame for simulation and plotting
thetas_states_ML = thetas_states_df
TF = thetas_states_ML == "free"
thetas_states_ML[TF] = ML_results$par
for (i in 1:ncol(thetas_states_ML))
{ thetas_states_ML[,i] = as.numeric(thetas_states_ML[,i]) }
thetas_states_ML$time[2:length(thetas_states_ML$time)] = thetas_states_ML$time[2:length(thetas_states_ML$time)] + delay_val
# Plot the results
trajectory = simulate_SIR_changes(thetas_states_df=thetas_states_ML, times=1:length(data.points))
maxy = max(max(trajectory$I), max(data.points))
xvals = c(trajectory$time, times)
yvals = c(trajectory$I, data.points)
plot(x=xvals, y=yvals, pch=".", col="white", xlim=c(0, max(trajectory$time)), ylim=c(0, maxy), xlab="Day", ylab="Number of individuals")
lines(x=trajectory$time, y=trajectory$I, lwd=3, col="firebrick2")
points(times, data.points)
legend(x="topleft", legend=c("Active COVID-19 case count", 'ML-fitted projection of "I" (Infected)'), lty=c("blank", "solid"), lwd=c(1,3), pch=c(1, "."), col=c("black","firebrick2"), cex=0.8)
titletxt = paste0("ML fit, Model 4 (a 4-regime model)\nTotal log-likelihood = ", round(ML_results$value, 2))
title(titletxt)
# Save this model's parameters and log-likelihood
thetas_states_ML_model4 = thetas_states_ML
total_lnL_Model4 = ML_results$value
#######################################################
# AIC comparison of the 4 models
#######################################################
# Below, I will do AIC comparison in R, but I will
# also provide the same example in Excel and in
# Google Sheets.
#
# (NOTE: The exact results you get may differ, if you
# change the input data e.g. by adding days.)
#
# Excel: # http://phylo.wikidot.com/local--files/biosci220:quantitative-biology/Week6_AIC_with_SIRs_fit_to_NZ_Ministry_of_Health_data.xlsx
#
# Google Sheet:
# https://docs.google.com/spreadsheets/d/1WNk1la1TfoIdKcadxm0kUW5JNTTl0ZDImzDEQq0fdSQ/edit?usp=sharing
#
# Making the AIC table in R:
#
# First, let's get the number of free parameters
# for each of the 4 models. We get this by counting
# the number of "free" statement in the
# theta_states_df, or more simply by just thinking
# about it.
#
# Model 1 estimated 1 R0 parameter, and a starting "I" count
# That is 2 free parameters estimated by ML.
#
# Model 2 added another R0, Model 3 added another, and
# Model 4 added another. So they have 3, 4, and 5
# parameters.
# Let's store the number of free parameters
num_free_params = c(2, 3, 4, 5)
# Let's collect the maximized log-likelihoods from
# above
lnL = c(total_lnL_Model1, total_lnL_Model2, total_lnL_Model3, total_lnL_Model4)
# AIC is just -2 * (lnL - num_free_params)
AIC_vals = -2*(lnL - num_free_params)
# Get the delta AICs -- the difference from the best (lowest!) AIC
best_AIC = min(AIC_vals)
delta_AICs = AIC_vals - best_AIC
# Calculate the relative likelihoods
relative_likelihoods = exp(-0.5 * delta_AICs)
# And, calculate the AIC weights
AIC_weight = relative_likelihoods / sum(relative_likelihoods)
# AIC weights as percentages (rounded to 2 decimals)
percent_weight = round(AIC_weight * 100, digits=2)
# Make model names, and results table
model_names = c("Model 1", "Model 2", "Model 3", "Model 4")
AIC_results_matrix = cbind(lnL, num_free_params, AIC_vals, AIC_weight, percent_weight)
AIC_results_df = as.data.frame(AIC_results_matrix, stringsAsFactors=FALSE)
# Make better row & column names
row.names(AIC_results_df) = model_names
names(AIC_results_df) = c("lnL", "nparams", "AIC", "AIC_weight", "pct_weight")
# Look at the results table!
AIC_results_df
# A "real" results table would often also include the parameter values
parameter_estimates = matrix(0, nrow=4, ncol=5)
# It takes some work to extract the parameter estimates from the
# "thetas_states_ML_model" objects where we saved them
parameter_estimates[1,1:5] = c(thetas_states_ML_model1$I[1], thetas_states_ML_model1$R0[1], thetas_states_ML_model1$R0[1], thetas_states_ML_model1$R0[1], thetas_states_ML_model1$R0[1])
parameter_estimates[2,1:5] = c(thetas_states_ML_model2$I[1], thetas_states_ML_model2$R0[1], thetas_states_ML_model2$R0[1], thetas_states_ML_model2$R0[1], thetas_states_ML_model2$R0[2])
parameter_estimates[3,1:5] = c(thetas_states_ML_model3$I[1], thetas_states_ML_model3$R0[1], thetas_states_ML_model3$R0[1], thetas_states_ML_model3$R0[2], thetas_states_ML_model3$R0[3])
parameter_estimates[4,1:5] = c(thetas_states_ML_model4$I[1], thetas_states_ML_model4$R0[1], thetas_states_ML_model4$R0[2], thetas_states_ML_model4$R0[3], thetas_states_ML_model4$R0[4])
parameter_estimates
# Add them to the table, add some nice rounding
first_two_columns = cbind(round(AIC_results_df[,1],1), AIC_results_df[,2])
middle_columns = cbind(round(parameter_estimates, 2))
last_columns = round(AIC_results_df[,3:ncol(AIC_results_df)], 2)
AIC_results_df2 = cbind(first_two_columns, middle_columns, last_columns)
names(AIC_results_df2) = c("lnL", "nparams", "I_ini", "R0_1", "R0_2", "R0_3", "R0_4", "AIC", "AIC_weight", "pct_weight")
# Look at the final table
AIC_results_df2
#######################################################
# Congratulations! You have now calculated AIC model
# weights for your 4 models! You will see AIC tables
# like this in the results sections of many scientific
# papers.
#
# The main interpretation of the table is about the
# AIC weights or percent weights (which give the
# same information).
#
# Model 1 (R0=2.18 for the entire epidemic)
# is clearly a horrible fit. It gets 0.0% of
# the model weight, and the plot also shows a
# huge difference between the data and the best-fit
# curve. Yay! This shows that NZ's epidemic clearly
# deviates from constant exponential growth! We are
# nowhere near the millions of cases required for
# herd immunity to take over, so we're not just
# flattening the curve, we are squashing the curve!
#
# Model 2 (R0=3.5 until day 24, then R0=0.3 after the
# lockdown) has a fit that is 4918.75 AIC units
# better (meaning lower) than Model 1. An AIC difference of 10
# translates to a model having about 150 times more
# likelihood (exp(0.5*10) = about 150), so a difference
# of thousands of AIC units is huge. So this is
# a pretty good result, for adding just 1 more
# parameter!
#
# But was Model 2 the best fit? We can see that, in
# Model 3, adding a third R0, allowing R0 to decrease
# before the lockdown starts, improves the fit by
# another 451 AIC units. This is an increase in
# relative likelihood of exp(0.5*451)=8.6*10^97 times.
# Model 3 has an initial R0 of 4.85, then a switch to
# R0=2.53 on Day 16, then a switch to R0=0.53 on Day 24.
#
# Can we do even better? Model 4 further subdivides
# Days 16-20 from Days 20-24. This *does* increase
# the log-likelihood, from -207.74 (Model 3) to
# -207.38 (Model 4). However, does this improvement
# in log-likelihood outweigh the penalty for having
# an extra parameter? We know that the penalty term
# is equivalent to 1 unit of log-likelihood, so a
# log-likelihood improvement of only 0.36 units is
# not going to overcome the penalty.
#
# The AIC weights allow us to quantify this effect.
# Model 3 gets 65.49% of the model weight, and Model 4
# gets 34.51%. The research would conclude that, of the
# models considered, Model 3 is the best fit (or at
# least, the least-worst fit). Model 4 retains
# plausibility, because it is always possible that a
# more complex model is true, but we can say that Model 3
# is an equivalent match to the data, but is more
# parsimonious.
#
#######################################################
#######################################################
# What is the point of all of this?
#
# In biology, we almost never know the "true model". We
# are just trying to find better approximations of
# reality. AIC is one common choice to balance fitting
# the observations against parsimony.
#
# AIC can be used to address common questions like
# "which model is the best, given the observations I
# have", or "out of this set of models, which models
# are reasonable, given the dataset"?
#
# AIC can be used any time you can get a maximized
# log-likelihood, and a count of the number of free
# parameters, on each model. It is used so often that
# R has some standard functions to provide them for
# linear models provided by the lm() function.
# Note: these functions will not work
# outside of standard "base" R functions like lm() and
# glm().
#
# Here is a quick example of that:
data("iris") # A default R dataset, measurements of iris flowers
plot(iris$Sepal.Length, iris$Petal.Length)
# Run a simple linear model (lm), with Sepal Length
# predicting Petal Length
lm_result = lm(iris$Petal.Length~Sepal.Length, data=iris)
# The "logLik" function reports the lnL, and the
# number of free parameters (df=3)
# (df means "degrees of freedom", i.e. number of free
# parameters)
logLik(lm_result)
AIC(lm_result)
# You can see that the AIC formula is being used:
lnL = as.numeric(logLik(lm_result))
-2*(lnL - 3)
AIC(lm_result)
#
# MAIN POINT: The skills you have used for using
# Maximum Likelihood and AIC on SIR models can be
# used anywhere from quite complex process-based
# models used in epidemiology, ecology, phylogenetics
# etc., all the way to the basic linear models used
# in introductory statistics.
#
#######################################################
################################################
# Appendix: The connection between Ordinary
# Least Squares (OLS), i.e. linear regression,
# and Maximum Likelihood.
#
# This is something I wish I had been taught
# early in my statistics education, so I am
# mentioning it here.
#
# This Appendix material is "bonus", I will not
# examine you on it for my section of the course.
#
################################################
#
# In standard "least squares" analyses
# (OLS, Ordinary Least Squares), where
# we talk about the line of
# best fit being the one that minimizes
# the sum of squared residuals (or RSS, residual
# sum of squares), this is just another way
# of talking about maximizing the log-likelihood,
# *if* we make the assumption that the residuals
# (the difference between the predicted line and
# the observed response) all follow the same
# normally distribution (in other words, the
# residuals are independent, normally distributed
# with constant variance).
#
# In this situation, the log-likelihood for the
# line of best fit is derived from the log of
# the probability density function for a normal
# distribution.
#
# I find that this is rarely explained in R code,
# so I am putting various formulations here.
#
# Here is the logLik() function code:
getAnywhere("logLik.lm")
# This is the formula from logLik, with weights (w)
# set to 1s (meaning equal weighting of all residuals):
0.5 * - N * (log(2 * pi) + 1 - log(N) +
log(sum(lm_result$residuals^2)))
# Here are several simpler formulations that
# give the same result:
#
# Formula for the log-likelihood of a linear model
# with x predicting y
# Source:
# https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf
# ...page 2, equation (3)
#
# N = number of observations
#
# RSS = residual sum of squares =
# sum of (yi - (b0 + b1*xi))^2
# where
# yi = ith response (observation)
# b0 = intercept
# b1 = slope
# xi = ith predictor
#
# stdev_squared = variance = RSS / N
# stdev = standard deviation = sqrt(stdev_squared) = sqrt(variance)
N = length(lm_result$residuals)
RSS = sum(lm_result$residuals^2)
stdev_squared = RSS/N
stdev = sqrt(RSS/N)
-N/2 * log(2*pi) - N*log(stdev) - 1/(2*stdev_squared)*RSS
# The formula works the same if you just use RSS
-N/2 * log(2*pi) - N*log(sqrt(RSS/N)) - 1/(2*RSS/N)*RSS
# R's logLik() function has the above formula re-arranged
# so that RSS appears only once:
0.5 * - N * (log(2 * pi) + 1 - log(N) +
log(RSS))
logLik(lm_result)
# These formulations all give the same result!
#######################################################
# References
#
# Cosma Shalizi (2015). "Lecture 6: The Method of
# Maximum Likelihood for Simple Linear Regression."
# Lecture 6 from Modern Regression, 36-401, Fall 2015.
# https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf
#
#######################################################
#######################################################
# References
#
# My SIR code relies heavily on code from the fitR package
# (Funk et al. 2019).
#
# Adam, David (2020). "UK has enough intensive care units for
# coronavirus, expert predicts." New Scientist, March 25, 2020.
# https://www-newscientist-com.ezproxy.auckland.ac.nz/article/2238578-uk-has-enough-intensive-care-units-for-coronavirus-expert-predicts/
#
# Anderson, Charles (2020). "New Zealand coronavirus deaths during lockdown could
# be just 20, modelling suggests." The Guardian. 27 Mar 2020.
# https://www.theguardian.com/world/2020/mar/27/new-zealand-coronavirus-deaths-during-lockdown-could-be-just-20-modelling-suggests
#
# Dillane, Tom; Knox, Chris (2020), "Coronavirus Covid 19: New Zealand's
# best and worst prepared DHBs by elderly population and ICU beds."
# NZ Herald. April 5, 2020
# https://www.nzherald.co.nz/nz/news/article.cfm?c_id=1&objectid=12321390
#
# Fefferman, Nina (2020). "The role of applied math in real-time pandemic
# response: How basic disease models work." Date: 3:30 EDT Tuesday,
# March 31, 2020.  http://www.nimbios.org/webinars#fefferman 
#
# Funk, Sebastian (2019). "Introduction to model fitting in R." http://sbfnk.github.io/mfiidd/introduction.html
#
# Funk, Sebastian; Camacho, Anton; Johnson, Helen; Minter, Amanda; O’Reilly, Kathleen; Davies, Nicholas (2019). "Model fitting and inference for infectious disease dynamics." Centre for the Mathematical Modelling of Infectious Diseases (CMMID), London School of Hygiene & Tropical Medicine.
#
# James, Alex; Hendy, Shaun C.; Plank, Michael J.; Steyn, Nicholas (2020). "Suppression and Mitigation
# Strategies for Control of COVID-19 in New Zealand." 25 March 2020.
# https://cpb-ap-se2.wpmucdn.com/blogs.auckland.ac.nz/dist/d/75/files/2017/01/Supression-and-Mitigation-New-Zealand-TPM-006.pdf
#
# Newton, Kate (2020). "The man modelling NZ's Covid-19 spread
# from his kitchen table." Radio New Zealand. 27 March 2020.
# https://www.rnz.co.nz/news/in-depth/412744/the-man-modelling-nz-s-covid-19-spread-from-his-kitchen-table
#########################################################
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.