Why does this test do not follow a chi-square distribution?

Question

I am trying to run some simulations to implement a test that is supposed (according to the theory) to follow a chi-squared distribution with 1 degree of freedom. My setting is as follows:

The data $Z_i$, $i=1,\cdots,n$ comes from the normal distribution $\mathcal{N}(\alpha_i,\theta_0)$ so that $Z_i$ has mean $\alpha_i$ and all the observations share the same standard deviation $\theta_0$. Each $\alpha_i$ is unknown with density $g$. Therefore, from the researcher's point of view, the density is $$ f(z) = \int \mathcal{N}(\alpha_i,\theta_0)g(\alpha_i) d\alpha_i $$ Assuming that $\alpha_i$ is discrete with points $r=1,\cdots,R$, I have: $$ f(z) = \sum^R_r \mathcal{N}(\alpha_r,\theta_0)g(\alpha_r) $$ Since the score of this model $S_\theta$ with respect to $\theta$ (and evaluated at $\theta_0$) should have mean zero, the CLT tells us that $\sqrt{n}\mathbb{E}_n[S_\theta(Z)]$ will have a normal distribution with variance $W = \mathbb{E}[S^2_\theta]$. In this case, the score will be

$$ S_\theta(z) = \frac{1}{f(z)}\sum^R_rp_\theta(z,\alpha_r\theta_0)g(\alpha_r), $$ where $p_\theta(z,\alpha_r\theta_0)$ is the derivative of $\mathcal{N}(\alpha_r,\theta)$ with respect to $\theta$ and evaluated at $\theta_0$.

This implies that $$ t = n*\mathbb{E}_n[S_\theta(Z)]*\mathbb{E}_n[S_\theta(Z)]/\hat{W}, $$ where $\mathbb{E}_n[S_\theta(Z)] = \frac{1}{n}\sum^n_{i=1}S_\theta(z_i)$ and $\hat{W}$ is the sample counterpart of $W$. is a chi-squared distribution with 1 df.

I want to show this result in R. I estimate $g$ using the R-package deconvolveR. Then I compute the score, compute its mean, and its variance. However, when I compute the 95% and 90% quantiles based on 1000 Monte Carlo simulations, I am not getting the right quantiles of the chi-squared distribution. This is my code:

theta_null <- 1   # Theta under the H0
n = 1500          # Number of observations in the sample
B = 1000          # Number of Monte Carlo repetitions
J = 11            # Number of basis. 
c = 3            # Regularization parameter 

test_statistic = rep(NA,B)
for(b in 1:B){
set.seed(238923 + b) 

##############################Create the data###################################

theta_true <- theta_null  #True Theta (standard deviation)


#Create the alphas 

alphas <- c(runif(n = 500, -1.7, -.7), runif(n = 1000, .7, 2.7))


#Create the Z 

z <- rnorm(n, mean = alphas, sd = theta_true)


#Create vector of evaluation points 

tau <- seq(from = -4, to = 5, by = 0.2)

####Compute the  test statistic #####

estimate_nc <- deconv_my(tau = tau, X = z, family = "Normal", pDegree = J, c0 = c, theta_null=theta_null)  
g_est  <- estimate_nc$stats[, "g"]


prob_nc = matrix(NA, length(z), length(tau))

for(i in 1:length(z)){
 for (m in 1:length(tau)){
   prob_nc[i,m] = (1/(theta_null*(sqrt(2*pi))))*exp((-1/2)*(((z[i]-tau[m])/theta_null)^2))
 }
}

f_hat_nc = prob_nc%*%g_est


num_score_nc = matrix(NA, length(z), length(tau))
for(i in 1:length(z)){
 for (m in 1:length(tau)){
   num_score_nc[i,m] = (1/(sqrt(2*pi)*(theta_null^4)))*exp((-1/2)*(((z[i]-tau[m])/theta_null)^2))*((z[i]-tau[m])^2 - (theta_null^2))
 }
}

num_score_avg_nc = num_score_nc%*%g_est


score_nc = num_score_avg_nc/f_hat_nc

m_nc = score_nc

m_nc_mean = mean(m_nc)

W_nc = mean(m_nc^2)

#The test statistic is 

c_no_robust = (n*m_nc_mean*m_nc_mean)/W_nc
test_statistic[b] = c_no_robust }
size_test_5 = sum(as.numeric(test_statistic >= qchisq(.95, df=1)))/nrow(test_statistic)
size_test_10 = sum(as.numeric(test_statistic >= qchisq(.90, df=1)))/nrow(test_statistic)
size_test_5
0.080 
size_test_10
0.22

Note the deconvolveR requires some evaluation points to conduct the estimation, this is $\tau$. As you can see I am not getting 5% and 10%. This is not improving once I increase the sample size

Do you think there is something wrong with my code?

Answer 1

You have made a couple of errors in forming the score test statistic.

First, the numerator of the score test statistic involves $S_\theta$ and not $E(S_\theta)$. The score statistic $t$ as you have defined it is not even a random variable, it is just a constant.

Second, you have failed to account for the fact that $g()$ and the $\alpha_r$ need to be estimated from the data. In other words, your model contains a large number of "nuisance parameters" in addition to the parameter of interest $\theta$ that you are conducting a test about. In the presence of nuisance parameters, the denominator of the score statistic is the information for $\theta$ adjusted for all the nuisance parameters that need to be estimated, but you have not done that adjustment.

Your model is even more complex that would normally be the case because $g()$ is being estimated semi-parametrically, which makes the effective number of nuisance parameters unclear. There is no existing theory for how to adjust the score statistic for such models.