Covariate dependent Markov models? Plot state transition probability along gradient of covariate values fist post here, came from Stack overflow as it was suggested to me this is more appropriate for the kind of question.
So, data consists of 4 variable, id, x1 and x2, continuous variables which are correlated with y, a binary variable. 0 and 1 in the binary variable represent different states. Is it possible to use Markov chain models to calculate and plot state transition probability along the gradient of covariate values for each id and subsequently for the pooled data?
set.seed(1)
id =rep(1, 100)
x1 = rnorm(100)           
x2 = rnorm(100)
z = 1 + 2*x1 + 3*x2      
pr = 1/(1+exp(-z))        
y = rbinom(100,1,pr)
a<-data.frame(id,x1,x2, y)


set.seed(2)
id =rep(2, 100)
x1 = rnorm(100)           
x2 = rnorm(100)
z = 1 + 2*x1 + 3*x2      
pr = 1/(1+exp(-z))        
y = rbinom(100,1,pr)
b<-data.frame(id,x1,x2, y)


set.seed(3)
id =rep(3, 100)
x1 = rnorm(100)           
x2 = rnorm(100)
z = 1 + 2*x1 + 3*x2      
pr = 1/(1+exp(-z))        
y = rbinom(100,1,pr)
c<-data.frame(id,x1,x2, y)


d<-rbind(a,b,c)

The ideal output would be a figure of state transition probabilities as a function of a covariate as found on link and described in moveHMM package, however, my data doesn't have coordinates as input and I am really unsure on how to proceed.
 A: It sounds like what you need is a (non-hidden) Markov model with multinomial logistic regression on the transition probabilities. It seems that there are packages for this; see for example the package hesim, which has a vignette dedicated to "Markov models with multinomial logistic regression".
For the sake of completeness, you could do this with moveHMM, using the option "knownStates" to pass the vector of observed states to the function fitHMM. I'm showing the code for this below, building on your example. Because moveHMM expects columns for "x" and "y", I created dummy variables, but they do not affect the state classification because the states are known. Similarly, initial parameter values need to be specified, but they do not matter here.
# Load package   
library(moveHMM)

# Add dummy "x" and "y" columns    
d$x <- rnorm(nrow(d))
d$y <- rnorm(nrow(d))

# Prepare data for moveHMM
data <- prepData(d)
    
# Initial parameters
stepPar0 <- c(1, 1, 1, 1)
anglePar0 <- c(0, 0, 1, 1)

# Formula for transition probabilities
f <- ~ x1 + x2
    
# Fit HMM with knownStates option
hmm <- fitHMM(data = data, 
              nbStates = 2, 
              stepPar0 = stepPar0, 
              anglePar0 = anglePar0,
              formula = f, 
              knownStates = data$state + 1)

After fitting the model, you can plot the relationship between the transition probabilities and the covariates using plot(hmm), and the relationship between the stationary state probabilities and the covariates using plotStationary(hmm). Both functions can take the option plotCI = TRUE to obtain pointwise confidence intervals.
