sir.model
Class SEIRState

java.lang.Object
  extended bysir.model.EpiState
      extended bysir.model.SEIRState
All Implemented Interfaces:
org.omegahat.Simulation.MCMC.MCMCState, Uncmin_methods
Direct Known Subclasses:
SEIRwDState

public class SEIRState
extends EpiState
implements org.omegahat.Simulation.MCMC.MCMCState, Uncmin_methods

Subclass implementing all functionality for an multivariate SEIR model. This class is used to represent the state as a SEIR model as mentioned in the O'Neill and Becker article "Inference for an epidemic when susceptibility varies", P.D. O'Neill and Niels G. Becker, Biostatistics 2(1), pp.99-108 (2001) and the M. Hoehle, E. Jørgensen, and P.D. O'Neill article "Inference in disease transmission experiments using stochastic epidemic models. I.e. we have

Note: Earlier versions of this program contained the possibility to specify fixed times where the set of parameters switched due to control measures etc. All this functionality is now removed -- consult the reference "Estimating Parameters for Stochastic Epidemics" by M. Höhle and E. Jørgensen, Dina Research Report no. 102, 2002.

Copyright: GPL

Version:
16 October 2002
Author:
Michael Höhle <hoehle@stat.uni-muenchen.de>

Field Summary
 double deltaE
          Parameters of the gamma distributed incbuation time
 int deltaE_
          Index of inc time parameters in param vector
 double deltaI
          Parameters of the gamma distributed infectious time
 int deltaI_
          Index of inc time parameters in param vector
 Event[] E
          All the exposed states
static double Eaccept
           
 boolean expTransformParams
          Are we operating with log transformed parameters in the loglik function.
 Event[][] firstEventInUnits
          First event in all units - a data structure updated each timegenXYFuncs is called.
 double gammaE
          Parameters of the gamma distributed incbuation time
 int gammaE_
          Index of inc time parameters in param vector
 double gammaI
          Parameters of the gamma distributed infectious time
 int gammaI_
          Index of inc time parameters in param vector
 boolean justDiagnostic
          Just a foobar thing - only use diagnostic times in the likelihood
 Event[][] lastEventInUnits
          Global variable to speed up the lastEventInit method so we do not need to new each time
 boolean MeanVarianceParameters
          Re-parameterization using mean and variance
 boolean useOrigins
          Are we using exposure origins in the loglik calculations.
 
Fields inherited from class sir.model.EpiState
accept, beta, beta_, betaN, betaN_, data, eta, I, loglik, modelParamNames, nogen, noOfModelParams, prob, R, runs, theta
 
Constructor Summary
SEIRState()
          Constructor - new necessary structures
 
Method Summary
protected  void addExtraEventsHook(int baseNoOfEvents)
          Hook function to add additional events to the events array.
 void array2Params(double[] theta, boolean expTransformParams)
          Helper function, given an array of size noOfModelParams assign the values of this vector to the correct names.
 double betaHaz(Unit u, Unit o, Event[][] le)
          Method to compute \lambda_e^j (t-) when the origin of the event is known.
 double betaHazAll(Unit u, Event[][] le)
          Method to compute \lambda_{all}^j(t-|...) for a specific individual (located in u_j).
 double betaSurvLn(Unit u, double tstop)
          Compute log S_j(t), i.e. the survival of an individual located in unit u up to time t, when overall hazard exerted on the individual is \lambda_j(t).
 java.util.LinkedList Eloc()
          Count the number of exposures originating from the unit itself.
 java.util.LinkedList Enbr()
          Count the number of exposures originating from a neighbour.
protected  int extraNoOfEventsHook()
          Hook function returning how many additional events are put into the events array by subclasses etc.
 double f_to_minimize(double[] theta)
          In case the object is given to a numerical minimizer, this is the function to minimize.
 void f_to_minimizeHook(double[] theta)
          Hook function to move additional parameters in numerical optimizer.
 EpiState generateNext(boolean justTheta)
          Generate method for MCMC chain Modifcation happen on the instance itself to keep things simple.
 void genXYFuncs()
          Generate X, Y, XY, and time by sorting appropriate events.
 void gradient(double[] x, double[] g)
          Gradient function, in case one wants to do numeric optimization of the likelihood.
 void hessian(double[] x, double[][] h)
          Hessian function, in case one wants to do numeric optimization of the likelihood.
 boolean isValidHook()
          Hook function to check for a valid state.
 Event[][] lastEventInit(boolean newArray)
          Prepare the last event data structure.
 double loglik()
          Just a dummy director to the currently used loglik function depending on the currently used parametrisation.
protected  double loglikAddIndividualHook(int individual)
          Hook function for adding additional events to the loglikelihood.
 double loglikGammaDelta()
          Abstract method to calculate the loglikelihood of the state.
 double[] ML(double[][] iHessian)
          Perform maximum likelihood estimation by numerically maximizing the likelihood equation.
 void MLParamsPrintHook(double[] theta)
          Hook function called inside ML once the results are showns
 void MLParamsSetup(double[] theta0, double[] typicalParamSize)
          Function to be called when initializing the theta0 and typicalParamSize initialization in the ML method.
 void MLProfiles()
          Foo function to be specialized in SEIRWState
 double[] params2Array()
          Helper functions to put all named parameters in an array.
 double[] periodsEV()
          Small stats function to compute mean and variance for inf or inc periods of current configuraion
 void sampleOrigins()
          Get origins right once we know E(t), I(t), R(t), i.e can't pick an origin where I(t) = 0.
 EpiState sampleState(boolean sampleCensored)
          Make a valid instance of the parameters to start the MCMC chain with.
protected  void sampleStateEventsHook(int i, boolean sampleCensored)
          Hook function to new any additional events in subclasses.
protected  void sampleStateParamsHook()
          Hook function to initialize any additional parameters introduced in subclasses.
 java.lang.String showStats(int noOfIterations)
           
 java.lang.String toHeaderHook()
          Hook function for header writing
 java.lang.String toStringDiagnostic()
          Show how state params relate to data
 java.lang.String toStringHook()
          Hook function for SEIR specific output
 java.lang.String toStringLong()
          Current State information is shown on screen in long format
 void updateBeta()
          Gibbs updating of beta - we know the full conditional (see paper) This version handles the multiple unit setup.
 void updateBetaN()
          Gibbs updating of betaN - we know the full conditional (see paper) This is a new version correcting an mistake of the conditioning on o_j.
 void updateDeltaE()
          Gibbs updating of deltaI - we know the full conditional (see paper)
 void updateDeltaI()
          Gibbs updating of deltaI - we know the full conditional (see paper)
 void updateE()
          Update E using full conditionals.
 void updateGammaE()
          Update gammaE using Metropolis Hastings.
 void updateGammaI()
          Update gammaI using Metropolis Hastings.
 void updateHook(boolean justTheta)
          Hook method to update additional parameters using MCMC.
 void updateI()
          Update I using Metropolis Hastings.
 void updateMetroHastings(double[] theta, int paramIndex, double sigma)
          Random Walk Metropolis-Hastings updating of a parameter with gamma prior distribution.
 void updateO()
          Update the o vector using Gibbs Sampling.
 void updateR()
          Update R using Metropolis Hastings.
 
Methods inherited from class sir.model.EpiState
acceptProposal, analLog, fromString, generateNext, isIValid, setData, setProb, showStats, toHeader, toString, toStringEvents, Yintegrate
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, wait, wait, wait
 

Field Detail

gammaE

public double gammaE
Parameters of the gamma distributed incbuation time


deltaE

public double deltaE
Parameters of the gamma distributed incbuation time


gammaI

public double gammaI
Parameters of the gamma distributed infectious time


deltaI

public double deltaI
Parameters of the gamma distributed infectious time


gammaE_

public int gammaE_
Index of inc time parameters in param vector


deltaE_

public int deltaE_
Index of inc time parameters in param vector


gammaI_

public int gammaI_
Index of inc time parameters in param vector


deltaI_

public int deltaI_
Index of inc time parameters in param vector


MeanVarianceParameters

public boolean MeanVarianceParameters
Re-parameterization using mean and variance


justDiagnostic

public boolean justDiagnostic
Just a foobar thing - only use diagnostic times in the likelihood


expTransformParams

public boolean expTransformParams
Are we operating with log transformed parameters in the loglik function. This is necassry when calling the optimizer, but not good when computing the Fisher info matrix


E

public Event[] E
All the exposed states


Eaccept

public static double Eaccept

useOrigins

public boolean useOrigins
Are we using exposure origins in the loglik calculations. Always true except in ML calculations


firstEventInUnits

public Event[][] firstEventInUnits
First event in all units - a data structure updated each timegenXYFuncs is called. The purpose is to speed up the lastEventInit() method.


lastEventInUnits

public Event[][] lastEventInUnits
Global variable to speed up the lastEventInit method so we do not need to new each time

Constructor Detail

SEIRState

public SEIRState()
Constructor - new necessary structures

Method Detail

toStringDiagnostic

public java.lang.String toStringDiagnostic()
Show how state params relate to data


toStringLong

public java.lang.String toStringLong()
Current State information is shown on screen in long format

Overrides:
toStringLong in class EpiState

betaHaz

public double betaHaz(Unit u,
                      Unit o,
                      Event[][] le)
Method to compute \lambda_e^j (t-) when the origin of the event is known.

Parameters:
u - Unit, where the individual is located.
o - Unit, where the exposure is originating from (o is equal to u or one of its neighbours)
le - Array with last event occured in each section. I.e. state just before to the event of interest.
Returns:
\lambda_e^j(t-) = beta^I(u=o) beta_n^I(o\in N_4(u)) Y_o(t) zero otherwise

betaHazAll

public double betaHazAll(Unit u,
                         Event[][] le)
Method to compute \lambda_{all}^j(t-|...) for a specific individual (located in u_j). That is the overall hazard exerted on a single individual just prior to occurence of an event. Corresponds to the setup, where no locations are given.

Parameters:
u - Unit, where the individual is located.
le - Array with last event occured in each section. I.e. state previous to the event of interest.
Returns:
\lambda_{all}^j(t-|...) = beta*Y_u(t) + \sum_{u'\in N_4(u)} beta_n * Y_u'(t-)

lastEventInit

public Event[][] lastEventInit(boolean newArray)
Prepare the last event data structure. This is an optimized implementation exploiting that the genXYFuncs created the desired array and that we reuse an array instead of newing something every time. Newing is expensive, but can be necessary at times.

Parameters:
newArray - Should we new a new instance. Necessary if we have nested calls to lastEventInit as in loglik
Returns:
Array where each unit is set to its first event.

betaSurvLn

public double betaSurvLn(Unit u,
                         double tstop)
Compute log S_j(t), i.e. the survival of an individual located in unit u up to time t, when overall hazard exerted on the individual is \lambda_j(t).

Parameters:
u - The unit, where the individual is located.
Returns:
-\int_{0)^{t} \lambda_j(u) du (notice this is log S_j(t))

loglikGammaDelta

public double loglikGammaDelta()
Abstract method to calculate the loglikelihood of the state. Parameterisation as I~\Gamma(gammaI,deltaI). The equation is taken from the article. Now extended to a spatial setting and optimized slightly for speed. The global variable @see useOrigins is used to check whether we use them in the calculations. Note: It is required that genXYFuncs has been called. Internal non-documented options: useOrigins : Ignore the origin trick.


loglik

public double loglik()
Just a dummy director to the currently used loglik function depending on the currently used parametrisation.

Specified by:
loglik in class EpiState

sampleOrigins

public void sampleOrigins()
Get origins right once we know E(t), I(t), R(t), i.e can't pick an origin where I(t) = 0.


sampleStateParamsHook

protected void sampleStateParamsHook()
Hook function to initialize any additional parameters introduced in subclasses.


sampleStateEventsHook

protected void sampleStateEventsHook(int i,
                                     boolean sampleCensored)
Hook function to new any additional events in subclasses.

Parameters:
i - Number of the individual we are newing events for.

sampleState

public EpiState sampleState(boolean sampleCensored)
Make a valid instance of the parameters to start the MCMC chain with. If a proper prior is specified for the model parameters we begin in the mean of the prior. For I/R a state is samples for all parametric instances. Semantic is that if there is at least one missing R it's a transmission exp. Note: Currently no regime or individual spatial param setup implemented.

Specified by:
sampleState in class EpiState

updateMetroHastings

public void updateMetroHastings(double[] theta,
                                int paramIndex,
                                double sigma)
Random Walk Metropolis-Hastings updating of a parameter with gamma prior distribution. This is general method not necessarily efficient as the entire likelihood is calculated.

Parameters:
theta - Vector of all parameters.
paramIndex - of the parameter to update
sigma - How large a standard deviation to use in the proposal distribution Note: The value is updated automatically by calling the array2Param function Therefore no value is returned

updateBeta

public void updateBeta()
Gibbs updating of beta - we know the full conditional (see paper) This version handles the multiple unit setup.


updateBetaN

public void updateBetaN()
Gibbs updating of betaN - we know the full conditional (see paper) This is a new version correcting an mistake of the conditioning on o_j.


updateDeltaE

public void updateDeltaE()
Gibbs updating of deltaI - we know the full conditional (see paper)


updateDeltaI

public void updateDeltaI()
Gibbs updating of deltaI - we know the full conditional (see paper)


updateGammaE

public void updateGammaE()
Update gammaE using Metropolis Hastings. This is a new improved version only using the relevant parts of loglik and prior. I.e. \pi(\gamma_E|\delta_E,D) and equation about full condition in article. As noted we could also try to do this using adaptive rejection sampling (ARS). Note though mean/var representation does not work now, but it sucks anyway.


updateGammaI

public void updateGammaI()
Update gammaI using Metropolis Hastings. This new version includes only relevant parts of the likelihood when calculating the acceptance prob. Quite a speed gain! Note though mean/var representation does not work now, but it sucks anyway.


updateO

public void updateO()
Update the o vector using Gibbs Sampling. This implementation is not very speedy and uses established \lambda functions instead of doing the necessary integration itself.


updateE

public void updateE()
Update E using full conditionals. Currently, Metropolis Hastings is used, but at this point we would very much prefer to use some adaptive rejection scheme. Still, the formula is quicker to update so we update all at the same time. This version does not work for mean/var parametrization, just the natural. Note: This does not take Diagnostic events into account. Specialized version necessary!


updateI

public void updateI()
Update I using Metropolis Hastings.


updateR

public void updateR()
Update R using Metropolis Hastings.


updateHook

public void updateHook(boolean justTheta)
Hook method to update additional parameters using MCMC. To be specialized by subclass.

Parameters:
justTheta - no data just parameters are updated.

generateNext

public EpiState generateNext(boolean justTheta)
Generate method for MCMC chain Modifcation happen on the instance itself to keep things simple.

Specified by:
generateNext in class EpiState
Parameters:
justTheta - Boolean indicating whether to just update the beta,...,gamma_D variables
Returns:
this

f_to_minimizeHook

public void f_to_minimizeHook(double[] theta)
Hook function to move additional parameters in numerical optimizer.

Parameters:
theta - The theta array

periodsEV

public double[] periodsEV()
Small stats function to compute mean and variance for inf or inc periods of current configuraion

Returns:
4Dim array containing E_e, E_i, V_e, V_i

params2Array

public double[] params2Array()
Helper functions to put all named parameters in an array. Remember Fortran indexing Function currently does not work if incu period a constant.

Specified by:
params2Array in class EpiState
Returns:
An @see noOfModelParams + 1 long array with params

array2Params

public void array2Params(double[] theta,
                         boolean expTransformParams)
Helper function, given an array of size noOfModelParams assign the values of this vector to the correct names.

Parameters:
theta - - A vector of length noOfModelParams containing the values to assign
expTransformParams - - When using the minimizer it is a good idea to operate with log(param) value instead to avoid negative values.

f_to_minimize

public double f_to_minimize(double[] theta)
In case the object is given to a numerical minimizer, this is the function to minimize. We will only do maximization of the likelihood, i.e. parameters have to be transmoglyffed to account for different indexing of arrays etc. Thus, this will work as a wrapper to logflambdaI.

Specified by:
f_to_minimize in interface Uncmin_methods
Parameters:
theta - The parameter vector(remember Fortran indexing of arrayx (starting at 1))

gradient

public void gradient(double[] x,
                     double[] g)
Gradient function, in case one wants to do numeric optimization of the likelihood. Contains an EMPTY body.

Specified by:
gradient in interface Uncmin_methods
Parameters:
x - The parameter (remember Fortran indexing of arrayx (starting at 1))
g - The computed gradient.

hessian

public void hessian(double[] x,
                    double[][] h)
Hessian function, in case one wants to do numeric optimization of the likelihood. Contains an EMPTY body.

Specified by:
hessian in interface Uncmin_methods
Parameters:
x - The parameter (remember indexing of arrayx (starting at 1))

extraNoOfEventsHook

protected int extraNoOfEventsHook()
Hook function returning how many additional events are put into the events array by subclasses etc. To be rewritten by subclasses.

Returns:
Number of additional events added.

addExtraEventsHook

protected void addExtraEventsHook(int baseNoOfEvents)
Hook function to add additional events to the events array. To be rewritten by subclasses.

Parameters:
baseNoOfEvents - Number of E,I,R events in SEIRState class.
Returns:
Number of additional events added.

loglikAddIndividualHook

protected double loglikAddIndividualHook(int individual)
Hook function for adding additional events to the loglikelihood. It is called from within loglikGammaDelta. The recovery events are already added. To be specialized in subclasses.

Parameters:
individual - The i'th individual to add events for
Returns:
Likelihood of the additional events.

genXYFuncs

public void genXYFuncs()
Generate X, Y, XY, and time by sorting appropriate events. Note: New Version using Event classes. Before, we sorted by zipping. Now we use general sort on all events. Slower? Probably - use zipper again, which required I to be sorted, but then we have to sort it twice??

Specified by:
genXYFuncs in class EpiState

MLParamsSetup

public void MLParamsSetup(double[] theta0,
                          double[] typicalParamSize)
Function to be called when initializing the theta0 and typicalParamSize initialization in the ML method. To be specialized int subclasses

Parameters:
theta0 - Array of theta0's, i.e. the start values
typicalParamSize - Typical size

MLParamsPrintHook

public void MLParamsPrintHook(double[] theta)
Hook function called inside ML once the results are showns

Parameters:
theta - Array with the optimized valued. Index 1-4 are taken care of.

ML

public double[] ML(double[][] iHessian)
Perform maximum likelihood estimation by numerically maximizing the likelihood equation. In special cases some params can be derived analytically.

Parameters:
iHessian - Inverse of Hessian as array (has to be newed before!)

MLProfiles

public void MLProfiles()
Foo function to be specialized in SEIRWState


Eloc

public java.util.LinkedList Eloc()
Count the number of exposures originating from the unit itself.

Returns:
E_{loc}

Enbr

public java.util.LinkedList Enbr()
Count the number of exposures originating from a neighbour.

Returns:
E_{nbr}

isValidHook

public boolean isValidHook()
Hook function to check for a valid state. In our case it's a special case if its a disease transmission experiment. Then all event times have to be in >= 0

Overrides:
isValidHook in class EpiState
Returns:
Boolean indicating whether all event times are >=0.

toHeaderHook

public java.lang.String toHeaderHook()
Hook function for header writing

Overrides:
toHeaderHook in class EpiState

toStringHook

public java.lang.String toStringHook()
Hook function for SEIR specific output

Overrides:
toStringHook in class EpiState

showStats

public java.lang.String showStats(int noOfIterations)