A couple of additional points:<br><br>(1) mixed logit models like yours are not fit using REML. The default fit uses Laplace approximation, which apparently is related to, but not the same as REML.<br><br>(2) the insignificance is, of course, due to the large standard error, but the large standard error is due to the distribution of your data -- the fact that you have not new-deaccented data points.<br>
<br>(3) you can recode Info to Info == "given", which makes it easy to show that given is different from the other two levels (consistent with much other previous research). But when you try to show that new is different from accessible, there is barely any evidence for that (as also shown, by splittin the data into a subset without "given" cases).<br>
<br>Florian<br><br><div class="gmail_quote">On Mon, Dec 22, 2008 at 4:30 PM, Roger Levy <span dir="ltr"><<a href="mailto:rlevy@ucsd.edu">rlevy@ucsd.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Laura de Ruiter wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div class="Ih2E3d">
Dear R-users and -experts,<br>
<br>
I am performing a rather simple analysis on a small data set (pasted below this email) and keep getting a to me inexplicable result. Perhaps I am missing something here - it would be great if someone could point out to me what I am doing wrong.<br>
<br>
I want to test whether the factor "Info" (which has three levels: "new", "given", "accessible") is a significant predictor for the binary variable "DeaccYN". The random factor is "Subject". The distribution of the data looks as follows:<br>
<br>
----------------------------------------------------------------------------- <br></div><div><div></div><div class="Wj3C7c">
Info<br>
DeaccYN given new accessible<br>
no 25 42 21<br>
yes 11 0 1<br>
------------------------------------------------------------------------------ <br>
<br>
This is the model:<br>
<br>
---------------------------------------------------------------------------------------------------------- <br>
deacc.lmer = lmer (DeaccYN ~ Info + (1|Subject), data = dat, family = "binomial")<br>
----------------------------------------------------------------------------------------------------------------- <br>
<br>
However, given the distribution above, this outcome seems rather weird to me:<br>
<br>
--------------------------------------------------------------------------------------------------------- <br>
summary (deacc.lmer)<br>
Generalized linear mixed model fit using Laplace<br>
Formula: DeaccYN ~ Info + (1 | Subject)<br>
Data: dat<br>
Family: binomial(logit link)<br>
AIC BIC logLik deviance<br>
60.4 70.82 -26.2 52.4<br>
Random effects:<br>
Groups Name Variance Std.Dev.<br>
Subject (Intercept) 0.18797 0.43356<br>
number of obs: 100, groups: Subject, 21<br>
<br>
Estimated scale (compare to 1 ) 0.7316067<br>
<br>
Fixed effects:<br>
Estimate Std. Error z value Pr(>|z|)<br>
(Intercept) -0.8635 0.3795 -2.2754 0.0229 *<br>
Infonew -18.7451 2764.2445 -0.0068 0.9946<br>
Infoaccessible -2.2496 1.1186 -2.0110 0.0443 *<br>
---<br>
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1<br>
<br>
</div></div></blockquote>
> [...]<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
---------------------------------------------------------------------------------------------------- <br><div class="Ih2E3d">
<br>
Why should the difference between 25/11 and 21/1 be significant, but the difference between 25/11 and 42/0 not? Very odd to me seems the standard error of 2764!<br>
<br>
</div></blockquote>
> [...]<div class="Ih2E3d"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
I was wondering: Is it perhaps a problem for the model that there are no cases in the DeaccYN == "yes" category for Info == "given"? And if this <br>
</blockquote></div>
^^^^^<br>
I believe you mean "new" here.<div class="Ih2E3d"><br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
is the case, why?<br>
Am I overlooking something here?<br>
</blockquote>
<br></div>
Dear Laura,<br>
<br>
Independently of the issue that Florian is raising...you are right that the lack of is a problem for the model (to be precise, it's a problem for estimating the significance of the parameter estimate using the z value). The z value, which is the basis of the significance level reported in the lmer summary, is based on the Wald statistic, which is<br>
<br>
z = \hat{b} / StdError(\hat{b})<br>
<br>
where \hat{b} is the parameter estimate. However, for parameter estimates of large magnitudes (as is the case for "new" here), the standard error is inflated, which leads to a small wald statistic. This problem is discussed in a number of places, including Agresti (1996, 2002); here's a decent mention of it online:<br>
<br>
<a href="http://userwww.sfsu.edu/%7Eefc/classes/biol710/logistic/logisticreg.htm" target="_blank">http://userwww.sfsu.edu/~efc/classes/biol710/logistic/logisticreg.htm</a><br>
<br>
You can see this quite clearly if you add just one "yes/given" example at the end of your data frame, for example:<br>
<br>
101 132 new yes<br>
<br>
Then the model gives (being sure to correctly specify the levels of "Info"):<br>
<br>
> summary(deacc.lmer)<br>
Generalized linear mixed model fit by the Laplace approximation<br>
Formula: Deacc ~ Info + (1 | Subject)<br>
Data: dat<br>
AIC BIC logLik deviance<br>
69.78 80.24 -30.89 61.78<div class="Ih2E3d"><br>
Random effects:<br>
Groups Name Variance Std.Dev.<br></div>
Subject (Intercept) 0.34332 0.58594<br>
Number of obs: 101, groups: Subject, 21<div class="Ih2E3d"><br>
<br>
Fixed effects:<br>
Estimate Std. Error z value Pr(>|z|)<br></div>
(Intercept) -0.8996 0.3936 -2.286 0.02226 *<br>
Iaccessible -2.2808 1.1460 -1.990 0.04656 *<br>
Inew -3.0050 1.1395 -2.637 0.00836 **<br>
<br>
<br>
The next question, of course, is how to deal with the dataset you actually have. If you weren't using a mixed-effects model, you could use a likelihood-ratio test by comparing your full model with a simpler model; the likelihood-ratio test isn't susceptible to problems with large parameter estimates the way the Wald test is. For example:<br>
<br>
> deacc.glm = glm (Deacc ~ Info, data = dat, family = "binomial")<br>
> deacc.glm1 = glm (Deacc ~ Info1, data = dat, family = "binomial")<br>
> anova(deacc.glm,deacc.glm1,test="Chisq")<br>
Analysis of Deviance Table<br>
<br>
Model 1: Deacc ~ Info<br>
Model 2: Deacc ~ Info1<br>
Resid. Df Resid. Dev Df Deviance P(>|Chi|)<br>
1 97 52.452<br>
2 98 71.600 -1 -19.148 1.210e-05<br>
<br>
In principle, you could do this with your mixed-effects model, being sure to use ML instead of REML fitting:<br>
<br>
> deacc.lmer = lmer (Deacc ~ Info + (1 | Subject), data = dat, family = "binomial",REML=F)<br>
> deacc.lmer1 = lmer (Deacc ~ Info1 + (1 | Subject), data = dat, family = "binomial",REML=F)<br>
> anova(deacc.lmer,deacc.lmer1)<br>
Data: dat<br>
Models:<br>
deacc.lmer1: Deacc ~ Info1 + (1 | Subject)<br>
deacc.lmer: Deacc ~ Info + (1 | Subject)<br>
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)<br>
deacc.lmer1 3 77.600 85.416 -35.800<br>
deacc.lmer 4 60.400 70.821 -26.200 19.2 1 1.177e-05 ***<br>
<br>
There is an argument that the likelihood-ratio test is anti-conservative and hence inappropriate for comparing mixed-effects models differing only in fixed-effects structure. (See Pinheiro & Bates, 2000, around page 76 or so. The argument is made only for linear mixed-effects models and I'm not sure of the status for logit mixed-effects models.) That being said, it's not clear how strong the anti-conservativity is, and in your case it seems like you have such an exceedingly powerful effect that you might be safe in using the likelihood-ratio test here and just mentioning the potential anti-conservativity as a caveat.<br>
<br>
So the summary is: believe it or not, how to assess the significance of the parameter estimate for "new" for your model & dataset is a bit of an open question, but it seems pretty clear that the estimate is significantly non-zero.<br>
<br>
Best & hope this helps.<br>
<br>
Roger<br><font color="#888888">
<br>
<br>
<br>
-- <br>
<br>
Roger Levy Email: <a href="mailto:rlevy@ucsd.edu" target="_blank">rlevy@ucsd.edu</a><br>
Assistant Professor Phone: 858-534-7219<br>
Department of Linguistics Fax: 858-534-4789<br>
UC San Diego Web: <a href="http://ling.ucsd.edu/%7Erlevy" target="_blank">http://ling.ucsd.edu/~rlevy</a></font><div><div></div><div class="Wj3C7c"><br>
<br>
_______________________________________________<br>
R-lang mailing list<br>
<a href="mailto:R-lang@ling.ucsd.edu" target="_blank">R-lang@ling.ucsd.edu</a><br>
<a href="http://pidgin.ucsd.edu/mailman/listinfo/r-lang" target="_blank">http://pidgin.ucsd.edu/mailman/listinfo/r-lang</a><br>
</div></div></blockquote></div><br>