I have a data set which includes seed set data from treatment and control plants, I want to include mean flower number as a random effect (derived from daily counts of each plant throughout the course of an experiment)(In the model: scaled_mean_flower_number
). Seed set should, of course, be somewhat (or strongly) dependent on the number of flowers produced by an individual, and there is a good deal of variation between individuals. Individual plants are the experimental units (in the model: plant_ID
).
Note: there is NO correlation between mean flower number and treatment, thankfully.
Though flower number was taken each day, I am using the mean flower count for each plant, and the response (seed_count
) is not a repeated measure (it was, of course, measured once per plant at the end of the season), so nesting scaled_mean_flower_number
in plant_ID
seems incorrect. However, I think that I do need to include both, as scaled_mean_flower_number
is a function of plant_ID
, though I might be wrong about this.
My questions are:
1: Should I include both plant_ID
and scaled_mean_flower_number
, or only one or the other?
2: If I include both, what is the appropriate way to structure the random effect term?
3: Which (if any) of these are actually correctly constructed?
Thank you in advance for any help you can offer.
Models
mod1 = glmer(data=data, seed_count ~ nectar_treatment * scaled_mean_flower_number + (1 + scaled_mean_flower_number | plant_ID), family= "poisson")
mod2 = glmer(data=data, seed_count ~ nectar_treatment + (1|plant_ID/scaled_mean_flower_number), family= "poisson")
mod3 = glmer(data=data, seed_count ~ nectar_treatment + (1|scaled_mean_flower_number), family= "poisson")
mod4 = glmer(data=data, seed_count ~ nectar_treatment + (1|plant_ID/scaled_mean_flower_number), family= "poisson")
mod5 = glmer(data=data, seed_count ~ nectar_treatment + scaled_mean_flower_number + (1|plant_ID), family= "poisson")
Output (summary(x))
>mod1: Error: number of observations (=90) < number of random effects (=180) for term (1 + scaled_mean_flower_number | plant_ID); the random-effects parameters are probably unidentifiable
########################################################################
>mod2: Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: seed_count ~ nectar_treatment + (1 | plant_ID/scaled_mean_flower_number)
Data: data
AIC BIC logLik deviance df.resid
844.5 854.5 -418.2 836.5 86
Scaled residuals:
Min 1Q Median 3Q Max
-0.96554 -0.07419 0.03030 0.05863 0.06804
Random effects:
Groups Name Variance Std.Dev.
scaled_mean_flower_number:plant_ID (Intercept) 1.25 1.118
plant_ID (Intercept) 1.73 1.315
Number of obs: 90, groups: scaled_mean_flower_number:plant_ID, 90; plant_ID, 90
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7079 0.2662 10.173 <2e-16 ***
nectar_treatmenttreatment -0.1310 0.3753 -0.349 0.727
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
nctr_trtmnt -0.700
########################################################################
>mod3: Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: seed_count ~ nectar_treatment + (1 | scaled_mean_flower_number)
Data: data
AIC BIC logLik deviance df.resid
867.1 874.6 -430.6 861.1 87
Scaled residuals:
Min 1Q Median 3Q Max
-3.7063 -0.1027 0.0313 0.0582 3.7442
Random effects:
Groups Name Variance Std.Dev.
scaled_mean_flower_number (Intercept) 3.045 1.745
Number of obs: 90, groups: scaled_mean_flower_number, 88
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7071 0.2721 9.948 <2e-16 ***
nectar_treatmenttreatment -0.1162 0.3836 -0.303 0.762
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
nctr_trtmnt -0.700
#########################################################################
>mod4: Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: seed_count ~ nectar_treatment + (1 | plant_ID/scaled_mean_flower_number)
Data: data
AIC BIC logLik deviance df.resid
844.5 854.5 -418.2 836.5 86
Scaled residuals:
Min 1Q Median 3Q Max
-0.96554 -0.07419 0.03030 0.05863 0.06804
Random effects:
Groups Name Variance Std.Dev.
scaled_mean_flower_number:plant_ID (Intercept) 1.25 1.118
plant_ID (Intercept) 1.73 1.315
Number of obs: 90, groups: scaled_mean_flower_number:plant_ID, 90; plant_ID, 90
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7079 0.2662 10.173 <2e-16 ***
nectar_treatmenttreatment -0.1310 0.3753 -0.349 0.727
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
nctr_trtmnt -0.700
########################################################################
>mod5: Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: seed_count ~ nectar_treatment + scaled_mean_flower_number + (1 | plant_ID)
Data: data
AIC BIC logLik deviance df.resid
814.6 824.6 -403.3 806.6 86
Scaled residuals:
Min 1Q Median 3Q Max
-1.16794 -0.04656 0.03795 0.07938 0.12562
Random effects:
Groups Name Variance Std.Dev.
plant_ID (Intercept) 1.879 1.371
Number of obs: 90, groups: plant_ID, 90
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.7227 0.2143 12.705 < 2e-16 ***
nectar_treatmenttreatment -0.1901 0.3012 -0.631 0.528
scaled_mean_flower_number 0.8789 0.1468 5.988 2.13e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) nctr_t
nctr_trtmnt -0.693
scld_mn_fl_ -0.061 -0.043