Finding the right fit for a monotonically-increasing distribution

Question

From N=30 subjects, I have the following distribution of responses to a certain measure of musical expertise (higher=more musical) (see histogram below). I would like to look for correlations between this measure and others, or use it as a regressor in other analyses.

1) Since this distribution is by no means normal, is it correct that I will need to use non-parametric correlation coefficients such as Spearman's, as opposed to parametric ones like Pearson's?

2) What curve (theoretical distribution) would best fit this observed distribution?

3) In case of two or more "candidates" that appear to make good fits, is it even statistically-sound to just try them all ("fish") until some goodness-of-fit measure (which one?!) is maximised?

Answer 1

With that information I would use linear model with Pearson, because there is only some data points that seem to fit fairly good to a line. If the model is not linear, or close to linear, then use Spearman. The model needs to be as simple as possible, because you can over-fit easily with low amount of data. Finding a good fit is not a goal. It is better to find something simple, but that actually holds. Already saying that some parameters do have a correlation is better that giving a model that explains completely the variation, but will fail always in the future to explain new data points. With low goodness of fit, you can still find meaningful results of correlation with good confidence intervals.

Fishing for good fits is often over-fitting. You can always explain everything with a polynomial with degree lower by one to the number of data points.

Drey has a great point :D For us that is few data points. Nothing actually concrete can be said. If you have some prior information, you can improve models and use more complex equations to explain data, without being fishy.

Answer 2

Wall of code is incoming - (sorry, I'm too lazy to upload the figures ;-).

library(tibble)
library(dplyr)
library(tidyr)
library(ggplot2)
# The data

definedInterval <- (0:3)*20 + 24

theHistogramData <- tibble::data_frame(binnedScore = 
                                         factor(c("[24,44]", "(44,64]", "(64,84]", "(84, 104]"), 
                                                levels = c("[24,44]", "(44,64]", "(64,84]", "(84, 104]"), 
                                                ordered = T), 
                                       nparticipants = c(1, 6, 10, 13),
                                       middleScore = (0:3)*20 + 34)

ggplot(theHistogramData, aes(x = binnedScore, y = nparticipants)) +
  geom_bar(stat = "identity")

# to make things easier we are gonna consider the middle score
ggplot(theHistogramData, aes(x = middleScore, y = nparticipants)) +
  geom_bar(stat = "identity")

What kind of data could produce such a plot ?

First things first: what is the response and what are the predictors? At first, I was thinking about the distribution of scores based on number of participants - but this is not what the plot suggests at a second though What is show is the number of participants (y) based on the score (x) - i.e. number of participants is the response to the predictor score. What a funny idea ;-)

The choise a model depens on the background knowledge of how the scores are actually achieved

Lets start simple: Lets assume it is a discrete probability distribution and the random variable is the middle score. Lets sample middle scores as weighted by the probability extracted from participants

uhhhSamples <- data.frame(smpls = sample((0:3)*20 + 34, 
                                         prob = with(theHistogramData, nparticipants/sum(nparticipants)), 
                                         size = 10000, replace = T))
ggplot(uhhhSamples, aes(x = smpls)) + 
  geom_bar()

But that was cheating. Sampling from only four scores is boring. A little better is this almost perfect distirubtion

theoreticalScores <- data.frame(score = 24:104) %>%
  mutate(gr = findInterval(score, definedInterval)) %>%
  inner_join(dplyr::select(theHistogramData, nparticipants) %>%
               mutate(prob = nparticipants/sum(nparticipants), gr = 1:4))
uhhhSamples2 <- data.frame(smpls = sample(theoreticalScores$score, prob = theoreticalScores$prob, size = 10000, replace = T))

uhhhSamples2summary <- uhhhSamples2 %>%
  mutate(gr = findInterval(smpls, definedInterval)) %>%
  group_by(gr) %>%
  summarize(nc = n())

ggplot(uhhhSamples2summary, aes(x = gr, y = nc)) + 
  geom_bar(stat = "identity")

Looks like the same only scaled by participants - what a surprise.

If we say there is no upper bound to the scores and we expect more and more people having expert knowledge with higer scores we could fit a linear model. From the shape we see a left shape of a parabola, a polynom of degree 2. We could take any points from the binned score but the scores in the middle seems to be the most general choise.

lmres <- lm(nparticipants ~ poly(middleScore, 2), theHistogramData)
summary(lmres)

This model looks good with high R^2 and a low error on the given data and it also seems the squared term x^2 does not play a role

theHistogramData %>%
  mutate(lmfitted = {lmres$fitted.values}) %>%
  ggplot(aes(x = middleScore, y = nparticipants)) +
  geom_bar(stat = "identity") +
  geom_line(mapping = aes(y = lmfitted), col = "red")

set.seed(4315)
predictionData <- data.frame(middleScore = sample(24:104, size = 1000,     replace = T))

predictionData %>%
  mutate(lmpred = predict(lmres, .),
         intervalGroup = findInterval(middleScore, definedInterval)) %>%
  group_by(intervalGroup) %>%
  # take mean of the prediction for each group
  summarise(lmmean = mean(lmpred)) %>%
  ggplot(aes(x = intervalGroup, y = lmmean)) +
  geom_bar(stat = "identity")

Well, that was a big round that actually does not say much. It was more or less a complicated way to show what a linear model does anyway.

So what about a theoretical distribution? If we were to impose a probabilistic interpretation to the score we could use a beta distribution; e.g.

curve(dbeta(x, shape1 = 5, shape2 = 2), from = 0, to = 1, n = 100)

Here we need to scale the scores to the interval [0,1] ... assuming that 104 is the maximal number of points

probabilisticData <- theHistogramData %>%
  mutate(participants = nparticipants/sum(nparticipants),
         scores = middleScore/104) %>%
  select(participants, scores) %>%
  mutate(cumpart = cumsum(participants))

uhhhSamples2prob <- uhhhSamples2 %>%
  mutate(normsmpls = smpls/104,
         predbeta  = pbeta(normsmpls, 6, 2))

uhhhSamples2prob %>%
  ggplot(aes(x = normsmpls, y = predbeta)) +
  geom_line() +
  geom_point(data = probabilisticData, mapping = aes(x = scores, y = cumpart), color = "red")

The first guess does not fit, but we can optimize prarameters of the beta distribution in order to get a desired coverage

optimFun <- function(pr, x, y) {
  return(sum( (pbeta(x, pr[1], pr[2]) - y)^2 ))
}
optimFun(c(5,2), 
         x = probabilisticData$scores, 
         y = probabilisticData$cumpart)
# hm, close!
optimFun(c(4,2), 
         x = probabilisticData$scores, 
         y = probabilisticData$cumpart)
# hm, closer!

respar <- optim(par = c(5 , 2), fn = optimFun, x =     probabilisticData$scores, y = (probabilisticData$cumpart))$par
# hm, closest! (+/- epsilon)

uhhhSamples2prob <- uhhhSamples2 %>%
  mutate(normsmpls = smpls/104,
         predbeta  = pbeta(normsmpls, respar[1], respar[2]))

uhhhSamples2prob %>%
  ggplot(aes(x = normsmpls, y = predbeta)) +
  geom_line() +
  geom_point(data = probabilisticData, mapping = aes(x = scores, y =     cumpart), color = "red")
# what a nice fit!

What do we learn ? We can do magic! If we don't know anything and only have 4 times 2 numbers we can do anything we want! Awesome. We can superimpose many models and many fits and they all will tell an incredible story. Putting a little bit more code around everything is just the cherry on the top of a cake! (And the cake is a lie)

Now, srsly, raw data and description of your variables and how you got their values is of utter importance! Without it, it is not possible to derive sensible conclusions.