[R-sig-ME] Question on calculating the conflidence limits on the intraclass correlation in a nested binomial mixed model

Thu Sep 19 11:19:56 CEST 2013

Dear all,
Would any of you happen to have any further ideas on how to accomplish profiling intraclass correlations calculated from a lme4 glmer object?
So far I just profiled the random effect variances (see below, keeping the total variance constant), but in actual fact I would like to profile the intraclass correlation (ratios of random effect variances over total variance) directly - is there any way to achieve that, or would this require delving into the profile.R code?

Any advice or pointers would be much appreciated!

Cheers,
Tom

icc.conflims = function(formula,dataset,family,level=0.9) {  # function to calculate conf lims on ICC based on computing likelihood profiles of random effect variances
  options(contrasts=c("contr.sum", "contr.poly"))
  fit=eval(bquote(glmer(.(formula),.(family),data=dataset)))
  vcor=VarCorr(fit)
  vcomps=as.numeric(vcor)
  V.link=ifelse(family=='binomial(link=logit)'|family=='binomial',(pi^2)/3,ifelse(family=='binomial(link=probit)',1,0))[1]
  V.tot=sum(vcomps) + V.link
  ICCs=vcomps/V.tot # intraclass correlations
  names(ICCs)=names(vcor)

  pr=profile(fit,signames=F)
  ICCclims=as.data.frame((confint(pr,level=level)^2)/V.tot) # lower and upper confidence limits on the random effect intraclass correlations
  ICCclims=ICCclims[1:length(names(vcor)),]
  rownames(ICCclims)[1:length(names(vcor))]=names(vcor)

  return(ICCclims)
}

dataset<-read.table("http://www.kuleuven.be/bio/ento/temp/data.txt",header=TRUE)
dataset$colony=as.factor(dataset$colony) 
dataset$patriline=as.factor(dataset$patriline) 

icc.conflims (trait~1+colony+(1|patriline),dataset=dataset,family=binomial(link=logit),level=0.9)
# 	                 5 %      95 %
# patriline 0.006382299 0.3027519

-----Original Message-----
From: ONKELINX, Thierry [mailto:Thierry.ONKELINX at inbo.be] 
Sent: 17 September 2013 11:11
To: Tom Wenseleers; r-sig-mixed-models at r-project.org
Subject: RE: [R-sig-ME] Question on calculating the conflidence limits on the intraclass correlation in a nested binomial mixed model

Dear Tom,

A random effect with only two levels is not a good idea. You will get very unreliable variance estimates at best. Note the profiled confidence interval goes from 0 to Inf... Therefore I recommend fitting group as a random effect.

Furthermore you profile the variance of random effects, but not the total variance. I think you should profile both, since they are correlated. Or calculate a profiled Vtot on the profile of the sigmas.

Best regards,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Tom Wenseleers
Verzonden: maandag 16 september 2013 18:46
Aan: r-sig-mixed-models at r-project.org
Onderwerp: [R-sig-ME] Question on calculating the conflidence limits on the intraclass correlation in a nested binomial mixed model

Dear all,
I calculated the intraclass correlations at two nested levels in a binomial mixed model from the random effect variances, and tried to calculate the 95% confidence limits on them from the confidence limits on the variance components, obtained by profiling the likelihood:
(I have 2 groups and ca 10 subgroups per group)

  fit = glmer(trait ~ 1 + (1|group/subgroup), family = binomial(link="logit"),data=mydata)
  vcomps = as.numeric(VarCorr(fit))
  V.link = (pi^2)/3
  V.tot= vcomps [1] + vcomps [2] + V.link
  ICC.group = results[2]/V.tot
  ICC.subgroup = results[1]/V.tot
 icc=c(ICC.group,ICC.subgroup)

  pr=profile(fit,signames=F)
  prclims=(confint(pr,level=level)^2)/V.tot
  low=prclims[2:1,1] # upper and lower CIs on the ICCs
  up=prclims[2:1,2]
  results = round(data.frame(icc, low, up),3)
  row.names(results) = c("group", "subgroup")
  results
#                 icc        low    up
# group    0.000    0.000   Inf
# subgroup 0.095 0.006 0.303

My question is whether this approach would be OK? I only consider the uncertainty in the estimation of the among-group and among-subgroup variances here, keeping the total variance constant. Is that OK? And what is the implicit assumption in calculating the confidence limits this way - e.g. what sources of uncertainty are considered then? (I also came across this reference, https://stuiterproxy.kuleuven.be/doi/abs/10.1080/,DanaInfo=www.tandfonline.com+03610910903324834, about calculating the conf lims on the ICC using likelihood profiling, not sure if that is relevant here).

I also tried calculating the confidence limits using parametric bootstrapping (using bootMer), but there I ran into problems with non convergence in the fits, owing to the small size of my dataset. And then I tried nonparametric cluster bootstrapping, whereby I resampled over subgroups (resampling the same nr of subgroups pers groups as I had initially), but I am not sure that that is entirely correct either (e.g. it doesn't take into account the binomial scatter per subgroup cluster, and also resampling within clusters biases my estimate, so I can't do that). Or does anyone perhaps happen to have other ideas about how to calculate confidence limits on my estimated intraclass correlation?

Cheers,
Tom

_______________________________________________________________________________________

Prof. Tom Wenseleers
*      Lab. of Socioecology and Social Evolution
           Dept. of Biology
           Zoological Institute
           K.U.Leuven
           Naamsestraat 59, box 2466
           B-3000 Leuven
           Belgium
* +32 (0)16 32 39 64 / +32 (0)472 40 45 96
* tom.wenseleers at bio.kuleuven.be
http://bio.kuleuven.be/ento/wenseleers/twenseleers.htm

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * * Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.