2
$\begingroup$

Can anyone suggest how I parameterize a Poisson random-intercept model, with a natural cubic spline function? I've been using glmer for a while and am happy with how I'm specifying the main fixed effects and random intercept, but I get scale warnings related to my spline basis.

I'm modelling counts of 'incidents' as dependent variable, predicted by counts of observations in demographic categories (lets just use age for this post), with time period as a natural cubic spline of months with knots every 6 months, and a random intercept for organisations/clusters. I'm assuming that, as I'm using a log link function, I should log-transform my count predictors, and a simplified version is:

mod <- glmer(incidents ~ (1|org_code)
                  + log(age17)
                  + log(age29)
                  + log(age49)
                  + log(age69)
                  + log(age70)
                  + ns(re_month, knots=seq(6, 54, 6))
                  , data=sub
                  , family=poisson(link=log))

I'm afraid I'm unable to share my data, but my 'log(count)' variables are in the range 4 to 12 on log scale, and my ns() spline basis columns are in the range -1 to 1.

I've followed Ben Bolker's trouble shooting article: http://rpubs.com/bbolker/lme4trouble1 and don't have singularity problems, mismatched scaled and absolute gradients etc. I don't think scale is appropriate, as I've already transformed to log scale. Different optimiser give similar results, but still getting:

In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

My questions are:

  • Am I going about this all wrong?

  • Is there another way to parametrize the spline? I've looked at
    gamm4, but I'd rather do it within glmer, as it develops from
    glm models to start with, and it's a bit of leap.

  • Is it reasonable to ignore the warning if I'm otherwise satisfied
    with the convergence? Is it 'safe' to use the AIC and parameter
    estimates from this?

Edit at the request of @IWS

Here is some additional info, as I'm afraid I can't supply the data. Here is are the summary stats for the fixed effects. The random-intercept is for repeated measure (60 in most cases) in 138 clusters. Coded as a factor.

          vars    n    mean      sd median trimmed     mad  min   max range skew kurtosis    se
incidents    1 8114  666.69  362.51  600.0  625.77  306.90    2  2784  2782 1.35     2.71  4.02
age17        3 8114 1245.34 1050.87  985.0 1056.90  553.01   40  8997  8957 2.99    12.68 11.67
age29        4 8114 2260.12 1268.73 1928.0 2095.78  984.45  259  9367  9108 1.39     2.51 14.08
age49        5 8114 4405.80 2502.30 3719.0 4054.05 1848.80  638 15533 14895 1.31     1.60 27.78
age69        6 8114 7303.33 3772.67 6302.5 6881.60 3438.89 1678 19142 17464 0.90     0.15 41.88
re_month     7 8114   30.51   17.33   31.0   30.52   22.24    1    60    59 0.00    -1.21  0.19

re_month is an integer of month from the start of a five year period e.g. 1 = Apr-2011, 2 = May-2011 etc. that is used to construct the spline with knots at 6-month intervals. The ns()1...10 variables below are generated by the call to ns().

(Intercept)   1, 1, 1, 1, 1, 1
log(age17)    5.796058, 6.086775, 5.948035, 5.940171, 5.872118, 5.717028
log(age29)    6.555357, 6.432940, 6.327937, 6.340359, 6.595781, 6.597146
log(age49)    7.189168, 7.183871, 7.262629, 7.117206, 7.180070, 7.176255
log(age69)    7.738924, 7.874359, 7.910591, 7.930566, 7.976595, 7.754053
log(age70)    8.761237, 8.800265, 8.674539, 8.700348, 8.646290, 8.620832
ns()1         0.6355301, 0.5711692, 0.4779412, 0.3700073, 0.2615287, ...
ns()2         0.2615741, 0.3703704, 0.4791667, 0.5740741, 0.6412037, ...
ns()3         0.0007716049, 0.0061728395, 0.0208333333, 0.0493827160,...
ns()4         0, 0, 0, 0, 0, 0
ns()5         0, 0, 0, 0, 0, 0
ns()6         0, 0, 0, 0, 0, 0
ns()7         0, 0, 0, 0, 0, 0
ns()8         -0.0260514200, -0.0133383270, -0.0056271067, -0.0016672...
ns()9         0.0781542599, 0.0400149811, 0.0168813201, 0.0050018726,...
ns()10        -0.0521028399, -0.0266766540, -0.0112542134, -0.0033345...
$\endgroup$
  • 1
    $\begingroup$ For me it is quite difficult to asses whether log-transformation of the independent/predictor variables is necessary without some information about your data (sample size, maybe a few rows of data and some descriptives?). However, I've not heard about a necessity of log transforming continuous variables because of a log link function. AFAIK a (log) transformation of your predictors can be necessary to obtain a linear association with your outcome (on the log scale), but is not mandatory. So my question to you would be: do you also get this error when you use the predictors as they are? $\endgroup$ – IWS Sep 6 '17 at 11:05
  • $\begingroup$ Thanks @IWS. I'll add some summary statistics for the data set in to the question, and the first few rows of my model matrix from glFormula. There are 8114 rows from 138 clusters. Some rows are missing, due to opening/closing organisations and data errors, but they are roughly 60 monthly samples from 138 clusters. I've tried to use the un-transformed variables as you suggest, but I have an additional Model failed to converge with max|grad| = 10.7588 (tol = 0.001, component 1) warning. Also still get the previous warning after several restarts. $\endgroup$ – Mainard Sep 6 '17 at 11:47
  • $\begingroup$ Then I am at a loss. I've upvoted your question and hope someone else might be able to help you. $\endgroup$ – IWS Sep 6 '17 at 11:48
2
$\begingroup$

This is on the line between a statistical and a computational question, but I'll take a shot at it.

tl;dr the bottom line is that if you have achieved similar results with different optimizers, then you can trust that the numerics of your fit are OK, and you don't need to worry about these warnings. But you should check for overdispersion.

  • (agreeing with @IWS) log-transforming your predictors is not necessarily required - this would be done for scientific reasons (e.g. if you had a priori reason to believe that the underlying relationship was a power-law, $y=a x^b$, then you would want to fit a linear relationship $\log(y) = \log(a) + b \log(x)$; if you thought an exponential relationship $y=\exp(a+bx)$ was more reasonable, then you wouldn't log-transform), or (phenomenologically) to improve the linearity of the relationship. If your predictors as well as your response are counts, though, it does seem reasonable to log-transform.
  • if you want to try centering the ageXX parameters on the log scale (which might help a bit for interpretation, although since you have no interactions in your model it will only affect the intercept parameter), you can use log(ageXX/mean(ageX)). I agree that scaling after logging makes little sense.
  • your spline model seems fine. The major advantage of gamm4 is that you wouldn't have to pick the number of knots yourself. It might not be that hard to switch, but I don't think it's necessary. If your data are non-uniformly distributed along re_month you might want to try specifying just the number of knots and let ns() pick the knot locations itself.
  • you should really check for overdispersion, and (if necessary) use either (1) an observation-level random effect (2) lme4::glmer.nb or (3) glmmTMB

The "large eigenvalue" error is triggered if the ratio of the largest to the smallest eigenvalue of the model Hessian (all eigenvalues should be positive) is too large. If you want to investigate further which model components are causing your problem, you can try something like this example:

 library(lme4)
 gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)
 evd <- eigen(gm1@optinfo$derivs$Hessian,symmetric=TRUE)
 ## identify params; this works with a single scalar random effect,
 ## would be more complicated with fancier RE model
 cnames <- c("theta",names(fixef(gm1)))
 dimnames(evd$vectors) <- list(cnames,cnames)
 print(evd,digits=3)
 ## eigen() decomposition
 ## $values
 ## [1] 68.01 51.20 28.51 17.95  9.56

 ## $vectors
 ##               theta (Intercept) period2 period3 period4
 ## theta       0.93012      0.3581  0.0278  0.0513 -0.0567
 ## (Intercept) 0.36178     -0.8559 -0.1011 -0.2506  0.2520
 ## period2     0.04846     -0.3050  0.7562  0.4692 -0.3356
 ## period3     0.03953     -0.2022 -0.6434  0.6410 -0.3643
 ## period4     0.00932     -0.0723 -0.0560 -0.5510 -0.8294

Here you would look at the first eigenvector and hope that something jumped out at you.

$\endgroup$
  • $\begingroup$ Thanks @Ben Bolker for the comprehensive answer. The X values are all counts, so I'll stick with the log transformation. I've applied the suggested centering and it reduces the ratio of eigen values, but it is still very large. The values for theta and intercept are 3 orders of magnitude smaller than other parameters. Will also try an alternative by fitting quantile values for the count variables, rather than binned counts. That won't require the log transformation and will be on the same scale as the spline. The model is very over-dispersed, and I'll try the suggestions above. $\endgroup$ – Mainard Sep 13 '17 at 8:45

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.