[R-sig-ME] Understanding/plotting fixed effects estimates &standard errors
Rafael Maia
queirozrafaelmv at yahoo.com.br
Fri Jun 22 06:45:34 CEST 2012
hi dr Johnson,
many thanks for the help. I tried following your suggestion but changing the values does not change the predictions, so I might be doing something plain wrong (maybe because the predictor is factorial and its just using the model matrix and the estimated fixed effects to calculate the predictions? I'm not entirely sure). The results, as well as mm & vcov follow below:
>newdat=expand.grid(success=c(100,300,500,10000),fail=c(100,300,500,10000),factor=c('0',"1"))
>mm=model.matrix(terms(m1),newdat)
>estims=mm %*% fixef(m1)
> as.vector(estims)
[1] 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342
[9] 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342 1.438342
[17] 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149
[25] 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149 1.649149
>pvar1=diag(mm %*% tcrossprod(vcov(m1),mm))
> sqrt(pvar1)
[1] 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946
[8] 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946 0.5228946
[15] 0.5228946 0.5228946 0.5229200 0.5229200 0.5229200 0.5229200 0.5229200
[22] 0.5229200 0.5229200 0.5229200 0.5229200 0.5229200 0.5229200 0.5229200
[29] 0.5229200 0.5229200 0.5229200 0.5229200
>#simplifying a bit to print the model matrix
> newdat=expand.grid(success=c(100,500),fail=c(100,500),factor=c('0',"1"))
> newdat
success fail factor
1 100 100 0
2 500 100 0
3 100 500 0
4 500 500 0
5 100 100 1
6 500 100 1
7 100 500 1
8 500 500 1
> mm=model.matrix(terms(m1),newdat)
> mm
(Intercept) factor1
1 1 0
2 1 0
3 1 0
4 1 0
5 1 1
6 1 1
7 1 1
8 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$factor
[1] "contr.treatment"
> vcov(m1)
2 x 2 Matrix of class "dpoMatrix"
[,1] [,2]
[1,] 0.2734187339 -0.0008978599
[2,] -0.0008978599 0.0018223310
I really appreciate the attention. Thanks again,
Abraços,
Rafael Maia
---
webpage: http://gozips.uakron.edu/~rm72
"A little learning is a dangerous thing; drink deep, or taste not the Pierian spring." (A. Pope)
Graduate Student - Integrated Bioscience
University of Akron
http://gozips.uakron.edu/~shawkey/
On Jun 22, 2012, at 12:18 AM, Paul Johnson wrote:
> Can I make a wild guess/accusation below? Please?
>
> On Thu, Jun 21, 2012 at 10:17 PM, Rafael Maia
> <queirozrafaelmv at yahoo.com.br> wrote:
>> hi dr Duffy,
>>
>> many thanks for the attention and the reply. I followed the instructions under "Predictions and/or confidence (or prediction) intervals on predictions" of the wiki FAQ; however, unless I missed something, the results were nearly identical (and therefore confidence intervals wide & overlapping) to what I had by removing the intercept from the model:
>>
>>> m1=lmer(cbind(success,fail) ~ factor + (1|spp/variable), bb, family='binomial')
>>> summary(m1)
>> ...
>> Fixed effects:
>> Estimate Std. Error z value Pr(>|z|)
>> (Intercept) 1.43834 0.52289 2.751 0.00595 **
>> factor1 0.21081 0.04269 4.938 7.88e-07 ***
>>
>>> newdat=expand.grid(success=0,fail=0,factor=c('0',"1"))
>>> mm=model.matrix(terms(m1),newdat)
>>> estims=mm %*% fixef(m1)
>>> estims
>> [,1]
>> 1 1.438342
>> 2 1.649149
>>> pvar1=diag(mm %*% tcrossprod(vcov(m1),mm))
>
> The predicted value is not very different between the two levels of
> "factor" when success=0 and fail=0. This is a binomial model, the
> precision of the predictions is not "linear" or evenly spaced from
> left to right.
>
> Can you supply the same calculations for some more meaningful
> combinations of success and fail?
>
> When you do, let us see
>
> mm
>
> and
>
> vcov(m1)
>
> please.
>
>
>>> sqrt(pvar1)
>> [1] 0.5228946 0.5229200
>>
>> #compare to
>>
>>> m3=lmer(cbind(success,fail) ~ factor -1 + (1|spp/variable), bb, family='binomial')
>>> summary(m3)
>> ...
>> Fixed effects:
>> Estimate Std. Error z value Pr(>|z|)
>> factor0 1.4383 0.5229 2.751 0.00595 **
>> factor1 1.6491 0.5229 3.154 0.00161 **
>>
>> many thanks,
>>
>> Abraços,
>> Rafael Maia
>> ---
>> webpage: http://gozips.uakron.edu/~rm72
>> "A little learning is a dangerous thing; drink deep, or taste not the Pierian spring." (A. Pope)
>> Graduate Student - Integrated Bioscience
>> University of Akron
>> http://gozips.uakron.edu/~shawkey/
>>
>> On Jun 21, 2012, at 10:48 PM, David Duffy wrote:
>>
>>>
>>> This is nothing to with mixed modelling per se, but how your contrasts for the fixed effects are set up. If this was fixed effects only, you could fit -1 + factor (no intercept) to give estimates for each level of factor with standard errors that you could correctly interpret as you are trying to. The FAQ (http://glmm.wikidot.com/faq) shows plotting of confidence intervals for predictions.
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>
> --
> Paul E. Johnson
> Professor, Political Science Assoc. Director
> 1541 Lilac Lane, Room 504 Center for Research Methods
> University of Kansas University of Kansas
> http://pj.freefaculty.org http://quant.ku.edu
More information about the R-sig-mixed-models
mailing list