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$
10
  • 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
    Commented Feb 8, 2017 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
    Commented Feb 8, 2017 at 21:51
  • 1
    $\begingroup$ @Stefan, thank you very much for those links. They are really very helpful! $\endgroup$
    – JKO
    Commented Feb 9, 2017 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
    Commented Feb 9, 2017 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
    Commented Feb 9, 2017 at 19:51

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.