5
$\begingroup$

I have some data like this:

id    pop   var
1     593   51
2     592   31
3     346   20
4    1214   70
5    1063   66
6    1370   71

each id represents a school. There are around 1847 of them, of which 1833 are unique . pop is the number of pupils in the school. var is the number of pupils with a particular trait. Interest in centered on the ratio of var/pop . I have used a binomial regression:

fit.bn <- glm(cbind(var,pop-var)~1, data=dt, family=binomial(link=logit))

which produced the results:

             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.748516   0.003295  -834.2   <2e-16 ***

And to get the estimate for the proportion and it's confidence interval I used:

invlogit(coef(fit.bn))
invlogit(confint(fit.bn,level = 0.95))

where

invlogit <- function (x) { return(exp(x)/(1+exp(x))) }

which produced:

(Intercept) 
 0.06017054 

     2.5 %     97.5 % 
0.05980603 0.06053643

Is this a valid confidence interval for these data ? I was surprised that it was so narrow, but I guess this is because the total number observations is high (altogether, sum(var)= 98008 and sum(pop)=1628981). Nevertheless, I am still unsure about this, due to:

sd(var/pop)
0.01373888

In other words the standard deviation of these proportions is around 20 times larger than their confidence interval.

Am I doing something stupid ? It just doesn't seem quite right.

I wondered about using a variance components model to try to account for variability between schools and within schools. I did this with

require(lme4)
fit.lme4 <- glmer(cbind(cbind(var,pop-var)~1+(1|id), data=dt, family=binomial(link=logit))

This produced:

Random effects:
 Groups   Name        Variance Std.Dev.
 id      (Intercept) 0.024125 0.15532 
Number of obs: 1847, groups: id, 1814

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.776882   0.005423  -327.7   <2e-16 ***

So I computed another confidence interval by:

invlogit(-1.776882  + (c(-1.96,1.96)* 0.005423))
0.1433781 0.1460089

which is wider than the first one, but the estimate is very different. Is this a valid procedure ? How do I interpret this estimate and confidence interval as compared to the first one ? And how do I interpret the variance of the random effects ? Am I barking up the wrong tree ?

Update:

Here is a histogram of var/pop enter image description here

Improve this question
$\endgroup$
6
  • $\begingroup$ I assume you have multiple 'var's and 'pop's per school? If so the glmer model is correct. CI - remember that the confidence interval is inversely proportional to the n, so if you have lots of schools your confidence interval will get smaller. Also, it doesn't seem like you have very much variability in the proportions for the snippet of data you gave. So the CIs seem ok to me. $\endgroup$
    – Matt Albrecht
    Nov 4, 2012 at 3:55
  • $\begingroup$ @MattAlbrecht, no, the data is aggregated so there is just 1 var and pop per school.......... $\endgroup$
    – Joe King
    Nov 4, 2012 at 7:07
  • $\begingroup$ I get error messages with glmer when there is only one measure per id: {glmerx <- glmer(cbind(var, pop-var) ~ 1 + (1|id), data=df1, family=binomial)} "Number of levels of a grouping factor for the random effects is equal to n, the number of observations". You have some double up id numbers in there as indicated by the glmer output "Number of obs: 1847, groups: id, 1814" 1847-1814 = 33 double ups $\endgroup$
    – Matt Albrecht
    Nov 4, 2012 at 8:44
  • $\begingroup$ @MattAlbrecht sorry you are right. I was told the data are all for unique schools, so I don't know whey there are duplicate IDs. $\endgroup$
    – Joe King
    Nov 4, 2012 at 12:02
  • $\begingroup$ @MattAlbrecht I have updated the question to reflect the duplicate ids. $\endgroup$
    – Joe King
    Nov 4, 2012 at 12:45

0

Reset to default

Your Answer

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