Solved – Markov chain model likelihood ratio test

degrees of freedomlikelihood-ratiomarkov-process

Suppose I am using two Markov Chain Models, one with order $k=1$ and a second one with order $k=2$. I am "reducing" the higher order model to a $k=1$ model in order to have easier calculation possibilities.

I train each model on the same data and also calculate the log likelihoods on the same data. Now I want to determine the log likelihood ratio test in order to make a model selection, as they are nested.

To do so I need the LRT (which is straight forward) and the degrees of freedom. Currently, I am determining the df by calculating the difference between the parameters of the $k=2$ ($m^2(m-1)$) and the null model $k=1$ ($m(m-1)$).

The problem now is that the degrees of freedom are very, very high and so I come up with a high p value all the time, which says that I should stick with my null model. I am unsure, if this is the right way to do so. The second order model is much sparser, so do I really need to calculate the worst case number of parameters, or can I make any limitations to that?

Maybe someone can help me out with that.
Cheers!

Best Answer

I see no problem with your resolution: if the number of states is m, there are $m\times(m-1)$ free parameters when $k=1$ and $m\times m\times(m-1)$ free parameters when $k=2$. Unless your data is strongly dependent upon the two past states, the likelihood ratio test will favour $k=1$.

If you want to reduce the number of parameters for $k=2$, you have to do it "by hand", i.e. by introducing restrictions on those $m\times m\times(m-1)$ free parameters... Or use a variable length Markov chain.

Related Solutions

Solved – Power calculation for likelihood ratio test

You can do this using simulation.

Write a function that does your test and accepts the lambdas and sample size(s) as arguments (you have a good start above).

Now for a given set of lambdas and sample size(s) run the function a bunch of times (the replicate function in R is great for that). Then the power is just the proportion of times that you reject the null hypothesis, you can use the mean function to compute the proportion and prop.test to give a confidence interval on the power.

Here is some example code:

tmpfunc1 <- function(l1, l2=l1, n1=10, n2=n1) {
    x1 <- rpois(n1, l1)
    x2 <- rpois(n2, l2)
    m1 <- mean(x1)
    m2 <- mean(x2)
    m <- mean( c(x1,x2) )

    ll <- sum( dpois(x1, m1, log=TRUE) ) + sum( dpois(x2, m2, log=TRUE) ) - 
            sum( dpois(x1, m, log=TRUE) ) - sum( dpois(x2, m, log=TRUE) )
    pchisq(2*ll, 1, lower=FALSE)
}

# verify under null n=10

out1 <- replicate(10000, tmpfunc1(3))
mean(out1 <= 0.05)
hist(out1)
prop.test( sum(out1<=0.05), 10000 )$conf.int

# power for l1=3, l2=3.5, n1=n2=10
out2 <- replicate(10000, tmpfunc1(3,3.5))
mean(out2 <= 0.05)
hist(out2)

# power for l1=3, l2=3.5, n1=n2=50
out3 <- replicate(10000, tmpfunc1(3,3.5,n1=50))
mean(out3 <= 0.05)
hist(out3)

My results (your will differ with a different seed, but should be similar) showed a type I error rate (alpha) of 0.0496 (95% CI 0.0455-0.0541) which is close to 0.05, more precision can be obtained by increasing the 10000 in the replicate command. The powers I computed were: 9.86% and 28.6%. The histograms are not strictly necessary, but I like seeing the patterns.

Solved – Likelihood ratio test in R

In addition to @Henry's answer and to @PeterEllis' comment, here are a few comments about the R code itself:

The third argument of mlog1 is sdev. You therefore need sdev in out1 and out2 (not sd);
The likelihood ratio test is the logarithm of the ratio between two likelihoods (up to a multiplicative factor). Equivalently, it is a difference between two log-likelihoods (up to a multiplicative factor). It is not the ratio between two log-likelihoods.
k1=out1$estimate returns an estimate of $\mu$. k2=out2$estimate returns another estimate of $\mu$ (based on another starting value). Both aim to estimate the same thing; that's why you always get something very close to $1$. Neither of them returns minus a log-likelihood.
You forgot the multiplicative factor ($-2$).

Related Question