p-values near 1.000 for all tests in a glm() (count and binary data, Poisson and binomial distributions)

Question

I want to run multiple GLM's to answer a few questions about changes in diet between species and over time.

Questions:

1. Does frequency of occurrence (FO) of pieces eaten differ between species or year?
2. Does number of pieces eaten (num.eaten) differ between species or year?
3. Does mass of pieces eaten (mass.eaten) differ between species or year?

And then I also want to answer:

4. Does frequency of occurrence (FO) differ between sex of species AD in year 2000, 2001 or 2011?
5. Does number of pieces eaten (num.eaten) differ between sex of species AD in year 2000, 2001 or 2011?
6. Does mass of pieces eaten (mass.eaten) differ between sex of species AD in year 2000, 2001 or 2011? Then I want to answer these same questions for species CD.

Example of data

table <- "year species   age sex num.eaten mass.eaten FO
1  2000    AB adult   f           10          0.23  1
2  2000    AB adult   m            0            NA  0
3  2001    AB adult   f            0            NA  0
4  2001    CD adult   m            2          0.01  1
5  2001    CD adult   f            6          0.12  1
6  2011    AB adult   m            5          0.01  1
7  2011    AB adult   m            5          0.06  1
8  2011    CD adult   f           20          0.28  1
9  2011    CD adult   m           14          0.36  1
10 2011    CD adult   f           11          0.46  1"

df <- read.table(text=table, header = TRUE)
df

Attempt

For the first 3 questions, I assumed I should do a GLM with species and year as factors, with a species * year interaction term.

For frequency of occurrence (FO), data are binary so follow a binomial distribution, thus:

FO.glm <- glm(FO ~ species * year,
                        data = dat, family=binomial())
summary(FO)

For num.eaten, data are counts so follow a Poisson distribution, thus:

num.eaten.glm <- glm(num.eaten ~ species * year,
                        data = dat, family=poisson())
summary(num.eaten)

For mass.eaten, data are also follow a Poisson distribution, thus:

mass.eaten.glm <- glm(mass.eaten ~ species * year,
                        data = dat, family=poisson())
summary(mass.eaten)

In the output from these 3 GLM's, I am getting p-values that range from 0.992 to 1.000 for each factor and interactions. I know there is a significant difference between some of my data (maybe not in this reproduced example, but in my real data). So I'm trying to determine what error I am making here.

In terms of questions 4-6, I cannot figure out how to run a GLM that will assess only one species at a time, and will test the difference of FO (then num.eaten, then mass.eaten) between sexes in year 2000 (then 2001, then 2011).

I am relatively new to R, and also new to GLM's. Any help/clarification would be appreciated.

Answer 1

When I read this question:"Does frequency of occurrence (FO) of pieces eaten differ between species or year?", I do not immediately think of using a 2-factor interaction model, but rather will have thought of testing two different models, one with "year" and one with "species". Furthermore the description of the variable "frequency of occurrence" makes me think, not of a binomial data situation, but rather of a Poisson situation.

> FO.glm <- glm(FO ~  year,
+               data = df, family=poisson())
> summary(FO.glm)

Call:
glm(formula = FO ~ year, family = poisson(), data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.10336  -0.00316  -0.00316   0.34311   0.50085  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -100.46255  141.00713  -0.712    0.476
year           0.04996    0.07025   0.711    0.477

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3.5703  on 9  degrees of freedom
Residual deviance: 3.0469  on 8  degrees of freedom
AIC: 23.047

Number of Fisher Scoring iterations: 5

> FO.glm2 <- glm(FO ~  species,
+               data = df, family=poisson())
> summary(FO.glm2)

Call:
glm(formula = FO ~ species, family = poisson(), data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0955   0.0000   0.0000   0.3531   0.4708  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.5108     0.5774  -0.885    0.376
speciesCD     0.5108     0.7303   0.699    0.484

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3.5703  on 9  degrees of freedom
Residual deviance: 3.0650  on 8  degrees of freedom
AIC: 23.065

Number of Fisher Scoring iterations: 5

Notice that when aggregated, the counts of FO in cells is not binomial:

with(df, table(FO, year, species) )
, , species = AB

   year
FO  2000 2001 2011
  0    1    1    0
  1    1    0    2

, , species = CD

   year
FO  2000 2001 2011
  0    0    0    0
  1    0    2    3

So this basically produces an equivalent set of conclusions, i.e., no significance. The reason you were getting p-values near unity was that you were building a model with data that had extremely small values. In such situation yopu might consider exact tests, but these data would not be likely to give a materially "improved" result. Note also that the model is being set up with "year" as a linear term whereas you might want to use factor(year) if you suspected a non-linear relationship of FO w.r.t year.