R: Markov Modulated GLM Object

mmglm {HiddenMarkov}

R Documentation

Markov Modulated GLM Object

Description

These functions create Markov modulated generalised linear model objects. These functions are in development and may change, see “Under Development” below.

Usage

mmglm0(x, Pi, delta, family, link, beta, glmformula = formula(y~x1),
       sigma = NA, nonstat = TRUE, msg = TRUE)
mmglm1(y, Pi, delta, glmfamily, beta, Xdesign,
       sigma = NA, nonstat = TRUE, size = NA, msg = TRUE)
mmglmlong1(y, Pi, delta, glmfamily, beta, Xdesign, longitude, 
           sigma = NA, nonstat = TRUE, size = NA, msg = TRUE)

Arguments

`x`	a dataframe containing the observed variable (i.e. the response variable in the generalised linear model) and the covariate. The function `mmglm0` requires that the response variable be named `y` and the covariate `x1`. Alternatively, `x` could be specified as `NULL`, meaning that the data will be added later (e.g. simulated). See Details below for the binomial case. The functions `mmglm1` and `mmglmlong1` do not have these naming restrictions.
`y`	numeric vector, response variable. In the case of binomial, it is the number of successes (see argument `size`).
`Pi`	is the `m \times m` transition probability matrix of the hidden Markov chain.
`delta`	is the marginal probability distribution of the `m` hidden states at the first time point.
`family`	character string, the GLM family, one of `"gaussian"`, `"poisson"`, `"Gamma"` or `"binomial"`.
`link`	character string, the link function. If `family == "binomial"`, then one of `"logit"`, `"probit"` or `"cloglog"`; else one of `"identity"`, `"inverse"` or `"log"`.
`glmfamily`	a `family` object defining the glm family and link function. It is currently restricted to Gaussian, Poisson, Binomial or Gamma models with the standard link functions provided by `glm`.
`Xdesign`	a `nN \times p` design matrix, where `p` is the number of parameters in the linear predictor, `N` is the number of subjects (`N=1` in `mmglm1`), and `n` is the number of observations for each subject (assumed to be the same).
`beta`	a `p \times m` matrix containing parameter values, used as initial values during estimation. In the case of the simple regression model of `mmglm0`, `p=2`. In the case of `mmglm1` and `mmglmlong1`, `p` is the number of columns of `Xdesign`.
`glmformula`	the only model formula for `mmglm0` is `y~x1`. Note that the functions `mmglm1` and `mmglmlong1` do not have this restriction, however, in those cases, the model formula is currently implicitly defined through `Xdesign`.
`sigma`	if `family == "gaussian"`, then it is the variance; if `family == "Gamma"`, then it is `1/sqrt(shape)`. It is of length `m` for each Markov state.
`nonstat`	is logical, `TRUE` if the homogeneous Markov chain is assumed to be non-stationary, default.
`longitude`	a vector the same length as `y` identifying the subject for each observation. The observations must be grouped by subject, and ordered by “time” within subject.
`size`	is number of Bernoulli trials in each observation when the glm `family` is binomial. It is the same length as `y`.
`msg`	is logical, suppress messages about developmental status.

Details

This family of models is similar in nature to those of the class dthmm, in that both classes have the distribution of the observed variable being “modulated” by the changing hidden Markov state. They differ slightly in the mechanism. This family assumes that the mean of the observation distribution can be expressed as a linear model of other known variables, but it is the parameters in the linear predictor that are being modulated by the hidden Markov process, thus causing the changes in the observed means. The linear model is assumed to be a generalised linear model as described by McCullagh & Nelder (1989).

The function mmglm0 is a very simple trivial case where the linear predictor is of the form \beta_0 + \beta_1 x_1. The version mmglm1 does not have this limitation. The model formula for mmglm1 is defined implicitly through the structure of the specified design matrix. The model mmglmlong1 is similar to mmglm1 but can be applied to longitudinal observations. Models of the form given by mmglm1 are assumed to have one time series, and from a theoretical perspective, one would be interested in the asymptotic properties of the parameter estimates as the series length gets very large. In the longitudinal case (mmglmlong1), the series of observations per individual is probably very small (<10), and hence interest is in the asymptotic properties as the number of individuals becomes large. Note that in the longitudinal case, the number of observations per individual is assumed to be the same. The responses are assumed to be conditionally independent given the value of the Markov chain and the explanatory variables in the linear predictor.

If family == "binomial" then the response variable y is interpreted as the number of successes. The dataframe x must also contain a variable called size being the number of Bernoulli trials. This is different to the format used by the function glm where y would be a matrix with two columns containing the number of successes and failures, respectively. The different format here allows one to specify the number of Bernoulli trials only so that the number of successes or failures can be simulated later.

When the density function of the response variable is from the exponential family (Charnes et al, 1976, Eq. 2.1), the likelihood function (Charnes et al, 1976, Eq. 2.4) can be maximised by using iterative weighted least squares (Charnes et al, 1976, Eq. 1.1 and 1.2). This is the method used by the R function glm. In this Markov modulated version of the model, the third term of the complete data log-likelihood, as given in Harte (2006, Sec. 2.3), needs to be maximised. This is simply the sum of the individual log-likelihood contributions of the response variable weighted by the Markov state probabilities calculated in the E-step. This can also be maximised using iterative least squares by passing these additional weights (Markov state probabilities) into the glm function.

Value

A list object with class "mmglm0", containing the above arguments as named components.

Under Development

These functions are still being developed. In previous releases of the package (< 1.3), there was only one function called mmglm. This has been renamed to mmglm0. The most recent version is mmglm1 along with mmglmlong1 which has flexibility to include longitudinal data. Further development versions will be numbered sequentially. The name mmglm has been reserved for the final stable version, at which point the numbered versions will become deprecated.

References

Cited references are listed on the HiddenMarkov manual page.

Examples

#--------------------------------------------------------
#     Gaussian with identity link function
#         using mmglm0

delta <- c(0,1)

Pi <- matrix(c(0.8, 0.2,
               0.3, 0.7),
             byrow=TRUE, nrow=2)

beta <- matrix(c(0.1, -0.1,
                 1.0,  5.0),
               byrow=TRUE, nrow=2)

x <- mmglm0(NULL, Pi, delta, family="gaussian", link="identity",
            beta=beta, sigma=c(1, 2))

n <- 1000
x <- simulate(x, nsim=n, seed=10)

#   Increase maxiter below to achieve convergence
#   Has been restricted to minimise time of package checks
y <- BaumWelch(x, bwcontrol(maxiter=2))

w <- hist(residuals(y))
z <- seq(-3, 3, 0.01)
points(z, dnorm(z)*n*(w$breaks[2]-w$breaks[1]), col="red", type="l")
box()

print(summary(y))
print(logLik(y))


#--------------------------------------------------------
#    Gaussian with log link function
#         using mmglm1

n <- 1000

#   the range of x needs changing according to the glmfamily
x <- seq(-0.9, 1.5, length.out=n)

colour <- c("blue", "green", "red")
colnum <- rep(1:3, n/3+1)[1:n] - 1

data <- data.frame(x=x, colour=colour[colnum+1])

#   will simulate response variable, not required in formula
#   design matrix only depends on RHS of formula
glmformula <- formula( ~ x + I(x^2) + colour)
glmfamily <- gaussian(link="log")
Xdesign <- model.matrix(glmformula, data=data)

# --- Parameter Values and Simulation ---

Pi <- matrix(c(0.8, 0.2,
               0.3, 0.7),
             byrow=TRUE, nrow=2)

delta <- c(1, 0)

sd <- c(1.2, 1)

beta <- matrix(c(-1, -1.2,
                 -2, -1.8,
                  3,  2.8,
                  1,  0.8, 
                  2,  2.2), 
               ncol=ncol(Pi), nrow=ncol(Xdesign), byrow=TRUE)

y <- mmglm1(NULL, Pi, delta, glmfamily, beta, Xdesign, sigma=sd)

y <- simulate(y, seed=5)

# --- Estimation ---

#   Increase maxiter below to achieve convergence
#   Has been restricted to minimise time of package checks
tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=2))
print(summary(tmp))


#-------------------------------------------------
#    Binomial with logit link function
#         using mmglm1

#   n = series length
n <- 1000

#   the range of x need changing according to the glmfamily
x <- seq(-1, 1.5, length.out=n)

colour <- c("blue", "green", "red")
colnum <- rep(1:3, n/3+1)[1:n] - 1

data <- data.frame(x=x, colour=colour[colnum+1])

glmformula <- formula( ~ x + I(x^2) + colour)
glmfamily <- binomial(link="logit")
Xdesign <- model.matrix(glmformula, data=data)

# --- Parameter Values and Simulation ---

Pi <- matrix(c(0.8, 0.2,
               0.3, 0.7),
             byrow=TRUE, nrow=2)

delta <- c(1, 0)

beta <- matrix(c(-1, -1.2,
                 -2, -1.8,
                  3,  2.8,
                  1,  0.8, 
                  2,  2.2), 
               ncol=ncol(Pi), nrow=ncol(Xdesign), byrow=TRUE)

y <- mmglm1(NULL, Pi, delta, glmfamily, beta, Xdesign, sigma=sd,
            size=rep(100, n))

#   each element of y$y is the number of successes in 100 Bernoulli trials
y <- simulate(y, seed=5)


# --- Estimation ---

#   Increase maxiter below to achieve convergence
#   Has been restricted to minimise time of package checks
tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=2))
print(summary(tmp))


#-------------------------------------------------
#    Gaussian with log link function, longitudinal data
#         using mmglmlong1

#   n = series length for each subject
#   N = number of subjects
n <- 5
N <- 1000

#   the range of x need changing according to the glmfamily
x <- seq(-0.9, 1.5, length.out=n)

colour <- c("blue", "green", "red")
colnum <- rep(1:3, n/3+1)[1:n] - 1

data <- data.frame(x=x, colour=colour[colnum+1])

#   will simulate response variable, not required in formula
#   design matrix only depends on RHS of formula
glmformula <- formula( ~ x + I(x^2) + colour)
glmfamily <- gaussian(link="log")
Xdesign0 <- model.matrix(glmformula, data=data)

#    multiple subjects
Xdesign <- NULL
for (i in 1:N) Xdesign <- rbind(Xdesign, Xdesign0)

# --- Parameter Values and Simulation ---

Pi <- matrix(c(0.8, 0.2,
               0.3, 0.7),
             byrow=TRUE, nrow=2)

delta <- c(0.5, 0.5)

sd <- c(1.2, 1)

beta <- matrix(c(-1, -1.2,
                 -2, -1.8,
                  3,  2.8,
                  1,  0.8, 
                  2,  2.2), 
               ncol=ncol(Pi), nrow=ncol(Xdesign), byrow=TRUE)

y <- mmglmlong1(NULL, Pi, delta, glmfamily, beta, Xdesign, sigma=sd,
                longitude=rep(1:N, each=n))

y <- simulate(y, seed=5)

# --- Estimation ---

#    Note: the "Not run" blocks below are not run during package checks
#    as the makePSOCKcluster definition is specific to my network,
#    modify accordingly if you want parallel processing.

cl <- NULL
## Not run: 
if (require(parallel)){
    cl <- makePSOCKcluster(c("localhost", "horoeka.localdomain", 
                             "horoeka.localdomain", "localhost"))
}
## End(Not run)

#   Increase maxiter below to achieve convergence
#   Has been restricted to minimise time of package checks
tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=2),
                 PSOCKcluster=cl)

## Not run: 
if (!is.null(cl)){
    stopCluster(cl)
    rm(cl)
}
## End(Not run)

print(summary(tmp))


#-------------------------------------------------
#    Binomial with logit link function, longitudinal data
#         using mmglmlong1

#   n = series length for each subject
#   N = number of subjects
n <- 10
N <- 100

#   the range of x need changing according to the glmfamily
x <- seq(-1, 1.5, length.out=n)

colour <- c("blue", "green", "red")
colnum <- rep(1:3, n/3+1)[1:n] - 1

data <- data.frame(x=x, colour=colour[colnum+1])

glmformula <- formula( ~ x + I(x^2) + colour)
glmfamily <- binomial(link="logit")
Xdesign0 <- model.matrix(glmformula, data=data)

#    multiple subjects
Xdesign <- NULL
for (i in 1:N) Xdesign <- rbind(Xdesign, Xdesign0)

# --- Parameter Values and Simulation ---

Pi <- matrix(c(0.8, 0.2,
               0.3, 0.7),
             byrow=TRUE, nrow=2)

delta <- c(0.5, 0.5)

beta <- matrix(c(-1, -1.2,
                 -2, -1.8,
                  3,  2.8,
                  1,  0.8, 
                  2,  2.2), 
               ncol=ncol(Pi), nrow=ncol(Xdesign), byrow=TRUE)

y <- mmglmlong1(NULL, Pi, delta, glmfamily, beta, Xdesign, sigma=sd,
                longitude=rep(1:N, each=n), size=rep(200, N*n))

y <- simulate(y, seed=5)

# --- Estimation ---

#   Increase maxiter below to achieve convergence
#   Has been restricted to minimise time of package checks
tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=1))
print(summary(tmp))

[Package HiddenMarkov version 1.8-13 Index]