Nuisance Parameters

Consider a classic problem:

Example

We have data \(x_1,..,x_n\) from a normal distribution mean \(\mu\) and standard deviation \(\sigma\), and we want to find a \(90\%\) confidence interval for \(\mu\).

So our model has two unknown parameters, but our interest is in just one of them. In such a case the other parameter(s) are called nuisance parameters. How can we handle them?

Frequentist Solution

There is an extension of the likelihood ratio test that says that (asymptotically) \(-2\log \Lambda(X)\) still has a chi-square distribution, but now with a degree of freedom that is the difference between the total number of parameters and the number of nuisance parameters.

Example

Consider this data set:

cat(dta)

1.25 4.65 4.69 5.16 5.33 6.83 6.94 7.6 8.25 8.28 8.42 8.44 8.45 8.53 8.69 8.8 8.85 8.91 8.96 8.99 9.08 9.1 9.13 9.4 9.44 9.44 9.48 9.64 9.72 9.79 9.97 10.03 10.09 10.11 10.18 10.37 10.4 10.61 11.16 11.53 11.64 11.78 11.8 12.01 12.46 13.03 14.4 15.81 16.22 16.79

say we want to test \(H_0:\mu=11\) vs \(H_a: \mu \ne 11\). Now

ll <- function(par) -sum(log(dnorm(dta, par[1], par[2])))
mle <- optim(c(mean(dta), sd(dta)), ll)$par               
round(mle, 3)

## [1] 9.613 2.808

ll.H0 <- function(par) -sum(log(dnorm(dta, 11, par)))
mle.H0 <- optim(sd(dta), ll.H0)$par               
round(mle.H0, 3)

## [1] 3.132

lrt <- 2*(ll.H0(mle.H0)-ll(mle))
round(c(lrt, qchisq(1-0.05, 2-1), 2*(1-pchisq(lrt, 2-1))), 4)

## [1] 10.9232  3.8415  0.0019

This is actually (almost) the same as the standard test:

t.test(dta, mu=11)

## 
##  One Sample t-test
## 
## data:  dta
## t = -3.4589, df = 49, p-value = 0.001132
## alternative hypothesis: true mean is not equal to 11
## 95 percent confidence interval:
##   8.806539 10.418661
## sample estimates:
## mean of x 
##    9.6126

If we wanted to find a confidence interval we could now invert this test. This would lead to a method known as profile likelihood.

Bayesian Solution

The treatment of nuisance parameters in a Bayesian analysis is completely straight-forward:

• find the complete posterior distribution
  
• find the marginal distribution of the parameters of interest

Example

For the data above find a 95% credible interval for \(\mu\) using the prior \(\pi(\mu, \sigma) = I_{[0, 20]}(\mu)I_{1, 5}(\sigma)\)

like <- function(mu, sig) 
  prod(dnorm(dta, mu, sig))
joint <- function(sig, mu) 
  like(mu, sig)*
  ifelse(mu<0 | mu>20, 0, 1)*
  ifelse(sig<1 | sig>5, 0, 1)
mx <- double.integral(joint, c(1, 0), c(5, 20))
post.mu <- function(mu) {
  y <- 0*mu
  for(i in seq_along(mu))
    y[i] <- integrate(joint, 1, 5, mu=mu[i])$value/mx
  y
}
ggcurve(fun=post.mu, A=8, B=11)

To find the credible interval we need the values of mu where the posterior distribution function has values 0.025 and 0.975:

mu <- 8
repeat {
  mu <- mu+0.01
  if(integrate(post.mu, 0, mu)$value>0.025) break
}
low <- mu
repeat {
  mu <- mu+0.01
  if(integrate(post.mu, 0, mu)$value>0.975) break
}
c(low, mu)

## [1] 8.91 9.97