You may want to try the lsmeans package to make comparisons among treatments.

The lsmeans package does handle model objects of the the type "glmerMod", which are produced by the glmer.nb function in the lme4 package.  glmer.nb can be used for mixed-effects negative binomial regression.

If you are interested, let me know, and I will send you some sample code to try.

You may want to try the lsmeans package to make comparisons among treatments.

The lsmeans package does handle model objects of the the type "glmerMod", which are produced by the glmer.nb function in the lme4 package.  glmer.nb can be used for mixed-effects negative binomial regression.

If you are interested, let me know, and I will send you some sample code to try.

>I am actually more interested in the comparison between treatments within each level of I(walk/10):population.

Try hierarhical formula 

.....(I(walk/10):population)/treat.....

You may try to write down carefully your model and the comparisons you're interested in with the coefficients of the (fixed part of your) model. Then build a contrast matrix; the package multcomp can be of help. In common cases, lsmeans should do the job also.

Note that you seem to have a very complex model, with at least 35 parameters if I counted correctly, which means you should have a huge number of lizards, roughly 350, to hope to have a reliable fit of the data considering the model. Especially when going to third order interactions that are tricky to estimate correctly (they can be seen as corrections on corrections on main parameters...)

You may also consider using population as a random effect, not a fixed effect, depending on the exact question you're asking.

Thank you all very much for your answers. 

I figured out that I have made a mistake taking out the two-way interactions, while the three-way interaction is still in. If I then run the full model the three-way interaction is not significant anymore.

Still I think that it should be. If I analyse each of the six populations separately, I get significant differences between treatments in two out of 6 populations when I take the time walking in consideration (because the more they walk, the more they tongue flick anyway; so a significant I(walk/10):treatment interaction). This is for me a very interesting result, as it represents interpopulation differences in response towards an odour. My sample sizes for those two populations are 15 and 18 with each lizard observed 4 times (so at least 60 observations per population). So this should be ok. But as it is often not advised to split up your dataset, I wanted to put all populations in one model. Then I would expect that the interaction I(walk/10):treat:pop would become significant and this is like I said earlier not the case. In total I have approximately 100 lizards which sums up to about 400 observations, which also seems to me as a decent sample size. 

I suppose this is due to a not so big difference between treatments in "only" two out of six populations for the reason that Emmanuel Curis gave. Do you know if bootstrapping would be an option here? Or would it make things even more complex and is splitting my data in populations maybe in this case the better option?

Let me describe another point of view at Your data (I hope that I understand the experiment).

1 You can calculate intensity of flick ( count of flicks/duration of walk). So, variable "walk" walks out from model.

2 intensity of flicks may have Gauss distribution ( or close to Gauss)

3 ...And the model would be simplified : lmer (intesity ~ pop*treat + (1|trial)+(1|ind), data=dta)

Advantages: less variables, linear mixed model allow use test from package lmerTest (lsmeans(), difflsmeans()), plot means by treatment and population.

Hi Charlotte,

It looks like it is possible to fit your model using the glmer.nb function with the lme4 package as recommended above by Salvatore. Having fit the model as suggested by Olga above  you can then use the lsmeans package to compute model means and then compare them.

model fit:

require(lme4)

fit

Thank you for your reply,

Unfortunately I tried using glmer.nb but I have convergence problems. I think it is the same problem as described here, https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q2/024673.html, because the model worked relatively fine with glmmadmb.

Olga's idea is also a good idea, but then walk and tongue flicks should be very highly correlated, I think. I will need to check this first. 

Thank you all very much. I am learning a lot.

Charlotte

Dear Charlotte,

walk and tongue flicks are correlated, yes.

Let me explain:

Poisson distribution (You tried it first) expresses the probability of a given number of events occurring in a fixed interval.Lambda (parameter of Poisson distribution) is the average number of events per interval.

Your intervals (walks) are different. I propose to compare parameter lambda instead of solving head-on. 

You can get acquainted with negative binomial distribution (NBD). It has two parameters and you waste DFs on estimation of two parameters NBD+ additional variable "walks".

....so, may be as result head-on solving,Your model has convergence problems.

Hi Charlotte,

Fitting separately your different populations may be worth the additional work of further comparisons if you expect not only changes in average, but also changes in variances (some populations may be more variable than others for the number of flicks), or eventually changes in distributions between populations, that may be difficult to handle in a single model (unless there is a mean/variance relationship, and this relationship is captured by the model, like variance=mean for the Poisson model).

However, handling the different time of walk as correction factor for the number of flicks seems a key question, depending on what you expect.

If flick and walk are "independant" and flicking frequency is constant, so that flick number is just a matter of time, so that the longer the walk, the higher the flick number, then the approach given by Olga is probably the best and the simpler one. It also means that what is interesting is the flicking frequency, not the flicking count. It may be included in a Poisson model using an offset parameter. After log transformation, it would also be a particular case of using ln(walk length) as continuous covariate, as you did, assuming a slope of 1, but with increased complexity to interpret the model and may be more problems of convergence/too much parameters for the data).

If however you cannot assume this, either because the relationship between flicking and time is not linear, or because walk itself acts on the number of flicks [let's say for instance longer walk make the lizard more tired so he flicks less, or longer walk is because the lizard did not found the flicking stimulus so flick number should not depend in fine of the walk length], then Olga's approach does not work, but simply adding walk length as an additionnal parameter will not help very much for the interpretation.

So in fine, think carefully of what you expect for the relationship between walk length and flicking number, build your model accordingly and then only you can check your assumptions confidently...

After checking I think indeed that the problem here is that the distributions differ between populations. When doing analyses for the mainland and one of the islands, a poisson fit is the best choice for the model and when analysing the remainder of the islands a negative binomial fit is better. If I fit the mainland populations with a negative binomial model, my model becomes overdispersed. And if I put everything together in one model, a negative binomial is actually the better fit in general. This probably masks the effect on the mainland populations. So I think a seperate analyses is in the end the best choice.

I have the feeling that lizards flick at a certain rate untill they find the odour source. After this they stop walking and continue flicking at the odour source to get more information. So I think this doesn't allow for deviding the tongue flicks by walk, if I understood correctly.

The approach with the confidence intervals is interesting! Thank you for sharing, Niklas. 

Do you have any information in your data to distinguish lizards that find a source and those that don't? This may be a covariate that could be at the origin of your apparent overdispersion... That does not preclude the use of the ratio (I mean, it may lead to statistically sound methods), but you need to carefully think to which question it will exactly answer.

I don't quite understand the sentences « When doing analyses for the mainland a Poisson fit is the best choice » and « If I fit the mainland populations with a negative binomial model, my model becomes overdispersed », which may seem  contradictory. Do you mean that the overdispersion parameter is not equal to 0 (or 1), the value it should have for a Poisson law, but seems evident since you allow it to vary? Or that it is _significantly_ overdispersed, but in that case how do you then consider that Poisson is a better model?

I ran the same model twice for the data of the two mainland populations, but first with a poisson fit and then with a negative binomial fit:

a1

I don't think the value of the dispersion parameter allows you to say if the model fit best or worst the model. You should compare likelihoods for that, for instance.

In the contrary, if you have a very high dispersion parameter, it may well be because it fits better since the data needs it (or because the distribution is not binomial negative)...

At this point, difficult to say without graphs, data and so on...