Nothing
R/hmm.R
In msm: Multi-State Markov and Hidden Markov Models in Continuous Time
Defines functions
msm.form.initprobs
in.range
msm.form.hranges
msm.form.hcovconstraint
msm.form.hconstraint
print_hmmcat
print.hmodel
msm.econstr2hconstr
msm.misccovinits2hcovinits
msm.misccov2hcov
msm.emodel2hmodel
msm.form.icmodel
msm.form.hcmodel
msm.form.mvhmodel
msm.form.univhmodel
msm.form.hmodel
msm.check.hcovariates
msm.check.hmodel
print.hmmdist
print.hmmMVdist
### FUNCTIONS FOR HIDDEN MARKOV MODELS IN CONTINUOUS TIME
### WITH ARBITRARY RESPONSE DISTRIBUTION
#' @export
print.hmmMVdist <- function(x, ...)
{
cat(sprintf("Multivariate hidden Markov model with %d outcomes:\n", length(x)))
for (i in x) print(i)
}
#' @export
print.hmmdist <- function(x, ...)
{
cat("Hidden Markov model", x$label, "distribution\n\n")
pnames <- if(x$label=="categorical") paste("P(",seq(x$pars[1]),")",sep="") else names(x$pars)
pars <- if(x$label=="categorical") x$pars[3:(2+x$pars[1])] else x$pars
cat("Parameters: ", paste(paste(pnames, pars, sep=" = "), collapse=", "))
cat("\n")
}
msm.check.hmodel <- function(hmodel, nstates)
{
if (!is.list(hmodel)) stop("Hidden model should be a list")
if (length(hmodel) != nstates) stop("hmodel of length ", length(hmodel), ", expected ", nstates)
for (i in hmodel) {
if (!inherits(i, "hmmdist")) stop("hmodel should be a list of HMM distribution objects")
}
}
msm.check.hcovariates <- function(hcovariates, qmodel)
{
if (!is.list(hcovariates)) stop("hcovariates should be a list")
if (length(hcovariates) != qmodel$nstates)
stop("hcovariates of length ", length(hcovariates), ", expected ", qmodel$nstates)
for (i in hcovariates) {
if (!is.null(i) )
if (!inherits(i, "formula")) stop("hcovariates should be a list of formulae or NULLs")
}
}
msm.form.hmodel <- function(hmodel, hconstraint=NULL, initprobs=NULL, est.initprobs)
{
nst <- length(hmodel)
if (is.null(initprobs))
initprobs <- if (est.initprobs) rep(1/nst, nst) else c(1, rep(0, nst-1))
else {
if (!is.numeric(initprobs)) stop("initprobs should be numeric")
if (is.matrix(initprobs)) {
if (ncol(initprobs) != nst) stop("initprobs matrix has ", ncol(initprobs), " columns, should be number of states = ", nst)
if (est.initprobs) { warning("Not estimating initial state occupancy probabilities since supplied as a matrix") }
initprobs <- initprobs / rowSums(initprobs)
est.initprobs <- FALSE
}
else {
if (length(initprobs) != nst) stop("initprobs of length ", length(initprobs), ", should be ", nst)
initprobs <- initprobs / sum(initprobs)
if (est.initprobs && any(initprobs==1)) {
est.initprobs <- FALSE
warning("Not estimating initial state occupancy probabilities, since some are fixed to 1")
}
}
}
nipars <- if (est.initprobs) nst else 0
if (any(sapply(hmodel, inherits, "hmmMVdist"))){
hmod <- msm.form.mvhmodel(hmodel)
} else hmod <- msm.form.univhmodel(hmodel)
hmod <- c(list(hidden=TRUE, nstates=nst, fitted=FALSE, nipars=nipars,
initprobs=initprobs, est.initprobs=est.initprobs, ematrix=FALSE),
hmod)
class(hmod) <- "hmodel"
hmod
}
msm.form.univhmodel <- function(hmodel){
nst <- length(hmodel)
labels <- sapply(hmodel, function(x) x$label)
models <- match(labels, .msm.HMODELS)
pars <- lapply(hmodel, function(x) x$pars)
plabs <- lapply(hmodel, function(x) names(x$pars))
## where non-misclassified outcome for hmmIdent distribution is not specified, this is
## just the state
pars[labels=="identity"][sapply(pars[labels=="identity"], length) == 0] <- which(labels=="identity")
plabs[labels=="identity"] <- "which"
names(plabs) <- paste("state", 1:nst, sep=".")
npars <- sapply(pars, length)
parstate <- rep(1:nst, npars)
firstpar <- c(0, cumsum(npars)[-nst])
pars <- as.numeric(unlist(pars))
plabs <- unlist(plabs)
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
names(pars) <- plabs
list(models=models, labels=labels, npars=npars,
nout=rep(1, nst), mv=FALSE, parout=rep(1, sum(npars)),
totpars=sum(npars), pars=pars, plabs=plabs, parstate=parstate,
firstpar=firstpar, locpars=locpars)
}
### CHANGES FROM UNIVARIATE TO MULTIVARIATE OUTCOME HMMS
### models should be matrix instead of vector, new dim given by nout (same for labels)
### npars was vector with one for each state. now matrix
### pars was ragged array, mapped to states with parstate. should still be. (same for plabs)
### need new vector parout to map pars to outcomes.
### firstpar, used in C to point to first parameters for a univariate model. was one for each state. now needs to be a matrix
msm.form.mvhmodel <- function(hmodel){
nst <- length(hmodel)
for (i in seq_along(hmodel))
if (!inherits(hmodel[[i]], "hmmMVdist")) hmodel[[i]] <- list(hmodel[[i]])
nout <- sapply(hmodel, length)
models <- labels <- firstpar <- matrix(nrow=max(nout), ncol=nst)
npars <- matrix(0, nrow=max(nout), ncol=nst)
pars <- plabs <- parstate <- parout <- vector(nst, mode="list")
for (i in 1:nst){
labels[1:nout[i],i] <- sapply(hmodel[[i]], function(x) x$label)
models[1:nout[i],i] <- match(labels[1:nout[i],i], .msm.HMODELS)
pars[[i]] <- lapply(hmodel[[i]], function(x)x$pars)
npars[1:nout[i],i] <- sapply(pars[[i]], length)
pars[[i]][labels[,i]=="identity" & npars[,i]==0] <- i # TESTME - hmmIdent with no arg: par is the state
plabs[[i]][labels[,i]=="identity" & npars[,i]==0] <- "which"
plabs[[i]] <- lapply(hmodel[[i]], function(x)names(x$pars))
parstate[[i]] <- rep(i, sum(npars[,i]))
parout[[i]] <- rep(1:nout[i], npars[1:nout[i],i])
}
firstpar <- matrix(c(0, cumsum(npars)[-length(npars)]), nrow=max(nout), ncol=nst)
firstpar[npars==0] <- NA
pars <- unlist(pars); plabs <- unlist(plabs); parout <- unlist(parout); parstate <- unlist(parstate)
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
names(pars) <- plabs
list(models=models, labels=labels, npars=npars, nout=nout, mv=TRUE,
totpars=sum(npars), pars=pars, plabs=plabs, parstate=parstate, parout=parout,
firstpar=firstpar, locpars=locpars)
}
## NOTE removed whichcovh, whichcovh.orig
msm.form.hcmodel <- function(hmodel, mm, hcovinits, hconstraint)
{
nst <- hmodel$nstates
ncovs <- if (is.null(mm)) rep(0, nst) else sapply(mm, function(x) {ncol(x)-1})
ncovs2 <- rep(rep(0, nst), hmodel$npars)
ncovs2[hmodel$locpars] <- ncovs[hmodel$parstate[hmodel$locpars]]
coveffstate <- rep(1:nst, tapply(ncovs2, hmodel$parstate, sum))
if (is.null(hcovinits)){
coveffect <- rep(0, sum(ncovs2))
}
else {
if (!(sum(ncovs2) == length(unlist(hcovinits)))) {
warning("Initial values for hidden covariate effects do not match numbers of covariates, ignoring")
coveffect <- rep(0, sum(ncovs2))
} else coveffect <- unlist(hcovinits)
if (!is.numeric(coveffect)) {
warning("hcovinits should be numeric")
coveffect <- rep(0, sum(ncovs2))
}
}
covlabels <- lapply(mm, function(x) colnames(x)[-1])
covlabels <- unlist(covlabels[hmodel$parstate[hmodel$locpars]])
names(coveffect) <- covlabels
hcmod <- list(ncovs=ncovs2, coveffect=coveffect, covlabels=covlabels,
coveffstate=coveffstate, ncoveffs=length(coveffect))
hmodel <- c(hmodel, hcmod)
hmodel$plabs[hmodel$plabs=="hcov"] <- paste("hcov.",covlabels,sep="")
class(hmodel) <- "hmodel"
hmodel
}
### NOTE whichcovi removed
msm.form.icmodel <- function(hmodel, mm, icovinits) {
nst <- hmodel$nstates
nicovs <- ncol(mm) - 1
nicovs <- rep(nicovs, nst-1)
if (!is.matrix(hmodel$initprobs))
nicovs[hmodel$initprobs[-1] == 0] <- 0 # don't estimate cov effects on probs which are fixed to zero
if (is.null(icovinits))
icoveffect <- rep(0, sum(nicovs))
else {
icoveffect <- unlist(icovinits)
if (!(length(icoveffect) == sum(nicovs))) {
warning("Initial values for initial state covariate effects do not match numbers of covariates, ignoring")
icoveffect <- rep(0, sum(nicovs))
}
else if (!is.numeric(icoveffect)) {
warning("icovinits should be numeric")
icoveffect <- rep(0, sum(nicovs))
}
}
names(icoveffect) <- rep(colnames(mm)[-1], each=sum(nicovs>0))
icmod <- list(nicovs=nicovs, icoveffect=icoveffect, nicoveffs=length(icoveffect))
hmodel <- c(hmodel, icmod)
class(hmodel) <- "hmodel"
hmodel
}
## Convert old-style misclassification model specification to a new-style HMM with categorical response
msm.emodel2hmodel <- function(emodel, qmodel)
{
nst <- qmodel$nstates
if (emodel$misc) {
hidden <- TRUE
nepars <- rowSums(emodel$imatrix)
models <- ifelse(apply(emodel$ematrix,1,function(x)any(x==1)), 2, 1)
npars <- ifelse(models==1, 2 + nst, 1)
pars <- plabs <- vector(nst, mode="list")
parstate <- rep(1:nst, npars)
names(pars) <- names(plabs) <- paste("state", 1:nst, sep=".")
for (i in seq(nst)) {
if (models[i]==1) {
ppars <- emodel$ematrix[i,]
plab <- rep("p", nst)
plab[ppars==0] <- "p0"
plab[i] <- "pbase"
plabs[[i]] <- c("ncats", "basecat", plab)
## Baseline category is the probability of no misclassification (diagonal of ematrix)
pars[[i]] <- c(nst, i, ppars)
}
else {
plabs[[i]] <- "which"
pars[[i]] <- which(emodel$ematrix[i,]==1)
}
}
firstpar <- c(0, cumsum(npars)[-qmodel$nstates])
pars <- unlist(pars)
plabs <- unlist(plabs)
names(pars) <- plabs
labels <- .msm.HMODELS[models]
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
hmod <- list(hidden=TRUE, fitted=FALSE, nstates=nst, models=models, labels=labels,
ematrix=TRUE, ## remember if obtained from ematrix, since could change meaning of obstrue
nout=rep(1, nst), mv=FALSE,
npars=npars, totpars=sum(npars), locpars=locpars,
pars=pars, plabs=plabs, parstate=parstate, parout=rep(1,sum(npars)),
firstpar=firstpar, nipars=emodel$nipars, initprobs=emodel$initprobs, est.initprobs=emodel$est.initprobs)
hmod$constr <- msm.econstr2hconstr(emodel$constr, hmod)
}
else {
hmod <- list(hidden=FALSE, fitted=FALSE, ematrix=FALSE, models=rep(0, qmodel$nstates), npars=0, ndpars=0)
}
class(hmod) <- "hmodel"
hmod
}
msm.misccov2hcov <- function(misccovariates, emodel)
{
nst <- nrow(emodel$imatrix)
whichst <- rep(1:nst, rowSums(emodel$imatrix))
hcov <- vector(nst, mode="list")
for (i in 1:nst) {
if (!any(whichst==i)) hcov[[i]] <- ~ 1
else hcov[[i]] <- misccovariates
}
hcov
}
msm.misccovinits2hcovinits <- function(misccovinits, hcovariates, emodel, ecmodel)
{
nst <- nrow(emodel$imatrix)
whichst <- rep(1:nst, rowSums(emodel$imatrix))
hcovinits <- vector(nst, mode="list")
for (i in 1:nst) {
if (is.null(misccovinits)) hcovinits[[i]] <- rep(0, ecmodel$ncovs * rowSums(emodel$imatrix)[i])
else if (!any(whichst==i)) hcovinits[[i]] <- numeric(0)
else hcovinits[[i]] <- as.vector(t(sapply(misccovinits, function(x)x[whichst==i])))
}
hcovinits
}
msm.econstr2hconstr <- function(econstr, hmodel)
{
constr <- seq(length.out=hmodel$totpars)
for (i in unique(econstr)) {
constr[hmodel$plabs == "p"][econstr == i] <-
min(constr[hmodel$plabs == "p"][econstr == i])
}
match(constr, unique(constr))
}
#' @export
print.hmodel <- function(x, ...)
{
ci <- (x$fitted && x$foundse)
cols <- if (ci) c("Estimate","LCL","UCL") else ("")
if (!x$hidden)
cat("Non-hidden Markov model\n")
else {
cat("Hidden Markov model, ")
nst <- x$nstates
cat(nst, "states\n")
if (x$est.initprobs){
cat("Initial state occupancy probabilities: ")
if (x$nipars > 0) {
cat("\n")
print(x$initprobs)
cat("\n")
if (any(x$nicovs > 0)) {
cat("Covariates on log odds of initial states relative to state 1\n")
print(x$icoveffect)
}
cat("\n")
}
else cat(paste(x$initprobs, collapse=","), "\n\n")
}
for (i in 1:nst) {
if (x$mv){
cat("State", i, "\n")
for (j in 1:x$nout[i]){
cat("Outcome", j, "-", x$labels[j,i], "distribution\n")
if (x$labels[j,i]=="categorical")
pars <- print_hmmcat(x, i, j, mv=TRUE) ## TESTME
else {
inds <- x$parstate==i & x$parout==j
pars <- as.matrix(x$pars[inds])
if (ci) pars <- cbind(pars, matrix(x$ci[inds,], ncol=2))
dimnames(pars) <- list(x$plabs[inds], cols)
}
## covs not supported for the moment
print(pars)
cat("\n")
}
} else {
cat("State", i, "-", x$labels[i], "distribution\n")
cat("Parameters: \n")
if (x$label[i]=="categorical")
pars <- print_hmmcat(x, i)
else {
pars <- as.matrix(x$pars[x$parstate==i])
if (ci) pars <- cbind(pars, matrix(x$ci[x$parstate==i ,], ncol=2))
dimnames(pars) <- list(x$plabs[x$parstate==i], cols)
}
if (any(x$ncovs[x$parstate==i] > 0)){
coveffs <- as.matrix(x$coveffect[x$coveffstate==i])
if (ci) coveffs <- cbind(coveffs, matrix(x$covci[x$coveffstate==i,], ncol=2))
rownames(coveffs) <- x$covlabels[x$coveffstate==i]
pars <- rbind(pars, coveffs)
}
print(pars)
cat("\n")
}
}
}
}
print_hmmcat <- function(x, i, j=NULL, mv=FALSE)
{
inds <- if (mv) x$parstate==i & x$parout==j else x$parstate==i
pars <- x$pars[inds]
res <- matrix(pars[3:(2+pars[1])], ncol=1)
rownames(res) <- paste("P(",seq(length.out=pars[1]),")",sep="")
if (x$fitted && x$foundse) {
ci <- matrix(x$ci[inds,], ncol=2)
res <- cbind(res, ci[3:(2+pars[1]),])
colnames(res) <- c("Estimate","LCL","UCL")
}
else colnames(res) <- "prob"
res
}
msm.form.hconstraint <- function(constraint, hmodel)
{
constr <- seq(length=hmodel$totpars)
for (con in names(constraint)) {
if ( ! (con %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", con, "\" in hconstraint unknown")
if (con %in% hmodel$plabs) {
tc <- constraint[[con]]
np <- length(tc)
if (np != sum(hmodel$plabs==con))
stop("constraint for \"", con, "\" of length ",
np, ", should be ", sum(hmodel$plabs==con))
for (i in unique(tc))
constr[hmodel$plabs == con][tc == i] <-
min(constr[hmodel$plabs == con][tc == i])
}
}
match(constr, unique(constr))
}
msm.form.hcovconstraint <- function(constraint, hmodel)
{
constr <- seq(length.out=hmodel$ncoveffs)
for (con in names(constraint)) {
if ( ! (con %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", con, "\" in hconstraint unknown")
if (con %in% hmodel$covlabels) {
tc <- constraint[[con]]
np <- length(tc)
if (np != sum(hmodel$covlabels==con))
stop("constraint for \"", con, "\" of length ",
np, ", should be ", sum(hmodel$covlabels==con))
for (i in unique(tc))
constr[hmodel$covlabels == con][tc == i] <-
min(constr[hmodel$covlabels == con][tc == i])
}
}
match(constr, unique(constr))
}
## User-supplied range constraints on HMM parameters.
msm.form.hranges <- function(ranges, hmodel)
{
hranges <- .msm.PARRANGES[hmodel$plabs,]
if (hmodel$ncoveffs>0) {
hcovranges <- matrix(rep(c(-Inf, Inf), hmodel$ncoveffs), nrow=hmodel$ncoveffs, byrow=TRUE)
rownames(hcovranges) <- hmodel$covlabels
hranges <- rbind(hranges, hcovranges)
}
if (!is.null(ranges)) {
if (!is.list(ranges)) stop("expected \"hranges\" to be a list")
for (i in names(ranges)) {
if ( ! (i %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", i, "\" in \"hranges\" unknown")
ran.default <- hranges[rownames(hranges)==i,,drop=FALSE]
ran.user <- do.call("cbind", ranges[[i]])
for (j in seq_len(nrow(ran.user))) {
if (ran.user[j,1] < ran.default[j,1]) {
warning("User-supplied lower bound of ",ran.user[j,1]," for ", i,
" less than theoretical minimum of ", ran.default[j,1], ", ignoring")
ran.user[j,1] <- ran.default[j,1]
}
if (ran.user[j,2] > ran.default[j,2]) {
warning("User-supplied upper bound of ",ran.user[j,2]," for ", i,
" less than theoretical maximum of ", ran.default[j,2], ", ignoring")
ran.user[j,2] <- ran.default[j,2]
}
hranges[which(rownames(hranges)==i)[j],] <- ran.user[j,]
}
}
}
## ideally should be strict here if estimated, but allow inits on Inf boundary if not estimated.
for (i in seq_along(hmodel$pars)){
if (!in.range(hmodel$pars[i], hranges[i,], strict=FALSE))
stop("Initial value ", hmodel$pars[i], " of parameter \"", hmodel$plabs[i], "\" outside allowed range ",
"[", paste(hranges[i,], collapse=","), "]")
}
hranges
}
in.range <- function(x, interval, strict=FALSE) {
if (!is.numeric(interval) || length(interval)!=2) stop("interval should be a numeric vector of length 2")
if (!is.numeric(x)) stop("x should be numeric")
if (strict)
( (x > interval[1]) & (x < interval[2]) )
else
( (x >= interval[1]) & (x <= interval[2]) )
}
msm.form.initprobs <- function(hmodel, initprobs, mf){
npts <- attr(mf,"npts")
if (!is.null(hmodel$phase.states) && is.null(initprobs)) {
hmodel$initprobs <- matrix(0, nrow=npts, ncol=hmodel$nstates)
initstate <- mf$"(state)"[!duplicated(mf$"(subject)")]
hmodel$initprobs[cbind(1:npts, match(initstate, hmodel$pars))] <- 1
if (hmodel$est.initprobs)
warning("Not estimating initial state occupancy probabilities: assuming everyone starts at first phase")
hmodel$est.initprobs <- FALSE
}
if (!hmodel$est.initprobs) {
if (is.matrix(hmodel$initprobs)) {
if (nrow(hmodel$initprobs) != npts)
stop("initial state occupancy probability should have ", npts, " (number of subjects) rows if supplied as a matrix, found ",nrow(hmodel$initprobs))
}
}
hmodel
}
Try the msm package in your browser
Any scripts or data that you put into this service are public.
msm documentation built on Oct. 5, 2024, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions msm.form.initprobs in.range msm.form.hranges msm.form.hcovconstraint msm.form.hconstraint print_hmmcat print.hmodel msm.econstr2hconstr msm.misccovinits2hcovinits msm.misccov2hcov msm.emodel2hmodel msm.form.icmodel msm.form.hcmodel msm.form.mvhmodel msm.form.univhmodel msm.form.hmodel msm.check.hcovariates msm.check.hmodel print.hmmdist print.hmmMVdist
### FUNCTIONS FOR HIDDEN MARKOV MODELS IN CONTINUOUS TIME
### WITH ARBITRARY RESPONSE DISTRIBUTION
#' @export
print.hmmMVdist <- function(x, ...)
{
cat(sprintf("Multivariate hidden Markov model with %d outcomes:\n", length(x)))
for (i in x) print(i)
}
#' @export
print.hmmdist <- function(x, ...)
{
cat("Hidden Markov model", x$label, "distribution\n\n")
pnames <- if(x$label=="categorical") paste("P(",seq(x$pars[1]),")",sep="") else names(x$pars)
pars <- if(x$label=="categorical") x$pars[3:(2+x$pars[1])] else x$pars
cat("Parameters: ", paste(paste(pnames, pars, sep=" = "), collapse=", "))
cat("\n")
}
msm.check.hmodel <- function(hmodel, nstates)
{
if (!is.list(hmodel)) stop("Hidden model should be a list")
if (length(hmodel) != nstates) stop("hmodel of length ", length(hmodel), ", expected ", nstates)
for (i in hmodel) {
if (!inherits(i, "hmmdist")) stop("hmodel should be a list of HMM distribution objects")
}
}
msm.check.hcovariates <- function(hcovariates, qmodel)
{
if (!is.list(hcovariates)) stop("hcovariates should be a list")
if (length(hcovariates) != qmodel$nstates)
stop("hcovariates of length ", length(hcovariates), ", expected ", qmodel$nstates)
for (i in hcovariates) {
if (!is.null(i) )
if (!inherits(i, "formula")) stop("hcovariates should be a list of formulae or NULLs")
}
}
msm.form.hmodel <- function(hmodel, hconstraint=NULL, initprobs=NULL, est.initprobs)
{
nst <- length(hmodel)
if (is.null(initprobs))
initprobs <- if (est.initprobs) rep(1/nst, nst) else c(1, rep(0, nst-1))
else {
if (!is.numeric(initprobs)) stop("initprobs should be numeric")
if (is.matrix(initprobs)) {
if (ncol(initprobs) != nst) stop("initprobs matrix has ", ncol(initprobs), " columns, should be number of states = ", nst)
if (est.initprobs) { warning("Not estimating initial state occupancy probabilities since supplied as a matrix") }
initprobs <- initprobs / rowSums(initprobs)
est.initprobs <- FALSE
}
else {
if (length(initprobs) != nst) stop("initprobs of length ", length(initprobs), ", should be ", nst)
initprobs <- initprobs / sum(initprobs)
if (est.initprobs && any(initprobs==1)) {
est.initprobs <- FALSE
warning("Not estimating initial state occupancy probabilities, since some are fixed to 1")
}
}
}
nipars <- if (est.initprobs) nst else 0
if (any(sapply(hmodel, inherits, "hmmMVdist"))){
hmod <- msm.form.mvhmodel(hmodel)
} else hmod <- msm.form.univhmodel(hmodel)
hmod <- c(list(hidden=TRUE, nstates=nst, fitted=FALSE, nipars=nipars,
initprobs=initprobs, est.initprobs=est.initprobs, ematrix=FALSE),
hmod)
class(hmod) <- "hmodel"
hmod
}
msm.form.univhmodel <- function(hmodel){
nst <- length(hmodel)
labels <- sapply(hmodel, function(x) x$label)
models <- match(labels, .msm.HMODELS)
pars <- lapply(hmodel, function(x) x$pars)
plabs <- lapply(hmodel, function(x) names(x$pars))
## where non-misclassified outcome for hmmIdent distribution is not specified, this is
## just the state
pars[labels=="identity"][sapply(pars[labels=="identity"], length) == 0] <- which(labels=="identity")
plabs[labels=="identity"] <- "which"
names(plabs) <- paste("state", 1:nst, sep=".")
npars <- sapply(pars, length)
parstate <- rep(1:nst, npars)
firstpar <- c(0, cumsum(npars)[-nst])
pars <- as.numeric(unlist(pars))
plabs <- unlist(plabs)
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
names(pars) <- plabs
list(models=models, labels=labels, npars=npars,
nout=rep(1, nst), mv=FALSE, parout=rep(1, sum(npars)),
totpars=sum(npars), pars=pars, plabs=plabs, parstate=parstate,
firstpar=firstpar, locpars=locpars)
}
### CHANGES FROM UNIVARIATE TO MULTIVARIATE OUTCOME HMMS
### models should be matrix instead of vector, new dim given by nout (same for labels)
### npars was vector with one for each state. now matrix
### pars was ragged array, mapped to states with parstate. should still be. (same for plabs)
### need new vector parout to map pars to outcomes.
### firstpar, used in C to point to first parameters for a univariate model. was one for each state. now needs to be a matrix
msm.form.mvhmodel <- function(hmodel){
nst <- length(hmodel)
for (i in seq_along(hmodel))
if (!inherits(hmodel[[i]], "hmmMVdist")) hmodel[[i]] <- list(hmodel[[i]])
nout <- sapply(hmodel, length)
models <- labels <- firstpar <- matrix(nrow=max(nout), ncol=nst)
npars <- matrix(0, nrow=max(nout), ncol=nst)
pars <- plabs <- parstate <- parout <- vector(nst, mode="list")
for (i in 1:nst){
labels[1:nout[i],i] <- sapply(hmodel[[i]], function(x) x$label)
models[1:nout[i],i] <- match(labels[1:nout[i],i], .msm.HMODELS)
pars[[i]] <- lapply(hmodel[[i]], function(x)x$pars)
npars[1:nout[i],i] <- sapply(pars[[i]], length)
pars[[i]][labels[,i]=="identity" & npars[,i]==0] <- i # TESTME - hmmIdent with no arg: par is the state
plabs[[i]][labels[,i]=="identity" & npars[,i]==0] <- "which"
plabs[[i]] <- lapply(hmodel[[i]], function(x)names(x$pars))
parstate[[i]] <- rep(i, sum(npars[,i]))
parout[[i]] <- rep(1:nout[i], npars[1:nout[i],i])
}
firstpar <- matrix(c(0, cumsum(npars)[-length(npars)]), nrow=max(nout), ncol=nst)
firstpar[npars==0] <- NA
pars <- unlist(pars); plabs <- unlist(plabs); parout <- unlist(parout); parstate <- unlist(parstate)
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
names(pars) <- plabs
list(models=models, labels=labels, npars=npars, nout=nout, mv=TRUE,
totpars=sum(npars), pars=pars, plabs=plabs, parstate=parstate, parout=parout,
firstpar=firstpar, locpars=locpars)
}
## NOTE removed whichcovh, whichcovh.orig
msm.form.hcmodel <- function(hmodel, mm, hcovinits, hconstraint)
{
nst <- hmodel$nstates
ncovs <- if (is.null(mm)) rep(0, nst) else sapply(mm, function(x) {ncol(x)-1})
ncovs2 <- rep(rep(0, nst), hmodel$npars)
ncovs2[hmodel$locpars] <- ncovs[hmodel$parstate[hmodel$locpars]]
coveffstate <- rep(1:nst, tapply(ncovs2, hmodel$parstate, sum))
if (is.null(hcovinits)){
coveffect <- rep(0, sum(ncovs2))
}
else {
if (!(sum(ncovs2) == length(unlist(hcovinits)))) {
warning("Initial values for hidden covariate effects do not match numbers of covariates, ignoring")
coveffect <- rep(0, sum(ncovs2))
} else coveffect <- unlist(hcovinits)
if (!is.numeric(coveffect)) {
warning("hcovinits should be numeric")
coveffect <- rep(0, sum(ncovs2))
}
}
covlabels <- lapply(mm, function(x) colnames(x)[-1])
covlabels <- unlist(covlabels[hmodel$parstate[hmodel$locpars]])
names(coveffect) <- covlabels
hcmod <- list(ncovs=ncovs2, coveffect=coveffect, covlabels=covlabels,
coveffstate=coveffstate, ncoveffs=length(coveffect))
hmodel <- c(hmodel, hcmod)
hmodel$plabs[hmodel$plabs=="hcov"] <- paste("hcov.",covlabels,sep="")
class(hmodel) <- "hmodel"
hmodel
}
### NOTE whichcovi removed
msm.form.icmodel <- function(hmodel, mm, icovinits) {
nst <- hmodel$nstates
nicovs <- ncol(mm) - 1
nicovs <- rep(nicovs, nst-1)
if (!is.matrix(hmodel$initprobs))
nicovs[hmodel$initprobs[-1] == 0] <- 0 # don't estimate cov effects on probs which are fixed to zero
if (is.null(icovinits))
icoveffect <- rep(0, sum(nicovs))
else {
icoveffect <- unlist(icovinits)
if (!(length(icoveffect) == sum(nicovs))) {
warning("Initial values for initial state covariate effects do not match numbers of covariates, ignoring")
icoveffect <- rep(0, sum(nicovs))
}
else if (!is.numeric(icoveffect)) {
warning("icovinits should be numeric")
icoveffect <- rep(0, sum(nicovs))
}
}
names(icoveffect) <- rep(colnames(mm)[-1], each=sum(nicovs>0))
icmod <- list(nicovs=nicovs, icoveffect=icoveffect, nicoveffs=length(icoveffect))
hmodel <- c(hmodel, icmod)
class(hmodel) <- "hmodel"
hmodel
}
## Convert old-style misclassification model specification to a new-style HMM with categorical response
msm.emodel2hmodel <- function(emodel, qmodel)
{
nst <- qmodel$nstates
if (emodel$misc) {
hidden <- TRUE
nepars <- rowSums(emodel$imatrix)
models <- ifelse(apply(emodel$ematrix,1,function(x)any(x==1)), 2, 1)
npars <- ifelse(models==1, 2 + nst, 1)
pars <- plabs <- vector(nst, mode="list")
parstate <- rep(1:nst, npars)
names(pars) <- names(plabs) <- paste("state", 1:nst, sep=".")
for (i in seq(nst)) {
if (models[i]==1) {
ppars <- emodel$ematrix[i,]
plab <- rep("p", nst)
plab[ppars==0] <- "p0"
plab[i] <- "pbase"
plabs[[i]] <- c("ncats", "basecat", plab)
## Baseline category is the probability of no misclassification (diagonal of ematrix)
pars[[i]] <- c(nst, i, ppars)
}
else {
plabs[[i]] <- "which"
pars[[i]] <- which(emodel$ematrix[i,]==1)
}
}
firstpar <- c(0, cumsum(npars)[-qmodel$nstates])
pars <- unlist(pars)
plabs <- unlist(plabs)
names(pars) <- plabs
labels <- .msm.HMODELS[models]
locpars <- which(plabs == rep(.msm.LOCPARS[labels], npars))
hmod <- list(hidden=TRUE, fitted=FALSE, nstates=nst, models=models, labels=labels,
ematrix=TRUE, ## remember if obtained from ematrix, since could change meaning of obstrue
nout=rep(1, nst), mv=FALSE,
npars=npars, totpars=sum(npars), locpars=locpars,
pars=pars, plabs=plabs, parstate=parstate, parout=rep(1,sum(npars)),
firstpar=firstpar, nipars=emodel$nipars, initprobs=emodel$initprobs, est.initprobs=emodel$est.initprobs)
hmod$constr <- msm.econstr2hconstr(emodel$constr, hmod)
}
else {
hmod <- list(hidden=FALSE, fitted=FALSE, ematrix=FALSE, models=rep(0, qmodel$nstates), npars=0, ndpars=0)
}
class(hmod) <- "hmodel"
hmod
}
msm.misccov2hcov <- function(misccovariates, emodel)
{
nst <- nrow(emodel$imatrix)
whichst <- rep(1:nst, rowSums(emodel$imatrix))
hcov <- vector(nst, mode="list")
for (i in 1:nst) {
if (!any(whichst==i)) hcov[[i]] <- ~ 1
else hcov[[i]] <- misccovariates
}
hcov
}
msm.misccovinits2hcovinits <- function(misccovinits, hcovariates, emodel, ecmodel)
{
nst <- nrow(emodel$imatrix)
whichst <- rep(1:nst, rowSums(emodel$imatrix))
hcovinits <- vector(nst, mode="list")
for (i in 1:nst) {
if (is.null(misccovinits)) hcovinits[[i]] <- rep(0, ecmodel$ncovs * rowSums(emodel$imatrix)[i])
else if (!any(whichst==i)) hcovinits[[i]] <- numeric(0)
else hcovinits[[i]] <- as.vector(t(sapply(misccovinits, function(x)x[whichst==i])))
}
hcovinits
}
msm.econstr2hconstr <- function(econstr, hmodel)
{
constr <- seq(length.out=hmodel$totpars)
for (i in unique(econstr)) {
constr[hmodel$plabs == "p"][econstr == i] <-
min(constr[hmodel$plabs == "p"][econstr == i])
}
match(constr, unique(constr))
}
#' @export
print.hmodel <- function(x, ...)
{
ci <- (x$fitted && x$foundse)
cols <- if (ci) c("Estimate","LCL","UCL") else ("")
if (!x$hidden)
cat("Non-hidden Markov model\n")
else {
cat("Hidden Markov model, ")
nst <- x$nstates
cat(nst, "states\n")
if (x$est.initprobs){
cat("Initial state occupancy probabilities: ")
if (x$nipars > 0) {
cat("\n")
print(x$initprobs)
cat("\n")
if (any(x$nicovs > 0)) {
cat("Covariates on log odds of initial states relative to state 1\n")
print(x$icoveffect)
}
cat("\n")
}
else cat(paste(x$initprobs, collapse=","), "\n\n")
}
for (i in 1:nst) {
if (x$mv){
cat("State", i, "\n")
for (j in 1:x$nout[i]){
cat("Outcome", j, "-", x$labels[j,i], "distribution\n")
if (x$labels[j,i]=="categorical")
pars <- print_hmmcat(x, i, j, mv=TRUE) ## TESTME
else {
inds <- x$parstate==i & x$parout==j
pars <- as.matrix(x$pars[inds])
if (ci) pars <- cbind(pars, matrix(x$ci[inds,], ncol=2))
dimnames(pars) <- list(x$plabs[inds], cols)
}
## covs not supported for the moment
print(pars)
cat("\n")
}
} else {
cat("State", i, "-", x$labels[i], "distribution\n")
cat("Parameters: \n")
if (x$label[i]=="categorical")
pars <- print_hmmcat(x, i)
else {
pars <- as.matrix(x$pars[x$parstate==i])
if (ci) pars <- cbind(pars, matrix(x$ci[x$parstate==i ,], ncol=2))
dimnames(pars) <- list(x$plabs[x$parstate==i], cols)
}
if (any(x$ncovs[x$parstate==i] > 0)){
coveffs <- as.matrix(x$coveffect[x$coveffstate==i])
if (ci) coveffs <- cbind(coveffs, matrix(x$covci[x$coveffstate==i,], ncol=2))
rownames(coveffs) <- x$covlabels[x$coveffstate==i]
pars <- rbind(pars, coveffs)
}
print(pars)
cat("\n")
}
}
}
}
print_hmmcat <- function(x, i, j=NULL, mv=FALSE)
{
inds <- if (mv) x$parstate==i & x$parout==j else x$parstate==i
pars <- x$pars[inds]
res <- matrix(pars[3:(2+pars[1])], ncol=1)
rownames(res) <- paste("P(",seq(length.out=pars[1]),")",sep="")
if (x$fitted && x$foundse) {
ci <- matrix(x$ci[inds,], ncol=2)
res <- cbind(res, ci[3:(2+pars[1]),])
colnames(res) <- c("Estimate","LCL","UCL")
}
else colnames(res) <- "prob"
res
}
msm.form.hconstraint <- function(constraint, hmodel)
{
constr <- seq(length=hmodel$totpars)
for (con in names(constraint)) {
if ( ! (con %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", con, "\" in hconstraint unknown")
if (con %in% hmodel$plabs) {
tc <- constraint[[con]]
np <- length(tc)
if (np != sum(hmodel$plabs==con))
stop("constraint for \"", con, "\" of length ",
np, ", should be ", sum(hmodel$plabs==con))
for (i in unique(tc))
constr[hmodel$plabs == con][tc == i] <-
min(constr[hmodel$plabs == con][tc == i])
}
}
match(constr, unique(constr))
}
msm.form.hcovconstraint <- function(constraint, hmodel)
{
constr <- seq(length.out=hmodel$ncoveffs)
for (con in names(constraint)) {
if ( ! (con %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", con, "\" in hconstraint unknown")
if (con %in% hmodel$covlabels) {
tc <- constraint[[con]]
np <- length(tc)
if (np != sum(hmodel$covlabels==con))
stop("constraint for \"", con, "\" of length ",
np, ", should be ", sum(hmodel$covlabels==con))
for (i in unique(tc))
constr[hmodel$covlabels == con][tc == i] <-
min(constr[hmodel$covlabels == con][tc == i])
}
}
match(constr, unique(constr))
}
## User-supplied range constraints on HMM parameters.
msm.form.hranges <- function(ranges, hmodel)
{
hranges <- .msm.PARRANGES[hmodel$plabs,]
if (hmodel$ncoveffs>0) {
hcovranges <- matrix(rep(c(-Inf, Inf), hmodel$ncoveffs), nrow=hmodel$ncoveffs, byrow=TRUE)
rownames(hcovranges) <- hmodel$covlabels
hranges <- rbind(hranges, hcovranges)
}
if (!is.null(ranges)) {
if (!is.list(ranges)) stop("expected \"hranges\" to be a list")
for (i in names(ranges)) {
if ( ! (i %in% c(hmodel$plabs, hmodel$covlabels)))
stop("parameter \"", i, "\" in \"hranges\" unknown")
ran.default <- hranges[rownames(hranges)==i,,drop=FALSE]
ran.user <- do.call("cbind", ranges[[i]])
for (j in seq_len(nrow(ran.user))) {
if (ran.user[j,1] < ran.default[j,1]) {
warning("User-supplied lower bound of ",ran.user[j,1]," for ", i,
" less than theoretical minimum of ", ran.default[j,1], ", ignoring")
ran.user[j,1] <- ran.default[j,1]
}
if (ran.user[j,2] > ran.default[j,2]) {
warning("User-supplied upper bound of ",ran.user[j,2]," for ", i,
" less than theoretical maximum of ", ran.default[j,2], ", ignoring")
ran.user[j,2] <- ran.default[j,2]
}
hranges[which(rownames(hranges)==i)[j],] <- ran.user[j,]
}
}
}
## ideally should be strict here if estimated, but allow inits on Inf boundary if not estimated.
for (i in seq_along(hmodel$pars)){
if (!in.range(hmodel$pars[i], hranges[i,], strict=FALSE))
stop("Initial value ", hmodel$pars[i], " of parameter \"", hmodel$plabs[i], "\" outside allowed range ",
"[", paste(hranges[i,], collapse=","), "]")
}
hranges
}
in.range <- function(x, interval, strict=FALSE) {
if (!is.numeric(interval) || length(interval)!=2) stop("interval should be a numeric vector of length 2")
if (!is.numeric(x)) stop("x should be numeric")
if (strict)
( (x > interval[1]) & (x < interval[2]) )
else
( (x >= interval[1]) & (x <= interval[2]) )
}
msm.form.initprobs <- function(hmodel, initprobs, mf){
npts <- attr(mf,"npts")
if (!is.null(hmodel$phase.states) && is.null(initprobs)) {
hmodel$initprobs <- matrix(0, nrow=npts, ncol=hmodel$nstates)
initstate <- mf$"(state)"[!duplicated(mf$"(subject)")]
hmodel$initprobs[cbind(1:npts, match(initstate, hmodel$pars))] <- 1
if (hmodel$est.initprobs)
warning("Not estimating initial state occupancy probabilities: assuming everyone starts at first phase")
hmodel$est.initprobs <- FALSE
}
if (!hmodel$est.initprobs) {
if (is.matrix(hmodel$initprobs)) {
if (nrow(hmodel$initprobs) != npts)
stop("initial state occupancy probability should have ", npts, " (number of subjects) rows if supplied as a matrix, found ",nrow(hmodel$initprobs))
}
}
hmodel
}
Try the msm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.