[R-sig-ME] Mixed Effects Logistic Regression Model and finding p-values

[Charla Patterson]

Hello List Members,

I am a masters student that is relatively new to R and am currently working on the analysis of my project. For my project, I collected data on individual birds every year (sex, migratory or resident) along with environmental measurements from their breeding grounds (average min. winter temperature, total winter precipitation, ect.). I am interested in finding out whether the decision to migrate in this population is affected by environmental cures.  

Here is a sample of my data from the original .csv file:

bnum	sex 	ystat 	year	wintering site	breeding site  direct.  minall  totprecall
39	1	1	1991.92	Burlington	Wye marsh	S	-4.5	98.75
39	1	1	1991.92	Burlington	Wye marsh	N	-4.5	98.75
75	0	1	1991.92	Aurora	        Wye marsh	        -4.5	98.75
78	1	0	1991.92	Wye Marsh	Wye Marsh		-4.5	98.75
79	0	0	1991.92	Wye Marsh	Wye Marsh		-4.5	98.75
80	1	0	1991.92	Wye Marsh	Wye Marsh		-4.5	98.75
961	0	1	1991.92	Burlington	Wye Marsh	S	-4.5	98.75
961	0	1	1991.92	Burlington	Wye Marsh	N	-4.5	98.75
74	1	0	1992.93	Wye marsh	Wye marsh		-3.9	125.65
75	0	0	1992.93	Wye marsh	Wye marsh		-3.9	125.65
78	1	0	1992.93	Wye Marsh	Wye Marsh		-3.9	125.65


Where bnum is the bird ID, sex (Female=1, male=0) and ystat is the decision to migrate (migrant=1, resident=0), minall is average minimum winter temperature for the year and totprecall is total winter precipitation for the year.

I wanted to use a mixed effects logistic regression model where ystat is the response (binary), the fixed effects would be sex, minall and totprecall  and the random effects would be bnum and year.

my proposed model and the corresponding results are:

> yr<- as.factor(year)
> mod3<- glmer(ystat~minall + totprecall+ (1 |yr)+(1 | bnum), data=birds,family=binomial(link='logit'))
> summary(mod3) 

Generalized linear mixed model fit by the Laplace approximation 
Formula: ystat ~ minall + totprecall + (1 | yr) + (1 | bnum) 
   Data: birds 
   AIC   BIC logLik deviance
 857.3 881.9 -423.6    847.3
Random effects:
 Groups Name        Variance Std.Dev.
 bnum   (Intercept) 20.0576  4.4786  
 yr     (Intercept)  2.2188  1.4896  
Number of obs: 1014, groups: bnum, 209; yr, 19

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.423079   2.088999  -0.681   0.4957  
minall      -0.419361   0.212170  -1.976   0.0481 *
totprecall  -0.006249   0.022842  -0.274   0.7844  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
           (Intr) minall
minall     -0.046       
totprecall -0.904  0.408


Although I think my model is okay, I was struggling with a couple things:

1. When I use year as a fixed effect in my model instead of as a random effect, the p-values for all other parameters are much lower. I don't understand why and it makes me apprehensive about the model and p-values I have. Are the p-values I reported here reliable?

2. I read that a good alternative is to use mcmcglmm to get p-values and HPDintervals, but I'm finding it tricky as well.

I tried to use mcmcsamp along with the pvals.fnc(mod3, nsim=10000) function, but it turns out that this doesn't work with lme4 package. I read that mcmcglmm is a good alternative, so I read the R notes but found that I can't get it to work. Here is my code and the associated warning messages:

m <- MCMCglmm(ystat ~ minall, random = ~bnum+yr, data = birds, family = "categorical")

Error in buildZ(rmodel.terms[r], data = data, formZ = TRUE, nginverse = names(ginverse)) : 
  missing values in random predictors

I have no idea what this means? I am really stuck and I feel like I'm missing something really important here.

3. This may be a very basic question, but could someone provide directions or a reference as to how to plot the logistic curve representing the relationship between minall and ystat? I used plot(minall,ystat) along with some other more complicated codes; however, I have not been able to produce a plot that shows the logistic curve.


Any guidance on any of these issues would be greatly appreciated! Thank you in advance for your time,

Charla