I think you can do so much better things with Bayesian inference than hypothesis testing. I would start with the estimation of the parameters of two independent normal distributions under a flat prior (see Gelman et al. 2003, Bayesian data analysis.) and check the assumptions using posterior predictive checks (Gelman et al. 2003; chapter 6). If the model is able to replicate the important aspects of your data, this approach is a good way to answer your hypothesis. If this model does not work for your data, i.e. important aspects are not explained well, you can use regression models like generalized linear models that are more appropriate for your data, or a t-distribution instead of a normal distribution... I think nonparametric hypothesis tests are not a good choice.

Here is the code for estimation of the parameters of two independent normal distributions and for model checking by posterior predictive checks:

## Calculate posterior distributions for mu1 and mu2; flat prior

## y1 is sample 1; y2 is sample 2

n1 = length(y1)

n2 = length(y2)

y1.mean = mean(y1)

y2.mean = mean(y2)

s21 = var(y1)

s22 = var(y2)

n.sim = 1000

sigma21 = 1/rgamma(n.sim, shape = (n1 - 1)/2, rate = (n1 - 1)/2 * s21)

sigma22 = 1/rgamma(n.sim, shape = (n2 - 1)/2, rate = (n2 - 1)/2 * s22)

mu1 = rnorm(n.sim, mean = y1.mean, sd = sqrt(sigma21/n1))

mu2 = rnorm(n.sim, mean = y2.mean, sd = sqrt(sigma22/n2))

## check your model using posterior predictive checks

## replicate the data

y1.rep = matrix(ncol = n1, nrow = n.sim)

y2.rep = matrix(ncol = n2, nrow = n.sim)

for(ii in 1:n.sim){

y1.rep[ii,] = rnorm(n1, mean = mu1[ii], sd = sqrt(sigma21[ii]))

y2.rep[ii,] = rnorm(n2, mean = mu2[ii], sd = sqrt(sigma22[ii]))

}

## display the replicated data, maybe the y-axis has to be adjusted using ylim = c(0, ??)

par(mfrow = c(2,2))

hist(y1, main = '', freq = FALSE, xlab = 'y1')

for(ii in 1:100){

points(density(y1.rep[ii,]), type = 'l', col = '#00000025')

}

hist(y2, main = '', freq = FALSE, xlab = 'y1')

for(ii in 1:100){

points(density(y2.rep[ii,]), type = 'l', col = '#00000025')

}

plot(y1.rep[1,], pch = 20, col = '#00000025', ylim = range(y1))

for(ii in 2:100){

points(y1.rep[ii,], pch = 20, col = '#00000025')

}

points(y1, pch = 20, col = '#0085FF')

plot(y2.rep[1,], pch = 20, col = '#00000025', ylim = range(y2))

for(ii in 2:100){

points(y2.rep[ii,], pch = 20, col = '#00000025')

}

points(y2, pch = 20, col = '#0085FF')

## if the important aspects of your data are replicated, proceed by analyzing the difference

## credible interval of the difference

quantile((mu1 - mu2), c(0.025, 0.975))

plot(mu1, mu2, pch = 20, col = '#00000050'); abline(0:1)

hist(mu1-mu2)

As a practicing Bayesian I suggest that if your interests are in a fully non-parametric hypothesis confirmation, you are better off avoiding Bayes methods entirely and simply confirm your hypothesis using permutation methods. The method mentioned above is very parametric, it is essentially the posterior sampled Bayes version of a two-variance t-test, but nothing in your initial question suggests that you are doing a 1-dimensional test of means, and if you have a hypothesis comparing proportions, variances, tails, and such, you are likely not going to perform the appropriate test.

If you are doing a general U test, what is the U statistic you are dealing with?

Is your sample size so small that can get away with modelling the sample density with a parametric form (perhaps a long tailed t, a bimodal density if both samples are definite mixtures, perhaps gammas for positive numbers,...)

If your sample size is large, can a kernel density method approximate the observed samples? If you have multivariate data the warning is that there's probably no helping from kernel non-parametrics if dimension is larger than 3.

How many simultaneous comparisons are you trying to perform? What is to be gained by reporting a Bayes factor, when permutation p-values against the null are likely more commonly accepted?

If your permutation scores are insignificant, the Bayes factors will be a wash, but as such I'm not sure the time spent developing a generalized procedure is going to benefit you more than either piecing together more data, or figuring out whether your test statistic can be modeled with a growth model more specific than non-parametrics.

Dear Alan, I do agreed with you .

Paul, thank you very much for your code. It works really good with some test data that i added but i can see the problem here that indeed my data for each variable are not normally distributed since they are ordinal. I can see though that the histogram being produced indicates the difference between mu1-mu2 but not sure if i can publish the results with this method, because of the violations of normality.

Alan, my sample size is about 100-110(currently the survey is in progress). I want to compared between independent groups between two conditions. That's why the appropriate test is a Mann-Whitney U test. The problem here is that the p value is not significant which could indicate that maybe the null is true. Frequentist statistics won't allow me to accept the null, so i need an alternative and Bayesian statistics allow me to do that.

I already used Bayes factor for some interval data that i had and it was extremely useful. For ordinal data i don't know what procedure i could use to perform an analysis of variance between two independent groups but i know it has to support the null. So i need to see how does my prior compare to the posterior at d=0 but for ordinal.

That's why since i found this paper i wondered if it could be ideal for this. There has to be a procedure for bayesian statistics that compared two groups and gives an estimate of how different or not they are. Do you think this is possible?

With the above post I wanted to say that it is quite easy to verify model assumptions by using posterior predictive checks. Clearly, ordinal data are not normally distributed, thus in your situation the above code is not suitable. Actually, the Mann-Whitney U-test is just as inappropriate for ordinal data. In the case of ordinal data regression models for ordinal responses like multinomial logistic regression (function MCMCmnl in the package MCMCpack) or ordered probit regression (function MCMCoprobit in the package MCMCpack) are common. See also the task view about Bayesian inference on CRAN: http://cran.r-project.org/web/views/Bayesian.html

I don't care much about tests of normality or the homogeneity of variances (Kologorov-Smirnoff-test, Levene-test, etc...) before doing the test. If your model (in my opinion, all tests can be seen as models) is able to replicate your data well, you can use it for inference. If not, adjust the model or use a more complex model. Unfortunately, almost no one reports the fit of a model. For me, the most obvious reason for that is that most people have no idea how to check their models. Or maybe they are afraid of observing a serious defect of their model because they don’t know how to fix it. That's why I preach the following two points.:

1. go beyond simple t-tests and AN(C)OVAs and

2. check the fit of your model.

Posterior predictive checks are the easiest way to deal with point 2.

I find it quite disturbing to mix Bayesian inference with hypothesis testing, because I use Bayesian inference to get rid of major disadvantages that come along with hypothesis testing. But I think this would be a longer discussion about the philosophy of Bayesian inference.

Paul, thank you very much for your suggestions. Indeed i can see your reasoning. I had my mind fixed in an analysis of variance test that didn't even consider that i could see it as a regression and work like that.

Ordinal probit model seems to be the way to go. Just that i am 100% certain can you verify this? I have an ordinal 7 point likert scale as the predicted variable and the predictor is a 10 point ordinal likert scale. An alternative predictor that i also have is a categorical variable with 3 values.

Michael,

I think Paul has gotten you on a better track now. Having noted that survey data fits into an ordered model, hopefully logit or probit gives a good fit to your survey scores.

Sometimes surveys are much more bimodal than a fixed layering system would admit, if it turns out respondents never give a 3 and only give 1's and 5's you might be better off collapsing columns. In fact, you might just want to collapse a few of your rows together (as perhaps, grouping answers "2-3-4" together on a 5 point scale so you have 1, 2-4, 5 as rows) and see if you have a signficant p-value on a test for independence after doing this. If that happens, you've got a "chi-squared test is a test of sample size" issue going on.

I think though, using the ordered logit/probit model, it will give you a picture of a regression effect of the "treatment" (your indicator that indicates the 2 populations), and you'll be able to say something like "Group B has between 5% and 10% more likelihood to report a 4-5 all other things considered" type statement, which is more quantitatively useful than a straight rejection/acceptance. MCMCpack seems to have a standard Bayes version for ordered logit that should make it easy to do.

Alan

The difference between the probit and the logistic model should be small. You should also try the function polr in the MASS package which implements the proportional odds cumulative logistic and probit model. These models are a bit easier to interpret. In all cases you should plot your data to get a feeling of the relations between your response and the possible predictors.

In order to understand the interpretation of regression models with ordinal responses, I would recommend you to read the following lectures. Another source of inspiration is the classic book 'Categorical Data Analysis' by A. Agresti.

http://sites.stat.psu.edu/~jls/stat544/lectures/lec19.pdf

http://sites.stat.psu.edu/~jls/stat544/lectures/lec20.pdf

http://sites.stat.psu.edu/~jls/stat544/lectures/lec21.pdf

I agree with Alan that it is sometimes necessary to collapse some categories of your ordinal variables, both the response and the predictors. Theory and experience should be the main reference for this.

Paul

Thank you both for your help! I really appreciate it.

Instead of using the grouping categorical variable i can use my 10 point scale which more accurate as a predictor.

A Spearman test between the 7 and the 10 gave me rho = .117 and p = 0.037.

This is a small effect but since i achieved significance i could run a regular ordinal logistic regression and hopefully get something useful. But,

the reason which i want to try Bayesian MCMC testing is that as far as i can tell it doesn't have significance testing which i consider non useful. We all know of the shortcomings of a p value. I do perceive Bayesian as more accurate as well but harder than traditional statistics(probably because of the lack of good textbooks, by good i mean as easy as the books you can find for traditional statistical testing).

Anyway,

I made a code that performs an ordinal probit:

require(MCMCpack)

perception=as.matrix(read.table("/home/micdhack/Dropbox/Bayesian/data.txt")[,1])

reveal=as.matrix(read.table("/home/micdhack/Dropbox/Bayesian/data.txt")[,2])

mcmc

Hi, maybe the attached document helps you to understand the model and the interpretation of the parameters.

Paul

Thank you very much Paul. Indeed it looks alot similar to the regular ordinal logistic regression with a couple of minor differences. Basically by using the values produced i can calculate the estimated group membership for each group(similar to figure 5.4.6. although this was calculated for two levels). I can make a table that will have all the estimated group membership based on each gamma level. I am though confused what to with the two missing extremeties gamma1 and gamma10.

I also have another question. What about goodness of fit, test of parallel lines and p values. P values probably have no place in Bayesian statistics but what about the other two?

http://www.restore.ac.uk/srme/www/fac/soc/wie/research-new/srme/modules/mod5/4/index.html

Don't be confused about the missing gammas, I don't think you need them for the correct interpretation of the effect of your predictor.

Goodness of fit tests can be done using posterior predictive checks: Here is a little guidance for that:

1. Write a function that calculates the probability for the levels of the response for a given value of the predictor. The formula for this can be found in any textbook about ordinal regression models.

2. Take a realization of your parameters and calculate these probabilities. Use them to generate fake data of your response for each observation using a multinomial distribution. Repeat this step with (about n.sim = 200) different realizations of your parameters. As a result you should have n.sim fake data of your response.

3. Now compare the distributions of your fake response data with the one you observed. If there is no large discrepancy and most aspects (for example the frequency distribution of the response) are explained well, then there is no reason to mistrust the results of your model. Otherwise, it could be a hint that an assumption is not met and you should consider changes about your model or alternatives.

If the missing gammas are a problem in step 1, try the polr function or bayespolr function of the arm package and simulate data using the sim function of the arm package. Gelman and Hill (2007) explain in detail how to check your model assumptions using fake data simulation. (I recommend every scientist who is not an expert in statistics to read this book!)

Instead of computing p-values one usually displays the results in an appropriate way so that the variation along the estimates and therefore the effect of the predictor is visible.

Sorry that I can't help you any further, but I'm quite busy at the moment.