[prev in list] [next in list] [prev in thread] [next in thread]
List: r-sig-mixed-models
Subject: Re: [R-sig-ME] Interpretation of lmer output in R
From: Julia Sommerfeld <Julia.Sommerfeld () utas ! edu ! au>
Date: 2011-02-28 8:17:15
Message-ID: AANLkTi=ChOR626hedN=UoGHBskjY1MX_EaKsiGA+2v1e () mail ! gmail ! com
[Download RAW message or body]
Dear Douglas and list,
Again thank you for the answers. I think I'm getting the "story" now - I was
definitely missing some background knowledge required to understand the
transition from logistic regression to the output of my lmer.
Keen to learn more...
Cheers
Julia
2011/2/24 Douglas Bates <bates@stat.wisc.edu>
> On Thu, Feb 24, 2011 at 6:33 AM, Julia Sommerfeld
> <Julia.Sommerfeld@utas.edu.au> wrote:
> > ... 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 *
>
> By the way, it is interesting to see that the p-values for the
> likelihood ratio test is slightly different from the p-value for the
> z-statistic, which is based on fitting only one model. In this
> circumstance I would be more inclined to believe the result of the
> likelihood ratio test.
>
> > > 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 *
>
> > 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).
>
> Yes.
>
> > 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?
>
> Yes.
>
> > plogis(1.6898789)
> >
> > plogis(1.6898789+0.0007075*second_column?)
> >
> > above plogis based on foraging=1, but what about "Depth"?
>
> The second_column is Depth. Check with
>
> model.matrix(fm)
>
> or, preferably,
>
> head(model.matrix(fm), n = 15)
>
> because you have 2249 rows in that model matrix.
>
> > I think I'm still getting someting wrong here...
> >
> > Cheers,
> >
> > Julia
> >
> >
> >
> >
> >
> > 2011/2/23 Douglas Bates <bates@stat.wisc.edu>
> > >
> > > On Wed, Feb 23, 2011 at 6:17 AM, Julia Sommerfeld
> > > <Julia.Sommerfeld@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/spatial_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@gmx.de
Julia.Sommerfeld@utas.edu.au
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic