week_06_Discussion_example_option2.R

week_06_Discussion_example_option2.R

#######################################################
# DISCUSSION ASSIGNMENT: FIT R0s FOR A COUNTRY OF INTEREST
# 
# For the final part of your assignment, you will 
# pretend you are an epidemiologist writing a report
# on the spread of Covid-19 in a particular country, 
# and commenting on what you can tell about the effectiveness
# of interventions at this point.
#
#
# You will use ML to fit three (3) or more SIR models 
# to some country's  Covid-19 active cases. You will 
# then compare these model fits with AIC.
#
# You will have to choose what "regimes" to use -- 
# i.e., the points at which R0 might change due to 
# the citizenry doing physical distancing, lockdown, 
# and/or release of these.
#
# This means you will have to decide how many days
# after the first case to start each regime.  You 
# will do this by changing "time" as before. These 
# decisions are made by you, the researcher, based on 
# looking at the chart of case counts, and/or based
# on using google etc. to research when a particular 
# country went into lockdown.
# 
# You might also manually change "delay_val," which 
# represents "how long after a change in R0, do you
# start to see it show up in the active cases data?"
#
# I found that a "delay_val" of 14 worked well, but
# you are welcome to change it if you want.
#
# There are 2 ways to do this assignment:
#
# Option #1. IF YOU HAVE R ON YOUR COMPUTER, AND HAVE
#            RELIABLE INTERNET.  Use country case count 
#            data provided by the European Union.  
# 
# Option #2: IF YOU ARE USING AN ONLINE VERSION OF R, 
#            OR HAVE WEAK INTERNET, you might not 
#            be able to access the online EU data
#            from your R session. Instead, update 
#            the New Zealand data 
#            by hand so that it extends to the present 
#            day, and then fit your models to that.
#            You can get the updated data from the 
#            "Total and active cases" chart at:
#
# Hancock, Faith (2020). Covid-19 in NZ. NZ Herald.
# https://www.newsroom.co.nz/2020/03/27/1101216/covid-19-in-nz-the-numbers
#
#            We just want the "active cases" count (hover your
#            mouse over the blue line to see the numbers).
#
# I am providing example scripts for both options. These 
# scripts will include all of the required functions. 
#
# You will just have to:
#
# 1. Copy the whole example script to a text editor where 
#    you can make/save changes.
#
# 2. Copy-paste the whole script into R 
#    and run it once,
#    just to make sure it works for you. If you hit some
#    problem with Option #1, go to the Option #2 script and
#    try it.
#
# 3. Modify the script as follows: 
#
#    For option #1, change the country you want to study
#    from something other than New Zealand. Change the 
#    "time" options for each model as desired.
#
#    For option #2, edit the New Zealand active cases data
#    in the script to make it up-to-date. Change the "time"
#    options for each model as desired.
#
#    You are free to experiment with other changes if you
#    like, of course. (E.g., changing plot titles, line
#    colors, other features of the model, whatever.)
#
# 4. Get the log-likelihoods under each model.
#
# 5. Calculate AIC, relative likelihoods, and AIC weights
#    for each model. You can do it in R, or modify the 
#    Excel or Google-Sheet examples I have provided.
#
# 6. In a Discussion post:
# 
#    * include a screenshot of your AIC table, including
#      at least these columns:
#      - Model name
#      - log-likelihood
#      - number of free parameters
#      - AIC
#      - AIC model weights
#
#    * include a screenshot that plots your AIC-best 
#      model against the data for your country
#
#    * Summarize your results and critically evaluate:
#     
#      E.g., "I fit 4 SIR models 
#      to the active-case data from country X. The models
#      had 1, 2, 3, or 4 regimes with different R0 values.
#      The AIC-best model was Model 3, gaining 55% of the
#      AIC weight. Under this best model, R0 was 6.1 for
#      the first phase of the epidemic, indicating extremely
#      rapid spread. However, after the country imposed a 
#      lockdown on day 14, R0 dropped to 0.1.  However, 
#      political disputes in the country led to the lockdown
#      being lifted after 7 days, causing a second wave of 
#      cases with R0 estimated as 4.5. This wave appears to
#      be ongoing.
#
#      Based on these model results, country X clearly needs
#      to return to lockdown. However, visual inspection of
#      the plot of the best-fit SIR model against the data
#      suggests that the available models were not a 
#      perfect fit, perhaps because of a change in case
#      definitions on Day 23 that made it appear as if 
#      1000 cases occurred suddenly on that day. A more 
#      sophisticated model would account for this change
#      in case definition, as well as XYZ.
#
#      References (if any)
#
# 7. Comment on at least one other student's summary, 
#    pointing out what you found interesting, making
#    suggestions, or asking questions. As always, stay
#    professional.
#      
#######################################################



#######################################################
# Example script for Option #2, Week 6 Discussion assignment
#######################################################

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 your data table:
head(dat)  # look at first 6 rows
tail(dat)  # look at last 6 rows
dim(dat)   # look at the number of rows & columns

# Plot the data:
country_name = "New Zealand"
dat$Date = as.Date(dat$Date, format="%Y-%m-%d")
plot(x=dat$Date, y=dat$Active_cases, xlab="Date", ylab="Active cases", col="firebrick2")
titletxt = paste0("Active COVID-19 cases in ", country_name, "\n(NZ Ministry of Health)")
title(titletxt)



#######################################################
# LOAD OUR FUNCTIONS FOR SIR SIMULATION AND INFERENCE
#######################################################

library(deSolve)

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

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

#######################################################
# 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()


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

#######################################################
# FINISHED LOADING OUR FUNCTIONS FOR SIR SIMULATION AND INFERENCE
#######################################################





#######################################################
# Model 1: A 1-regime model
#######################################################
times = 1:nrow(dat)      # Number of days since first case
time = c(1, nrow(dat)-1) # The 2nd entry is the 2nd-to-last day,
                         # this forces most of the data to
                         # use the first R0
R0 = c(3.0, 2.5)
D_inf = c(14.0, 14.0)
S = c(5000000, 0) # putting in the country's population size
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=times)
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, EUDC active cases from: ", country_name, "\nModel 1 (a 1-regime model) max lnL = ", 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
#######################################################
times = 1:nrow(dat)      # Number of days since first case
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) # putting in the country's population size
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=times)
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, EUDC active cases from: ", country_name, "\nModel 2 (a 2-regime model) max lnL = ", round(ML_results$value, 2))

# 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
#######################################################
times = 1:nrow(dat)      # Number of days since first case
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) # putting in the country's population size
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=times)
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, EUDC active cases from: ", country_name, "\nModel 3 (a 3-regime model) max lnL = ", 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
#######################################################
times = 1:nrow(dat)      # Number of days since first case
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) # putting in the country's population size
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=times)
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, EUDC active cases from: ", country_name, "\nModel 4 (a 4-regime model) max lnL = ", 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



# Calculate AIC and AIC model weights by
# 
# (Option A) Re-purposing the R code from 
# the Week 6 exercise
# 
# (Option B) Editing the example Excel or Google Sheet files:
#
#
# Example spreadsheet calculating AIC and AIC model weights:
#
# Excel file on Canvas: 
# BIOSCI 220 -> Files -> Module 2 -> Week 6 -> Week6_AIC_with_SIRs_fit_to_NZ_Ministry_of_Health_data.xlsx
# 
# On public web: 
# http://phylo.wikidot.com/local--files/biosci220:quantitative-biology/Week6_AIC_with_SIRs_fit_to_NZ_Ministry_of_Health_data.xlsx
# 
# On Google Sheets: 
# (You will have to copy-paste to your own Google Sheet, this one is edit-locked.)
# https://docs.google.com/spreadsheets/d/1WNk1la1TfoIdKcadxm0kUW5JNTTl0ZDImzDEQq0fdSQ/edit?usp=sharing
#