How to correctly plot step length gamma distributions using the converged means and standard deviations of a fitted model summary?

[matt-needle]

I am working on recreating the default plots of the function plot.momentuHMM, so I can limit the axes and customize the look more to my needs. For the step parameter of my model, I took the mean and sd for each state from the model summary. I then calculated the shape and scale of each gamma distribution using the respective formulae:

shape = mean^2 / sd^2
scale = sd^2 / mean

When plotting the distributions, they have a similar form as those given by the default plot, but the scale of the y-axis is very different. Have I got something fundamentally wrong here, or is momentuHMM calculating the distributions differently under the hood?

library(momentuHMM)
library(ggplot2)

modelFit <- readRDS("fitHMM.rds")
colours.states <- c("#E69F00", "#56B4E9", "#009E73")

x <- seq(0, 2000, length=1000)

meanEn <- modelFit$mle$step[1,1]   # 10.87739
sdEn <- modelFit$mle$step[2,1]     # 8.775794

meanEx <- modelFit$mle$step[1,2]   # 60.38591
sdEx <- modelFit$mle$step[2,2]     # 43.38034

meanTr <- modelFit$mle$step[1,3]   # 304.8606
sdTr <- modelFit$mle$step[2,3]     # 211.0313

sh <- function(mean, sd) { return(mean^2 / sd^2) }
sc <- function(mean, sd) { return(sd^2 / mean) }

y_en <- dgamma(x, shape = sh(meanEn,sdEn), scale = sc(meanEn,sdEn))
y_ex <- dgamma(x, shape = sh(meanEx,sdEx), scale = sc(meanEx,sdEx))
y_tr <- dgamma(x, shape = sh(meanTr,sdTr), scale = sc(meanTr,sdTr))

ggplot() +
  geom_line(aes(x=x, y=y_en), colour=colours.states[[1]], size=1) +
  geom_line(aes(x=x, y=y_ex), colour=colours.states[[2]], size=1) +
  geom_line(aes(x=x, y=y_tr), colour=colours.states[[3]], size=1) +
  theme_minimal() +
  ylab("density") +
  xlab("step")

[bmcclintock]

If you check out the details in the help file (?plot.momentuHMM), you'll find the following:

"The state-dependent densities are weighted by the frequency of each state in the most probable state sequence (decoded with the function viterbi). For example, if the most probable state sequence indicates that one third of observations correspond to the first state, and two thirds to the second state, the plots of the densities in the first state are weighted by a factor 1/3, and in the second state by a factor 2/3."

So if you weight the densities by the state frequencies from the Viterbi algorithm, your plot should match that produced by plot.momentuHMM.

[matt-needle]

Many thanks. That has helped me immensely. Apologies for overlooking that in the documentation!

my updated code to calculate and plot the weighted densities

library(momentuHMM)
library(ggplot2)

modelFit <- readRDS("fitHMM.rds")

# calculate frequencies of states
v <- viterbi(modelFit)
stateFreq <- table(v) / length(v)

# default plot.momentuhmm colours
colours.states <- c("#E69F00", "#56B4E9", "#009E73")

# generate sequence for x axis of density functions
x <- seq(0, 2000, length=1000)

# get converged mean and sd for each state 
meanEn <- modelFit$mle$step[1,1]  # 10.87739
sdEn <- modelFit$mle$step[2,1]    # 8.775794

meanEx <- modelFit$mle$step[1,2]    # 60.38591
sdEx <- modelFit$mle$step[2,2]      # 43.38034

meanTr <- modelFit$mle$step[1,3]   # 304.8606 
sdTr <- modelFit$mle$step[2,3]     # 211.0313

# calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}

# get density functions of the distributions
y_en <- dgamma(x, shape=sh(meanEn,sdEn),  scale=sc(meanEn,sdEn)) * stateFreq[[1]]
y_ex <- dgamma(x, shape=sh(meanEx,sdEx),  scale=sc(meanEx,sdEx)) * stateFreq[[2]]
y_tr <- dgamma(x, shape=sh(meanTr,sdTr),  scale=sc(meanTr,sdTr)) * stateFreq[[3]]
# sum densities to get total
y_tot <- y_en + y_ex + y_tr

# combine densities in a single dataframe for more convenient plotting
df.y_en <- data.frame(dens=y_en, state="encamped", x=x)
df.y_ex <- data.frame(dens=y_ex,  state="exploratory", x=x)
df.y_tr <- data.frame(dens=y_tr,  state="transit", x=x)
df.y_tot <- data.frame(dens=y_tot,  state="total", x=x)

cmb <- rbind(df.y_en,df.y_ex,df.y_tr,df.y_tot)

# reorder factor levels so "total" appears bottom of the legend
cmb$state <- factor(cmb$state, levels=c("encamped","exploratory","transit","total"))

# plot distributions
ggplot() +
  geom_line(data=cmb,aes(x=x,y=dens,colour=state,linetype=state), size=1) +
  scale_colour_manual(values=c(colours.states,"#000000")) +
  scale_linetype_manual(values=c("solid","solid","solid","dashed")) +
  scale_y_continuous(limits=c(0,0.021)) +
  scale_x_continuous(limits=c(0,2000)) +
  theme_minimal() +
  theme(legend.title = element_blank()) +
  ylab("density") +
  xlab("step")