2
$\begingroup$

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" )

Boxplot of 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?

$\endgroup$

1 Answer 1

2
$\begingroup$

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).

A giant thank you to Professor Stasinopoulos for clearing up this problem!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.