[R-sig-ME] Model specification/family for a continuous/proportional, response with many zeros
Michael Lawson
mrm|500 @end|ng |rom york@@c@uk
Thu May 20 11:34:16 CEST 2021
Dear Thierry,
Sorry, I didn't realise I was sending everything in HTML mode. I have
enabled plain text mode so hopefully it has worked this time (code
below). Thanks for the suggestion and link - comparing the
observed-expected density looks a lot better than just comparing the
means like I did. I'll have a go!
Thanks,
Mike
bins <- seq(from = 0, to = 1, by = 0.01)
sim_seq <- c()
for(i in bins){
a <- mean(unlist(lapply(sim_beta_glmm, function(x){ length(which(x>i
& x <= i+0.01))/length(x)}), use.names = FALSE))
sim_seq <- c(sim_seq,a)
}
real_seq <- c()
for(i in bins){
a <- print(mean(unlist(lapply(glmm_zone_data$prop_time, function(x){
length(which(x>i & x <= i+0.01))/length(x)}), use.names = FALSE)))
real_seq <- c(real_seq,a)
}
barplot(rbind(sim_seq,real_seq),col=c("green","red"),beside = TRUE,
legend.text = c("simulated","real"))
cbind(bins,sim_seq,real_seq)
bins sim_seq real_seq
[1,] 0.00 0.2199795081967 0.232240437
[2,] 0.01 0.1021612021858 0.166666667
[3,] 0.02 0.0750964480874 0.103825137
[4,] 0.03 0.0592128415301 0.073770492
[5,] 0.04 0.0478316939891 0.038251366
[6,] 0.05 0.0397267759563 0.030054645
[7,] 0.06 0.0331961748634 0.016393443
[8,] 0.07 0.0279322404372 0.013661202
[9,] 0.08 0.0238013661202 0.024590164
[10,] 0.09 0.0201101092896 0.010928962
On Thu, 20 May 2021 at 10:08, Thierry Onkelinx <thierry.onkelinx using inbo.be> wrote:
>
> Dear Mike,
>
> Your code is completely unreadable. Did you send it in plain text? Sending HTML email tends to mangle code.
>
> I check the distribution by running a large number of simulations to get the range of potential values. Then I divide the observed density (from the data) by the expected density (from the simulations). See https://inlatools.netlify.app/articles/distribution.html#distribution-checks-1 for an example. This is based on INLA models and doesn't handle beta distributions (yet). But you can get an idea of how you can do this yourself.
>
> Best regards,
>
> Thierry
>
> ir. Thierry Onkelinx
> Statisticus / Statistician
>
> Vlaamse Overheid / Government of Flanders
> INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
> Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
> thierry.onkelinx using inbo.be
> Havenlaan 88 bus 73, 1000 Brussel
> 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
> ///////////////////////////////////////////////////////////////////////////////////////////
>
>
>
>
> Op do 20 mei 2021 om 10:46 schreef Michael Lawson via R-sig-mixed-models <r-sig-mixed-models using r-project.org>:
>>
>> Hi again,
>>
>> I wanted to look a little more at the rest of the simulated values
>> (non-zeros), so I took the mean proportion of simulated values in set bins
>> of 0.01 and compared this to the real values.
>>
>>
>> *bins <- seq(from = 0, to = 1, by = 0.01)*
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> *sim_seq <- c()for(i in bins){ a <- mean(unlist(lapply(sim_beta_glmm,
>> function(x){ length(which(x>i & x <= i+0.01))/length(x)}), use.names =
>> FALSE)) sim_seq <- c(sim_seq,a)}real_seq <- c()for(i in bins){ a <-
>> print(mean(unlist(lapply(glmm_zone_data$prop_time, function(x){
>> length(which(x>i & x <= i+0.01))/length(x)}), use.names = FALSE)))
>> real_seq <- c(real_seq,a)}*
>>
>>
>> *barplot(rbind(sim_seq,real_seq),col=c("green","red"),beside = TRUE,
>> legend.text = c("simulated","real"))*
>>
>>
>> *cbind(bins,sim_seq,real_seq)*
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> *bins sim_seq real_seq [1,] 0.00 0.2199795081967 0.232240437
>> [2,] 0.01 0.1021612021858 0.166666667 [3,] 0.02 0.0750964480874
>> 0.103825137 [4,] 0.03 0.0592128415301 0.073770492 [5,] 0.04
>> 0.0478316939891 0.038251366 [6,] 0.05 0.0397267759563 0.030054645 [7,]
>> 0.06 0.0331961748634 0.016393443 [8,] 0.07 0.0279322404372 0.013661202
>> [9,] 0.08 0.0238013661202 0.024590164 [10,] 0.09 0.0201101092896
>> 0.010928962*
>>
>> As you can see in the output above (and the plot attached), the real data
>> is quite a bit more right-skewed compared to the simulated values. Does
>> this look like a good enough fit or will I have to try a different model?
>>
>> Many thanks,
>> Mike
>>
>> On Wed, 19 May 2021 at 12:38, Michael Lawson <mrml500 using york.ac.uk> wrote:
>>
>> > Dear Alain,
>> >
>> > Thank you for the suggestion. I tried simulating the 0's as you suggested
>> > in the following way. The output is reassuring - with the actual number of
>> > zeros in the middle of the distribution.
>> >
>> > sim_beta_glmm <- simulate(beta_mod, nsim = 10000)
>> > sim_zeros <- unlist(lapply(sim_beta_glmm, function(x){
>> > length(which(x==0))/length(x)}), use.names = FALSE)
>> > hist(sim_zeros, breaks = c(100))
>> > abline(v = plogis(-1.1804), col = "red")
>> >
>> > All the best,
>> > Mike
>> >
>> >
>> > On Wed, 19 May 2021 at 11:19, Highland Statistics Ltd via
>> > R-sig-mixed-models <r-sig-mixed-models using r-project.org> wrote:
>> >
>> >>
>> >> > Today's Topics:
>> >> >
>> >> > 1. Re: Model specification/family for a continuous/proportional
>> >> > response with many zeros (Michael Lawson)
>> >> >
>> >> > ----------------------------------------------------------------------
>> >> >
>> >> > Message: 1
>> >> > Date: Wed, 19 May 2021 09:31:46 +0100
>> >> > From: Michael Lawson <mrml500 using york.ac.uk>
>> >> > To: Thierry Onkelinx <thierry.onkelinx using inbo.be>
>> >> > Cc: r-sig-mixed-models <r-sig-mixed-models using r-project.org>
>> >> > Subject: Re: [R-sig-ME] Model specification/family for a
>> >> > continuous/proportional response with many zeros
>> >> > Message-ID:
>> >> > <
>> >> CACtWw1GAnDJ1d8CWh_DODP6jPqfsyptc-RV0gw0_5z9bWEwubw using mail.gmail.com>
>> >> > Content-Type: text/plain; charset="utf-8"
>> >> >
>> >> > Hi Thierry,
>> >> >
>> >> > I understood after your previous message that a high dispersion
>> >> parameter
>> >> > in beta models does not signify overdispersion. Both answers addressed
>> >> this
>> >> > misconception in those two links I provided.
>> >> >
>> >> > I suppose the real answer I am now looking for is - how do I assess the
>> >> > validity and fit of the model? Are there specifics or assumptions I must
>> >> > take into account with zero-inflated beta glmms?
>> >>
>> >> Just simulate 1000 data sets from your model, and see whether the
>> >> simulated data sets are comparable to your original data. One of the
>> >> things you can look at is whether the simulated data sets contain
>> >> similar number of zeros as in the observed data. I'm not sure whether
>> >> the DHARMA package can do this for you. If not..it is easy to program.
>> >>
>> >> See also Figure 8 in (sorry for self-citing):
>> >>
>> >>
>> >> https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12577
>> >>
>> >> Alain
>> >>
>> >>
>> >> > Thanks for your help,
>> >> > Mike
>> >> >
>> >> > On Tue, 18 May 2021 at 13:57, Thierry Onkelinx <
>> >> thierry.onkelinx using inbo.be>
>> >> > wrote:
>> >> >
>> >> >> Dear Mike,
>> >> >>
>> >> >> I think you misread my reply. I never stated that there's something
>> >> wrong
>> >> >> with the model. The only "problem" I highlighted was your misconception
>> >> >> about the "high overdispersion". In this case, a high parameter value
>> >> >> indicates a low variance, which is what we mostly want to see.
>> >> >>
>> >> >> Best regards,
>> >> >>
>> >> >> ir. Thierry Onkelinx
>> >> >> Statisticus / Statistician
>> >> >>
>> >> >> Vlaamse Overheid / Government of Flanders
>> >> >> INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
>> >> AND
>> >> >> FOREST
>> >> >> Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
>> >> >> thierry.onkelinx using inbo.be
>> >> >> Havenlaan 88 bus 73, 1000 Brussel
>> >> >> 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
>> >> >>
>> >> >>
>> >> ///////////////////////////////////////////////////////////////////////////////////////////
>> >> >>
>> >> >> <https://www.inbo.be>
>> >> >>
>> >> >>
>> >> >> Op di 18 mei 2021 om 12:49 schreef Michael Lawson <mrml500 using york.ac.uk
>> >> >:
>> >> >>
>> >> >>> Dear Thierry,
>> >> >>>
>> >> >>> Thanks for the help. So if the dispersion parameter in this model
>> >> doesn't
>> >> >>> fit with the beta distribution, are there any alternative approaches
>> >> I can
>> >> >>> use?
>> >> >>>
>> >> >>> I can't seem to find much information on this elsewhere other than
>> >> these
>> >> >>> two threads:
>> >> >>> https://stats.stackexchange.com/a/451453/233414
>> >> >>> https://stats.stackexchange.com/a/466951/233414
>> >> >>>
>> >> >>> All the best,
>> >> >>> Mike
>> >> >>>
>> >> >>> On Tue, 18 May 2021 at 08:12, Thierry Onkelinx <
>> >> thierry.onkelinx using inbo.be>
>> >> >>> wrote:
>> >> >>>
>> >> >>>> Dear Mike,
>> >> >>>>
>> >> >>>> The zero-inflation is specified on the logit scale. plogis(-1.18) =
>> >> >>>> 0.235 23.5% zero seems reasonable when reading your story. (Didn't
>> >> look at
>> >> >>>> the data).
>> >> >>>>
>> >> >>>> You need to look at the definition for the "over"dispersion
>> >> parameter.
>> >> >>>> For a beta distribution is \phi with Var(y) = \mu * (1 - \mu) /
>> >> (\phi + 1)
>> >> >>>> (see ?glmmTMB::beta_family) Hence a large value of \phi implies a low
>> >> >>>> variance.
>> >> >>>>
>> >> >>>> Best regards,
>> >> >>>>
>> >> >>>> ir. Thierry Onkelinx
>> >> >>>> Statisticus / Statistician
>> >> >>>>
>> >> >>>> Vlaamse Overheid / Government of Flanders
>> >> >>>> INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR
>> >> NATURE
>> >> >>>> AND FOREST
>> >> >>>> Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
>> >> >>>> thierry.onkelinx using inbo.be
>> >> >>>> Havenlaan 88 bus 73, 1000 Brussel
>> >> >>>> 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
>> >> >>>>
>> >> >>>>
>> >> ///////////////////////////////////////////////////////////////////////////////////////////
>> >> >>>>
>> >> >>>> <https://www.inbo.be>
>> >> >>>>
>> >> >>>>
>> >> >>>> Op ma 17 mei 2021 om 15:45 schreef Michael Lawson <
>> >> mrml500 using york.ac.uk>:
>> >> >>>>
>> >> >>>>> Hi Thierry,
>> >> >>>>>
>> >> >>>>> Thank you for your advice and speedy response.
>> >> >>>>>
>> >> >>>>> Most of the data is closer to the lower bound (0). e.g. the mean
>> >> time
>> >> >>>>> for group A in zone A = 15.1 seconds and group A in zone B = 3.8
>> >> seconds.
>> >> >>>>> However there are a very small number of outliers near the upper
>> >> bound, the
>> >> >>>>> largest being 294 out of the 300 seconds (see the attached file if
>> >> you want
>> >> >>>>> to look at the data).
>> >> >>>>>
>> >> >>>>> I have taken a stab at running a Zero-inflated Beta GLMM using
>> >> glmmTMB
>> >> >>>>> in R like so:
>> >> >>>>>
>> >> >>>>> betta_mod <- glmmTMB(prop_time ~ group*zone + (1|id),
>> >> >>>>> family = beta_family(),
>> >> >>>>> ziformula=~1,
>> >> >>>>> data = glmm_zone_data)
>> >> >>>>>
>> >> >>>>> summary(beta_mod)
>> >> >>>>>
>> >> >>>>> *Family: beta ( logit )*
>> >> >>>>>
>> >> >>>>>
>> >> >>>>>
>> >> >>>>>
>> >> >>>>>
>> >> >>>>>
>> >> >>>>> *Formula: prop_time ~ group * zone + (1 | id)Zero
>> >> inflation:
>> >> >>>>> ~1Data: glmm_zone_data AIC BIC logLik deviance
>> >> >>>>> df.resid -763.6 -736.3 388.8 -777.6 359Random
>> >> >>>>> effects:Conditional model: Groups Name Variance Std.Dev. id
>> >> >>>>> (Intercept) 2.386e-09 4.885e-05Number of obs: 366, groups: id,
>> >> >>>>> 14Overdispersion parameter for beta family (): 13.1Conditional
>> >> model:
>> >> >>>>> Estimate Std. Error z value Pr(>|z|) (Intercept)
>> >> >>>>> -2.7685 0.1031 -26.844 < 2e-16 ***groupB -0.4455
>> >> >>>>> 0.1498 -2.975 0.002932 **zonezone_B -0.4179 0.1524
>> >> -2.741
>> >> >>>>> 0.006124 **groupB:zonezone_B 0.8443 0.2190 3.855 0.000116
>> >> >>>>> ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
>> >> >>>>> 1Zero-inflation model: Estimate Std. Error z value
>> >> Pr(>|z|)
>> >> >>>>> (Intercept) -1.1804 0.1233 -9.575 <2e-16 ***---Signif.
>> >> codes: 0
>> >> >>>>> ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1*
>> >> >>>>>
>> >> >>>>> Does this look like the correct way of specifying the model? I am a
>> >> >>>>> little confused about specifying and interpreting the zero-inflation
>> >> >>>>> component - I have only just begun reading about this.
>> >> >>>>>
>> >> >>>>> I noticed that the dispersion parameter is quite high at 13.1. I'm
>> >> not
>> >> >>>>> sure if this matters for beta models?. I tried running DHARMa
>> >> >>>>> simulateResiduals on the model output and got significant
>> >> deviations in the
>> >> >>>>> dispersion (<2.2e-16) and KS tests. e.g.
>> >> >>>>>
>> >> >>>>> DHARMa::testDispersion(beta_mod)
>> >> >>>>>
>> >> >>>>> *DHARMa nonparametric dispersion test via sd of residuals fitted vs.
>> >> >>>>> simulated*
>> >> >>>>>
>> >> >>>>> *data: simulationOutput*
>> >> >>>>> *ratioObsSim = 1.3612, p-value < 2.2e-16*
>> >> >>>>> *alternative hypothesis: two.sided*
>> >> >>>>>
>> >> >>>>>
>> >> >>>>>
>> >> >>>>> Many thanks,
>> >> >>>>> Mike
>> >> >>>>>
>> >> >>>>> On Mon, 17 May 2021 at 13:22, Thierry Onkelinx <
>> >> >>>>> thierry.onkelinx using inbo.be> wrote:
>> >> >>>>>
>> >> >>>>>> Dear Michael,
>> >> >>>>>>
>> >> >>>>>> Your data has bounds (lower bound at 0 and upper bound at 300) and
>> >> you
>> >> >>>>>> have a lot of data close to a boundary. In such a case, a
>> >> continuous
>> >> >>>>>> distribution which ignores those bound is not a good idea. If the
>> >> time
>> >> >>>>>> spent outside of both zones is limited, then a long time in zone A
>> >> excludes
>> >> >>>>>> a long time in zone B by definition. Then I'd look towards a
>> >> multinomial
>> >> >>>>>> distribution. If the time spent outside both zones is dominant,
>> >> then you
>> >> >>>>>> can use a zero-inflated beta as you suggested. A zero-inflated
>> >> gamma might
>> >> >>>>>> be OK if the data is not too close to the upper boundary. If you
>> >> are
>> >> >>>>>> considering zero-inflated beta vs zero-inflated gamma, then you
>> >> should
>> >> >>>>>> choose zero-inflated beta IMHO.
>> >> >>>>>>
>> >> >>>>>> Best regards,
>> >> >>>>>>
>> >> >>>>>> ir. Thierry Onkelinx
>> >> >>>>>> Statisticus / Statistician
>> >> >>>>>>
>> >> >>>>>> Vlaamse Overheid / Government of Flanders
>> >> >>>>>> INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR
>> >> NATURE
>> >> >>>>>> AND FOREST
>> >> >>>>>> Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality
>> >> Assurance
>> >> >>>>>> thierry.onkelinx using inbo.be
>> >> >>>>>> Havenlaan 88 bus 73, 1000 Brussel
>> >> >>>>>> 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
>> >> >>>>>>
>> >> >>>>>>
>> >> ///////////////////////////////////////////////////////////////////////////////////////////
>> >> >>>>>>
>> >> >>>>>> <https://www.inbo.be>
>> >> >>>>>>
>> >> >>>>>>
>> >> >>>>>> Op ma 17 mei 2021 om 13:52 schreef Michael Lawson via
>> >> >>>>>> R-sig-mixed-models <r-sig-mixed-models using r-project.org>:
>> >> >>>>>>
>> >> >>>>>>> Hello,
>> >> >>>>>>>
>> >> >>>>>>> I am new to GLMMs and have a dataset where I have two distinct
>> >> groups
>> >> >>>>>>> (A
>> >> >>>>>>> and B) of 7 individuals each. The data consists of repeated
>> >> >>>>>>> measurements of
>> >> >>>>>>> each individual where the amount of time spent at either zone_A or
>> >> >>>>>>> zone_B
>> >> >>>>>>> is recorded (out of a total time of 300s/observation period). For
>> >> >>>>>>> most of
>> >> >>>>>>> the time period the individuals are in neither zone.
>> >> >>>>>>>
>> >> >>>>>>> I want to test if group A and group B spend more time in zone A
>> >> >>>>>>> compared to
>> >> >>>>>>> zone B (and vice versa).
>> >> >>>>>>>
>> >> >>>>>>> Speaking to someone else, they said I should use a Binomial GLMM
>> >> using
>> >> >>>>>>> cbind. i.e.
>> >> >>>>>>> cbind(time_at_zone_A, time_at_zone_B) ~ group + (1| id).
>> >> >>>>>>>
>> >> >>>>>>> However, the response variable is continuous (albeit with an upper
>> >> >>>>>>> bound of
>> >> >>>>>>> 300 seconds per observation period), so I'm not sure if this is
>> >> >>>>>>> appropriate?
>> >> >>>>>>>
>> >> >>>>>>> Should I convert the response into a proportion and use something
>> >> >>>>>>> like a
>> >> >>>>>>> Beta GLMM or else use a continuous (Gamma) GLMM? e.g. something
>> >> like:
>> >> >>>>>>> prop_time ~ zone*group + (1|id)
>> >> >>>>>>>
>> >> >>>>>>> The data is quite heavily right-skewed and contains a lot of 0's,
>> >> so
>> >> >>>>>>> reading around it also looks like I may need to convert these
>> >> into a
>> >> >>>>>>> zero-inflated/hurdle model?
>> >> >>>>>>>
>> >> >>>>>>> Thank you for any suggestions,
>> >> >>>>>>> Mike
>> >> >>>>>>>
>> >> >>>>>>> [[alternative HTML version deleted]]
>> >> >>>>>>>
>> >> >>>>>>> _______________________________________________
>> >> >>>>>>> R-sig-mixed-models using r-project.org mailing list
>> >> >>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >> >>>>>>>
>> >> > [[alternative HTML version deleted]]
>> >> >
>> >> >
>> >> >
>> >> >
>> >> > ------------------------------
>> >> >
>> >> > Subject: Digest Footer
>> >> >
>> >> > _______________________________________________
>> >> > R-sig-mixed-models mailing list
>> >> > R-sig-mixed-models using r-project.org
>> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >> >
>> >> >
>> >> > ------------------------------
>> >> >
>> >> > End of R-sig-mixed-models Digest, Vol 173, Issue 19
>> >> > ***************************************************
>> >>
>> >> --
>> >> Dr. Alain F. Zuur
>> >> Highland Statistics Ltd.
>> >> 9 St Clair Wynd
>> >> AB41 6DZ Newburgh, UK
>> >> Email: highstat using highstat.com
>> >> URL: www.highstat.com
>> >>
>> >> _______________________________________________
>> >> R-sig-mixed-models using r-project.org mailing list
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >>
>> >
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list