[R-sig-ME] Curved residuals vs fitted plot

Wed Mar 12 12:46:16 CET 2014

Hi Victoria,

I get quite a different residual plot from your model - code pasted below - which looks much better. There still appears to be a weak residual trend but these plots can be deceptive especially if there's a mass of points in one part of the graph. The apparent trend might be an artefact of skew at either end of the probability scale. I believe there's no requirement for these residuals to be normal, or even symmetrical, and they're more likely to be skewed close to predicted probabilities of 0 or 1.  E.g. you have several rows with zero responses, which inevitably produce a lump of negative residuals. (This residual plot should also be viewed with caution for the reason I mentioned earlier about standardising the residuals using the binomial variance function rather than the logitnormal-binomial variance function (does it exist, anyone?).)

Probably a better way to assess the model is to plot on the same graph the data and the model predictions, with prediction intervals. E.g. plot prediction intervals on boxplot(Aphid.rm.prop ~ Treatment*Distance.m, data=dframe1).

> > Would you recommend I abandon the two-vector response and try for binary? If I use 0, 1s it can't be overdispersed, but has less power.

If you expand the data from binomial to binary, and still fit the same random effects, you are fitting the same model, with the same power, except that the observation-level overdispersion random effect becomes a more conventional "group" random effect. The only difference will be to make the model fit more slowly. E.g. "+ (1|Site)" implies the the same model in both of these data sets:

Collapsed:
1 2 Site1
2 1 Site2

Expanded:
0 Site1
1 Site1
0 Site1
1 Site2
1 Site2
0 Site2

Best wishes,
Paul

sapply(list(Field,Treatment,Field.section,Transect,Distance.m,Aphid.rm.success,Aphid.remain.fail,fffPlot,Bait.ID),function(x)length(unique(x)))
### Bait.ID has a duplicated level so has n-1 levels - remake with n levels:
Bait.ID<-paste("ID",1:length(Bait.ID),sep="")
dframe1<-data.frame(Field,Treatment,Field.section,Transect,Distance.m,Aphid.rm.success,Aphid.remain.fail,fffPlot,Bait.ID)
rm(Field,Treatment,Field.section,Transect,Distance.m,Aphid.rm.success,Aphid.remain.fail,fffPlot,Bait.ID) # removes ambiguity about source of objects
head(dframe1)
dframe1$n <- dframe1$Aphid.rm.success + dframe1$Aphid.remain.fail
dframe1$Aphid.rm.prop <- dframe1$Aphid.rm.success / dframe1$n

plot(Aphid.rm.prop ~ jitter(Distance.m,factor=0.2), data=dframe1,col=as.numeric(Treatment),pch=as.numeric(Treatment))
boxplot(Aphid.rm.prop ~ Treatment, data=dframe1)
boxplot(Aphid.rm.prop ~ Treatment*Distance.m, data=dframe1)

### I changed Field from random to fixed, because although it's conceptually a random effect, two levels are too few to get a reasonable variance estimate
model3<-glmer(cbind(Aphid.rm.success,Aphid.remain.fail) ~ Treatment*Distance.m+Field+(1|Field.section/Transect/fffPlot/Bait.ID), # Bait.ID to remove overdispersion
binomial(link=logit),data=dframe1) # note I kept fffPlot smallest scale in (better residuals)
summary(model3)
plot(model3) # issue? Or just messy ecological data?

### look at random effect residuals
dotplot(ranef(model3,condVar=TRUE))$Bait.ID
dotplot(ranef(model3,condVar=TRUE))$fffPlot
dotplot(ranef(model3,condVar=TRUE))$Transect  # close to zero variance at this level
dotplot(ranef(model3,condVar=TRUE))$Field.section  # close to zero variance at this level
### these look sufficiently close to normal, but should do a QQ plot, at least on the factors that vary
qqmath(ranef(model3,condVar=TRUE))$Bait.ID
qqmath(ranef(model3,condVar=TRUE))$fffPlot

# shift overdispersion random effect from fitted values to residuals
Fitted <- plogis(qlogis(fitted(model3)) - ranef(model3)$Bait.ID[[1]])
Resid <- (dframe1$Aphid.rm.prop - Fitted) / sqrt(Fitted * (1 - Fitted)/dframe1$n)
# plot residuals differentiating between treatments and distances
plot(Fitted, Resid, col=dframe1$Treatment, pch=21, bg=grey(pnorm(scale(dframe1$Distance.m))))
abline(h=0)

On 11 Mar 2014, at 15:57, Victoria Wickens <V.J.Wickens at pgr.reading.ac.uk> wrote:

> Apologies, the link to the plot did not work. Here is the functioning link to illustrate the bizzare plot  http://imgur.com/rx2Oc3s
>
> Would you suggest this plot illustrates a bad model fit?
>
> Would you recommend I abandon the two-vector response and try for binary? If I use 0, 1s it can't be overdispersed, but has less power. However, if the first model assumptions cannot be upheld because of the residuals, binary might be my only option?
>
> Thank once again for taking the time to read this,
>
> Victoria
> ______________
> Victoria Wickens
> PhD student
> Room GU08
> Centre for Agri-Environmental Research (CAER)
> School of Agriculture, Policy & Development
> University of Reading
> Reading RG6 7BE
> ____________
>
> Email: v.j.wickens at pgr.reading.ac.uk
> http://www.reading.ac.uk/caer/staff_students.html
> https://www.facebook.com/UniRdgAPD
>
> ________________________________________
> From: r-sig-mixed-models-bounces at r-project.org <r-sig-mixed-models-bounces at r-project.org> on behalf of Victoria Wickens <V.J.Wickens at pgr.reading.ac.uk>
> Sent: 10 March 2014 14:29
> To: Paul Johnson; r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Curved residuals vs fitted plot
>
> Hi Paul,
>
> Thank you for reading my email and for providing clear explanation of why I have the curvature. I found your response and the link very helpful. Sadly I am having trouble getting the R book which is not available at my library, but I do have a copy of the Mixed Effects Models and Extensions in Ecology with R by Zuur.
>
> I now understand the fitted vs resid plots produce curved scatter trends due to the observational random effect and have implemented the code you provided to shift the effect from the fitted values to the residuals to create a new plot (see here http://imgur.com/edit).
>
> My plot is not nearly as well behaved as your example though, with a steep tail of residuals with high fitted values. Does this mean my model is a bad fit? I'm not sure of the next best move to take.
>
> I have tried the following (all with the observational random effect, Bait.ID, to remove overdispersion but with similar resid vs fitted plots as seen on the web-link):
>
> 1. glmer
> 2. glmer + polynomial regression (logDistance factor)
> 3. glmer + quadratic term (Distance^2)
> 4. gamm (limited success)
>
> Is there another way to remove overdispersion or account for it? Or am I going about this wrong?
>
> Thank you once again for taking time to help me. I really appreciate it,
>
> Victoria
> ______________
> Victoria Wickens
> PhD student
> Room GU08
> Centre for Agri-Environmental Research (CAER)
> School of Agriculture, Policy & Development
> University of Reading
> Reading RG6 7BE
> ____________
>
> Email: v.j.wickens at pgr.reading.ac.uk
> http://www.reading.ac.uk/caer/staff_students.html
> https://www.facebook.com/UniRdgAPD
>
> ________________________________________
> From: Paul Johnson <paul.johnson at glasgow.ac.uk>
> Sent: 08 March 2014 01:26
> To: Victoria Wickens; r-sig-mixed-models at r-project.org
> Subject: RE: Curved residuals vs fitted plot
>
> Hi Victoria,
>
> Positive correlation between the residuals and fitted values is expected whenever you have a row-level random effect to mop up overdispersion. It happens because overdispersion implies that the low responses are lower than predicted and the high responses are higher than predicted, even allowing for binomial sampling error. The gap between a low predicted proportion and a still lower observed proportion is bridged partly by a negative residual and partly by a negative overdispersion random effect (which can be thought of as an additional 'free' residual to the constrained residual due to binomial sampling error). The same happens at the highest predicted proportions, except that both residual and random effect will be positive, hence the positive correlation between the predicted (fitted) values, which contain the overdispersion random effects, and the residuals.
>
> This problem is dealt with, for binomial GLMMs, in chapter 7 of Alain Zuur's Beginner's guide to GLM and GLMM with R. I don't have this book handy and I can't remember his solution. One solution is to shift the overdispersion random effect from the fitted values to the residuals, which makes sense given that they are really just add-on residuals - see the example code below.
>
> This thread discusses the same problem for lognormal-Poisson GLMMs:
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2013q3/020819.html
>
> Best wishes,
> Paul
>
> # simulate binomial data with a positive trend and overdispersion
> obs <- factor(1:100)
> x <- seq(-1, 0.98, 0.02)
> set.seed(12345678)
> y <- rbinom(100, 20, plogis(x + rnorm(100, 0, 1)))
> par(mfrow=c(3, 1))
> plot(y/20 ~ x)
> # fit logitnormal-binomial GLMM
> fit <- glmer(cbind(y, 20-y) ~ I(1:100) + (1|obs), family="binomial")
> summary(fit)
> # standard residual vs fitted plot
> plot(fitted(fit), resid(fit))
> abline(h=0)
> # shift overdispersion random effect from fitted values to residuals
> Fitted <- plogis(qlogis(fitted(fit)) - ranef(fit)$obs[[1]])
> Resid <- (y/20 - Fitted) / sqrt(Fitted * (1 - Fitted)/20)
> plot(Fitted, Resid)
> abline(h=0)
> # NB THIS GETS RID OF THE TREND BUT ISN'T RIGHT!
> # This is because I've standardised the variance using the
> # variance function for the binomial, p(1-p)/n, but the distribution
> # underlying the residuals is now a logitnormal-binomial
> # so should incorporate the overdispersion variance. I don't
> # know how to do this.
>
>
> ________________________________________
> From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] On Behalf Of Victoria Wickens [V.J.Wickens at pgr.reading.ac.uk]
> Sent: 07 March 2014 19:20
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Curved residuals vs fitted plot
>
> Dear mixed-models group,
>
> Thank you for taking time to read this plea for help.
>
> I wondered if you could give me any advice regarding the curved residuals vs fitted plot with my model:
>
> model3<-glmer(y~Treatment*Distance.m+(1|Field/Field.section/Transect/fffPlot/Bait.ID),
> binomial(link=logit),data=dframe1)
>
> My aim is to explain aphid removal from bait trays in response to treatment (presence of absence of a flower strip) and distance (five distances, or plots, moving away from the margins; 1, 5, 10, 15, and 20m). (More experimental details after message in case helpful).
>
> I am using R version 3.0.2 and the latest version of lme4. The resid vs fitted figure is not great, is this because I need to look at a logit fitted scale? (See here to view the figure: <http://imgur.com/CYe5CxC> http://imgur.com/CYe5CxC)
>
> Any help at all would be more than amazing,
>
> Victoria
>
> Experimental setup
> There are two fields overall, each contains one treatment (Mit) and one control (Non). Looking at one treatment, in one field section, there are three transects parallel to each other and perpendicular to the margin treatment (T1, T2 and T3). There are five plots along each transect at specific distances from the treatment margins (1, 5, 10, 15, 20m). At each plot, three replica bait trays were placed (a, b and c). Each tray contained 20 aphids glued to foam. In total, I had two farm fields at the largest scale and between them there were 180 bait trays overall (90 per field).
>
>
> Script attached, also listed below
>
>
> library(lme4)
>
> Field<-c("L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,"L" ,
> "B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B")
>
> Treatment<-c("Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,
> "Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,
> "Mit" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,
> "Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Non" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit" ,"Mit")
>
> Field.section<-c("D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"D" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"C" ,"B" ,"B" ,
> "B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"B" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A" ,"A")
>
> Transect<-c("T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T10" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T11" ,"T12" ,"T12" ,"T12" ,"T12" ,"T12" ,"T12" ,"T12" ,
> "T12" ,"T12" ,"T12" ,"T12" ,"T12" ,"T12" ,"T12" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T7" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T8" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,
> "T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T9" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T4" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T5" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,"T6" ,
> "T6" ,"T6" ,"T6" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T1" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T2" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3" ,"T3")
>
> Distance.m<-c(1 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10
> ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20 ,1 ,1 ,1 ,5 ,5 ,5 ,10 ,10 ,10 ,15 ,15 ,15 ,20 ,20 ,20)
>
> Aphid.rm.success<-c(3 ,7 ,19 ,19 ,2 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,0 ,0 ,8 ,16 ,2 ,0 ,2 ,1 ,1 ,3 ,0 ,0 ,9 ,10 ,15 ,18 ,1 ,0 ,14 ,10 ,4 ,1 ,3 ,12 ,2 ,4 ,0 ,20 ,17 ,19 ,11 ,19 ,5 ,0 ,0 ,3 ,0 ,0 ,0 ,3 ,11 ,0 ,19 ,13 ,13 ,10 ,9 ,5 ,19 ,1 ,2 ,0 ,12 ,2 ,0 ,1 ,7 ,15 ,19 ,6 ,19 ,18 ,17 ,1 ,0 ,0 ,1 ,1 ,9 ,18 ,11 ,19 ,1 ,3 ,0 ,0 ,1 ,1 ,0 ,1 ,1 ,3 ,0 ,0 ,0 ,4 ,1 ,0 ,0 ,2 ,1 ,6 ,0 ,0 ,0 ,1 ,0 ,0 ,1 ,1 ,0 ,0 ,2 ,1 ,11 ,1 ,0 ,3 ,3 ,0 ,1 ,3 ,0 ,0 ,1 ,0 ,0 ,8 ,1 ,3 ,2 ,1 ,3 ,0 ,
> 1 ,2 ,0 ,1 ,1 ,0 ,0 ,0 ,8 ,0 ,0 ,0 ,7 ,1 ,0 ,1 ,1 ,0 ,1 ,0 ,1 ,0 ,1 ,5 ,0 ,20 ,0 ,17 ,0 ,10 ,1 ,0 ,5 ,0 ,0 ,0 ,0 ,0)
>
> Aphid.remain.fail<-c(17 ,13 ,1 ,1 ,18 ,20 ,20 ,20 ,20 ,20 ,20 ,19 ,19 ,19 ,17 ,19 ,19 ,19 ,20 ,20 ,12 ,4 ,18 ,20 ,18 ,19 ,19 ,17 ,20 ,20 ,10 ,9 ,5 ,2 ,19 ,20 ,6 ,10 ,16 ,19 ,17 ,8 ,18 ,16 ,20 ,0 ,3 ,1 ,9 ,2 ,15 ,20 ,20 ,17 ,20 ,20 ,20 ,17 ,9 ,20 ,1 ,7 ,7 ,10 ,11 ,15 ,1 ,19 ,18 ,20 ,8 ,18 ,20 ,19 ,13 ,5 ,1 ,14 ,0 ,1 ,3 ,19 ,20 ,20 ,19 ,18 ,11 ,2 ,9 ,0 ,19 ,18 ,20 ,20 ,19 ,19 ,20 ,19 ,19 ,16 ,20 ,20 ,20 ,16 ,19 ,20 ,20 ,18 ,19 ,14 ,20 ,20 ,20 ,19 ,20 ,20 ,19 ,19 ,20 ,20 ,18 ,19 ,
> 9 ,19 ,19 ,18 ,17 ,20 ,19 ,17 ,20 ,20 ,19 ,20 ,17 ,12 ,19 ,17 ,18 ,19 ,17 ,20 ,19 ,18 ,20 ,18 ,20 ,20 ,19 ,20 ,12 ,20 ,20 ,20 ,13 ,19 ,20 ,19 ,19 ,20 ,20 ,20 ,19 ,20 ,19 ,15 ,20 ,0 ,20 ,3 ,20 ,10 ,19 ,20 ,15 ,20 ,20 ,20 ,20 ,20)
>
> fffPlot<-c("D T1 1" ,"D T1 1" ,"D T1 1" ,"D T1 1" ,"D T1 5" ,"D T1 5" ,"D T1 5" ,"D T1 10" ,"D T1 10" ,"D T1 10" ,"D T1 15" ,"D T1 15" ,"D T1 15" ,"D T1 20" ,"D T1 20" ,"D T1 20" ,"D T2 1" ,"D T2 1" ,"D T2 1" ,"D T2 5" ,"D T2 5" ,"D T2 5" ,"D T2 10" ,"D T2 10" ,"D T2 10" ,"D T2 15" ,"D T2 15" ,"D T2 15" ,"D T2 20" ,"D T2 20" ,"D T2 20" ,"D T3 1" ,"D T3 1" ,"D T3 1" ,"D T3 5" ,"D T3 5" ,"D T3 5" ,"D T3 10" ,"D T3 10" ,"D T3 15" ,"D T3 15" ,"D T3 15" ,"D T3 20" ,"D T3 20" ,"D T3 20" ,
> "C T1 1" ,"C T1 1" ,"C T1 1" ,"C T1 5" ,"C T1 5" ,"C T1 5" ,"C T1 10" ,"C T1 10" ,"C T1 10" ,"C T1 15" ,"C T1 15" ,"C T1 15" ,"C T1 20" ,"C T1 20" ,"C T1 20" ,"C T2 1" ,"C T2 1" ,"C T2 1" ,"C T2 5" ,"C T2 5" ,"C T2 5" ,"C T2 10" ,"C T2 10" ,"C T2 10" ,"C T2 15" ,"C T2 15" ,"C T2 15" ,"C T2 20" ,"C T2 20" ,"C T2 20" ,"C T3 1" ,"C T3 1" ,"C T3 1" ,"C T3 5" ,"C T3 5" ,"C T3 5" ,"C T3 10" ,"C T3 10" ,"C T3 10" ,"C T3 15" ,"C T3 15" ,"C T3 15" ,"C T3 20" ,"C T3 20" ,"C T3 20" ,"B T1 1" ,
> "B T1 1" ,"B T1 1" ,"B T1 5" ,"B T1 5" ,"B T1 5" ,"B T1 10" ,"B T1 10" ,"B T1 10" ,"B T1 15" ,"B T1 15" ,"B T1 15" ,"B T1 20" ,"B T1 20" ,"B T1 20" ,"B T2 1" ,"B T2 1" ,"B T2 1" ,"B T2 5" ,"B T2 5" ,"B T2 5" ,"B T2 10" ,"B T2 10" ,"B T2 10" ,"B T2 15" ,"B T2 15" ,"B T2 15" ,"B T2 20" ,"B T2 20" ,"B T2 20" ,"B T3 1" ,"B T3 1" ,"B T3 1" ,"B T3 5" ,"B T3 5" ,"B T3 5" ,"B T3 10" ,"B T3 10" ,"B T3 10" ,"B T3 15" ,"B T3 15" ,"B T3 15" ,"B T3 20" ,"B T3 20" ,"B T3 20" ,"A T1 1" ,"A T1 1" ,"A T1 1" ,
> "A T1 5" ,"A T1 5" ,"A T1 5" ,"A T1 10" ,"A T1 10" ,"A T1 10" ,"A T1 15" ,"A T1 15" ,"A T1 15" ,"A T1 20" ,"A T1 20" ,"A T1 20" ,"A T2 1" ,"A T2 1" ,"A T2 1" ,"A T2 5" ,"A T2 5" ,"A T2 5" ,"A T2 10" ,"A T2 10" ,"A T2 10" ,"A T2 15" ,"A T2 15" ,"A T2 15" ,"A T2 20" ,"A T2 20" ,"A T2 20" ,"A T3 1" ,"A T3 1" ,"A T3 1" ,"A T3 5" ,"A T3 5" ,"A T3 5" ,"A T3 10" ,"A T3 10" ,"A T3 10" ,"A T3 15" ,"A T3 15" ,"A T3 15" ,"A T3 20" ,"A T3 20" ,"A T3 20")
>
> Bait.ID<-c("ID1" ,"ID2" ,"ID3" ,"ID3" ,"ID4" ,"ID5" ,"ID6" ,"ID7" ,"ID8" ,"ID9" ,"ID10" ,"ID11" ,"ID12" ,"ID13" ,"ID14" ,"ID15" ,"ID16" ,"ID17" ,"ID18" ,"ID19" ,"ID20" ,"ID21" ,"ID22" ,"ID23" ,"ID24" ,"ID25" ,"ID26" ,"ID27" ,"ID28" ,"ID29" ,"ID30" ,"ID31" ,"ID32" ,"ID33" ,"ID34" ,"ID35" ,"ID36" ,"ID37" ,"ID38" ,"ID39" ,"ID40" ,"ID41" ,"ID42" ,"ID43" ,"ID44" ,"ID45" ,"ID46" ,"ID47" ,"ID48" ,"ID49" ,"ID50" ,"ID51" ,"ID52" ,"ID53" ,"ID54" ,"ID55" ,"ID56" ,"ID57" ,"ID58" ,"ID59" ,"ID60" ,"ID61" ,"ID62" ,
> "ID63" ,"ID64" ,"ID65" ,"ID66" ,"ID67" ,"ID68" ,"ID69" ,"ID70" ,"ID71" ,"ID72" ,"ID73" ,"ID74" ,"ID75" ,"ID76" ,"ID77" ,"ID78" ,"ID79" ,"ID80" ,"ID81" ,"ID82" ,"ID83" ,"ID84" ,"ID85" ,"ID86" ,"ID87" ,"ID88" ,"ID89" ,"ID90" ,"ID91" ,"ID92" ,"ID93" ,"ID94" ,"ID95" ,"ID96" ,"ID97" ,"ID98" ,"ID99" ,"ID100" ,"ID101" ,"ID102" ,"ID103" ,"ID104" ,"ID105" ,"ID106" ,"ID107" ,"ID108" ,"ID109" ,"ID110" ,"ID111" ,"ID112" ,"ID113" ,"ID114" ,"ID115" ,"ID116" ,"ID117" ,"ID118" ,"ID119" ,"ID120" ,"ID121" ,"ID122" ,
> "ID123" ,"ID124" ,"ID125" ,"ID126" ,"ID127" ,"ID128" ,"ID129" ,"ID130" ,"ID131" ,"ID132" ,"ID133" ,"ID134" ,"ID135" ,"ID136" ,"ID137" ,"ID138" ,"ID139" ,"ID140" ,"ID141" ,"ID142" ,"ID143" ,"ID144" ,"ID145" ,"ID146" ,"ID147" ,"ID148" ,"ID149" ,"ID150" ,"ID151" ,"ID152" ,"ID153" ,"ID154" ,"ID155" ,"ID156" ,"ID157" ,"ID158" ,"ID159" ,"ID160" ,"ID161" ,"ID162" ,"ID163" ,"ID164" ,"ID165" ,"ID166" ,"ID167" ,"ID168" ,"ID169" ,"ID170" ,"ID171" ,"ID172" ,"ID173" ,"ID174" ,"ID175" ,"ID176" ,"ID177" ,"ID178" ,
> "ID179")
>
>
> dframe1<-data.frame(Field,Treatment,Field.section,Transect,Distance.m,Aphid.rm.success,Aphid.remain.fail,fffPlot,Bait.ID)
>
> y<-cbind(Aphid.rm.success,Aphid.remain.fail)
>
> model3<-glmer(y~Treatment*Distance.m+(1|Field/Field.section/Transect/fffPlot/Bait.ID), # Bait.ID to remove overdispersion
> binomial(link=logit),data=dframe1) # note I kept fffPlot smallest scale in (better residuals)
> summary(model3)
>
> plot(model3) # issue? Or just messy ecological data?
>
> ______________
> Victoria Wickens
> PhD student
> Room GU08
> Centre for Agri-Environmental Research (CAER)
> School of Agriculture, Policy & Development
> University of Reading
> Reading RG6 7BE
> ____________
>
> Email: v.j.wickens at pgr.reading.ac.uk<mailto:v.j.wickens at pgr.reading.ac.uk>
> http://www.reading.ac.uk/caer/staff_students.html
> https://www.facebook.com/UniRdgAPD
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models