DepmixS4 is one of many package to fit hidden Markov models (HMM). In this example, we will demostrate how to simulate and fit covariate-dependent transition HMM using depmixS4.
tempdat <- data.frame(time=rep(1:24,100),y=NA)
head(tempdat)## time y
## 1 1 NA
## 2 2 NA
## 3 3 NA
## 4 4 NA
## 5 5 NA
## 6 6 NACreate a depmix object
## Loading required package: nnet## Loading required package: MASS## Loading required package: Rsolnp## Initial state probabilties model
## pr1 pr2
## 0.5 0.5
##
## Transition model for state (component) 1
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 0
## cos(2 * pi * time/24) 0 0
## sin(2 * pi * time/24) 0 0
## Probalities at zero values of the covariates.
## 0.5 0.5
##
## Transition model for state (component) 2
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 0
## cos(2 * pi * time/24) 0 0
## sin(2 * pi * time/24) 0 0
## Probalities at zero values of the covariates.
## 0.5 0.5
##
##
## Response parameters
## Resp 1 : gaussian
## Re1.(Intercept) Re1.sd
## St1 0 1
## St2 0 1Set HMM parameters
tempmod2 <- setpars(tempmod, c(0.1, 0.9 # Initial prob parameters
, 0 , 2 , 0 , 2 , 0 , 2 , 0 , -2 , 0 , -2 , 0 , -2 # Sinusoidal parameters for state switch (mlogit)
, 2, 1 , 0 , 1 # emission distribution parameters
)
)
summary(tempmod2)## Initial state probabilties model
## pr1 pr2
## 0.1 0.9
##
## Transition model for state (component) 1
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 2
## cos(2 * pi * time/24) 0 2
## sin(2 * pi * time/24) 0 2
## Probalities at zero values of the covariates.
## 0.1192029 0.8807971
##
## Transition model for state (component) 2
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 -2
## cos(2 * pi * time/24) 0 -2
## sin(2 * pi * time/24) 0 -2
## Probalities at zero values of the covariates.
## 0.8807971 0.1192029
##
##
## Response parameters
## Resp 1 : gaussian
## Re1.(Intercept) Re1.sd
## St1 2 1
## St2 0 1Simulate hidden states and response data and put it back into a data frame. Strongly recommend simulating the response data using R, especially for user defined emission distrubtions.
source("simulate.R")
simdat <- simhmm(tempmod2,nsim=1)
simdf <- data.frame(y=simdat@response[[1]][[1]]@y
, states = simdat@states
, time = tempdat$time)
head(simdf)## y states time
## 1 2.1675800 1 1
## 2 -0.1524507 2 2
## 3 2.7712699 1 3
## 4 1.2789240 2 4
## 5 2.4461855 1 5
## 6 0.2599461 2 6Create a new depmix object and fit the simulated data above
source("simulate.R")
simdf <- data.frame(y=simdat@response[[1]][[1]]@y
, states = simdat@states
, time = tempdat$time)
mod <- (depmix(y~1 # response
, nstates = 2 # number of states
, transition=~cos(2*pi*time/24)+sin(2*pi*time/24) # transition structure
, family= gaussian()
, data=simdf
)
)
set.seed(1)
fitmod <- fit(mod,verbose=FALSE)## converged at iteration 111 with logLik: -3894.125summary(fitmod)## Initial state probabilties model
## pr1 pr2
## 0 1
##
## Transition model for state (component) 1
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 2.278506
## cos(2 * pi * time/24) 0 2.653767
## sin(2 * pi * time/24) 0 2.087787
## Probalities at zero values of the covariates.
## 0.09291877 0.9070812
##
## Transition model for state (component) 2
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 -1.813430
## cos(2 * pi * time/24) 0 -1.793069
## sin(2 * pi * time/24) 0 -1.835081
## Probalities at zero values of the covariates.
## 0.8597759 0.1402241
##
##
## Response parameters
## Resp 1 : gaussian
## Re1.(Intercept) Re1.sd
## St1 0.000 1.012
## St2 1.987 0.999source("Weibull_setup.R")
simdf2 <- data.frame(y=exp(simdat@response[[1]][[1]]@y)
, states = simdat@states
, time = tempdat$time)
rModels <- list()
rModels[[1]] <- list()
rModels[[1]][[1]] <- weibull(simdf2$y, pstart = c(0.5,0.5),data=simdf2)
rModels[[2]] <- list()
rModels[[2]][[1]] <- weibull(simdf2$y, pstart = c(1.5,1.5),data=simdf2)
trstart <- rep(1/2,4)
transition <- list()
transition[[1]] <- transInit(~ cos(2*pi*time/24)+ sin(2*pi*time/24), nst = 2, data=simdf2,
pstart = c(1/2,1/2,0,1/2,0,1/2))
transition[[2]] <- transInit(~ cos(2*pi*time/24)+ sin(2*pi*time/24), nst = 2,data=simdf2,
pstart = c(1/2,1/2,0,1/2,0,1/2))
inMod <- transInit(~ 1, ns = 2, pstart = rep(1/2, 2),
data = data.frame(1))
mod2 <- makeDepmix(response = rModels, transition = transition,
prior=inMod,homogeneous = FALSE)
fitmod2 <- fit(mod2, verbose = FALSE, emc=em.control(rand=FALSE, maxit=460))## converged at iteration 21 with logLik: -6337.316summary(fitmod2)## Initial state probabilties model
## St1 St2
## (Intercept) 0 -10.629
##
## Transition model for state (component) 1
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 1.440191
## cos(2 * pi * time/24) 0 1.551058
## sin(2 * pi * time/24) 0 1.548372
## Probalities at zero values of the covariates.
## 0.1915158 0.8084842
##
## Transition model for state (component) 2
## Model of type multinomial (mlogit), formula: ~cos(2 * pi * time/24) + sin(2 * pi * time/24)
## Coefficients:
## St1 St2
## (Intercept) 0 -2.755347
## cos(2 * pi * time/24) 0 -3.044262
## sin(2 * pi * time/24) 0 -2.218542
## Probalities at zero values of the covariates.
## 0.9402146 0.05978538
##
##
## Response parameters
## Re1.shape Re1.scale
## St1 0.984 11.488
## St2 1.165 1.600