[R-sig-ME] Fwd: Interpretation of lmer output in R
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Thu Feb 24 13:56:44 CET 2011
Dear Julia,
Plotting the mean probability of foraging versus depth will make thing more clear.
Depth <- seq(0, 50, length = 101)
plot(Depth, plogis(-1.6898789+0.0007075*Depth), type = "l")
Note that the difference in probability between Depth = 0 and Depth = 50 is quite small and thus not relevant. So you have a variable (depth) in the model which improves the model, but has no biological relevance. Such thing can happen when you have a lot of data.
I would keep depth in the model, and conclude that is it have a neglectable effect. This is a case were rejecting the STATISICAL null-hypothesis allows to accept the BIOLOGICAL null-hypotheses. Which is stronger than not being able to reject the BIOLOGICAL null-hypothesis.
Best regards,
Thierry
----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
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
> Julia Sommerfeld
> Verzonden: donderdag 24 februari 2011 13:34
> Aan: r-sig-mixed-models at r-project.org
> Onderwerp: [R-sig-ME] Fwd: Interpretation of lmer output in R
>
> ---------- Forwarded message ----------
> From: Julia Sommerfeld <Julia.Sommerfeld at utas.edu.au>
> Date: 2011/2/24
> Subject: Re: [R-sig-ME] Interpretation of lmer output in R
> To: Douglas Bates <bates at stat.wisc.edu>
>
>
> ... it's iluminating the confusion in my brain every day a
> little bit more!
>
> Yes, the confusion comes from my coding 0,1 both parameter. I
> hope you don't mind that I keep asking....
>
> 1. Maybe it helps if I have a look at a different data set
> where only the response variable has a coding of 0,1.
>
> I want to test if foraging (0,1 where 0=means "NO foraging"
> and 1=means "foraging" ) is influenced by water depth (Depth)
> (in meters), i.e. does my model including "Depth" fit better
> than without?
>
> My model:
>
> >fm<-lmer(Foraging~Depth+(1|Bird),family=binomial)
> >fm1<-lmer(Foraging~1+(1|Bird),family=binomial)
> >anova(fm,fm1)
>
> Data:
> Models:
> fd1: Foraging ~ 1 + (1 | Bird)
> fd: Foraging ~ Depth + (1 | Bird)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> fd1 2 1909.0 1920.5 -952.52
> fd 3 1904.8 1922.0 -949.40 6.2385 1 0.0125 *
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>
> > summary(fm)
> Generalized linear mixed model fit by the Laplace approximation
> Formula: Foraging ~ Depth + (1 | Bird)
> AIC BIC logLik deviance
> 1905 1922 -949.4 1899
> Random effects:
> Groups Name Variance Std.Dev.
> Bird (Intercept) 0.29505 0.54319
> Number of obs: 2249, groups: Bird, 23
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -1.6898789 0.1463908 -11.544 <2e-16 ***
> Depth 0.0007075 0.0002884 2.453 0.0142 *
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Correlation of Fixed Effects:
> (Intr)
> Depth 0.244
>
>
> >From this ouput I conclude that the model including "Depth" (fm)
> >explains
> the observed data significanlty better than without Depth (fm1).
>
> 2. But, if I want to know HOW Depth is influencing/related to
> the foraging event (1). I.e., do birds preferably forage in
> shallow or deep waters? How can I get this information out of
> the summary? Would plogis() still be the right approach in this case?
>
> plogis(1.6898789)
>
> plogis(1.6898789+0.0007075*second_column?)
>
> above plogis based on foraging=1, but what about "Depth"?
>
> I think I'm still getting someting wrong here...
>
> Cheers,
>
> Julia
>
>
>
>
>
> 2011/2/23 Douglas Bates <bates at stat.wisc.edu>
>
> > On Wed, Feb 23, 2011 at 6:17 AM, Julia Sommerfeld
> > <Julia.Sommerfeld at utas.edu.au> wrote:
> > > Dear Douglas and list,
> >
> > > Again, thanks for the thorough answers.
> >
> > > Regarding my last email: I'm mixing up my own data! Sorry
> for that.
> >
> > > 1. plogis(x)
> >
> > > SameSite=1 means that birds did use the SAME nest site,
> hence: they
> > > show site-fidelity.
> >
> > > SameSite=0 means that birds did CHANGE nest site, hence: no
> > site-fidelity.
> >
> > > This means that site-fidelity of "bird" who failed breeding
> > > corresponds
> > to a
> > > probability of 44% and that bird who was successful
> corresponds to a
> > > probability of 72% (see below).
> >
> > Yes. When you code the response as 0 and 1 then the model provides
> > the probability that the response will be a 1.
> >
> > > Thus, with "plogis" of the Intercept and/or Intercept+Fixed Effect
> > Estimate
> > > I can calculate the probability of my observed event
> (R-help: plogis
> > gives
> > > the distribution function; being the inverse logit-function)?
> >
> > Yes. It happens that the logit function (see
> > http://en.wikipedia.org/wiki/Logit) is the quantile
> function for the
> > standard logistic distribution (the R function qlogis) and
> its inverse
> > is the cumulative probability function for the standard
> logistic (the
> > R function plogis). Because I know that the qlogis and plogis
> > functions are carefully written to handle difficult cases
> gracefully I
> > use them instead of writing out the expressions explicitly.
> >
> > > That is: plogis of the intercept estimate is based upon
> BreedSuc1=0,
> > > SameSite=1 and plogis of the intercept+BreedSuc1 Estimate
> is based
> > > upon BreedSuc1=1, SameSite=1?
> >
> > Yes. Another way of looking at it would be to extract the model
> > matrix for your fitted model, using the function model.matrix. You
> > will see that the first column is all 1's and the second
> column will
> > be 0's and 1's according to that bird's breeding success.
> Technically
> > the fixed-effects part of the predicted log-odds for a
> particular bird
> > is (Intercept) * first_column + BreedSuc1 * second_column
> but, because
> > every row is either 1 0 or 1 1, the only possible values of the
> > predicted log-odds are (Intercept) and (Intercept) + BreedSuc1
> >
> > > And the "result" will always depend upon my coding (BreedSuc1=1
> > > means successful, SameSite=1 means site-fidelity). This
> is still a
> > > bit confusing...
> >
> > I would consider the two factors separately. Regarding the
> SameSite
> > encoding just think of it as you are eventually modeling the
> > probability of an "event" for which there are two possible outcomes.
> > You can call them Failure/Success or Off/On or No/Yes or
> whatever you
> > want but when you get to a numerical coding as 0/1 the
> model will be
> > for the probability of getting a 1.
> >
> > As for the breeding success factor, you again need to
> characterize it
> > according to some numerical encoding to be able to create the model
> > matrix. You happened to choose a 0/1 encoding so in the model the
> > coefficient for that term is added to the intercept when
> there is a 1
> > for that factor and not added when there is a 0.
> >
> > > My question was:
> > > What if both codings change? I.e. BreedSuc1=0 means
> successful and
> > > SameSite=0 means site-fidelity? Would it then still be
> 1-p instead of p?
> >
> > If change the encoding of the response, then your coefficients will
> > all flip signs. What was previously a log-odds of site
> fidelity of,
> > say, 1.2 will now become a log-odds of site switching of -1.2. You
> > can check that
> >
> > plogis(-x) = 1 - plogis(x)
> >
> > If you change the encoding of BreedSuc then the (Intercept)
> > coefficient will be the log-odds of a '1' outcome (whether
> that means
> > site fidelity or site switching) for a bird who was successful and
> > (Intercept) + BreedFail1 will be the log-odds of a 1 outcome for a
> > bird who was successful.
> >
> > I think part of the confusion for you is that you are trying to
> > combine the encoding of site-fidelity and breeding success. It is
> > best to keep them distinct because site-fidelity is the
> response and
> > breeding success is a covariate.
> >
> > > 2. Z-value:
> > >
> > > I assume that z-value is the same as z-score? I've found the
> > > following
> > link
> > > (
> >
> http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/spatia
> > l_statistics_toolbox/what_is_a_z_score_what_is_a_p_value.htm
> > ).
> > >
> >
> > > There, the z-score is defined as:
> > >
> > > "The Z score is a test of statistical significance that helps you
> > > decide whether or not to reject the null hypothesis".
> > >
> > > "Z scores are measures of standard deviation. For
> example, if a tool
> > returns
> > > a Z score of +2.5 it is interpreted as "+2.5 standard deviations
> > > away
> > from
> > > the mean". P-values are probabilities. Both statistics are
> > > associated
> > with
> > > the standard normal distribution. This distribution
> relates standard
> > > deviations with probabilities and allows significance and
> confidence
> > > to
> > be
> > > attached to Z scores and p-values".
> >
> > > Is this correct?
> >
> > Well, not really. The logic behind statistical hypothesis
> testing is
> > subtle and frequently misunderstood. When you try to simplify the
> > explanation, as is done above, it often ends up not quite accurate.
> >
> > The way I present it to students is that we have two
> competing claims,
> > which for us will boil down to a simpler model being
> adequate (called
> > the "null hypothesis" as in "nothing going on here, folks") or the
> > simpler model is not adequate and we need to extend it. The second
> > claim is called the alternative hypothesis. In your case the null
> > hypothesis is that breeding success doesn't influence the
> probability
> > of site fidelity and the alternative hypothesis is that it
> does. You
> > want to establish the alternative hypothesis but, because of the
> > variability in the data, you can't "prove" it directly.
> You have to
> > use an indirect method, which is to assume that it is not the case
> > (i.e. the null hypothesis is true) and demonstrate that it
> is unlikely
> > to obtain the data that you did, or something even more unusual, if
> > the null hypothesis were true.
> >
> > So that is the subtlety. You can't prove your point
> directly so you
> > set up the "straw man" argument (the null hypothesis) and try to
> > refute it. But nothing is certain. You can't prove that it is
> > impossible to get the data that you did if breeding success did not
> > influence site fidelity. The best you can do is say it is
> extremely
> > unlikely.
> >
> > Now we need to formalize this argument by settling on a "test
> > statistic", which is a quantity that can be calculated from the
> > observed data. Along with the test statistic we need a reference
> > distribution which is the distribution of the test
> statistic when the
> > null hypothesis is true. Formally, we are supposed to set
> up all the
> > rules about the test statistic and the reference
> distribution before
> > we see the data - sort of like the "pistols at dawn, march
> ten paces
> > then turn and fire" rules of a duel. In practice we
> sometimes look at
> > the data before we decide exactly what hypotheses we will
> test. This
> > is the "data snooping" that Ben wrote about.
> >
> > Anyway, when we have set up the rules and we have the data
> we evaluate
> > the test statistic and then determine the probability of
> seeing that
> > value or something even more extreme when the null
> hypothesis is true.
> > This is what we call the p-value. If that probability is
> very small
> > then either we encountered a rare case or the null hypothesis is
> > false. That's why you, as the researcher, want the p-value to be
> > small, when you are trying to establish the alternative.
> >
> > (As an aside, this is where most people lose track of the argument.
> > They think the p-value is related to the probability of one of the
> > hypotheses being true. That is not the case. The p-value is the
> > probability of seeing the data that we did, or something
> more unusual,
> > if the null hypothesis - the "no change" assumption - were true.)
> >
> > Now we get to the point of how do we calculate the p-value.
> We need
> > to come up with a way of characterizing "unusual" data.
> What I prefer
> > to do is to see how good the fit is without allowing for
> differences
> > due to breeding success and how good it is when allowing for
> > differences due to breeding success. The measure of "how
> good is the
> > fit" is called the deviance. So we fit the model with the
> BreedSuc1
> > term and without it and compare the deviance values. The reference
> > distribution for this difference is generally taken to be a
> > chi-squared distribution. Technically all we can say is
> that when the
> > sample sizes get very large, the difference in the deviance values
> > tends to a chi-squared distribution. But I would claim
> that it is a
> > reasonable way of judging the difference in quality of fit
> for finite
> > samples too.
> >
> > An alternative approach to judging if a coefficient is
> "significantly
> > different" from zero is to take the value of the coefficient and
> > divide by its standard error. If everything is behaving well, this
> > ratio would have a standard Gaussian (or "normal")
> distribution when
> > the null hypothesis is true. We write a standard normal random
> > variable as Z so this ratio is called a z-statistic.
> >
> > For me the problem with using a z-statistic is that you are
> comparing
> > two models (with and without the term of interest) by
> fitting only one
> > of them and then extrapolating to decide what the other fit
> should be
> > like. This was a necessary short-cut when fitting a model might
> > involve weeks of hand calculation. But if fitting a model involves
> > only a few seconds of computer time, why not fit both and do a
> > comparison of the actual difference in the quality of the fits?
> >
> > So, why should we even both printing the z-statistics? I consider
> > them to be a guide but if I actually want to do the hypothesis test
> > involving a simple model versus a more complex model I fit both
> > models.
> >
> > I hope this is more illuminating than confusing.
> >
>
>
>
> --
> Julia Sommerfeld - PhD Candidate
> Institute for Marine and Antarctic Studies University of
> Tasmania Private Bag 129, Hobart TAS 7001
>
> Phone: +61 458 247 348
> Email: julia.somma at gmx.de
> Julia.Sommerfeld at utas.edu.au
>
> [[alternative HTML version deleted]]
>
>
More information about the R-sig-mixed-models
mailing list