2
$\begingroup$

I am using a parametric bootstrap/monte carlo hypothesis testing method to generate the null distribution of the log likelihood ratio statistic. However, I am worried I might be doing it wrong because about half of likelihood ratio statistics calculated on my simulated data are zero, meaning I never get any p-values between 1 and 0.5.

My null hypothesis is that the data comes from an i.i.d. normal distribution, and my alternate hypothesis is that the data comes from a convex combination of two multivariate normal distributions parameterized by $\lambda$ as follows: $\sigma_{1}^{2}(\lambda A_{n}+(1-\lambda)\mathbb{I}_{n})$, where $A_{n}$ is the variance-covariance matrix describing the correlation between any two observations $x_{i}$ and $x_{j}$ and $\mathbb{I}_{n}$ is the identity matrix, and $\lambda$ is restricted to $[0,1]$.

For my given data set, I estimate the variance under the null hypothesis according to $(1/(n-1))\sum_i(x_{i}-x_{mean})^{2}$. Then I generate 999 datasets each of size n from a normal distribution with the variance I just estimated. Then for each of these datasets I fit the $\lambda$ value by maximum likelihood and compare it to the maximum likelihood value of the nested hypothesis which fixes $\lambda=0$. Then having calculated the ratio on 999 simulated datasets, I can calculate the the p-value by counting how many log likelihood ratio values are greater than or equal to that on the original dataset.

The problem is, more than half of the log likelihood ratios are exactly zero. Is this the result of some sort of mistake I'm making? Or is it possible for the distribution of log likelihood ratios to have such a weird distribution? And in this case are the p-values reliable?

EDIT 1: (I'm actively editing this and will remove this warning when I finish.): Thanks all for your speculation, here is an example of $A_{n}$ derived from the spatial structure of ten samples, n=10.

        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
 [1,] 1.0000 0.9416 0.7284 0.7284 0.5834 0.5834 0.5834 0.0014 0.0000 0.0000
 [2,] 0.9416 1.0000 0.7284 0.7284 0.5834 0.5834 0.5834 0.0014 0.0000 0.0000
 [3,] 0.7284 0.7284 1.0000 0.9213 0.5834 0.5834 0.5834 0.0014 0.0000 0.0000
 [4,] 0.7284 0.7284 0.9213 1.0000 0.5834 0.5834 0.5834 0.0014 0.0000 0.0000
 [5,] 0.5834 0.5834 0.5834 0.5834 1.0000 0.9168 0.6186 0.0014 0.0000 0.0000
 [6,] 0.5834 0.5834 0.5834 0.5834 0.9168 1.0000 0.6186 0.0014 0.0000 0.0000
 [7,] 0.5834 0.5834 0.5834 0.5834 0.6186 0.6186 1.0000 0.0014 0.0000 0.0000
 [8,] 0.0014 0.0014 0.0014 0.0014 0.0014 0.0014 0.0014 1.0000 0.0000 0.0000
 [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.8664
[10,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8664 1.0000

While this "toy" example of $A_{n}$ above is given for demonstration purposes, to test the suggestion jbowman that I allow $\lambda$ to be negative, I'm going to use

A_n<-matrix(c(1,0.3,0.3,0.3,0.3,0.3,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,1,0.72,0.72,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.72,1,0.79,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.72,0.79,1,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.54,0.54,0.54,1,0.83,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.54,0.54,0.54,0.83,1,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,1,0.92,0.91,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,0.92,1,0.91,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,0.91,0.91,1,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,1,0.4,0.4,0.4,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,1,0.6,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,0.6,1,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,0.6,0.9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.43,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0,0,0,0,0,0,0,0,0,0,0,0,0,0.43,1,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,1,1,0.76,0.76,0.76,0.76,0.76,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,1,1,0.76,0.76,0.76,0.76,0.76,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,1,0.99,0.89,0.89,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.99,1,0.89,0.89,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.89,0.89,1,0.97,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.89,0.89,0.97,1,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.79,0.79,0.79,0.79,1,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,1,0.86,0.86,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,0.86,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,0.86,1,1),25,25)

where $A_{n}$ is (25 x 25) just like in jbowman's example below. Thus I will run the following code:

A_n<-matrix(c(1,0.3,0.3,0.3,0.3,0.3,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,1,0.72,0.72,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.72,1,0.79,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.72,0.79,1,0.54,0.54,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.54,0.54,0.54,1,0.83,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.3,0.54,0.54,0.54,0.83,1,0.15,0.15,0.15,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,1,0.92,0.91,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,0.92,1,0.91,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.15,0.15,0.15,0.15,0.15,0.15,0.91,0.91,1,0.08,0.08,0.08,0.08,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,1,0.4,0.4,0.4,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,1,0.6,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,0.6,1,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.4,0.6,0.9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.43,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0,0,0,0,0,0,0,0,0,0,0,0,0,0.43,1,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0.18,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,1,1,0.76,0.76,0.76,0.76,0.76,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,1,1,0.76,0.76,0.76,0.76,0.76,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,1,0.99,0.89,0.89,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.99,1,0.89,0.89,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.89,0.89,1,0.97,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.89,0.89,0.97,1,0.79,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.76,0.76,0.79,0.79,0.79,0.79,1,0.71,0.71,0.71,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,1,0.86,0.86,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,0.86,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0.18,0.18,0.71,0.71,0.71,0.71,0.71,0.71,0.71,0.86,1,1),25,25)
I_n <- diag(1,25)
func <- function(lambda, x) {
    S <- lambda*A_n + (1-lambda)*I_n
    retval <- mvtnorm::dmvnorm(x, mean=rep(0,25), sigma=S, log=TRUE)
}
mle <- rep(0,1000)
mle2<- rep(0,1000)
for (i in seq_along(mle)) {
    x <- rnorm(25)
    res <- optimize(func, lower=0, upper=1, maximum=TRUE, x=x)
    mle[i] <- res$maximum
    # I think -0.14 is the correct "true" lower limit:
    res2<- optimize(func, lower=-0.14, upper=1, maximum=TRUE, x=x)
    mle2[i]<- res2$maximum
}
summary(mle)
summary(mle2)
hist(mle)
hist(mle2)

Notice that for my type of $A_{n}$ I can replicate the results shown in this answer for my problem where $\lambda$ is constrained in $[0,1]$. However, one potentially strange thing is that when I allow $\lambda$ to take on negative values, despite the data generating process using $\lambda$ equals zero, the actual mode of the estimate is the limiting negative value (-0.14 for this specific $A_{n}$ matrix), though the mean and median are both close to zero.

summary(mle)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#0.0000408 0.0000614 0.0000763 0.1007287 0.1293151 0.9183230 
summary(mle2)
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#-0.13996 -0.12253 -0.03151  0.04491  0.12745  0.91830

enter image description here enter image description here

Improve this question
$\endgroup$
9
  • 6
    $\begingroup$ It was recently discussed/lamented in chat that investigators have decided to use Wilcoxon tests with sample sizes too small to give $p<0.05$ under any circumstances. That isn't the same as what you're doing, but, yes, some p-values can be impossible, with that serving as one example. $\endgroup$
    – Dave
    Commented Jul 5 at 16:29
  • 1
    $\begingroup$ Any non-randomized test statistic based on realizations of random variables having at most countable supports will have at most a countable number of possible values. Thus, any p-value computed from such a test statistic cannot have any possible value in the interval $[0,1].$ This is very apparent in nonparametric tests of small datasets (where the random values in question are ranks, signs, or limited to the data values). In your case, because you generate data from a continuous unbounded distribution, I am inclined to suspect an error in the calculation. $\endgroup$
    – whuber ♦
    Commented Jul 5 at 17:51
  • 1
    $\begingroup$ I suspect that this is because the true value of the parameter $\lambda$ is on the boundary of the parameter space, which, as is well-known, messes up a lot of maximum likelihood estimator properties. In your case, if $A_n$ differs substantially from the identity matrix - where "substantially" depends on the sample size - I could well imagine the MLE actually is zero in a lot of cases. Try a) finding something close to the smallest (most negative) $\lambda$ for which the covariance matrix is PD, and b) relaxing the restriction on the lower bound of $\lambda$ in the obvious way... $\endgroup$
    – jbowman
    Commented Jul 5 at 19:30
  • 1
    $\begingroup$ Then the true value is in the interior of the parameter space, and the log-likelihood ratios will likely be more-or-less centered around 0, which is what they would converge to as $n \to \infty$. $\endgroup$
    – jbowman
    Commented Jul 5 at 19:31
  • 2
    $\begingroup$ The LR test remains valid but the statistic no longer asymptotically has a chi-squared distribution. Typically it's a mixture of a $\chi^2(0)$ distribution and what you would expect. $\endgroup$
    – whuber ♦
    Commented Jul 6 at 15:07

2 Answers 2

Reset to default
4
$\begingroup$

As others have mentioned, yes there are cases where some p-values are impossible. But I don't think that is the case here.

It is also possible that there is an error in your calculations. But it is also possible that you are just not generating enough bootstrap samples to see the small p-values.

If the true p-value for a given $\lambda$ is $0.001$, then there is a $0.63$ probability of seeing 0 out of 1,000 samples be less extreme then your observed data. If you rerun your analysis using 9,999 or 99,999 bootstrap samples you will probably see some additional non-zero but very close to 0 p-values. Or you could just increase your bootstrap samples to about 3,700 and report your 0 p-values as "p < 0.001" (with 95% confidence).

Improve this answer
$\endgroup$
3
$\begingroup$

If $A_n$ differs substantially from the identity matrix - where "substantially" depends upon the sample size - the MLE of $\lambda$ will be driven to zero quite often, and your results will follow.

Conceptually, there are $n(n-1)/2$ distinct off-diagonal entries in $\lambda A_n + (1-\lambda) \mathbb{I}_n$, all of which, under the null hypothesis and data generating process, are equal to zero. Even for modest $n$, that's a lot of parameters that are equal to zero and can be estimated as such simply by setting the single parameter $\hat{\lambda} = 0$. With $n = 25$, we have 300 distinct off-diagonal entries in the covariance matrix; that's a lot of "pressure" on the MLE towards $\hat{\lambda} = 0$.

We can illustrate this effect with a simple experiment. We generate 25 standard Gaussian variates, so $\mathbb{I_n}$ is the true covariance matrix. We skip the estimation of $\sigma^2$, as it is equal to one, and by constructing $A_n$ so that the diagonal elements are all equal to one, the diagonal entries of the covariance matrix equal one regardless of $\lambda$. We then find the MLE of $\lambda$ using a standard one-dimensional optimization routine, and repeat the experiment 1,000 times:

library(mvtnorm)
A_n <- diag(1,25) + 0.01
I_n <- diag(1,25)

func <- function(lambda, x) {
    S <- lambda*A_n + (1-lambda)*I_n
    retval <- dmvnorm(x, mean=rep(0,25), sigma=S, log=TRUE)
}

mle <- rep(0,1000)
for (i in seq_along(mle)) {
    x <- rnorm(25)
    res <- optimize(func, lower=0, upper=1, maximum=TRUE, x=x)
    mle[i] <- res$maximum
}

Now for the results:

> summary(mle)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000661 0.0000661 0.0000661 0.0791413 0.0481997 0.9999581

Evidently, 0.000061 is as small a positive number as optimize can return; the default tolerance for testing convergence on my hardware is 0.00012207. As we can see, the MLE is essentially equal to zero over half the time:

> mean(mle < 0.0000662)
[1] 0.695

In fact, we can reduce the off-diagonal elements of $A_n$ by several orders of magnitude while getting the same result.

Now let's redesign $A_n$ so that most of the off-diagonal elements are zero, i.e., it's closer to the identity matrix:

A_n <- diag(1,25) 
for (i in 1:5) A_n[i,i+1] <- A_n[i+1,i] <- 0.1

which results in:

> summary(mle)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000661 0.0000661 0.2166391 0.4853882 0.9999339 0.9999339

with histogram:

Histogram of MLE of <span class=$\lambda$" />

Clearly our results are critically dependent on $A_n$, and for $A_n$ with even a small number of non-zero off-diagonal elements, the tendency will be to push $\hat{\lambda}$ to zero or, perhaps, one.

Now for some simulation-based investigation of unbiasedness. To do this, we set $\lambda = 0.4$, in the interior of the parameter space, and consider three sample sizes: $20, 100, 1000$. We run the simulation 10,000 times for the $n=20$ case, 1,000 times for the $n = 100$ and $n=1000$ cases, making lunch while the last one is executing. Our $A_n$ matrix is constructed as:

A_n <- diag(0.8,N) + 0.2 
I_n <- diag(1,N)
mu <- rep(0,N)
sigma <- 0.4*A_n + 0.6*I_n # lambda = 0.4

A summary version of the results is that in none of the three cases could a t-test reject the null hypothesis that $\lambda = 0.4$ at the 95% level of confidence. As our sample sizes increased, the distribution of $\hat{\lambda}$ slowly became unimodal:

enter image description here enter image description here enter image description here

In a limited and informal way, bias does not really seem to be an issue; the real issue is the bimodal sampling distribution (with the modes at the extremes) of $\hat{\lambda}$ even for sample sizes in the hundreds.

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ wow. this looks very promising thank you. I think your analysis is completely correct, but I'm confused at the step of "so what do I do about it?" Do I just carry on doing the normal likelihood ratio test? Or maybe you have a better way I can study this problem? The two things I would be most interested in are hypothesis testing and an unbiassed estimate of lambda. So far my approaches have been, generate the null statistic with parametric bootstrap to test significance, and to estimate the true lambda value... $\endgroup$
    – A Friendly Fish
    Commented Jul 5 at 22:35
  • $\begingroup$ ... my approach has been identify high correlation points and assume their correlation is actually one, so that any variation in these points is assumed to be due to a nugget effect, and then after measuring the nugget, conditional on that, finding the maximum likelihood estimate of $\sigma^2$ for $A_{n}$ which explains the full dataset. (this approach stems from the comments of stats.stackexchange.com/q/650189/398015 ) $\endgroup$
    – A Friendly Fish
    Commented Jul 5 at 22:38
  • $\begingroup$ 1. I don't think the normal LR test will work here, as the true value / null hypothesis is on the boundary of the parameter space, which typically messes up the asymptotic distribution of the statistic. 2. I also wouldn't assume the correlation is as high as1; what will happen if you have a collection of three points with pairwise assumed correlations = 1 whose values don't support that? 3. Does this mean you are, in effect, estimating $A_n$ (albeit in a rough sort of way) in a pre-processing step? ... I'll have to think about this a little more... and review your previous post. $\endgroup$
    – jbowman
    Commented Jul 5 at 22:43
  • 1
    $\begingroup$ WRT "pressure" - because it appeared to me that the data generating process was assuming all those values were zero, so, in that case, the true values are zero, and the MLE of $\lambda$ will consequently tend to prefer $\lambda$ values that make the estimated values close to zero, i.e., small $\lambda$s. $\endgroup$
    – jbowman
    Commented Jul 6 at 0:06
  • 1
    $\begingroup$ Correct. Simulation $\neq$ proof, of course, but it's indicative of something. $\endgroup$
    – jbowman
    Commented Jul 6 at 22:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.