[R-sig-ME] Binomial model variance and repeatability estimates with MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Oct 21 22:05:10 CEST 2014
Hi,
The residual variance of a binary response cannot be estimated, so use
prior1 = list(R = list(V = 1, fix=1),
G = list(G1 = list(V = 1, nu = 0.002)))
In this example it is more efficient to aggregate success/failures of
an individual into a multi-trial binomial response and use:
prior2 = list(R = list(V = 1, nu=0.002))
sim.mcmc2<-MCMCglmm(cbind(Fail,Success)~1,
family="multinomial2", prior=prior2,
nitt = 260000, thin = 200, burnin = 60000,
verbose=FALSE,data=ind.data)
sim.mcmc2$VCV/(sim.mcmc2$VCV+pi^2/3)
Cheers,
Jarrod
Quoting Ned Dochtermann <ned.dochtermann at gmail.com> on Tue, 21 Oct
2014 14:44:19 -0500:
> List,
> A colleague recently asked me some questions about issues with
> estimating repeatabilities for binomial/binary data. In trying to
> help out I started playing around with simulating some data and have
> been struggling to get decent effective sample sizes even when
> parameters are being well estimated. I suspect the problem is partly
> the priors (inverse-gamma) but have had issues with other priors too
> (priors as specified in MCMCglmm are still quite unclear to me in
> general, however).
>
> While quite often repeatability gets estimated very well despite
> effective sample sizes of only a few dozen I also am getting poorly
> estimated repeatabilities (~0) with a population level repeatability
> of <.4 and sample sizes less than n=100 and reps=2, suggesting the
> inverse-gama being informative around zero. Is the poor mixing
> potentially due to using multinomial2 rather than categorical? I
> don't want to go the categorical route because repeatabilities are
> of particular interest.
>
> I'm much less comfortable with generalized rather than standard
> linear models so I wonder if there are any glaring issues with the
> coding. Any thoughts on what might be going on?
>
> Thanks for any help.
> Ned
>
>
>
> #make function to calculate inverse logit:
> inv.log<- function(x){
> exp(x)/(1+exp(x))}
>
> #specify some start conditions and bookkeeping stuff
>
> #specify data distributions:
> repeatability=0.3 #Change this for other repeatabilities (sampled
> population's repeatability, not sample's realized latent
> repeatability)
> sigma2.e=6 #arbitrary, set repeatability above
> sigma2.a=(repeatability*sigma2.e+repeatability*((pi^2)/3))/(1-repeatability)
> B0=0 #could do different fixed effects with some more work...
>
> #(n=number of individuals, i=number of replications per individual)
> counter=0; n=150; i=4
> ind.data=matrix(0,n*i,3) #make place to store data
> tru.lat=matrix(0,n*i,2) #store sampled underlying scores
>
> for(j in 1:n){
> alpha.j=rnorm(1,0,sigma2.a)
> for(rep in 1:i){
> counter=counter+1
> ind.data[counter,1]=j
> eij=rnorm(1,0,sigma2.e)
> pij=inv.log((B0+alpha.j+eij))
> ind.data[counter,3]=rbinom(1,1,pij)
> ind.data[counter,2]=1-ind.data[counter,3]
> tru.lat[counter,1]=alpha.j
> tru.lat[counter,2]=eij
> }}
>
> ind.data=as.data.frame(ind.data)
> names(ind.data)=c("Ind", "Fail", "Success")
>
> prior1 = list(R = list(V = 1, nu=0.002),
> G = list(G1 = list(V = 1, nu = 0.002)))
> sim.mcmc<-MCMCglmm(cbind(Fail,Success)~1, random= ~Ind,
> family="multinomial2", prior=prior1,
> nitt = 260000, thin = 200, burnin = 60000,
> verbose=FALSE,data=ind.data)
>
> #estimated repeatabilities
> posterior.mode(sim.mcmc$VCV[,1]/(sim.mcmc$VCV[,1]+sim.mcmc$VCV[,2]+((pi^2)/3)))
> HPDinterval(sim.mcmc$VCV[,1]/(sim.mcmc$VCV[,1]+sim.mcmc$VCV[,2]+((pi^2)/3)))
> #effective sample sizes, tells you if model is mixing well
> summary(sim.mcmc)[6]$Gcovariances[4]
> summary(sim.mcmc)[8]$Rcovariances[4]
>
> #sample repeatability, not sure this is quite right
> tru.a2=var(tru.lat[1:(n/i)*i,1])
> tru.e2=var(tru.lat[,2])
> tru.a2/(tru.a2+tru.e2+((pi^2)/3))
>
>
> --
> Ned A. Dochtermann
> Assistant Professor / Department of Biological Sciences
> NORTH DAKOTA STATE UNIVERSITY
> p: 701.231.7353 / f: 701.231.7149 / www.ndsu.edu
>
> https://sites.google.com/site/neddochtermann/
> ned.dochtermann at ndsu.edu
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list