Sensitivity analysis of continuous-time Markov models

Matlab code for the analyses described in M. Spencer (unpublished), Sensitivity analysis of Markov models for communities of competing sessile organisms.  Tested with Matlab release 14 service pack 1 for Linux.


To get the derivative dwn of the stationary distribution with respect to element i,j of the continuous-time rate matrix Q:

dw1 =Qstat_deriv(Q,i,j,1);

To get the derivative of the matrix P(t) = e^(Qt) with respect to element i,j of Q (e is matrix exponential):

dPqij=Pmatrixderivs(Q,i,j,t);

To get the directional derivative of the stationary distribution of P(t) in the direction given by dPqij:

dw1_P=Pstat_deriv(expm(Q*t),1,dPqij);

dw1_P and dw1 should be identical up to rounding error.  This answers questions about the effects of changing a single direct interaction rate.

To get the directional derivatives of the stationary distribution of P(t) corresponding to scenarios that maintain the column-stochastic nature of P(t), as discussed in Caswell, H., 2001, Matrix Population Models, second edition, Sinauer Associates, Sunderland, MA, page 253:

[Pdwnspec,Pdwnscalespec,Pdwnuni,Pdwnscaleuni,Pdwnprop,Pdwnscaleprop]=Pstat_Vderivs(expm(Q*t),i,j,1);

In order, 
	Pdwnspec and Pdwnscalespec are specific compensation in element j,j
	Pdwnuni, Pdwnscaleuni are uniform compensation in all elements of column j
	Pdwnprop, Pdwnscaleprop are proportional compensation in all elements of column j
	
(in each case, the ...scale... versions are explicitly scaled so the stationary distribution always sums to 1, but this happens automatically through appropriate scaling of the direction matrix).  All of these answer questions about the effects of changing a large number of direct interaction rates simultaneously.

To get the changes in direct interactions corresponding to a specified change V in P(t):

[Vc,dPqij]=invert_Pmatrixderivs(Q,V,t);

Vc is the corresponding change in direct interactions.
dPqij is just a check that it worked properly: it redoes the transformation on Vc that invert_Pmatrixderivs undoes, so it should be the same as V.

Matthew Spencer (matts@mathstat.dal.ca), 20 December 2005