首页 >
> 详细

Lab 6: Parameter Estimation

Epid 633

Scenario & Setup

A flu-like illness has been spreading through Stardew Valley County. You have been asked to develop a predictive model for the disease, which the Stardew Valley health authorities will use to estimate several important parameters that can help in designing an intervention:

the infectious period of the disease

the reporting fraction (i.e. what fraction of cases are reported)

the transmission rate and/or basic reproduction number (R0)

In addition, they are hoping that your model will be able to successfully forecast the number of cases there will be in the upcoming weeks.

Data set

Here’s the data set that we’ll be working with–weekly data on the number of currently ill individuals:

times = c(0, 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98)

data = c(97, 271, 860, 1995, 4419, 6549, 6321, 4763, 2571, 1385, 615, 302, 159, 72, 34)

1) SIR Model Setup

Start by setting up an SIR model, using the equations:

In this case, we’ll take S, I, and R to be the fraction of the population that’s susceptible, infectious, and recovered. The measurement equation we’ll be using is:

y=kNI

where N is the total population size and k is the reporting fraction. To break down this equation, note that I is the fraction of the population currently infected, so NI is the total number of actual people currently infected. Then y=kNI represents the total number of observed cases currently infected.

Code up the model above, noting that:

Stardew Valley County has a total population of N=200000. We will treat this as a fixed (i.e. known, not estimated) parameter, so you can just code this as a number (i.e. it’s not part of the parameter vector that you’ll use for estimation).

At t=0, approximiately 200 people were infected, and no one had recovered yet. This means that the initial conditions were that 99.9% of the population was susceptible, and 0.1% was infected (i.e. x0 = c(0.999,0.001,0) for S, I, and R respectively).

For starting parameters, we can take β=0.4. So far, the infectious period appears to be roughly 4 days, so we can start with a guess of γ=1/4=0.25. For the reporting fraction, let’s start by supposing we have perfect reporting (i.e. every case is reported), so that k=1.

Now, plot your model output (y) from 0 to 100 days, and plot the data on the same plot (use a line for your model and points/dots for the data). How well does the model match the data using our initial guesses for the parameters? Do you notice any issues with the model fit currently?

2) Setting up the likelihood

Next, we need to make a function that will evaluate the goodness-of-fit for the model—in this case, the negative log likelihood. We will pass this function to the optimizer, optim, which will use our function to find the best-fit parameter estimates (i.e. the parameter values that minimize our function).

Code up a likelihood function for the SIR model, using supposing that your data is taken from a Poisson distribution with rate parameter y (I’ve included the Poisson likelihood form in the template code below, but you can try deriving it yourself if you want! It’s just like the Gaussian likelihood we did on Monday, just use a Poisson distribution instead of a normal distribution). Here’s the overall structure that the likelihood function will take (note that this is just a template! You need to adapt this for your own code, depending on how you name things, etc.):

CalcNLL = function(params, ODEfunction, x0, times, data){

# The function will calculate the negative log likelihood---to do that, we need:

# - the parameters (this is what optim will fiddle with to find the best fit)

#

# - the ODE function for the SIR model you wrote in part 1). Optim will use

# to simulate the model so it can compare the model to the data

#

# - to simulate the model with SIRode, we also need the initial conditions (x0)

# and the time points (times)

#

# - lastly, we need to pass the data into NLL, so we can compare the model and

# data

# Since we don't want to include negative values for the parameters, start by taking the absolute value of whatever the current parameter values are (note this is a bit of a silly trick, ask me more if you're confused)

params = abs(params)

# Simulate model

xcurr = ode(x0, times, ODEfunction, params) # replace this with whatever ODE code you're using to run your model!

# Measurement equation

y = k*N*I # replace "k", "N", and "I" with whatever they are in your code, e.g. params[1], xcurr[,3], etc.

# Negative Log Likelihood (NLL)

NLL = sum(y) - sum(data*log(y)) # Poisson ML

# note this is a slightly shortened version--there's an additive constant term missing but it shouldn't alter the optimal estimates. Alternatively, can do:

# NLL = -sum(log(dpois(round(data),round(y)))) # the round is b/c Poisson is for (integer) count data

# return(NLL)

}

Evaluate your likelihood function for a couple of values—try your default parameters, and then try taking β=0.42 and β=0.38 (leaving the other parameters equal to the defaults). Based on this, if you were optim, which direction would you move β to improve the fit? Does this match with what you saw in your initial model fit in Part 1?

3) Parameter Estimation

Okay, time to run the optimizer! Estimate your model parameters from the data set provided. While you’re at it, also pull the final, best-fit negative log likelihood value. The function optim takes the following inputs:

starting parameters the first input for optim will be the starting values for the model parameters

cost function - in our case, this is our likelihood-calculator we made above

arguments for the cost function - all the extra arguments that the cost function takes (i.e. everything else the cost function needs, like data,times,your ODE function, etc.)

So the structure is: optim(starting parameters, cost function, arguments for cost function). Here’s a template for how to set up your optimizer:

### First, run optim ###

results = optim(initialparams, fn = CalcNLL, ODEfunction = SIRode, times = times, data = data, x0 = x0) # Ask me if you're confused about why we have 'times = times', 'data = data', etc.!

### Pull out the parameter estimates ###

paramests = results$par

### Pull out the final cost function value (-log likelihood) ###

negLL = results$value

Plot your model and data on the same plot (like you did in Part 1). How does your model fit look now? How did your parameters change from your initial guess? What are the values you estimated for each of your parameters (β, γ, and k)? What implications might your estimate of k have for intervention efforts?

4) Forecasting

Let’s consider the case where you are attempting to fit and forecast an ongoing epidemic (i.e. with incomplete data). Truncate your data to only include the first five data points, then re-fit the model parameters.

How do your parameter estimates change?

Simulate your model for the full time (0 to 100 days) and plot your model predictions along with the complete data set. How well did you predict the dynamics of the epidemic?

Now re-do the forecast, but suppose that we only have the first four data points of the epidemic—how do your estimates and predictions change?

5) Model Comparison and the Akaike Information Criterion (AIC)

Lastly, let’s compare the SIR model we’ve been using to a slightly more complex model, the SEIR model:

Here, E represents the fraction of the population who are infected but not yet transmissible (i.e. latently infected).

Estimate the parameters using this model and the complete data set, supposing that:

The initial conditions for S, E, I, and R are: x0 = c(0.9985,0.001,0.0015,0)

The initial parameter guesses are params = c('beta'=0.4,'alpha' = 1, 'gamma'=0.25, 'kappa'=1)

You can use the same measurement equation as before (y=kNI). Does the fit look better/different than what you saw with the SIR model in Part 3? How are the parameter estimates different with this model compared to the SIR model?

Calculate the AIC for both models and compare. Which model has a lower (better) AIC? How different are the two AICs? Recall that the AIC is given by:

2(−LL+nP)

where −LL is the negative log likelihood and nP is the number of parameters to be estimated (i.e. 3 for the SIR model and 4 for the SEIR model).

2π) Extra problems (just for fun—these are not required and don’t need to be turned in!)

Try out alternative cost functions, e.g. using dpois for the Poisson likelihood, or using least squares. Do you notice differences in fit? What about in how long it takes to generate the parameter estimates?

Try making a heat map of the likelihood as a function of β and γ (you could try out a fancier plotting package like ggplot2 or plotly for this). How does the shape of the likelihood change as you reduce the amount of data you have?

联系我们

- QQ：99515681
- 邮箱：99515681@qq.com
- 工作时间：8:00-23:00
- 微信：codinghelp

- Data留学生作业代做、代写r程序语言作业、代做r课程设计作业、代写boos 2019-12-13
- 代写cmt307留学生作业、代做machine Learning作业、代写p 2019-12-13
- Comp2050作业代做、Computer Science作业代写、代做ja 2019-12-13
- 代写compsci 671D作业、代做modeling留学生作业、代写r程序 2019-12-13
- 代做gui留学生作业、代做programming课程作业、Java程序语言作 2019-12-13
- 代做framework课程作业、Java程序设计作业代写、Java实验作业代 2019-12-13
- Ae2ace留学生作业代写、代做stock Trading作业、Java编程 2019-12-13
- Eee101留学生作业代做、Programming作业代写、C++编程语言作 2019-12-13
- Csci 340作业代做、代写java程序语言作业、代做java实验作业代写 2019-12-12
- Data课程作业代做、代写nbershade作业、代做r课程设计作业、代写r 2019-12-12
- 代写csci 1100作业、Program课程作业代做、Python语言作业 2019-12-12
- Data留学生作业代做、代写sql实验作业、Sql编程有作业调试、Pseud 2019-12-12
- 代做g6077留学生作业、System课程作业代写、代做web编程语言作业、 2019-12-12
- 代写comp529作业、代做analysis留学生作业、代写java语言作业 2019-12-12
- Ce235留学生作业代写、Program课程作业代写、C/C++程序语言作业 2019-12-12
- 代写system留学生作业、代做python语言作业、代写java，C/C+ 2019-12-12
- 代写ma705留学生作业、代写python程序语言作业、代写python实验 2019-12-11
- Stat 3312作业代做、R语言作业代写、代做r编程设计作业、代写sas 2019-12-11
- Comp201作业代做、代写software Engineering作业、J 2019-12-11
- Statistics 3022作业代做、代写data留学生作业、R编程设计作 2019-12-11