[Depmix-commits] r434 - in pkg/depmixS4: . inst inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 4 16:39:00 CEST 2010
Author: ingmarvisser
Date: 2010-08-04 16:39:00 +0200 (Wed, 04 Aug 2010)
New Revision: 434
Removed:
pkg/depmixS4/inst/CITATION
pkg/depmixS4/inst/doc/dpx4Rev.Rnw
pkg/depmixS4/inst/doc/graphs/
Modified:
pkg/depmixS4/DESCRIPTION
pkg/depmixS4/NEWS
pkg/depmixS4/README
Log:
Updated news/readme etc
Modified: pkg/depmixS4/DESCRIPTION
===================================================================
--- pkg/depmixS4/DESCRIPTION 2010-08-04 14:26:44 UTC (rev 433)
+++ pkg/depmixS4/DESCRIPTION 2010-08-04 14:39:00 UTC (rev 434)
@@ -1,7 +1,7 @@
Package: depmixS4
Version: 1.0-0
-Date: 2010-06-20
-Title: Dependent Mixture Models in S4 - Hidden Markov Models
+Date: 2010-08-04
+Title: Dependent Mixture Models - Hidden Markov Models of GLMs and Other Distributions in S4
Author: Ingmar Visser <i.visser at uva.nl>, Maarten Speekenbrink <m.speekenbrink at ucl.ac.uk>
Maintainer: Ingmar Visser <i.visser at uva.nl>
Depends: R (>= 2.11.1), stats, nnet, methods, MASS, Rsolnp
Modified: pkg/depmixS4/NEWS
===================================================================
--- pkg/depmixS4/NEWS 2010-08-04 14:26:44 UTC (rev 433)
+++ pkg/depmixS4/NEWS 2010-08-04 14:39:00 UTC (rev 434)
@@ -1,9 +1,9 @@
Changes in depmixS4 version 1.0-0
- o
+ o added a vignette to the package and upped the version number 1.0-0 to
+ celebrate publication in the Journal of Statistical Software.
-
Changes in depmixS4 version 0.9-1
o fixed a bug in setting the lower and upper bounds for GLMresponse
Modified: pkg/depmixS4/README
===================================================================
--- pkg/depmixS4/README 2010-08-04 14:26:44 UTC (rev 433)
+++ pkg/depmixS4/README 2010-08-04 14:39:00 UTC (rev 434)
@@ -10,6 +10,7 @@
the initial state probabilities are also modeled as multinomial logistics
(or multinomials with identity link, which is the default when no
covariates are present) with the possibility of including covariates.
+
Optimization is by default done using the EM algorithm. When linear
constraints are included, package Rsolnp is used for optimization (there is
also support for using Rdonlp2 as optimizer, see details below). New
Deleted: pkg/depmixS4/inst/CITATION
===================================================================
--- pkg/depmixS4/inst/CITATION 2010-08-04 14:26:44 UTC (rev 433)
+++ pkg/depmixS4/inst/CITATION 2010-08-04 14:39:00 UTC (rev 434)
@@ -1,12 +0,0 @@
-citHeader("To cite package depmixS4 in publications, use the following:")
-
-citEntry(entry="Article",
- title = "{depmixS4}: An R-package for Hidden Markov Models",
- author = personList(person(first="Ingmar", last="Visser"),
- person(first="Maarten", last="Speekenbrink")),
- journal = "Journal of Statistical Software",
- year = "2010",
- url = "http://cran.r-project.org",
- textVersion = paste("Ingmar Visser & Maarten Speekenbrink (2010).
- depmixS4: An R-package for Hidden Markov Models")
-)
Deleted: pkg/depmixS4/inst/doc/dpx4Rev.Rnw
===================================================================
--- pkg/depmixS4/inst/doc/dpx4Rev.Rnw 2010-08-04 14:26:44 UTC (rev 433)
+++ pkg/depmixS4/inst/doc/dpx4Rev.Rnw 2010-08-04 14:39:00 UTC (rev 434)
@@ -1,1077 +0,0 @@
-% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
-
-\documentclass[nojss]{jss}
-\usepackage{amsmath}
-
-%\usepackage[]{amsmath, amsfonts, amstext, amsthm}
-%\usepackage{amssymb}
-%\usepackage[]{graphics}
-%\usepackage{epsfig}
-%\usepackage{epstopdf}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%% almost as usual
-\author{Ingmar Visser\\University of Amsterdam \And
- Maarten Speekenbrink\\University College London}
-
-\title{\pkg{depmixS4} : An \proglang{R}-package for Hidden Markov Models}
-
-% \thanks{This article was originally published in the {\em Journal of
-% Statistical Software}, volume ??, etc. Please refer to that paper
-% when citing this package.}
-
-
-%% for pretty printing and a nice hypersummary also set:
-\Plainauthor{Ingmar Visser, Maarten Speekenbrink} %% comma-separated
-\Plaintitle{depmixS4: An R-package for Hidden Markov Models} %% without formatting
-
-\Shorttitle{depmixS4: Hidden Markov Models} %% a short title (if necessary)
-
-%% an abstract and keywords
-\Abstract{
-
- \pkg{depmixS4} implements a general framework for defining and
- estimating dependent mixture models in the \proglang{R}
- programming language \citep{R2010}. This includes standard Markov
- models, latent/hidden Markov models, and latent class and finite
- mixture distribution models. The models can be fitted on mixed
- multivariate data with distributions from the \code{glm} family,
- the (logistic) multinomial, or the multivariate normal
- distribution. Other distributions can be added easily, and an
- example is provided with the {\em exgaus} distribution.
- Parameters are estimated by the EM algorithm or, when (linear)
- constraints are imposed on the parameters, by direct numerical
- optimization with the \pkg{Rsolnp} or \pkg{Rdonlp2} routines.
-
-}
-
-%\VignetteIndexEntry{hidden Markov model}
-
-\Keywords{hidden Markov model, dependent mixture model, mixture
-model, constraints}
-
-\Plainkeywords{hidden Markov model, dependent mixture model, mixture
-model, constraints} %% without formatting
-%% at least one keyword must be supplied
-
-%% publication information
-%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
-%% \Volume{13}
-%% \Issue{9}
-%% \Month{September}
-%% \Year{2004}
-%% \Submitdate{2004-09-29}
-%% \Acceptdate{2004-09-29}
-
-%% The address of (at least) one author should be given
-%% in the following format:
-\Address{
-Ingmar Visser\\
-Department of Psychology\\
-University of Amsterdam\\
-Roetersstraat 15\\
-1018 WB, Amsterdam\\
-The Netherlands\\
-E-mail: \email{i.visser at uva.nl} \\
-URL: \url{http://www.ingmar.org/}
-}
-
-%% It is also possible to add a telephone and fax number
-%% before the e-mail in the following format:
-%% Telephone: +43/1/31336-5053
-%% Fax: +43/1/31336-734
-
-\newcommand{\vc}{\mathbf}
-\newcommand{\mat}{\mathbf}
-\newcommand{\greekv}[1]{\mbox{\boldmath$\mathrm{#1}$}}
-\newcommand{\greekm}[1]{\mbox{\boldmath$#1$}}
-
-%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%\batchmode
-
-\SweaveOpts{echo=TRUE}
-\usepackage{a4wide}
-
-%% \usepackage{Sweave}
-
-\begin{document}
-
-
-%set width of figures produced by Sweave
-\setkeys{Gin}{width=0.9\textwidth}
-
-\maketitle
-
-<<echo=FALSE, results=HIDE>>=
-.Options$digits<-3
-library(depmixS4)
-@
-
-
-NOTE: This article was originally published in the {\em Journal of
-Statistical Software}, volume ??, etc. Please refer to that paper
-when citing this package.
-
-
-%% include your article here, just as usual
-%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.
-
-\section{Introduction}
-%% Note: If there is markup in \(sub)section, then it has to be escape as above.
-
-Markov and latent Markov models are frequently used in the social
-sciences, in different areas and applications. In psychology, they
-are used for modelling learning processes (see \citealp{Wickens1982},
-for an overview, and e.g. \citealp{Schmittmann2006}, for a recent
-application). In economics, latent Markov models are so-called regime
-switching models (see e.g., \citealp{Kim1994} and
-\citealp{Ghysels1994}). Further applications include speech
-recognition \citep{Rabiner1989}, EEG analysis \citep{Rainer2000}, and
-genetics \citep{Krogh1998}. In these latter areas of application,
-latent Markov models are usually referred to as hidden Markov models.
-See for example \citet{Fruhwirth2006} for an overview of hidden Markov
-models with extensions. Further examples of applications can be found
-in e.g.\ \citet[][chapter~1]{Cappe2005}. A more gentle introduction
-into hidden Markov models with applications is the book by
-\citet{Zucchini2009}.
-
-The \pkg{depmixS4} package was motivated by the fact that while Markov
-models are used commonly in the social sciences, no comprehensive
-package was available for fitting such models. Existing software for
-estimating Markovian models include Panmark \citep{Pol1996}, and for
-latent class models Latent Gold \citep{Vermunt2003}. These programs
-lack a number of important features, besides not being freely
-available. There are currently some packages in \proglang{R} that
-handle hidden Markov models but they lack a number of features that we
-needed in our research. In particular, \pkg{depmixS4} was designed to
-meet the following goals:
-
-\begin{enumerate}
-
- \item to be able to estimate parameters subject to general
- linear (in)equality constraints
-
- \item to be able to fit transition models with covariates, i.e.,
- to have time-dependent transition matrices
-
- \item to be able to include covariates in the prior or initial
- state probabilities
-
- \item to be easily extensible, in particular, to allow users to
- easily add new uni- or multivariate response distributions and
- new transition models, e.g., continuous time observation models
-
-\end{enumerate}
-
-Although \pkg{depmixS4} was designed to deal with longitudinal or time
-series data, for say $T>100$, it can also handle the limit case when
-$T=1$. In this case, there are no time dependencies between observed
-data and the model reduces to a finite mixture or latent class model.
-While there are specialized packages to deal with mixture data, as far
-as we know these don't allow the inclusion of covariates on the prior
-probabilities of class membership. The possibility to estimate the
-effects of covariates on prior and transition probabilities is a
-distinguishing feature of \pkg{depmixS4}. In the next section, we
-provide an outline of the model and likelihood equations.
-
-
-\section{The dependent mixture model}
-
-The data considered here have the general form $\vc{O}_{1:T}=
-(O_{1}^{1}, \ldots, O_{1}^{m}$, $O_{2}^{1}, \ldots, O_{2}^{m}$,
-\ldots, $O_{T}^{1}, \ldots, O_{T}^{m})$ for an $m$-variate time series
-of length $T$. In the following, we use $\vc{O}_{t}$ as shorthand for
-$O_{t}^{1}, \ldots, O_{t}^{m}$. As an example, consider a time series
-of responses generated by a single participant in a psychological
-response time experiment. The data consists of three variables:
-response time, response accuracy, and a covariate which is a pay-off
-variable reflecting the relative reward for speeded and/or accurate
-responding. These variables are measured on 168, 134 and 137
-occasions respectively (the first series of 168 trials is plotted in
-Figure~\ref{fig:speed}). These data are more fully described in
-\citet{Dutilh2010}, and in the next section a number of example models
-for these data is described.
-
-\begin{figure}[htbp]
-\begin{center}
-<<fig=TRUE,echo=FALSE,height=5,width=7,eps=FALSE>>=
-data(speed)
-plot(as.ts(speed[1:168,]),main="Speed-accuracy trade-off")
-@
- \caption{Response times (rt), accuracy (corr) and pay-off
- values (Pacc) for the first series of responses in dataset
- \code{speed}.}
- \label{fig:speed}
-\end{center}
-\end{figure}
-
-The latent Markov model is usually associated with data of this type,
-in particular for multinomially distributed responses. However,
-commonly employed estimation procedures \citep[e.g.,][]{Pol1996}, are
-not suitable for long time series due to underflow problems. In
-contrast, the hidden Markov model is typically only used for `long'
-univariate time series \citep[][chapter~1]{Cappe2005}. We use the
-term "dependent mixture model" because one of the authors (Ingmar
-Visser) thought it was time for a new name to relate these
-models\footnote{Only later did I find out that \citet{Leroux1992}
-already coined the term dependent mixture models in an application
-with hidden Markov mixtures of Poisson count data.}.
-
-The fundamental assumption of a dependent mixture model is that at any
-time point, the observations are distributed as a mixture with $n$
-components (or states), and that time-dependencies between the
-observations are due to time-dependencies between the mixture
-components (i.e., transition probabilities between the components).
-These latter dependencies are assumed to follow a first-order Markov
-process. In the models we are considering here, the mixture
-distributions, the initial mixture probabilities and transition
-probabilities can all depend on covariates $\vc{z}_t$.
-
-In a dependent mixture model, the joint likelihood of observations
-$\vc{O}_{1:T}$ and latent states $\vc{S}_{1:T} = (S_1,\ldots,S_T)$,
-given model parameters $\greekv{\theta}$ and covariates $\vc{z}_{1:T}
-= (\vc{z}_1,\ldots,\vc{z}_T)$, can be written as:
-\begin{equation}
- \Prob(\vc{O}_{1:T},\vc{S}_{1:T}|\greekv{\theta},\vc{z}_{1:T}) =
- \pi_{i}(\vc{z}_1) \vc{b}_{S_t}(\vc{O}_1|\vc{z}_{1})
- \prod_{t=1}^{T-1} a_{ij}(\vc{z}_{t}) \vc{b}_{S_t}(\vc{O}_{t+1}|\vc{z}_{t+1}),
-\end{equation}
-where we have the following elements:
-\begin{enumerate}
-
- \item $S_{t}$ is an element of $\mathcal{S}=\{1\ldots n\}$, a set
- of $n$ latent classes or states.
-
- \item $\pi_{i}(\vc{z}_1) = \Prob(S_1 = i|\vc{z}_1)$, giving the
- probability of class/state $i$ at time $t=1$ with covariate
- $\vc{z}_1$.
-
- \item $a_{ij}(\vc{z}_t) = \Prob(S_{t+1}=j|S_{t}=i,\vc{z}_t)$,
- provides the probability of a transition from state $i$ to state
- $j$ with covariate $\vc{z}_t$,
-
- \item $\vc{b}_{S_t}$ is a vector of observation densities
- $b_{j}^k(\vc{z}_t) = \Prob(O_{t}^k|S_t = j, \vc{z}_t)$ that
- provide the conditional densities of observations $O_{t}^k$
- associated with latent class/state $j$ and covariate $\vc{z}_t$,
- $j=1, \ldots, n$, $k=1, \ldots, m$.
-
-\end{enumerate}
-
-For the example data above, $b_j^k$ could be a Gaussian distribution
-function for the response time variable, and a Bernoulli distribution
-for the accuracy variable. In the models considered here,
-both the transition probability functions $a_{ij}$ and the initial
-state probability functions $\greekv{\pi}$ may depend on covariates as
-well as the response distributions $b_{j}^{k}$.
-
-
-\subsection{Likelihood}
-
-To obtain maximum likelihood estimates of the model parameters, we
-need the marginal likelihood of the observations. For hidden Markov
-models, this marginal (log-)likelihood is usually computed by the
-so-called forward-backward algorithm \citep{Baum1966,Rabiner1989}, or
-rather by the forward part of this algorithm. \cite{Lystig2002}
-changed the forward algorithm in such a way as to allow computing the
-gradients of the log-likelihood at the same time. They start by
-rewriting the likelihood as follows (for ease of exposition the
-dependence on the model parameters and covariates is dropped here):
-\begin{equation}
- L_{T} = \Prob(\vc{O}_{1:T}) = \prod_{t=1}^{T}
- \Prob(\vc{O}_{t}|\vc{O}_{1:(t-1)}),
- \label{condLike}
-\end{equation}
-where $\Prob(\vc{O}_{1}|\vc{O}_{0}):=\Prob(\vc{O}_{1})$. Note that for a
-simple, i.e.\ observed, Markov chain these probabilities reduce to
-$\Prob(\vc{O}_{t}|\vc{O}_{1},\ldots,
-\vc{O}_{t-1})=\Prob(\vc{O}_{t}|\vc{O}_{t-1})$.
-The log-likelihood can now be expressed as:
-\begin{equation}
- l_{T} = \sum_{t=1}^{T} \log[\Prob(\vc{O}_{t}|\vc{O}_{1:(t-1)}].
- \label{eq:condLogl}
-\end{equation}
-
-To compute the log-likelihood, \cite{Lystig2002} define the following
-(forward) recursion:
-\begin{align}
- \phi_{1}(j) &:= \Prob(\vc{O}_{1}, S_{1}=j) = \pi_{j} \vc{b}_{j}(\vc{O}_{1})
- \label{eq:fwd1} \\
- \phi_{t}(j) &:= \Prob(\vc{O}_{t}, S_{t}=j|\vc{O}_{1:(t-1)}) %\\
- = \sum_{i=1}^{N} [\phi_{t-1}(i)a_{ij} \vc{b}_{j}(\vc{O}_{t})] \times
-(\Phi_{t-1})^{-1},
- \label{eq:fwdt}
-\end{align}
-where $\Phi_{t}=\sum_{i=1}^{N} \phi_{t}(i)$. Combining
-$\Phi_{t}=\Prob(\vc{O}_{t}|\vc{O}_{1:(t-1)})$, and
-equation~(\ref{eq:condLogl}) gives the following expression for the
-log-likelihood:
-\begin{equation}
- l_{T} = \sum_{t=1}^{T} \log \Phi_{t}.
- \label{eq:logl}
-\end{equation}
-
-
-\subsection{Parameter estimation}
-
-Parameters are estimated in \pkg{depmixS4} using the EM algorithm or
-through the use of a general Newton-Raphson optimizer. In the EM
-algorithm, parameters are estimated by iteratively maximising the
-expected joint log-likelihood of the parameters given the observations and
-states. Let $\greekv{\theta} = (\greekv{\theta}_1,
-\greekv{\theta}_2,\greekv{\theta}_3)$ be the general parameter vector
-consisting of three subvectors with parameters for the prior model,
-transition model, and response models respectively. The joint
-log-likelihood can be written as:
-\begin{equation}
-\log \Prob(\vc{O}_{1:T}, \vc{S}_{1:T}|\vc{z}_{1:T},\greekv{\theta}) = \log
-\Prob(S_1|\vc{z}_{1},\greekv{\theta}_1)
-+ \sum_{t=2}^{T} \log \Prob(S_t|S_{t-1},\vc{z}_{t-1},\greekv{\theta}_2)
-+ \sum_{t=1}^{T} \log \Prob(\vc{O}_t|S_t,\vc{z}_t,\greekv{\theta}_3)
-\end{equation}
-This likelihood depends on the unobserved states $\vc{S}_{1:T}$. In the
-Expectation step, we replace these with their expected values given a set of
-(initial) parameters $\greekv{\theta}' = (\greekv{\theta}'_1,
-\greekv{\theta}'_2,\greekv{\theta}'_3)$ and observations $\vc{O}_{1:T}$.
-The expected log-likelihood:
-\begin{equation}
-Q(\greekv{\theta},\greekv{\theta}') = E_{\greekv{\theta}'}
-(\log \Prob(\vc{O}_{1:T},\vc{S}_{1:T}|\vc{O}_{1:T},\vc{z}_{1:T},\greekv{\theta})),
-\end{equation}
-can be written as:
-\begin{multline}
-\label{eq:Q}
-Q(\greekv{\theta},\greekv{\theta}') =
-\sum_{j=1}^n \gamma_1(j) \log \Prob(S_1=j|\vc{z}_1,\greekv{\theta}_1) \\
-+ \sum_{t=2}^T \sum_{j=1}^n \sum_{k=1}^n \xi_t(j,k) \log \Prob(S_t = k|S_{t-1}
-= j,\vc{z}_{t-1},\greekv{\theta}_2) \\
- + \sum_{t=1}^T \sum_{j=1}^n \sum_{k=1}^m \gamma_t(j)
-\log \Prob(O^k_t|S_t=j,\vc{z}_t,\greekv{\theta}_3),
-\end{multline}
-where the expected values $\xi_t(j,k) = P(S_t = k, S_{t-1} =
-j|\vc{O}_{1:T}, \vc{z}_{1:T},\greekv{\theta}')$ and $\gamma_t(j) =
-P(S_t = j|\vc{O}_{1:T}, \vc{z}_{1:T},\greekv{\theta}')$ can be
-computed effectively by the forward-backward algorithm \citep[see
-e.g.,][]{Rabiner1989}. The Maximisation step consists of the
-maximisation of (\ref{eq:Q}) for $\greekv{\theta}$. As the right hand
-side of (\ref{eq:Q}) consists of three separate parts, we can maximise
-separately for $\greekv{\theta}_1$, $\greekv{\theta}_2$ and
-$\greekv{\theta}_3$. In common models, maximisation for
-$\greekv{\theta}_1$ and $\greekv{\theta}_2$ is performed by the
-\code{nnet.default} routine in the \pkg{nnet} package
-\citep{Venables2002}, and maximisation for $\greekv{\theta}_3$ by the
-standard \code{glm} routine. Note that for the latter maximisation,
-the expected values $\gamma_t(j)$ are used as prior weights of the
-observations $O^k_t$.
-
-The EM algorithm however has some drawbacks. First, it can be slow to
-converge towards the end of optimization. Second, applying
-constraints to parameters can be problematic; in particular, EM can
-lead to wrong parameter estimates when applying constraints. Hence,
-in \pkg{depmixS4}, EM is used by default in unconstrained models, but
-otherwise, direct optimization is used. Two options are available for
-direct optimization using package \pkg{Rsolnp}
-\citep{Ghalanos2010,Ye1987}, or \pkg{Rdonlp2}
-\citep{Tamura2009,Spellucci2002}. Both packages can handle general
-linear (in)equality constraints (and optionally also non-linear
-constraints).
-
-
-\section[Using depmixS4]{Using \pkg{depmixS4}}
-
-Two steps are involved in using \pkg{depmixS4} which are illustrated
-below with examples:
-\begin{enumerate}
-
- \item model specification with function \code{depmix} (or with
- \code{mix} for latent class and finite mixture models, see example
- below on adding covariates to prior probabilities)
-
- \item model fitting with function \code{fit}
-\end{enumerate}
-We have separated the stages of model specification and model fitting
-because fitting large models can be fairly time-consuming and it is
-hence useful to be able to check the model specification before
-actually fitting the model.
-
-\subsection[Example data: speed]{Example data: \code{speed}}
-
-Throughout this article a data set called \code{speed} is used. As
-already indicated in the Introduction, it consists of three time
-series with three variables: response time, accuracy, and a covariate,
-Pacc, which defines the relative pay-off for speeded versus accurate
-responding. Before describing some of the models that are fitted to
-these data, we provide a brief sketch of the reasons for gathering
-these data in the first place.
-
-Response times are a very common dependent variable in psychological
-experiments and hence form the basis for inference about many
-psychological processes. A potential threat to such inference based
-on response times is formed by the speed-accuracy trade-off: different
-participants in an experiment may respond differently to typical
-instructions to `respond as fast and accurate as possible'. A popular
-model which takes the speed-accuracy trade-off into account is the
-diffusion model \citep{Ratcliff1978}, which has proven to provide
-accurate descriptions of response times in many different settings.
-
-One potential problem with the diffusion model is that it predicts a
-continuous trade-off between speed and accuracy of responding, i.e.,
-when participants are pressed to respond faster and faster, the
-diffusion model predicts that this would lead to a gradual decrease in
-accuracy. The \code{speed} data set that we analyze below was
-gathered to test this hypothesis versus the alternative hypothesis
-stating that there is a sudden transition from slow and accurate
-responding to fast responding at chance level. At each trial of the
-experiment, the participant is shown the current setting of the
-relative reward for speed versus accuracy. The bottom panel of
-figure~\ref{fig:speed} shows the values of this variable. The
-experiment was designed to investigate what would happen when this
-reward variable changes from reward for accuracy only to reward for
-speed only. The \code{speed} data that we analyse here are from
-participant A in Experiment 1 in \citet{Dutilh2010}, who provide a
-complete description of the experiment and the relevant theoretical
-background.
-
-The central question regarding this data is whether it is indeed best
-described by two modes of responding rather than a single mode of
-responding with a continuous trade-off between speed and accuracy.
-The hallmark of a discontinuity between slow versus speeded responding
-is that switching between the two modes is asymmetric \citep[see
-e.g.][for a theoretical underpinning of this claim]{Maas1992}. The
-\code{fit} help page of \pkg{depmixS4} provides a number of examples
-in which the asymmetry of the switching process is tested; those
-examples and other candidate models are discussed at length in
-\citet{Visser2009b}.
-
-
-\subsection{A simple model}
-
-A dependent mixture model is defined by the number of states and the
-initial state, state transition, and response distribution functions.
-A dependent mixture model can be created with the
-\code{depmix}-function as follows:
-<<>>=
-library(depmixS4)
-data(speed)
-set.seed(1)
-mod <- depmix(response=rt~1, data=speed, nstates=2,
- trst=runif(4))
-@
-
-The first line of code loads the \pkg{depmixS4} package and
-\code{data(speed)} loads the \code{speed} data set. The line
-\code{set.seed(1)} is necessary to get starting values that will
-result in the right model, see more on starting values below.
-
-The call to \code{depmix} specifies the model with a number of
-arguments. The \code{response} argument is used to specify the
-response variable, possibly with covariates, in the familiar format
-using formulae such as in \code{lm} or \code{glm} models. The second
-argument, \code{data=speed}, provides the \code{data.frame} in which
-the variables from \code{response} are interpreted. Third, the number
-of states of the model is given by \code{nstates=2}.
-
-
-\paragraph{Starting values} Note also that start values for the
-transition parameters are provided in this call by the \code{trstart}
-argument. Starting values for parameters can be provided using three
-arguments: \code{instart} for the parameters relating to the prior
-probabilities, \code{trstart}, relating the transition models, and
-\code{respstart} for the parameters of the response models. Note that
-the starting values for the initial and transition models as well as
-multinomial logit response models are interpreted as {\em
-probabilities}, and internally converted to multinomial logit
-parameters (if a logit link function is used). Start values can also
-be generated randomly within the EM algorithm by generating random
-uniform values for the values of $\gamma_{t}$ in (\ref{eq:Q}) at
-iteration 0. Automatic generation of starting values is not available
-for constrained models (see below).
-
-\paragraph{Fitting the model and printing results} The \code{depmix}
-function returns an object of class \code{depmix} which contains the
-model specification, and not a fitted model. Hence, the model needs
-to be fitted by calling \code{fit}:
-<<>>=
-fm <- fit(mod)
-@
-
-The \code{fit}-function returns an object of class
-\code{depmix.fitted} which extends the \code{depmix} class, adding
-convergence information (and information about constraints if these
-were applied, see below). The class has \code{print} and
-\code{summary} methods to see the results. The \code{print} method
-provides information on convergence, the log-likelihood and
-the AIC and BIC values:
-<<>>=
-fm
-@
-
-These statistics can also be extracted using \code{logLik}, \code{AIC}
-and \code{BIC}, respectively. By comparison, a 1-state model for
-these data, i.e. assuming there is no mixture, has a log-likelihood of
-$-305.33$, and 614.66, and 622.83 for the AIC and BIC respectively.
-Hence, the 2-state model fits the data much better than the 1-state
-model. Note that the 1-state model can be specified using \code{mod <-
-depmix(rt~1, data=speed, nstates=1)}, although this model is trivial
-as it will simply return the mean and the sd of the rt variable.
-
-The \code{summary} method of \code{fit}ted models provides the
-parameter estimates, first for the prior probabilities model, second
-for the transition models, and third for the response models.
-<<>>=
-summary(fm)
-@
-
-Since no further arguments were specified, the initial state, state
-transition and response distributions were set to their defaults
-(multinomial distributions for the first two, and a Gaussian
-distribution for the response). The resulting model indicates two
-well-separated states, one with slow and the second with fast
-responses. The transition probabilities indicate rather stable
-states, i.e., the probability of remaining in either of the states is
-around 0.9. The initial state probability estimates indicate that
-state 1 is the starting state for the process, with a negligible
-probability of starting in state 2.
-
-\subsection{Covariates on transition parameters}
-
-By default, the transition probabilities and the initial state
-probabilities are parameterized using a multinomial model with an
-identity link function. Using a multinomial logistic model allows one
-to include covariates on the initial state and transition
-probabilities. In this case, each row of the transition matrix is
-parameterized by a baseline category logistic multinomial, meaning
-that the parameter for the base category is fixed at zero
-\citep[see][p.\ 267 ff., for multinomial logistic models and various
-parameterizations]{Agresti2002}. The default baseline category is the
-first state. Hence, for example, for a 3-state model, the initial
-state probability model would have three parameters of which the first
-is fixed at zero and the other two are freely estimated.
-\citet{Chung2007} discuss a related latent transition model for
-repeated measurement data ($T=2$) using logistic regression on the
-transition parameters; they rely on Bayesian methods of estimation.
-Covariates on the transition probabilities can be specified using a
-one-sided formula as in the following example:
-<<>>=
-set.seed(1)
-mod <- depmix(rt~1, data=speed, nstates=2, family=gaussian(),
- transition=~scale(Pacc), instart=runif(2))
-fm <-fit(mod,verbose=FALSE)
-@
-
-Note the use of \code{verbose=FALSE} to suppress printing of
-information on the iterations of the fitting process. Applying
-\code{summary} to the \code{fit}ted object gives (only transition
-models printed here by using argument \code{which}):
-<<>>=
-summary(fm, which="transition")
-@
-
-The summary provides all parameters of the model, also the (redundant)
-zeroes for the base-line category in the multinomial model. The
-summary also prints the transition probabilities at the zero value of
-the covariate. Note that scaling of the covariate is useful in this
-regard as it makes interpretation of these intercept probabilities
-easier.
-
-\subsection{Multivariate data}
-
-Multivariate data can be modelled by providing a list of formulae as
-well as a list of family objects for the distributions of the various
-responses. In above examples we have only used the response times
-which were modelled as a Gaussian distribution. The accuracy
-variable in the \code{speed} data can be modelled with a multinomial
-by specifying the following:
-<<>>=
-set.seed(1)
-mod <- depmix(list(rt~1,corr~1), data=speed, nstates=2,
- family=list(gaussian(),multinomial("identity")),
- transition=~scale(Pacc),instart=runif(2))
-fm <- fit(mod,verbose=FALSE)
-@
-
-This provides the following fitted model parameters (only the
-response parameters are given here):
-<<>>=
-summary(fm,which="response")
-@
-
-As can be seen, state 1 has fast response times and accuracy is
-approximately at chance level (.474), whereas state 2 corresponds with
-slower responding at higher accuracy levels (.904).
-
-Note that by specifying multivariate observations in terms of a list,
-the variables are considered conditionally independent (given the
-states). Conditionally \emph{dependent} variables must be handled as
-a single element in the list. Effectively, this means specifying a
-multivariate response model. Currently, \pkg{depmixS4} has one
-multivariate response model which is for multivariate normal
-variables.
-
-
-\subsection{Fixing and constraining parameters}
-
-Using package \pkg{Rsolnp} \citet{Ghalanos2010} or \pkg{Rdonlp2} by
-\citet{Tamura2009}, parameters may be fitted subject to general linear
-(in-)equality constraints. Constraining and fixing parameters is done
-using the \code{conpat} argument to the \code{fit}-function, which
-specifies for each parameter in the model whether it's fixed (0) or
-free (1 or higher). Equality constraints can be imposed by giving two
-parameters the same number in the \code{conpat} vector. When only
-fixed values are required, the \code{fixed} argument can be used
-instead of \code{conpat}, with zeroes for fixed parameters and other
-values (ones e.g.) for non-fixed parameters. Fitting the models
-subject to these constraints is handled by the optimization routine
-\code{solnp} or, optionally, by \code{donlp2}. To be able to
-construct the \code{conpat} and/or \code{fixed} vectors one needs the
-correct ordering of parameters which is briefly discussed next before
-proceeding with an example.
-
-\paragraph{Parameter numbering} When using the \code{conpat} and
-\code{fixed} arguments, complete parameter vectors should be supplied,
-i.e., these vectors should have length equal to the number of parameters of
-the model, which can be obtained by calling \code{npar(object)}. Note
-that this is not the same as the degrees of freedom used e.g.\ in the
-\code{logLik} function because \code{npar} also counts the baseline
-category zeroes from the multinomial logistic models. Parameters are
-numbered in the following order:
-\begin{enumerate}
- \item the prior model parameters
- \item the parameters for the transition models
- \item the response model parameters per state (and subsequently
- per response in the case of multivariate time series)
-\end{enumerate}
-
-To see the ordering of parameters use the following:
-<<term=FALSE>>=
-setpars(mod, value=1:npar(mod))
-@
-
-To see which parameters are fixed (by default only baseline parameters
-in the multinomial logistic models for the transition models and the
-initial state probabilities model):
-<<term=FALSE>>=
-setpars(mod, getpars(mod,which="fixed"))
-@
-When fitting constraints it is useful to have good starting values
-for the parameters and hence we first fit the following model without
-constraints:
-<<>>=
-trst <- c(0.9,0.1,0,0,0.1,0.9,0,0)
-mod <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,
- nstates=2,family=list(gaussian(),multinomial("identity")),
- trstart=trst,inst=c(0.99,0.01))
-fm1 <- fit(mod,verbose=FALSE)
-@
-
-After this, we use the fitted values from this model to constrain the
-regression coefficients on the transition matrix (parameters number~6
-and~10):
-<<term=FALSE>>=
-pars <- c(unlist(getpars(fm1)))
-pars[6] <- pars[10] <- 11
-pars[1]=0
-pars[2]=1
-pars[13] <- pars[14] <- 0.5
-fm1 <- setpars(mod,pars)
-# start with fixed and free parameters
-conpat <- c(0,0,rep(c(0,1),4),1,1,0,0,1,1,1,1)
-# constrain the beta's on the transition parameters to be equal
-conpat[6] <- conpat[10] <- 2
-fm2 <- fit(fm1,equal=conpat)
-@
-
-Using \code{summary} on the fitted model shows that the regression
-coefficients are now estimated at the same value of 12.66. Setting
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 434
More information about the depmix-commits
mailing list