# Interpreting predict() with gamlss() and beta-binomial family=BB

Note after answer is posted: The issue here was actually about formatting the data of the dependent variables in the formula when using (beta-)binomial families. Not about predict(). See the answer.

The output of predict(...,type="response"), when I'm using gamlss(...,family=BB), does not match the proportion in the data. But when I exponentiate the parameter values (instead of using predict()), the result does match the proportion in the data. What am I not understanding about the output of predict.gamlss()?

Here's a minimal example.

First, I generate some data. Each "point" has a total count, N, and a number of successes S ($$0 \le S \le N$$). N comes from a negative-binomial; S comes from a beta-binomial relative to N.

library("gamlss")
set.seed(47405)
N = as.integer( rnbinom( n=100 , mu=250 , size=500 ) )
S = rBB( n=length(N) , mu=0.8 , sigma=0.1 , bd=N )
theData = data.frame( S=S , N=N )
str(theData)

'data.frame':   100 obs. of  2 variables:
$$S: int 163 172 245 172 203 168 240 168 179 215 ...$$ N: int  234 249 257 246 225 220 265 198 238 231 ...


Here's a boxplot of the data proportions, S/N:

dataProp = sum(theData$S)/sum(theData$N)
boxplot( S/N , data=theData , notch=TRUE , ylab="S/N" )


Notice, above, that the overall proportion in the data, S/N, is 0.7917529. In my mind, this value is what predict(...,type="response") should (approximately) match.

But here's what actually is produced:

gamlssInfo = gamlss( cbind( S , N ) ~ 1 ,
family = BB , # beta-binomial
data = theData )
gamlssSumm = summary( gamlssInfo )
gamlssSumm

              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -0.2411631 0.01549786 -15.56105 3.131276e-28
(Intercept) -5.6019584 0.23435204 -23.90403 1.122643e-42


Above, the first intercept is mu and the second intercept is sigma.

Finally, the output of predict() vs the result of exponentiating mu, compared with the data:

c( "predict()" = predict( gamlssInfo , type="response" )[1] ,
"exponential()" = exp( gamlssSumm[1,1] ) ,
"data" = dataProp )

    predict() exponential()          data
0.4399998     0.7857135     0.7917529


Notice, above, that the output of predict() is way less than the data proportion, but exponentiating the estimate matches the data pretty well. What is predict() doing in this case?

I emailed the primary creator of the gamlss package, Professor Mikis Stasinopoulos, and he was very kind to provide an extensive explanation and correction to my code.
The solution is simple: Where I had specified gamlss( cbind( S , N ) ~ 1 , ... it should be gamlss( cbind( S , N-S ) ~ 1 , .... That is, the count data should be specified as successes and failures, not as successes and totals. Then the output of predict.gamlss() matches the data proportion.
Professor Stasinopoulos also verified that the output of gamlss matches the output of glm when using the binomial family (instead of the beta-binomial family).