model check for negative binomial model

CORRECTION


Dear Alessandra,


My apologies, I msread the output you pasted in your previous messages. I
did not realize that the parameter estimates that you were providing
correspond to a negative binomial model, I assumed that they came from a
binomial (logit) model. As Profr. Bolker mentioned, you cannot fit a
negative binomial model to your data, it has to be a binomial (logit model)
in which you model the number of eggs that hatch out of the total number of
eggs in the nest. What the model predicts (after back transforming the odds
ratio) is the probability, p, of an egg hatching, which can then be
translated into the expected number of eggs that hatch simply by
multiplying it by the number of eggs in the nest, np.


I hope this helps and, again, my apologies for the confusion and for
messing things yp.


Salvador


En Lun, 17 Febrero, 2020 en 20:46, yo <salvadorsanchezcolon at prodigy.net.mx>
escribi?:

Para: bielli.alessandra at gmail.com
Cc: bbolker at gmail.com; r-sig-mixed-models at r-project.org

Cara Alessandra,


If I may interject into your conversation, the key to your question lies
in the parameter estimates:for the fixed effects:


                               Estimate           Std. Error             z
value        Pr(>|z|)
(Intercept)               -1.38864             0.08227
-16.879      < 2e-16 ***
Relocation..Y.N.Y       0.32105              0.09152                3.508
     0.000452 ***
SPL                          0.22218               0.08793
2.527       0.011508 *


As I understand your design, you have two independent variables (factors):
Species and treatment, each with two levels: L vs. G and Y vs N.


Thus, the logit model that you obtain (omitting the random effects) is:


logit(p) = log(p/(1-p) = -1.38864 + 0.32105*Treatment + 0.22218*Species

Treatment and species are two-level factors which are coded as 0 or 1 and
the model has to have parameter estimates for each of those levels. By
design, in generalized linear models parameters are estimated by taking one
of the levels of each factor as the reference level and assigned a value of
0, and the parameters for the other levels are expressed in relation to the
reference level. Thus, in your case, treatment N and species G are
designated as the reference levels for their corresponding factors and,
hence, their parameter values are both 0 (and not listed in the output).
Then, treatment Y has a positive effect (0.32) on the odds ratio of
hatching/not-hatching compared to treatment N; and species L also has a
positive effect (0.22 times) on the odds ratio of hatching vs not-hatching.


The intercept then (-1.38864) denotes the odds ratio for the combination
of treatment N and species G; that is, when both factors are at their
reference level of 0. This translates into a probability of hatching (p):


p =  exp(-1.38864)/(1+exp(-1.38864) = 0.1996


For treatment Y and species G, the model becomes:


p =  exp(-1.38864 + 0.32105)/(1+exp(-1.38864 + 0.32105) = 0.25569


as treatment Y has a positive effect on the odds ratio of hatching.


For treatment N and species L, the model becomes:


p =  exp(-1.38864 + 0.22218)/(1+exp(-1.38864 + 0.22218) = 0.2375


as species L has a positive effect on the odds ratio of hatching.


Finally, for the combination of treament Y and species L, the model
becomes:


p =  exp(-1.38864 + 0.32105 + 0.22218)/(1+exp(-1.38864 + 0.32105 + 0.22218)
= 0.3002


as both levels Y and L have positive effects on the odds ratio of hatching.


I hope this helps.


Best regards,


Salvador






En Lun, 17 Febrero, 2020 en 18:2, Alessandra Bielli <
bielli.alessandra at gmail.com> escribi?:

Para: Ben Bolker
Cc: r-sig-mixed-models at r-project.org
Dear Ben

I am trying to make a prediction for the combination of species (L or G)
and treatment (control/experiment).

I am still confused about the prediction values. I would like to present
results as a success rate for a nest, to say that treatment
increases/decreases success by ...%. But the value I have is the
probability that 1 egg in the nest succeeds, correct? I am not sure how to
use these predictions.

Thanks for your help!

Alessandra

On Mon, Feb 17, 2020 at 2:15 PM Ben Bolker <bbolker at gmail.com> wrote:

That's correct. There are some delicate issues about prediction:

* do you want to use the original (potentially unbalanced) data for
prediction? (That's what you're doing here).
* or, do you want to make predictions for a "typical" nest and week
combination, in which case you would use


pframe <- with(your_data,
expand.grid(Relocation..Y.N.=unique(Relocation..Y.N.),
SP=unique(SP))
predict(m.unhatched,type="response",re.form=NA,newdata=pframe))

You could also use expand.grid() to generate a balanced design (i.e.
all combinations of weeks and nests), which would give yet another

answer.

There are a lot of packages designed for doing these kinds of
post-fitting manipulations (e.g. 'margins', 'emmeans'), you might find
them useful ...



On 2020-02-17 1:48 p.m., Alessandra Bielli wrote:

Dear Thierry

Thanks for your reply.

I read a bit about the prediction for a binomial model with
success/failures and I have a couple of questions.

If I use the predict function with the model you recommended, I obtain
log.odds or probabilities if I use "type=response":

tapply(predict(m.unhatched,type="response"),list(main$SP,main$Relocation..Y.N.),mean)

N Y
G 0.7314196 0.6414554
L 0.6983576 0.6003087

Are these probabilities of success (i.e. hatched) in one nest?

Thanks,

Alessandra

On Mon, Feb 17, 2020 at 7:18 AM Thierry Onkelinx
<thierry.onkelinx at inbo.be <mailto:thierry.onkelinx at inbo.be>> wrote:

Dear Alessandra,

Since you have both the number hatched and the total clutch size you
can calculate the number of successes and failures. That is
sufficient for a binomial distribution.

glmer(cbind(Hatched, Unhatched) ~ Relocation..Y.N. + SP + (1 |
Beach_ID) + (1 | Week), family = binmial)

A negative binomial or Poisson allow predictions larger than the
offset. Which is nonsense given that the number hatched cannot
surpass the total clutch size.

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR
NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be <mailto:thierry.onkelinx at inbo.be>
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be <http://www.inbo.be>

///////////////////////////////////////////////////////////////////////////////////////////

To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given
body of data. ~ John Tukey

///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op wo 12 feb. 2020 om 18:42 schreef Alessandra Bielli
<bielli.alessandra at gmail.com <mailto:bielli.alessandra at gmail.com>>:

Dear Ben

Thanks for your quick response.

Yes, emergence success is usually between 60 and 80% or higher.
I am not sure how to use a binomial, if my data are counts?

Can you explain why the approximation doesn't work well if
success gets
much above 50%? Does it make sense, then, to have "unhatched" as
dependent
variable, so that I predict mortality (usually below 50%) using
a nb with
offset(log(total clutch)) ?

summary(m.emerged)

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
Family: Negative Binomial(2.2104) ( log )
Formula: Unhatched ~ Relocation..Y.N. + SP +
offset(log(Total_Clutch)) +
(1 | Beach_ID) + (1 | Week)
Data: main

AIC BIC logLik deviance df.resid
5439.4 5466.0 -2713.7 5427.4 614

Scaled residuals:
Min 1Q Median 3Q Max
-1.4383 -0.7242 -0.2287 0.4866 4.0531

Random effects:
Groups Name Variance Std.Dev.
Week (Intercept) 0.003092 0.0556
Beach_ID (Intercept) 0.025894 0.1609
Number of obs: 620, groups: Week, 31; Beach_ID, 8

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.38864 0.08227 -16.879 < 2e-16 ***
Relocation..Y.N.Y 0.32105 0.09152 3.508 0.000452 ***
SPL 0.22218 0.08793 2.527 0.011508 *
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Correlation of Fixed Effects:
(Intr) R..Y.N
Rlct..Y.N.Y -0.143
SPL -0.540 -0.038

Thanks,

Alessandra

On Tue, Feb 11, 2020 at 7:29 PM Ben Bolker <bbolker at gmail.com
<mailto:bbolker at gmail.com>> wrote:

Short answer: if emergence success gets much above 50%, then

the

approximation you're making (Poisson + offset for binomial, or

NB +

offset for negative binomial) doesn't work well. You might

try a

beta-binomial (with glmmTMB) or a binomial + an

observation-level random

effect.

(On the other hand, your intercept is -0.3, which

corresponds to a

baseline emergence of 0.42 - not *very* high (but some beaches

and years

will be well above that ...)

Beyond that, are there any obvious patterns of mis-fit in the
predicted values ... ?

On 2020-02-11 8:09 p.m., Alessandra Bielli wrote:

Dear list

I am fitting a poisson model to estimate the effect of a

treatment on

emergence success of hatchlings. To estimate emergence

success, I use

number of emerged and an offset(log(total clutch).

However, overdispersion was detected:

overdisp_fun(m.emerged) #overdispersion detected

chisq ratio rdf p
3490.300836 5.684529 614.000000 0.000000

Therefore, I switched to a negative binomial. I know

overdispersion is

not

relevant for nb models, but the model plots don't look too

good. I also

tried to fit a poisson model with OLRE, but still the plots

don't look

good.
How do I know if my model is good enough, and what can I do

to improve

it?

summary(m.emerged)

Generalized linear mixed model fit by maximum likelihood

(Laplace

Approximation) ['glmerMod']
Family: Negative Binomial(7.604) ( log )
Formula: Hatched ~ Relocation..Y.N. + SP +

offset(log(Total_Clutch)) + (1

|Beach_ID) + (1 | Year)
Data: main

AIC BIC logLik deviance df.resid
6015.6 6042.2 -3001.8 6003.6 614

Scaled residuals:
Min 1Q Median 3Q Max
-2.6427 -0.3790 0.1790 0.5242 1.6583

Random effects:
Groups Name Variance Std.Dev.
Beach_ID (Intercept) 0.004438 0.06662
Year (Intercept) 0.001640 0.04050
Number of obs: 620, groups: Beach_ID, 8; Year, 5

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.29915 0.04055 -7.377 1.62e-13 ***
Relocation..Y.N.Y -0.16402 0.05052 -3.247 0.00117 **
SPL -0.08311 0.04365 -1.904 0.05689 .
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ?

1

Correlation of Fixed Effects:
(Intr) R..Y.N
Rlct..Y.N.Y -0.114
SPL -0.497 -0.054


Thanks for your help,

Alessandra

<mailto:R-sig-mixed-models at r-project.org> mailing list

<mailto:R-sig-mixed-models at r-project.org> mailing list

[[alternative HTML version deleted]]