mmglm {HiddenMarkov} | R Documentation |
These functions create Markov modulated generalised linear model objects. These functions are in development and may change, see “Under Development” below.
mmglm0(x, Pi, delta, family, link, beta, glmformula = formula(y~x1), sigma = NA, nonstat = 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)
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*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)*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*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. |
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 beta0 + beta1*x1. 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.
A list
object with class "mmglm0"
, containing the above arguments as named components.
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.
The functions mmglm
and mmglmlong1
currently have no methods for the generic functions residuals
and logLik
.
#-------------------------------------------------------- # 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)) x <- simulate(x, nsim=5000, seed=10) y <- BaumWelch(x) hist(residuals(y)) 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 --- tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, tol=1e-05)) print(summary(tmp)) #------------------------------------------------- # Binomial with logit link function # using mmglm1 # n = series length for each subject n <- 5000 # 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)) y <- simulate(y, seed=5) # --- Estimation --- tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=100)) 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 --- # library(snow) # cl <- makeSOCKcluster(c("localhost", "localhost", # "horoeka.localdomain", "horoeka.localdomain")) cl <- NULL tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, tol=1e-03), SNOWcluster=cl) # stopCluster(cl) 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 --- tmp <- BaumWelch(y, bwcontrol(posdiff=FALSE, maxiter=500)) print(summary(tmp))