Solved – use binomial distribution for beta distribution when using lme4 package (glmer function)

beta distributionbinomial distributionlme4-nlmemeta-analysismixed model

I am conducting a meta-analysis dealing with the proportional data. I would like to use "proportional mean scores (0-100%; mean scores divided by the maximum test score" as a response variable and weight using inverse variance (1/variance). I also would like to use multilevel modeling because some of the response variables are nested in each article. I am thinking the simplest command for the glmer function will be like this:

res <- glmer (proportional_mean_score ~ moderator1 * moderator2 +(1|article_number), family = binomial (link=logit), weights = 1/variance, data = dat))

I was told that I should be using the beta distribution instead of binomial, and weights should represent the total of total count for the proportion. However, I cannot find any package which let me use beta-distribution while specifying the weights and random effects. I first tried with the lmer function but the estimation values go beyond 100%, so I decided to use glmer function with binomial distribution. Can I just use the glmer function in the lme4 package and weight with variance instead of the total counts or maximum scores?

I was told that I do not have to stick with the lme4 package; instead, I was recommended to use the brains package. https://journal.r-project.org/archive/2018/RJ-2018-017/RJ-2018-017.pdf

Uing the brms package, I can use bata models while specifying random effects. One flip side would be that I have to learn how to interpret the results from the Bayesian perspective.

I guess this question does not seem to attract that much attention. However, in case, some one who got to my question wondering the same thing, I would like to share what I learned.

P.S.
No offence to the lme4 package. I like lme4 a lot too.

Best Answer

The beta mixed effects model is available in the GLMMadaptive; for an example, check here (the purpose of the example is to show how you can define a user-specified family object; in your case you could directly use mixed_model(..., family = beta.fam())).

To account for the differences in the total counts you could perhaps consider adding an offset instead of a weight.

Related Solutions

R – Which R Package to Use for Calculating Component Parameters in Mixture Models?

Try mixdist

Here's an example:

library(mixdist)  

#Build data vector "x" as a mixture of data from 3 Normal Distributions  
x1 <- rnorm(1000, mean=0, sd=2.0)  
x2 <- rnorm(500, mean=9, sd=1.5)  
x3 <- rnorm(300, mean=13, sd=1.0)  
x <- c(x1, x2, x3)  

#Plot a histogram (you'll play around with the value for "breaks" as    
#you zero-in on the fit).   Then build a data frame that has the  
#bucket midpoints and counts.  
breaks <- 30  
his <- hist(x, breaks=breaks)  
df <- data.frame(mid=his$mids, cou=his$counts)  
head(df)  

#The above Histogram shows 3 peaks that might be represented by 3 Normal  
#Distributions.  Guess at the 3 Means in Ascending Order, with a guess for  
#the associated 3 Sigmas and fit the distribution.  
guemea <- c(3, 11, 14)  
guesig <- c(1, 1, 1)  
guedis <- "norm"  
(fitpro <- mix(as.mixdata(df), mixparam(mu=guemea, sigma=guesig), dist=guedis))  

#Plot the results  
plot(fitpro, main="Fit a Probability Distribution")  
grid()  
legend("topright", lty=1, lwd=c(1, 1, 2), c("Original Distribution to be Fit", "Individual Fitted Distributions", "Fitted Distributions Combined"), col=c("blue", "red", rgb(0.2, 0.7, 0.2)), bg="white")  


===========================  


Parameters:  
      pi     mu  sigma  
1 0.5533 -0.565 1.9671  
2 0.2907  8.570 1.6169  
3 0.1561 12.725 0.9987  

Distribution:  
[1] "norm"  

Constraints:  
   conpi    conmu consigma   
  "NONE"   "NONE"   "NONE"

enter image description here

Binomial Count Data – Should glmer, lmer Be Used or Just Average It All?

Let's start with the first, best, model:

mixed <- glmer(cbind(Correct, Max-Correct) ~ Accent + Noise_Factor + (1 | Participant),
               data=noise.data, family="binomial")

where I've tidied the syntax and removed the extra reference to Accent. This seems reasonable as a starting point.

Should you ignore (by leaving out) Noise_Factor?

Surely not, if you don't have to. You put it in there for a reason I presume. Consequently the right hand side of the third model seems underspecified. In the event that it truly has no effect itself and it also does not moderate the effects of anything else then leaving it out would be reasonable, but I assume you are not in a position to assert something that strong. Certainly you don't have to.

Is the data set fit better by one of the first or second model?

That's very hard to say. The DV in the first and second models are not the same and the models are not nested so any analysis of deviance (or even an R-squared) comparison doesn't seem to make any sense.

Are you specifying the DV in the first model correctly?

Yes.

Recommendation

In sum, I'd recommend you start with model one as above. If you wonder about interactions then you can use the anova function to compare models with and without them because they are nested and have the same DV.

You might look also to see whether the effect of other covariates vary by subject. Currently subjects vary in the ratio of correct to incorrect responses they make on average, but they might possibly differ in more than that. This is quite straightforward to check in lme4.