How to compare centiles from different models?

I am comparing the centiles from different GAMLSS models. Which model has a better performance. Why?

For model 1, the number of cases below 0.4 centile is 0.5. Is this possible? The number of cases below (forexample) 0.4 centile should not be less than 0.4? What does it mean when the number of cases of 0.4 centile is more that 0.4?

********  Model 1 ********
% of cases below  0.4 centile is  0.5006258
% of cases below  10 centile is  9.803922
% of cases below  50 centile is  49.35336
% of cases below  90 centile is  90.2378
% of cases below  99.6 centile is  99.70797
********  Model 2 ********
% of cases below  0.4 centile is  0.458907
% of cases below  10 centile is  9.595327
% of cases below  50 centile is  48.47726
% of cases below  90 centile is  89.27826
% of cases below  99.6 centile is  99.70797
********  Model 3 ********
% of cases below  0.4 centile is  0.458907
% of cases below  10 centile is  9.386733
% of cases below  50 centile is  47.43429
% of cases below  90 centile is  88.77764
% of cases below  99.6 centile is  99.58281
********  Model 4 ********
% of cases below  0.4 centile is  0.5006258
% of cases below  10 centile is  9.762203
% of cases below  50 centile is  50.31289
% of cases below  90 centile is  90.15436
% of cases below  99.6 centile is  99.66625


• I'd rather compared the worm plots and Q statistics of the models. Jan 8 at 16:41
• How many data values do you have? The answers depend, roughly, on whether you have more or less than ten thousand. Your legend and the dense overplotting of symbols might be obscuring many thousands of points, making it impossible to guess accurately.
– whuber
Jan 8 at 16:42
• @ whube, thanks for the reply. I have around 3000 observations. Jan 8 at 16:44
• @ COOLSerdash, thanks for the reply. I have two predictors (age and hight) in this project and Q statitics of exproximately all models (by using: BCCGO and No distributions, such as m4 <- gamlss(y ~ log(W) + pb(log(Age)), sigma.fo =~pb(log(Age)), nu.fo=~log(Age), family = BCCGo(mu.link = "log"), data=DAT1.F)) ) ) identifies a region of age and weight when the model is inadequate. So, I decide to compare them based on their centiles. Jan 8 at 16:55

Some of the deviations you point to can be due to random variation alone. There are simple ways to test them.

When there are independent observations of a variable whose conditional distributions are continuous, then the chance that an observation is less than a quantile $$q$$ is -- by the very definition of quantile -- equal to $$q.$$ Thus, the chance the observation lies between quantiles $$q$$ and $$q^\prime$$ is $$|q-q^\prime|.$$

Your centiles partition the observations into six non-overlapping bins bounded by the quantiles $$(0, 0.4, 10, 50, 90, 99.6, 100)/100.$$ Accordingly, the likelihood (conditional on the model) of observing counts in this bins is Multinomial. Specifically, the probability of observing a count vector $$(n_1,n_2,\ldots, n_6)$$ (when the model is correct across the entire range of $$x$$ values in the data) is

$$\mathcal L(\mathbf n) = \binom{n}{n_1,n_2,\ldots,n_6} \left(\frac{0.4 - 0}{100}\right)^{n_1}\left(\frac{10 - 0.4}{100}\right)^{n_2}\cdots \left(\frac{100 - 99.6}{100}\right)^{n_6}$$

where $$n = n_1 + n_2 + \cdots n_6$$ is the total.

Assuming there are approximately $$n=3000$$ observations, the (natural) log likelihoods of these four models work out to $$(-17.1, -18.5, -21.2, -16.7)$$ in order. The largest likelihood of $$-16.7$$ identifies model (4) as having the highest likelihood. This does not mean it's a great model: to assess that, we should look at the details. A method that is easy to interpret is to plot the differences between the observed proportions $$x_i$$ and the proportions indicated by the centiles, $$q_i,$$ divided by their standard errors. The latter equal $$\sqrt{q_i(1-q_i)/n}.$$ Although these standardized residuals aren't truly independent, the expected counts in each bin are large enough that this makes no difference. Generally we hope most standardized residuals for any model will lie between $$-2$$ and $$2,$$ perhaps with a few of them slightly exceeding these thresholds.

Evidently models 4 and 1 have consistently near-zero residuals across the range of quantiles: they don't appreciably under- or over-estimate the conditional quantiles anywhere. Model 2 tends to have too many low counts (that is, it underestimates the middle quantiles) while Model 3 is significantly worse than any of the others.

Significance can be determined by the Likelihood Ratio (LR) test: the p-value is the tail area of a $$\chi^2(1)$$ distribution evaluated at twice the difference of log likelihoods. Again using $$n=3000,$$ I find p-values of $$(0.372, 0.057, 0.003, 1.000)$$ for these four models (relative to the best). Even accounting for the multiple testing, $$0.003$$ is so low we ought to reject model 3 as poorer than the best model. However, we shouldn't reject the other models outright: their differences could be due to chance alone.

Here is R code to produce the image. Most of it reads and formats the data and then initializes the plot. The last block of code does the calculations.

x <- scan(text = "% of cases below  0.4 centile is  0.5006258
% of cases below  10 centile is  9.803922
% of cases below  50 centile is  49.35336
% of cases below  90 centile is  90.2378
% of cases below  99.6 centile is  99.70797
% of cases below  0.4 centile is  0.458907
% of cases below  10 centile is  9.595327
% of cases below  50 centile is  48.47726
% of cases below  90 centile is  89.27826
% of cases below  99.6 centile is  99.70797
% of cases below  0.4 centile is  0.458907
% of cases below  10 centile is  9.386733
% of cases below  50 centile is  47.43429
% of cases below  90 centile is  88.77764
% of cases below  99.6 centile is  99.58281
% of cases below  0.4 centile is  0.5006258
% of cases below  10 centile is  9.762203
% of cases below  50 centile is  50.31289
% of cases below  90 centile is  90.15436
% of cases below  99.6 centile is  99.66625",
what = c(rep(character(), 7), numeric()))

n <- 3e3
p <- c(0.4, 10, 50, 90, 99.6) / 100
x <- matrix(as.numeric(matrix(x, 8)[8,]), length(p))

q <- apply(x, 2, function(x) pbinom(round(x/100 * n), n, p))

plot(0:1,c(-3,3), type = "n", xlab = "Quantile", ylab = "Z",
main = "Standardized Residuals")
abline(h = 0, lwd = 2, col = "gray")

sapply(1:ncol(q), function(j) {
r <- x[,j] / 100 - p
se <- sqrt(p * (1-p) / n)
lines(p, r/se, col = hsv((j-1)/6, .9, .5))
points(p, r/se, pch = as.character(j), col = hsv((j-1)/6, .9, .5))
})

• thank you so much for time and kindness. I just did not understand how you could make the plot of standardized residuals based on the centiles results. Could you please share the programming function or formula with me? Agian, thank you so much Jan 8 at 20:50
• great, I am very grateful for your advice Jan 9 at 0:47