0
$\begingroup$

I'm performing a meta-analysis of studies that compared treatment and control groups from pre-test to post-test using means of count data. Other authors have performed similar analyses using raw counts (Welsh & Farrington, 2008, Farington & Welsh (2013), and Jones (2005)). What I would like to do is drop variables into the model one-by-one to test whether they significantly moderate the treatment effects. I have created a glm log reg that models the random effects using Jones' (2005, p37) code:

id <- c(1:10)
a <- c(0.28, 2.26, 0.44, 1.22, 1.80, 1.45, 1.09, 2.00, 2.27, 2.11)
b <- c(0.15, 1.02, 0.16, 0.51, 0.76, 0.51, 0.13, 0.59, 0.57, 0.79)
c <- c(0.30, 2.21, 0.21, 1.05, 1.78, 1.26, 0.84, 1.86, 2.17, 2.58)
d <- c(0.42, 2.27, 0.28, 0.59, 0.91, 0.54, 0.51, 0.73, 0.90, 0.85)
mod1 <- c("Yes", "Yes", "No", "Yes", "No", "Yes", "Yes", "Yes", "No", "Yes")
mod2 <- c(2010, 2013, 2010, 2017, 2001, 2009, 2012, 2006, 2015, 2015)
bef <- c(a,c)
aft <- c(b,d)
n <- bef+aft 
treat <- scan(,"")
1: E 
2: C
treat<- gl(2,10,20,labels = treat)
study<- gl(10,1,20,labels = id)
model_effect_sizes <- glm(bef/n ~ treat + study,family=quasibinomial, weights=n)
summary(model_effect_sizes)

I'd like to add moderator variables into the model one by one but am unsure on how to code them. I have tried simply doubling the variable as follows so that it is the right length to be include in the model, but its coefficients are NAs. I have read here that its likely due to multicollinearity: 

mod1a <- as.factor(c(mod1, mod1))
moderator1a<-glm(bef/n ~ treat + study + mod1a,family=quasibinomial, weights=n)
summary(moderator1a)

I've also tried coding the second half of the vector with 0s and including the variable as a numerical vector which produces coefficients for the moderator:

mod1b <- c(2, 2, 1, 2, 1, 2, 2, 2, 1, 2, 0,0,0,0,0,0,0,0,0,0)
moderator1b<-glm(bef/n ~ treat + study + mod1b,family=quasibinomial, weights=n)
    summary(moderator1b)

But I'm uncertain if this produces reliable coefficients, and would also prefer to be able to include mod1 as a factor. Similarly, if I try this approach with with mod2 the intercept and the treatment variable produce almost identical coefficients:

mod2a <- c(2010, 2013, 2010, 2017, 2001, 2009, 2012, 2006, 2015, 2015,0,0,0,0,0,0,0,0,0,0)
moderator2a<-glm(bef/n ~ treat + study + mod2a,family=quasibinomial, weights=n)
        summary(moderator2a)

So the question comes down to: how can these moderators be included in this model to avoid singularity and to keep their original structure? Thanks for your help,

Jones, H. E. (2005). Measuring Effect Size in Area-Based Crime Prevention Research. Statistical Laboratory. Cambridge, UK, Cambridge University. Masters of Philosophy.

Improve this question
$\endgroup$
2
  • $\begingroup$ Hey.. you have only 10 observations and you have 12 variables, dim(model.matrix(~ treat + study + mod1a)), because the factors will be expanded. You cannot fit this model unless you group your factors or you have more observations $\endgroup$
    – StupidWolf
    Commented Mar 4, 2020 at 9:05
  • $\begingroup$ Great point @StupidWolf, and one that I hadn't thought of. How could we group the factors in this case, I'm not sure what you mean by this? $\endgroup$
    – m.tripley
    Commented Mar 4, 2020 at 21:25

1 Answer 1

Reset to default
0
$\begingroup$

Here are a few comments:

  • why not use the original count data?
  • I think you should use study as a random effect
  • if you switch "study" and "mod1a" in your model "moderator1a", you will get a result for "mod1a", but an NA for a level of "study". This is linked to the comment from @StupidWolf: you have too little observations for your number of variables.
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for your response. There are differing sample numbers for both treatment and control at each of the time points, and the count data does not reflect this and overestimates the pooled effect. Yes, part of the problem is that we want to model the random effects and need to keep 'study' in the model to do so. $\endgroup$
    – m.tripley
    Commented Mar 4, 2020 at 21:20
  • $\begingroup$ To specify random effects, you should use for instance glmer() from the package lme4. The syntax would look something like: glmer(bef/n ~ treat + (1|study), family = binomial, weights = n). glmer() does not accept quasibinomial family, and to use the binomial family, you would need to use either a binary response variable (0 or 1) or of the type cbind(success, failure) for each study (i.e., the count data). $\endgroup$
    – KrisAnathema
    Commented Mar 5, 2020 at 8:03

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.