An Alternative Interpretation of the \( p \)-value

The likelihood principle is often used as the basis for criticizing the use of a \( p \)-value in classical statistics. Reading the recent discussions in the March issue of Ecology, I had an alternative explanation of what a \( p \)-value is.

The \( p \)-value is defined as the conditional probability of observing something as extreme or more than the data if the “null'' hypothesis is true. In a one sample \( t \)-test problem, we are interested in testing whether the population mean is equal to a specific value \[ H_0: \mu=\mu_0. \]

The test is based on the central limit theorem, which states that the sampling distribution of a sample mean \( \bar{x} \) is \( N(\mu,\sigma^2/n) \), where \( \mu, \sigma \) are the population mean and standard deviation and \( n \) is the sample size. In a hypothesis testing, we compare the sample mean \( \bar{x} \) to the sampling distribution under the null hypothesis \( N(\mu_0,\sigma^2/n) \).

For simplicity, we assume that \( \sigma \) is known and use a one-sided \( p \)-value. The typical definition of the \( p \)-value is the shaded area in the figure, the probability of observing sample means as extreme as and more extreme than \( \bar{x} \), hence the criticisms of violating the likelihood principle because data not observed must be used to calculate the \( p \)-value.

The \( p \)-value can be interpreted as a probability, but it can also be interpreted as an indicator of the likelihood, the density value of \( \bar{x} \). That is, the tail area is a monotonic function of the likelihood. A small \( p \)-value is related to a small likelihood value and vice versa. A density value cannot be easily understood, while a probability is scaled and easy to understand. Alternatively, we can also measure the evidence against the null hypothesis using the distance between \( \bar{x} \) and \( \mu_0 \): \( d=|\bar{x} - \mu_0| \), and \( d \) is also a monotonic function of the likelihood. Because \( \bar{x} = \pm d \) share the same likelihood, the probabilistic interpretation of \( p \)-value must be 2 times the shaded area. As long as we know the rule of translating a likelihood value to a \( p \)-value, whether it is 1 or 2 times the shaded area is irrelevant.

Using a \( p \)-value, we now interpret the evidence'' against the hypothesis in terms of a probability and a binary decision rule is now easily justifiable (no matter how arbitrary). I would argue now that the \( p \)-value itself does not violate the likelihood principle. It can be seen as an easy-to-understand indicator of the likelihood (of the observed sample mean being a random sample of the sampling distribution defined by the null hypothesis). It is the literal interpretation of the indicator that introduces confusion.

Many years ago, when I took my first statistics course, we must find the \( p \)-value in the standard normal distribution table in the back of the text. That is, we calculate \( z=\frac{\bar{x}-\mu_0}{\sigma/\sqrt{n}} \) and find a \( p \)-value in a table. Could it (using \( p \)-value) be simply a computational short cut for the likelihood? After all, a likelihood is a monotonic function of the \( p \)-value.

Using Bayesian to Cheat

A few days before the end of last semester, I stumbled upon a paper in the journal Ecological Indicators [Song and Guan(2013)]. The title indicates a Bayesian estimation method was used. The authors started with a discussion of the “environmental efficiency analysis” and introduced an indicator, which is some kind of mathematical function of input and output. There was no intuitive explanation on the indicator. The indicator was based on three input variables (population, “fixed capital formation,” and non renewable energy consumption), one desirable output variable (GDP), and one undesirable output (industrial SO2 emission). The calculation resulted in 17 environmental efficiency scores (9 cities in 2 years). The main objective of the paper is to explore factors affecting these environmental efficiency scores (EE), using regression. The potential factors are (1) per capita GDP (RGDP), (2) total import and export volume (IE), (3) the proportion of the “second industry” in GDP (GY), (4) the proportion of “the industry of the second industry” in GDP (GGY), and (5) the proportion of environmental spending in GDP (HZ). The authors explained that (1) is an indicator of economic scale, (2) is a measure of economic exchange with the outside world (the authors used the term “opening up,” a bad translation of a Chinese term), (3) and (4) are measures of industry structure, and (5) is the “government factor.” Not being an economist, I don’t want to comment on the choice of these potential factors, except that GDP is now used both as part of the environmental efficiency score and as a potential factor that will be used to explain the variance of the score.
The authors used a multiple regression approach, but regression coefficients were estimated using MCMC. I was expecting a discussion on the choice of prior distributions of these coefficients. But it was soon clear that there was no prior distribution. So, why did they use MCMC for a multiple regression problem? Based on the authors’ affiliation (School of Statistics and Mathematics), I assume that they know that there should be no substantive difference between using MCMC and using OLS. The authors presented the following regression model:



The estimated model coefficients may have revealed the answer:

coefficient     estimate standard error

      α            -0.1188 1.7260

      β1          -0.1564 0.1562

      β2            3.8130 1.5510

      β3           4.8050     8.0620

      β4          -3.3960     5.5560

      β5           0.9362 18.8900

All slopes, except β2 are statistically not different from 0! If the authors used the OLS, a typical regression model output would include a column of p-values, which will make the paper unpublishable. Using MCMC, the authors are able to present the estimated coefficients, standard error, and selected quantiles. Without the column of p-values, a busy reviewer may not be able to catch the problem. (But all my students in an introductory biostatistics class recognized the problem.)

Is this a successful story of cheating by using “Bayesian” statistics?

References
[Song and Guan(2013)] Malin Song and Youyi Guan. The environmental efficiency of Wanjian demonstration area: a Bayesian estimation approach. Ecological Indicators, 36:59–67, 2013. 

The line between Bayesian and classical statistics is very thin

A while ago, we wrote a paper on correcting re-transformation bias in regression analysis using Bayesian regression [Stow et al., 2006], where we used a Bayesian regression model implemented in WinBUGS and obtained the predictive posterior distribution of the response at given predictor values. The re-transformation bias was avoided by performing the re-transformation on MCMC samples from the posterior predictive distribution before calculating the predictive means in the original scale. This approach obviously works. But the added burden of using WinBUGS is a real obstacle for many. In Gelman and Hill [2007], the same Bayesian posterior sampling was done by using the sampling distributions of regression model parameters. For a simple linear regression problem of y = β₀ + β₁x + ε, where ε ~ N(0,σ²), the classical statistics provides sampling distributions of σ² and regression coefficients β₀ and β₁. The sampling distribution of σ² is a scaled inverse chi square distribution:

where n is the sample size, p is the number of coefficients, andis the estimated residual standard deviation. Conditional on σ², the joint distribution of (β ₀,β₁) is a multivariate normal:

where V _β = (X^T X)^-1 is the unscaled covariance matrix and (₀, ₁) are the estimated model coefficients. These results can be used to derive a Monte Carlo simulation algorithm for generating samples of model coefficients form their joint posterior distribution.

1. Extract necessary information from the fitted linear model object lm.obj:

summ <- summary (lm.obj)

coef <- summ$coef[,1:2]

sigma.hat <- summ$sigma
beta.hat <- coef[,1]
V.beta <- summ$cov.unscaled
n <- summ$df[1]+summ$df[2]
p <- summ$df[1]

2. Draw n.sim posterior samples:       

out<-matrix(0, nrow=n.sim, ncol=p+1)
for (i in 1:n.sim){
  chi2 <- rchisq(1, n-p)
  sigma <- sigma.hat * sqrt((n-p)/chi2)
  beta <- mvrnorm(1, beta.hat, V.beta*sigma^2)
  out[i,] <- c(beta, sigma)
}

3. Put the above two components into a function (e.g., lm.sims(lm.obj, n.sim) we can draw random samples from the joint posterior distribution of model coefficients, which can be used for drawing posterior predictive samples:

       coef.post <- lm.sims(lm.obj=my.lm, n.sim=1000)

       y.tilde <- rnorm(n.sim, coef.post[,1]+coef.post[,2]*x.tilde, 

                                  coef.post[,3])

where x.tilde is the predictor value at which the predictive distribution y.tilde is evaluated.

In other words, we simplified the Bayesian operation into two lines of R code in step 3. In Section 9.2.2.1 of Qian [2010], I compared the simple Monte Carlo simulation approach to the Bayesian regression (using MCMC).  

Drawing random samples of σ in step 2 is definitely a Bayesian step, as the classical sampling distribution is defined for .  In classical statistics, is a random variable, and in Bayesian statistics, it is fixed.  In any case, summarizes information about σ in the data.  The connection between a Bayesian approach and the corresponding classical approach is the shared likelihood function. If we use “non-informative” flat priors, the Bayesian posterior distribution is proportional to the likelihood function. Consequently, if we can find quick solutions from the classical repertoire of models, we should use them.  The line that separates a Bayesian from a classical approach is the definition of probability. The interpretation of probability can be subjective (e.g., the chance of a candidate wins the election) or objective (e.g., the chance of a head in a coin toss). In most scientific endeavors, we are more likely to use the former interpretation (e.g., the likelihood that a model is true). 

Whether we call ourselves Bayesian or not should not be based on whether we use MCMC (WinBUGS, OpenBUGS, JAGS, Stan), rather based on how we interpret the results. For most standard statistical modeling problems, we should use standard classical model results. MCMC gives us a tool for complicated computations. The use of MCMC is not equivalent to being Bayesian.

References

    A. Gelman and J. Hill. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, New York, 2007.

    S.S. Qian. Environmental and Ecological Statistics with R. Chapman and Hall/CRC Press, 2010.

    C.A. Stow, K.H. Reckhow, and S.S. Qian. A Bayesian approach to retransformation bias in transformed regression. Ecology, 87 (6): 1472–1477, 2006.

Environmental and Ecological Statistics with R Webpage

The webpage for Environmental and Ecological Statistics with R was removed by Duke.  I was not told about the move (nor was I notified when they disabled my Duke email). I am in the process of moving the page to Toledo.  In the meanwhile, I have uploaded the data folder and two script files here.

To Binomial or To Poisson? That is the Question

Modeling pathogenic protozoa such as Cryptosporidium oocysts in water is a problem of modeling count data. An important aspect of the problem is to understand the method detection probability. The currently used method for detecting and quantifying Cryptosporidium oocysts in a water sample has an average recovery rate of about 40%. That is, the probability of detecting an oocyst in a water sample is about 0.4 when the oocyst is present. The method used for quantifying this recovery rate is to have labs analyze water samples spiked with known number of oocysts. In the recent US EPA study, 50 certified labs were given several spiked samples and their recoveries were calculated as the number of oocysts identified divided by the number of oocysts spiked. The number of oocysts spiked is typically about 100. The counting process has a standard deviation of 2 to 3 oocysts. Given the relatively accurate total counts, we generally use a binomial model:

where the index ij represents the ith observation from lab j. A simple multilevel model is then to model the detection probability as lab specific:

and imposing a common prior on the logits of lab means:

This model assumes that detection probability varies by lab and lab-specific detection probabilities are exchangeable.

There are two problems with the model. The first problem is a numerical one. The likelihood of a binomial is of the form ∏ p^y(1 - p)^n-y. The likelihood function can be numerically challenging to evaluate under certain situations. This is, obviously, a problem that can be addressed in programming. The other problem is that the total (N_ij) is not observed exactly. When using binomial, this uncertainty is ignored. Because ignored uncertainty will not vanish, it will be reflected in the inflated uncertainty level in the estimated model parameters. In this case, the only parameter we have is the binomial mean p.

Section 1.2.5 of Agresti [2002] discusses the relationship between multinomial and Poisson. Simplifying the discussion to a binomial situation, we have the following.

1.

We use Y ₁ to denote the vector of detected number of oocysts and Y ₂ the number of missed oocysts. The vector of total number of oocysts is N = Y ₁ + Y ₂.

2.

We model Y ₁ and Y ₂ as two independent Poisson random variates with means λ₁ and λ₂, respectively.

3.

If Y ₁ and Y ₂ are Poisson random variables, the sum of Y ₁ and Y ₂ is also a Poisson random variable with mean λ₁ + λ₂.

4.

In the binomial model the sum of Y ₁ and Y ₂ is given, the joint distribution of Y ₁ and Y ₂ must then be conditional on the sum:

which is shown [Agresti, 2002] to be a binomial distribution characterized by π₁ = λ₁ ∕ (λ₁ + λ₂) and π₂ = λ₂ ∕ (λ₁ + λ₂) = 1 - π₁.

In other words, the binomial distribution is the conditional distribution of two independent Poisson variates. Or, we can model the binomial count data as two independent Poisson random variables and derive the binomial parameter.

Since we know the total (N), we can parameterize N into the Poisson mean. For example, we can parameterize the Poisson mean of the detected number of oocysts as a product of the spiked total (N) and detection probability (p₁):

Likewise, the model for λ₂ is Np₂. N is used as an offset in both models. Given that p₁ + p₂ = 1, the estimated binomial mean is

Only one Poisson model is needed. Modeling using Poisson with the total counts as the offset is the same as modeling using binomial.

By using a Poisson model, we can add a multiplicative error term (e.g., LN(0,σ²)) to account for the uncertainty in the number of spiked oocysts, and/or other sources of over-dispersion.

References

   A. Agresti. Categorical Data Analysis. Wiley, 2002.

Stan versus JAGS

Stan was born on August 30, 2012.  It was built to speed up multilevel generalized linear models.  Instead of improving the sampler used in BUGS, Stan used a better sampler -- the Hamiltonian Monte Carlo.  The HMC converges faster.  I tried Stan in the last two weeks and did a small comparison using a change point model.  The model was discussed in Environmental and Ecological Statistics with R (section 6.1.2).  In the book, lilacs first bloom dates from 4 sites were used to fit 4 hockey stick models, one for each site.  I pooled data from the 4 sites and formed a simple multilevel model, where the 4 model parameters were assumed exchangeable.  For example:

[; \beta_{0j} \sim N(\mu_0, \sigma_0) ;]

Implementing the model in BUGS (or JAGS) is fairly straight forward:


## multilevel hockey stick model
## 
model{
for (i in 1:n){
   y[i] ~ dnorm(y.hat[i], tau.y[stid[i]])
   y.hat[i] <- beta0[stid[i]]+(beta1[stid[i]]-delta[stid[i]]*step(x[i]-phi[stid[i]]))*(x[i]-phi[stid[i]])
}
for (j in 1:n.st){
   beta0[j] ~ dnorm(mu0, tau[1])
            beta1[j] ~ dnorm(mu1, tau[2])
   delta[j] ~ dnorm(mu2, tau[3])
       phi[j] ~ dnorm(mu3,tau[4])T(L[j],U[j])
   tau.y[j] <- pow(sigma.y[j], -2)
     sigma.y[j] ~ dunif(0, 10)
        } 
mu0 ~ dnorm(0, 0.0001)
mu1 ~ dnorm(0, 0.0001)
mu2 ~ dnorm(0, 0.0001)
mu3 ~ dnorm(0, 0.0001)

for (k in 1:4){
   tau[k] <- pow(sigma[k], -2)    
   sigma[k] ~ dunif(0, 100)
}
}

But the model converges very slowly. After a warm-up of 100,000 iterations, the next 100,000 iterations (thin=150) are still jumpy. 

The Stan code looks like this:

 /*  multilevel hockey stick model without group predictors */ 

data{

  int n; 

  int nst;

  real y[n]; 

  int stid[n]; // station ID

  real x[n]; // standardized year

  real minYr;

  real maxYr;

}

parameters{

  real beta0[nst]; // intercept at change point

  real beta1[nst];

  real delta[nst];

  real phi[nst];

  real mu0;

  real mu1;

  real mu2;

  real mu3;

  real sigmay;

  real sigma[4];

}

model{

    mu0 ~ normal(0, 100);

    mu1 ~ normal(0, 100);

        mu2 ~ normal(0, 100);

        mu3 ~ normal(0, 100);

        beta0 ~ normal(mu0, sigma[1]);// vector

        beta1 ~ normal(mu1, sigma[2]);// vector

        phi ~ normal(mu2,sigma[3]);   // vector

        delta ~ normal(mu3,sigma[4]); // vector

for (i in 1:n){

   y[i]~normal(beta0[stid[i]]+(beta1[stid[i]] - delta[stid[i]] * int_step(x[i]-phi[stid[i]]))*(x[i]-phi[stid[i]]), sigmay);

}

}

Two big differences between Stan and BUGS/JAGS: one is that we need to declare data and parameters, and the second is that the order of the code matters.  The first part is easy to adjust as we need to supply data and I always provide initials for all parameters.  The second is somewhat hard for me as I always write BUGS code (no ordering is needed) in the order of "likelihood -- mean function -- prior" and I must reverse the order in Stan.  But this is a small price to pay for a big reward in terms of speed of convergence and the possibility of doing something more.  

The Stan model output is somewhat different from the hockey stick model output.  For example, the estimated change point is 1975(3.5), 1983(4.9), 1974(8), and 1976(6.7) in Table 6.1 of the book, and the multilevel model estimates are 1976, 1976, 1974, and 1978 (with a more or less consistent sd of 6.6).  The multilevel estimated slopes of the first line segment ([;\beta_1;]) are much closer to 0 than the same from the hockey stick model.  The R package rstan is also very easy to use:

bugs.in <- function(infile=USlilac, n.chains=3){

  y <- infile$FirstBloom

  temp <- !is.na(y) ## Stan cannot take NA in y

  infile <- infile[temp,]

  y <- y[temp]

  xx <- infile$Year

  x <- (xx-mean(xx))/sd(xx)

  stid <- as.numeric(ordered(infile$STID))

  n <- length(y)

  n.st <- max(stid)

  L <- tapply(x, stid, min)

  U <- tapply(x, stid, max)

  inits <- list()

  bugs.dat <- list(n=n, nst=n.st, y=y, x=x, stid=stid,              

                   minYr=min(L), maxYr=max(U))

  for (i in 1:n.chains)

    inits[[i]] <- list(beta0=rnorm(n.st,130,2),

                       beta1=rnorm(n.st, 0, 0.1),

                       delta=runif(1, 0, 3), 

                       phi=runif(n.st, max(L), min(U)),

                       sigmay=runif(1, 0, 3), 

                       sigma=runif(2, 0, 1),

                       mu0 = rnorm(1), mu1 = rnorm(1))

  parameters <- c("beta0","beta1","delta", "phi", "mu0", 

                  "mu1","mu2","mu3","sigma", "sigmay")

  return(list(para=parameters, data=bugs.dat, inits=inits,

              meanYR=mean(xx), sdYR=sd(xx), n.chains=n.chains,

              model="

 /*  multilevel hockey stick model */ 

data{

  int n; 

  int nst;

  real y[n]; 

  int stid[n]; // station ID

  real x[n]; // standardized year

  real minYr;

  real maxYr;

}

parameters{

  real beta0[nst]; // intercept at change point

  real beta1[nst];

  real delta[nst];

  real phi[nst];

  real mu0;

  real mu1;

  real mu2;

  real mu3;

  real sigmay;

  real sigma[4];

}

model{

    mu0 ~ normal(0, 100);

    mu1 ~ normal(0, 100);

        mu2 ~ normal(0, 100);

        mu3 ~ normal(0, 100);

        beta0 ~ normal(mu0, sigma[1]);

        beta1 ~ normal(mu1, sigma[2]);

        phi ~ normal(mu2,sigma[3]);

        delta ~ normal(mu3,sigma[4]);

for (i in 1:n){

   y[i]~normal(beta0[stid[i]]+

(beta1[stid[i]] - delta[stid[i]] * int_step(x[i]-phi[stid[i]]))*(x[i]-phi[stid[i]]), sigmay);

}

}

"

  ))

}

keep4 <- USlilac$STID==354147 | USlilac$STID==456974 | 

         USlilac$STID==426357 | USlilac$STID==456624

input.to.bugs <- bugs.in(infile=USlilac[keep4,])

n.chains <- input.to.bugs$n.chains

n.iter <- 100000

n.keep <- 2000

thin <- max(c(1,floor((n.iter)*n.chains/n.keep)))

fit <- stan(model_code=input.to.bugs$model,

            data=input.to.bugs$data, init=input.to.bugs$inits,

            pars=input.to.bugs$para,

            iter=n.iter, thin=thin, chains=n.chains)

print(fit)

quartz()

plot(fit)

My experience with Stan is a pleasant one, especially the level of help provided by the Stan team.  The last question I submitted to them was at 10:30 pm and I got an answer back around 11:00 pm the same day!