0
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ What is the clustering variable? Is scaled_mean_flower_number a factor with discrete levels? I'm having trouble understanding why you need a random effect. $\endgroup$ – Andrew M Feb 8 '17 at 20:59
  • 1
    $\begingroup$ I would suggest providing a working example that explains your data and reproduces the error/warning (see here). Here are a few links that may help:1,2,3, and 4 for some worked examples (inlc. some that give the "converge" warning. $\endgroup$ – Stefan Feb 8 '17 at 21:51
  • 1
    $\begingroup$ @Stefan, thank you very much for those links. They are really very helpful! $\endgroup$ – JKO Feb 9 '17 at 16:26
  • 1
    $\begingroup$ I still can't tell for sure, because the data are not described in enough detail. But the proposed models look inappropriate--you only have one observation per plant_ID. I'm surprised any of them are converging. (I think that the identification, hence convergence is possible due to the mean-variance assumption of the Poisson model.) Just because scaled_mean_flower_number is not fixed by design does not mean that it can't be used as a fixed effect. "Random" covariates can be "fixed" effects with no problem whatsoever (as long as the model is not mis-specified). $\endgroup$ – Andrew M Feb 9 '17 at 19:44
  • 1
    $\begingroup$ On "fixed" vs "random" regressors (aka predictors aka covariates): stats.stackexchange.com/questions/139957/fixed-regressors?rq=1 . (NB: the distinction between a "random effect", which is a parameter that we treat as coming from a probability distribution, and would have been less confusing if we only called it a "variance component") $\endgroup$ – Andrew M Feb 9 '17 at 19:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.