Maximum Likelihood Ratio with restriction

by user1329307   Last Updated October 10, 2019 02:19 AM

I'm trying to code the likelihood ratio test about two categories for observations from a multinomial experiment with twelve categories; e.g.,

$H_0: p_1 \leq p2$ versus $H_1: p_1>p_2$.

As a test of concepts, I attempted to compare the power of the likelihood ratio test and a naive use of a binomial proportions test for indepentent samples (ignoring the negative correlation between $p_1$ and $p_2$).

library(extraDistr)

mnom.ll<-function(x,theta,size,neg=FALSE){
  p<-theta
  ll<-dmnom(x=x,size=size,prob=p,log=TRUE)
  if(sum(p)!=1){
    ll<- -Inf
  }
  ifelse(!neg,ll,-ll)
}

mnom.like.unrestricted<-function(x,theta,size,neg=FALSE){
  p<-theta
  like<-dmnom(x=x,size=size,prob=p)
  if(sum(p)!=1){
    like<-0
  }
  ifelse(!neg,like,-like)
}

mnom.like.restricted<-function(x,theta,size,neg=FALSE){
  p<-theta
  like<-dmnom(x=x,size=size,prob=p)
  if(p[1]!=p[2]|sum(p)!=1){
    like<- 0 # penalty
  }
  ifelse(!neg,like,-like)
}

size=1483
binomial.test<-rep(NA,5000)
lr.test<-rep(NA,5000)
for(i in 1:5000){
  survey.dat<-c(rmnom(n=1,size=size,prob=c(0.29,0.26,0.16,0.04,0.03,0.03,0.01,0.02,0.02,0.01,0.01,0.12))) #power
  binomial.summary<-binom.test(x = survey.dat[c(1:2)], n = rep(sum(survey.dat),2),alternative = "two.sided")
  binomial.test[i]<-binomial.summary$p.value

  mle.unrestricted<-optim(par = survey.dat/size,
                          fn = mnom.like.unrestricted,
                          x=survey.dat,
                          size=size,
                          neg=TRUE)
  mle.restricted<-optim(par =  survey.dat/size,
                        fn = mnom.like.restricted,
                        x=survey.dat,
                        size=size,
                        neg=TRUE)
  mle.restricted

  ll.unrestricted<-mnom.ll(x=survey.dat,theta=mle.unrestricted$par,size = size,neg = FALSE)
  ll.restricted<-mnom.ll(x=survey.dat,theta=mle.restricted$par,size = size,neg = FALSE)
  lr.test[i]<-1-pchisq(-2*(ll.restricted - ll.unrestricted),df=1)
}
length(which(binomial.test<0.05))/5000
length(which(lr.test<0.05))/5000

The issue is that the penalty in the restricted case isn't working even though it's there and I get the same mles in both the restricted and unrestricted case, which match the estimates based on the sample. The issue is that the restricted mles have a convergence code of 0 and the unrestricted have a convergence code of 1. Is it an issue with my penalty? Is there a quick hack around this issue?



Related Questions


Updated June 03, 2015 22:08 PM

Updated March 04, 2017 19:19 PM

Updated May 18, 2018 16:19 PM

Updated August 27, 2018 09:19 AM