Welcome State space notation Functions > kal.fil.mul.pom
Last update Jul 10, 2004

DESCRIPTION:

Extended Kalman filter for multinomial data assuming pom link and linear Gaussian latent process. It is a dynamic extension of the proportional odds model of McCullagh (1980). The specific model is
Observation equation:
(yt | ht)
~
Mn(nt,pt),
Link function:




logit(g1t)
:
logit(gk-1,t)




=
FtTqt,
System equation:
qt
=
Gtqt-1 + wt,
    wt
~
Np(0,Wt),
Initial prior:
q0
~
Np(m0,C0),
where git is the cumulated probabilities. Note that FtT and qt have a special structure. This structure is determined by the interpretation of the proportional odds model as a discretization of an underlying logistic variable. Please refer to Section 5.4 of Klein (2003) orsee the example below.

USAGE:

      kal.fil.mul.pom(ssm)

REQUIRED ARGUMENTS:

ssm:
Object of class ssm (state space model).
The following attributes must be specified in ssm:
Yt:
n × k matrix containing the multinomial observations arranged row-wise.
nt:
1 × n vector containing the number of trials at each time-point.
Should be nt <- apply(Yt,1,sum).
Ft:
Function returning a p × (k-1) design matrix.
Gt:
Function returning a p × p evolution transfer matrix.
Wt:
Function returning a p × p evolution variance matrix.
m0:
p × 1 vector (prior mean).
C0:
p × p matrix (prior variance).

OPTIONAL ARGUMENTS:

The following attributes of the ssm object might be needed:
Xt:
Matrix containing the covariates needed in Ft, Gt and/or Wt. This matrix must have n rows, such that row t contains the covariates needed for generating the system matrices at time t.
psi:
Parameter vector must be supplied if it is needed in the calculation of the system matrices.
The following attributes of the ssm object are optional:
fam:
Tacitly assumed to have value multinomial.
link:
Tacitly assumed to have value pom.
m.start:
Used to specify the Taylor expansion points. These are not meant to be set by the user. Applying the Kalman smoother on a ssm object that has been filtered will set the m.start attribute to the smoothed states. This is needed in the iterated extended Kalman filter and smoother.

VALUE:

Returns an object of class ssm with the same attributes as the object in the call, but with attribute filtered updated with:
mt:
Posterior means mt given the observations up to and including time t. These are arranged row-wise in a matrix.
Ct:
Posterior variances Ct given the observations up to and including time t. These are arranged in a 3-dimensional array such that Ct is placed in filtered$Ct[,,t].
llh:
Log-likelihood of of the approximating Gaussian state space model.
Rt:
Prior variances for the states arranged in 3-dimensional array such that Rt is placed in filtered$Rt[,,t]. The reason for returning these is that they are needed in the Kalman smoother.

SIDE EFFECTS:

None

DETAILS:

The filter works by sequentially approximating the multinomial observation density with a Gaussian density based on a Taylor expansion around the values supplied in the attribute m.start of the ssm object. If m.start is not supplied, the Taylor expansion will be around the one-step forecast mean of the observation.

REFERENCES:

Durbin & Koopman, (2002), Time Series Analysis by State Space Models, Oxford Statistical Science Series.
Fahrmeir & Tutz, (1994), Multivariate Statistical Modelling Based on Generalized Linear Models, Springer Series in Statistics.
Klein, (2003), State Space Models for Exponential Family Data, Ph.D. Thesis, Department of Statistics, University of Southern Denmark.
McCullagh, (1980), Regression models for ordinal data (with discussion) Journal of the Royal Statistical Society, Series B, Vol. 42, 109-142.

EXAMPLES:

# Application of the iterated extended Kalman filter and smoother
# to fit a proportional odds model
# (example from Chapter 5)

# Observations
Yt <- matrix(c(0, 0, 1, 7, 8, 8,19, 8, 1,
               6, 9,12,11, 7, 6, 1, 0, 0,
               1, 1, 6, 8,23, 7, 5, 1, 0,
               0, 0, 0, 1, 3, 7,14,16,11), 
             ncol=9, byrow=T)

# Number of trials
nt <- rep(52,4)

# Design matrix
Ft <- function(i,xt,phi)
{       
  result <- rep(0,11)
  if (((i-1)%%4+1)==1)  result <- t(cbind(diag(8),0,0,0))
  if (((i-1)%%4+1)==2)  result <- t(cbind(diag(8),-1,0,0))
  if (((i-1)%%4+1)==3)  result <- t(cbind(diag(8),0,-1,0))
  if (((i-1)%%4+1)==4)  result <- t(cbind(diag(8),0,0,-1))      
  result
}

# Evolution transfer matrix is the identity
Gt <- function(i, xt,phi) diag(11)

# Zero evolution noise
Wt <- function(i,xt,phi) diag(rep(0,11))        

# Specify the state space model object
cheese.ss <- ssm(Ft=Ft,
Gt=Gt,
Wt=Wt,
m0=c(-5.5,-4.4,-3.3,-2,-1,0,1.5,3,-3.4,-1.7,1.6)/2,
C0=diag(rep(1000,11)),
Yt=Yt,
nt=nt,
fam="multinomial", 
link="pom")

# Apply the iterated extended Kalman filter and smoother
cheese.ss <- ieks(cheese.ss, eps=10^-9)

# Estimates and standard errors
cbind(cheese.ss$smoothed$m.tilde[4,],
      sqrt(diag(cheese.ss$smoothed$C.tilde[,,4])))

########################################################


# Specify state space model 
ss <- ssm(Ft = function(i,x,phi)  
                       {t(cbind(diag(5),-1))},
          Gt = function(i,x,phi)  
                       {diag(6)},
          Wt = function(i,x,phi)  
                       {diag(c(rep(phi[1],5),phi[2]))},
          phi = c(0.001,0.1),
          m0 = c(-2,-1,0,1,2,0),
          C0 = diag(1,6),
          nt = floor(runif(50,min=20,max=50)),
          fam = "multinomial",
          link = "pom")

# Simulate observations
ss <- simulate.ssm(ss, n=50)

# Plot simulated latent process
plot(ss$simulated$theta[,6], type="l", ylim=range(ss$simulated$theta))
for (i in 1:5)
        lines(ss$simulated$theta[,i], lty=2)
title("Simulated latent process")

# Plot cumulated observations
tmp <- t(apply(ss$Yt,1,cumsum))/ss$nt
plot(tmp[,6], ylim=c(0,1), type="l")
for(i in 1:5)
        lines(tmp[,i])
title("Cumulated observations")

# Apply the extended Kalman filter
ss <- kal.fil.mul.pom(ss)

# Plot the filtered states.
plot(ss$filtered$mt[,6], xlab="time", 
     ylim=range(ss$filtered$mt), type="l")
for(i in 2:6) lines(ss$filtered$mt[,i])
title("Filtered latent process")

# Apply the iterated extended Kalman filter and smoother
ss <- ieks(ss)

# Smoothed trajectory of beta_t
upper <- ss$smoothed$m.tilde[,6] + sqrt(ss$smoothed$C.tilde[6,6,])
lower <- ss$smoothed$m.tilde[,6] - sqrt(ss$smoothed$C.tilde[6,6,])
plot(ss$smoothed$m.tilde[,6], type="l", lty=1, ylim=range(upper,lower))
lines(ss$simulated$theta[,6], lty=2)
lines(upper, lty=3)
lines(lower, lty=3)
title("Smoothed trajectory of beta_t")

# Plot the smoothed states.
plot(ss$smoothed$m.tilde[,6], xlab="time", 
     ylim=range(ss$smoothed$m.tilde), type="l")
for(i in 1:5) lines(ss$smoothed$m.tilde[,i], lty=2)
title("Smoothed latent process")