Overview

In Part 1, we generated data from a variety of ODE models. Here in Part 2, we outline the entire fitting process. The code here is, in part, a “brain dump” - so we will look to flesh the details out a bit more later on.

Outline

In this session we cover …

A. Generating Data from a Differential Equation Model
B. Model Fitting: Integrating The Model Approach
C. Model Fitting: Differentiating the Data Approach
D. Coupled Damped Linear Oscillator Example

Prelim - Loading libraries used in this script.

library(deSolve)  #for the differential equation solvers
library(ggplot2)   #for plotting
library(nlme)     #for fitting mixed effets models
library(zoo)       #for making time-series embedded data
library(reshape2)  #for reshaping data
library(tseriesChaos) #for embedd function 

Setting a seed becasue we are generating data.

set.seed(1234)

A. Generating Data from a Differential Equation Model

Let’s start with a very simple model - the constant (linear) growth model. In differential form, the model is …
\[\frac{dx}{dt} = a \]

We can generate some data from this model

1. Write the differential equation as a function.

#diffeq model as function
model.constant = function(times, y, parms) {
  ## times = current times
  ## y = current state vector: in this case the 0th derivative, c(x)
  ## parms = vector of model parameters.
  with(as.list(c(y,parms)),{
  derivs = c( 
    dx = a
    )
  ## Return
  list(derivs)
})
}

2. Specify vectors with the initial values, model parameters, and times over which to integrate.

#Initial Values
initvalues.1 <- c(x = 0.1)
#Model Parameters
parms.1 <- c(a = 0.2)
#Times over which to integrate
times <- seq(0, 100, 1)

3. Solve the differential equation using numerical methods.

This is done using the ode() function, which is itself a wrapper of the lsoda() function (and a collection of similar functions for different types of integrations (e.g., Runge-Katta), some of which work better with particular types of models (e.g., stiff and non-stiff systems)).

out <- ode(y = initvalues.1, times = times, func = model.constant, parms = parms.1) 

This produces an object of class deSolve that can then be examined, summarized, used for diagnosis, and plotted.

head(out)
##      time   x
## [1,]    0 0.1
## [2,]    1 0.3
## [3,]    2 0.5
## [4,]    3 0.7
## [5,]    4 0.9
## [6,]    5 1.1
summary(out)
##                  x
## Min.      0.100000
## 1st Qu.   5.100000
## Median   10.100000
## Mean     10.100000
## 3rd Qu.  15.100000
## Max.     20.100000
## N       101.000000
## sd        5.860034
diagnostics(out)
## 
## --------------------
## lsoda return code
## --------------------
## 
##   return code (idid) =  2 
##   Integration was successful.
## 
## --------------------
## INTEGER values
## --------------------
## 
##   1 The return code : 2 
##   2 The number of steps taken for the problem so far: 102 
##   3 The number of function evaluations for the problem so far: 103 
##   5 The method order last used (successfully): 2 
##   6 The order of the method to be attempted on the next step: 2 
##   7 If return flag =-4,-5: the largest component in error vector 0 
##   8 The length of the real work array actually required: 36 
##   9 The length of the integer work array actually required: 21 
##  14 The number of Jacobian evaluations and LU decompositions so far: 0 
##  15 The method indicator for the last succesful step,
##            1=adams (nonstiff), 2= bdf (stiff): 1 
##  16 The current method indicator to be attempted on the next step,
##            1=adams (nonstiff), 2= bdf (stiff): 1 
##  
## --------------------
## RSTATE values
## --------------------
## 
##   1 The step size in t last used (successfully): 1 
##   2 The step size to be attempted on the next step: 1 
##   3 The current value of the independent variable which the solver has reached: 100.002 
##   4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0 
##   5 The value of t at the time of the last method switch, if any: 0 
## 

These diagnostics output may be useful for checking the kind of integration method used (ode() chooses on its own), and if there were problems. (I have not run enough of these yet to know how, specifically, to use the diagnostics)

4. Plot the output trajectory (solution).

plot(out, main="Constant Growth", lwd=2, xlab="Time", ylab="x")

Given our proclivity for plotting using ggplot2 function, we can also make the plot there, after conversion of the deSolve object into a data frame.

outplot=as.data.frame(out)

ggplot(data=outplot, aes(x=time, y=x)) +
  geom_point() + 
  geom_line() +
  xlab("Time") +
  ylab("x") +
  ggtitle("Constant Growth")

Expansion to Multiple Scenarios (Persons)

We can expand by plotting multiple scenarios (e.g., simualted trajectories for a set of selected or randomly chosen parameters (or starting values)).

First, we generate parameter values for a hypothetical sample of 30 persons.

#Creating a list of ids and parameter values 
outlist <- list()
plist <- cbind(id = 1:30,
               a = runif(30,min=0.1,max=.4))

Then, we use the ode() function to generate longitudinal data for those 30 persons

#Loop ode() through the matrix of parameter values and create a data set
#create individaul ids and parameters
plist <- cbind(id = 1:30,
               a = rnorm(30,mean=.25,sd=0.1))
#create a list to hold ode() generated values
outlist <- list()
#Loop through the matrix of parameter values
for(i in 1:nrow(plist)){
  outlist[[i]] <- ode(y = initvalues.1, times = times, func = model.constant, parms = plist[i,"a"])
}

#combne into data frame with original data
#creating data.frame holder 
outdata <- data.frame()
#converting plist and outlist into a combined data frame
for(i in 1:nrow(plist)){
  idata <- (cbind(id = plist[i,"id"], a = plist[i,"a"], outlist[[i]]))
  outdata <- rbind(outdata,idata)
}

And then plot the resulting longitudinal data.

#plotting
ggplot(data = outdata, aes(x = time, y = x, group = id, color=factor(id))) +
  #geom_point() + 
  geom_line() +
  xlab("Time") + 
  ylab("x") +
  scale_x_continuous(breaks=seq(0,100,by=10))