2
$\begingroup$

I am not asking for help in the coding unless that may resolve the issue. I am wondering if this is even possible from a statistics standpoint and if it is, how I go about resolving it, because all my attempts have failed and this is a fairly niche statistical question. I am not sure if my issue is in my syntax or with the ability of this experimental design to be analyzed by a Poisson MLM.

Thus the warnings and errors returned by the code below may help determine which direction this goes: Coding/syntax error(s) or issue between my experimental structure and the statistics themselves. The coding and syntax errors seem unlikely to me given that I have spent 40+ hours attempting to resolve this on my own with much more than the code and many corrections far beyond what has been attempted below. Individuals with statistical backgrounds would be more apt to help with this determination, hence posting here and not SO, or elsewhere because this question is about both application of statistics to the experimental structure and statistical computing as the issue may lie in the underlying math. Both are in the purview of the CV community.

I have a set of count data that is repeated measures. I am using lme4::GLMER with a log link function (Poisson GLMER). Here are my variables and their attributes:

Variables:

Age: Categorical, 2 levels (Adolescent, Adult)

Condition: Categorical, 2 levels (IAE, CTRL)

  • CTRL never had access to ethanol; Total Ethanol Consumed was set to zero so the NA-values for Total Ethanol Consumed would not just get dropped from the model.
  • IAE had access and could choose to drink ethanol

Concentration: Continuous, repeated measure

Total Ethanol Consumed: Continuous, nested in Condition

Y` = Total.Aversive, count, Continuous

Experimental Structure:

Experiment Structure

What I Have Tried:

I have tried a number of times to specify in the model that Total Ethanol Consumed is nested within Condition, however any time I attempt to perform my GLMER it does not work and I cannot use emmeans:emmeans or emtrends to get estimates or means from the CTRL group even when I specify the nesting. emmeans and emtrends seem to think the nesting is with age and not c.totaletoh. This implies an issue with the structure of the model, but from what I can tell, everything is properly structured. Does the issue lie with the ability of the variables in the model to carry the information and not be collinear?

###Variable Coding Adjustment###

#adjust Age to contrast coding
contrasts(mydetoh$Age)=contr.sum(2)
contrasts(mydetoh$Age)
#Returns:
#           [,1]
#Adolescent    1
#Adult        -1

#adjust Condition to contrast coding
contrasts(mydetoh$Condition)=contr.sum(2)
contrasts(mydetoh$Condition)
#Returns:
#     [,1]
#CTRL    1
#EtOH   -1


####Testing Adding c.totaletoh (Total Ethanol Consumed) nested within Condition to the fixed effect structure.
#Here c.conc is centered Concentration.
#Total Ethanol Consumed was centered after all CTRLs were assigned the value of zero

#Total Ethanol Consumed during IAE (c.totale) used as a "nested predictor variable" nested in Condition="EtOH" in the fixed effects structure)

NEavers <-glmer(Total.Aversive ~ Age*c.conc*(Condition/c.totaletoh) + (c.conc|RatID), data=mydetoh, family=poisson)
summary(NEavers)

# Model did not converge, used code below to extend # of iterations and start from where the previous model left off.
#write over old model, if it doesnt converge, it doesn't matter anyway, can just rerun if necessary
ss1 <- getME(NEavers,c("theta","fixef")) #gets values from old model and puts them into an object
NEavers <- update(NEavers,start=ss1,control=glmerControl(optCtrl=list(maxfun=2e6))) #restarts the model from leave-off point
summary(NEavers)
#CHECK TO SEE IF MODEL CONVERGED w/in GRADIENT (in this case it does, the warning goes away)

#####OUTPUT#####
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: Total.Aversive ~ Age * c.conc * (Condition/c.totaletoh) + (c.conc |      RatID)
   Data: mydetoh
Control: glmerControl(optCtrl = list(maxfun = 2e+06))

     AIC      BIC   logLik deviance df.resid 
  1878.2   1926.6   -924.1   1848.2      171 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.4261 -1.1942 -0.2747  0.7115  7.9505 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 RatID  (Intercept) 0.4311   0.6566        
        c.conc      9.3194   3.0528   -0.20
Number of obs: 186, groups:  RatID, 63

Fixed effects:
                                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            2.747463   0.097973  28.043  < 2e-16 ***
Age1                                   0.226499   0.097553   2.322   0.0202 *  
c.conc                                 2.735081   0.471254   5.804 6.48e-09 ***
Condition1                             0.453746   0.097774   4.641 3.47e-06 ***
Age1:c.conc                            0.174713   0.467186   0.374   0.7084    
ConditionEtOH:c.totaletoh             -0.024920   0.004534  -5.496 3.88e-08 ***
Age1:Condition1                       -0.175511   0.097549  -1.799   0.0720 .  
c.conc:Condition1                      0.292815   0.469522   0.624   0.5329    
Age1:ConditionEtOH:c.totaletoh         0.002011   0.004520   0.445   0.6564    
c.conc:ConditionEtOH:c.totaletoh      -0.002285   0.022907  -0.100   0.9206    
Age1:c.conc:Condition1                -0.469138   0.467162  -1.004   0.3153    
Age1:c.conc:ConditionEtOH:c.totaletoh -0.042898   0.022891  -1.874   0.0609 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Age1   c.conc Cndtn1 Ag1:c. CEOH:. Ag1:C1 c.c:C1 A1:CEO c.:CEO A1:.:C1
Age1         0.088                                                                       
c.conc      -0.212 -0.014                                                                
Condition1   0.172 -0.098 -0.035                                                         
Age1:c.conc -0.014 -0.209  0.077  0.024                                                  
CndtnEtOH:. -0.221 -0.214  0.022  0.223  0.051                                           
Age1:Cndtn1 -0.099  0.180  0.025  0.090 -0.039  0.213                                    
c.cnc:Cndt1 -0.035  0.024  0.139 -0.209 -0.091 -0.024 -0.015                             
Ag1:CnEOH:. -0.213 -0.228  0.054  0.213  0.028 -0.094  0.228 -0.054                      
c.cn:CEOH:.  0.022  0.049 -0.195 -0.024 -0.223 -0.236 -0.048  0.198  0.096               
Ag1:c.cn:C1  0.025 -0.039 -0.092 -0.014  0.154 -0.051 -0.209  0.079 -0.028  0.222        
A1:.:CEOH:.  0.050  0.026 -0.221 -0.050 -0.205  0.099 -0.026  0.222 -0.233 -0.120  0.205 
fit warnings:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
#####END OUTPUT#####

#Seems fine above until you try to look at the marginal means:
#Age x Control emmeans
NOTE: A nesting structure was detected in the fitted model:
    Age %in% Condition
NOTE: Results may be misleading due to involvement in interactions

emm.NEavers.axc
Condition = CTRL:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent nonEst    NA  NA        NA        NA
 Adult      nonEst    NA  NA        NA        NA

Condition = EtOH:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent   2.70 0.197 Inf      2.31      3.08
 Adult        1.89 0.156 Inf      1.59      2.20

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
#No estimate for controls


#Age x Control x Total Ethanol
emm.NEavers.axcxt <- emmeans(NEavers, ~ Age|Condition*c.totaletoh)
NOTE: A nesting structure was detected in the fitted model:
    Age %in% Condition
NOTE: Results may be misleading due to involvement in interactions

emm.NEavers.axcxt
Condition = CTRL, c.totaletoh = 9.88e-16:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent nonEst    NA  NA        NA        NA
 Adult      nonEst    NA  NA        NA        NA

Condition = EtOH, c.totaletoh = 9.88e-16:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent   2.70 0.197 Inf      2.31      3.08
 Adult        1.89 0.156 Inf      1.59      2.20

Results are given on the log (not the response) scale. 
Confidence level used: 0.95


#Age x Condition trends for c.totaletoh (which I wouldn't expect CTRLs to have an estimate for, or to have a 0)
emt.NEavers.axc.t <- emtrends(NEavers, ~ Age|Condition, var="c.totaletoh")
NOTE: A nesting structure was detected in the fitted model:
    Age %in% Condition
NOTE: Results may be misleading due to involvement in interactions

emt.NEavers.axc.t
Condition = CTRL:
 Age        c.totaletoh.trend     SE  df asymp.LCL asymp.UCL
 Adolescent            nonEst     NA  NA        NA        NA
 Adult                 nonEst     NA  NA        NA        NA

Condition = EtOH:
 Age        c.totaletoh.trend     SE  df asymp.LCL asymp.UCL
 Adolescent           -0.0229 0.0061 Inf   -0.0349   -0.0110
 Adult                -0.0269 0.0067 Inf   -0.0401   -0.0138

Trends are based on the log (transformed) scale 
Confidence level used: 0.95


#This happens even when each part of the model structure formula is spelled out:

##TEST: PRODUCES SAME RESULT AS THE NESTED MODEL STRUCTURE ABOVE.
NEAtest2 <-glmer(Total.Aversive ~ c.conc+
                                  Age+
                                  Condition+
                                  Condition:c.totaletoh+
                                  c.conc:Age+
                                  c.conc:Condition+
                                  c.conc:Condition:c.totaletoh+
                                  Age:Condition+
                                  Age:Condition:c.totaletoh+
                                  c.conc:Age:Condition+
                                  c.conc:Age:Condition:c.totaletoh
                                  + (c.conc|RatID), data=mydetoh, family=poisson)
summary(NEAtest2)
ss2 <- getME(NEAtest2,c("theta","fixef"))
NEAtest2 <- update(NEAtest2,start=ss2,control=glmerControl(optCtrl=list(maxfun=2e6)))
summary(NEAtest2)
emmeans(NEAtest2,~ Age*Condition, weights="show.levels")
#SAME ISSUE AS ABOVE

Checking to see if centering is the issue

#UNCENTERD Total Ethanol Consumed (totaletoh) used as a "nested predictor variable" nested in Condition="EtOH" in the fixed effects structure)

NEavers.uc <-glmer(Total.Aversive ~ Age*c.conc*(Condition/totaletoh) + (c.conc|RatID), data=mydetoh, family=poisson)
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0319191 (tol = 0.001, component 1)

summary(NEavers.uc)

####OUTPUT####
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: Total.Aversive ~ Age * c.conc * (Condition/totaletoh) + (c.conc |      RatID)
   Data: mydetoh

     AIC      BIC   logLik deviance df.resid 
  1878.2   1926.6   -924.1   1848.2      171 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.4256 -1.1949 -0.2747  0.7115  7.9508 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 RatID  (Intercept) 0.4307   0.6563        
        c.conc      9.3011   3.0498   -0.20
Number of obs: 186, groups:  RatID, 63

Fixed effects:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          3.098333   0.128137  24.180  < 2e-16 ***
Age1                                 0.198664   0.127969   1.552   0.1206    
c.conc                               2.767484   0.620158   4.463 8.10e-06 ***
Condition1                           0.102979   0.128073   0.804   0.4214    
Age1:c.conc                          0.779974   0.618956   1.260   0.2076    
ConditionEtOH:totaletoh             -0.024920   0.004532  -5.499 3.83e-08 ***
Age1:Condition1                     -0.147588   0.127973  -1.153   0.2488    
c.conc:Condition1                    0.259891   0.619365   0.420   0.6748    
Age1:ConditionEtOH:totaletoh         0.001994   0.004518   0.441   0.6590    
c.conc:ConditionEtOH:totaletoh      -0.002334   0.022889  -0.102   0.9188    
Age1:c.conc:Condition1              -1.071302   0.618931  -1.731   0.0835 .  
Age1:c.conc:ConditionEtOH:totaletoh -0.042886   0.022874  -1.875   0.0608 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Age1   c.conc Cndtn1 Ag1:c. CnEOH: Ag1:C1 c.c:C1 A1:CEO c.:CEO A1:.:C1
Age1         0.190                                                                       
c.conc      -0.199 -0.022                                                                
Condition1  -0.316 -0.195  0.058                                                         
Age1:c.conc -0.021 -0.198  0.186  0.027                                                  
CndtnEtOH:t -0.667 -0.116  0.138  0.668 -0.013                                           
Age1:Cndtn1 -0.196 -0.315  0.028  0.190  0.058  0.116                                    
c.cnc:Cndt1  0.058  0.028 -0.345 -0.198 -0.194 -0.140 -0.022                             
Ag1:CndEOH: -0.116 -0.670 -0.009  0.116  0.140 -0.094  0.671  0.010                      
c.cnc:CEOH:  0.133 -0.011 -0.668 -0.134 -0.106 -0.234  0.012  0.670  0.096               
Ag1:c.cn:C1  0.027  0.058 -0.195 -0.021 -0.344  0.013 -0.199  0.187 -0.141  0.105        
Ag1:.:CEOH: -0.011  0.134 -0.106  0.012 -0.675  0.099 -0.134  0.106 -0.231 -0.120  0.675 
fit warnings:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
convergence code: 0
Model failed to converge with max|grad| = 0.0319191 (tol = 0.001, component 1)
####END OUTPUT####

# Model did not converge, used code below to extend # of iterations as far as they can go and start from where the previous model left off.

ss1.1 <- getME(NEavers.uc,c("theta","fixef")) #gets values from old model and puts them into an object
NEavers.uc <- update(NEavers.uc,start=ss1.1,control=glmerControl(optCtrl=list(maxfun=2e9))) #restarts the model from leave-off point
#model still did not converge:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.00100564 (tol = 0.001, component 1)

summary(NEavers.uc)
####OUTPUT####
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: Total.Aversive ~ Age * c.conc * (Condition/totaletoh) + (c.conc |      RatID)
   Data: mydetoh
Control: glmerControl(optCtrl = list(maxfun = 2e+09))

     AIC      BIC   logLik deviance df.resid 
  1878.2   1926.6   -924.1   1848.2      171 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.4261 -1.1942 -0.2747  0.7114  7.9505 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 RatID  (Intercept) 0.4311   0.6566        
        c.conc      9.3196   3.0528   -0.20
Number of obs: 186, groups:  RatID, 63

Fixed effects:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          3.098093   0.128192  24.168  < 2e-16 ***
Age1                                 0.198213   0.128020   1.548   0.1215    
c.conc                               2.767154   0.620717   4.458 8.27e-06 ***
Condition1                           0.103117   0.128127   0.805   0.4209    
Age1:c.conc                          0.778274   0.619484   1.256   0.2090    
ConditionEtOH:totaletoh             -0.024920   0.004534  -5.496 3.88e-08 ***
Age1:Condition1                     -0.147239   0.128026  -1.150   0.2501    
c.conc:Condition1                    0.260672   0.619929   0.420   0.6741    
Age1:ConditionEtOH:totaletoh         0.002010   0.004520   0.445   0.6565    
c.conc:ConditionEtOH:totaletoh      -0.002283   0.022907  -0.100   0.9206    
Age1:c.conc:Condition1              -1.072659   0.619496  -1.732   0.0834 .  
Age1:c.conc:ConditionEtOH:totaletoh -0.042897   0.022892  -1.874   0.0609 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Age1   c.conc Cndtn1 Ag1:c. CnEOH: Ag1:C1 c.c:C1 A1:CEO c.:CEO A1:.:C1
Age1         0.190                                                                       
c.conc      -0.201 -0.023                                                                
Condition1  -0.316 -0.195  0.058                                                         
Age1:c.conc -0.022 -0.200  0.186  0.027                                                  
CndtnEtOH:t -0.667 -0.116  0.139  0.668 -0.013                                           
Age1:Cndtn1 -0.196 -0.315  0.029  0.190  0.058  0.116                                    
c.cnc:Cndt1  0.059  0.028 -0.345 -0.200 -0.194 -0.141 -0.023                             
Ag1:CndEOH: -0.116 -0.670 -0.009  0.116  0.142 -0.094  0.671  0.009                      
c.cnc:CEOH:  0.134 -0.011 -0.668 -0.135 -0.106 -0.236  0.011  0.670  0.096               
Ag1:c.cn:C1  0.028  0.058 -0.195 -0.021 -0.344  0.013 -0.201  0.187 -0.142  0.105        
Ag1:.:CEOH: -0.011  0.135 -0.106  0.011 -0.675  0.099 -0.135  0.106 -0.233 -0.120  0.675 
fit warnings:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
convergence code: 0
Model failed to converge with max|grad| = 0.00100564 (tol = 0.001, component 1)

#Convergence failure is unresolved. the centering of the totaletoh variable helped.
####END OUTPUT####

Edit based on rvl's comments:

#Tried to call table(NEavers), did not work; 
#tried as.data.frame.table(), error below. I've never used table function
as.data.frame.table(emm.NEavers.axc, row.names = NULL, stringsAsFactors = TRUE, 
                    sep = "", base = list(LETTERS))

Error in dimnames(x) <- dnx : 'dimnames' applied to non-array

Trying emmeans() call with covnest = TRUE

Using mean-centered totaletoh because uncentered totaletoh did not work. c.totaletoh has no 0 values once mean centered.

emm.NEavers.axc <- emmeans(NEavers, ~ Age|Condition, covnest = TRUE)
NOTE: A nesting structure was detected in the fitted model:
    Age %in% Condition
NOTE: Results may be misleading due to involvement in interactions

emm.NEavers.axc
Condition = CTRL:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent nonEst    NA  NA        NA        NA
 Adult      nonEst    NA  NA        NA        NA

Condition = EtOH:
 Age        emmean    SE  df asymp.LCL asymp.UCL
 Adolescent   2.70 0.197 Inf      2.31      3.08
 Adult        1.89 0.156 Inf      1.59      2.20

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

Can anyone explain to me why this will not work? Is it the statistic itself? Is it the nesting structure? Is it something else I am not seeing? (Edit: when I tried to see if it might work in the random effects structure, it resulted in a singularity due to 0 variance with the CTRL which was reported in a previous edit). My best guess is that the fixed effects structure did not carry the information from the nested variable correctly, but I am relatively new to nesting. Also, I am completely stuck.

If anyone knows why this doesn't work, I'd love an explanation. If anyone knows how to get it to work, please let me know. It needs to be a nested FIXED effect structure as that is the setup of the data.

$\endgroup$

closed as off-topic by mkt, user158565, Michael Chernick, Peter Flom Aug 6 at 15:41

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – user158565, Michael Chernick, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ in your figure both blocks are labeled "adult" - should one be "adolescent"? I don't think this question is specific to mixed models. I think you want to set "ethanol consumed" to zero for all control individuals, so that the "control vs treatment" comparison will contract treatment individuals who consumed no ethanol with control individuals ... $\endgroup$ – Ben Bolker Aug 6 at 3:00
  • 1
    $\begingroup$ First, it appears that the variable totaletoh is a numeric predictor. It doesn't make sense to have a continuous predictor as a nested effect. Notice that the emmeans results with that as a by variable just output results for the mean totaleoh. Second, emmeans detects nesting based on either model structure or data structure. In this case it appears that it is the data structure. I think via a table() call you'd verify that the levels of Age are entirely distinct for each Condition. Finally, it is a mystery how you can be showing emmeans results when that call errored. $\endgroup$ – rvl Aug 7 at 18:35
  • 1
    $\begingroup$ I don't understand this: "CTRL never had access to ethanol, so they could never drink, but this was set to zero anyway so the NA-values would not just get dropped." What is "this"? You're not saying CTRL was set to zero, so are you saying totaletoh was set to 0? $\endgroup$ – rvl Aug 7 at 18:50
  • 1
    $\begingroup$ You need a twiddle in the emmeans call: emmeans(NEavers, ~ Age | Condition) $\endgroup$ – rvl Aug 7 at 18:56
  • 1
    $\begingroup$ You might try adding covnest = TRUE to the emmeans() call. That includes covariates in the nesting-detection calculations. I had forgotten I had added this feature, and it kind of flies in the face of my "first" comment earlier. It won't work right though if there are totaletoh values of 0 in other than CTRL cases. $\endgroup$ – rvl Aug 8 at 14:41