Factors in nlme: singularity in backsolve error

Question

I am parameterizing exponential fits for some metabolic scaling models. I have done this in lmer already, without problem, using logged dependent and independent variables. However, I would now like to incorporate other parameters that aren't necessarily exponentially related to the dependent variable. Hence, I've turned to nlme (lme4::nlmer doesn't seem to handle fixed effects), but I don't have much experience with it. Apologies in advance for newbie mistakes.

With the code below, I am getting the following error. I'm guessing that it has something to do with the 'site' factor being misspecified:

Error in nlme.formula(dep ~ scaling_fun(alpha, beta, ind, site), data = scale_df,  :
  Singularity in backsolve at level 0, block 1

When I fit a simpler function that does not involve 'site', the model seems to work correctly.

Any thoughts would be greatly appreciated!

Thanks, Allie

# dput for data
# copy from http://pastebin.com/WNHhi2kZ (too large to include here)

> head(scale_df)
          dep        ind spp site
2  0.28069471 -0.0322841 157    A
3 -0.69719050 -1.2568901 183    A
4  0.29252012  0.1592420 246    A
5  0.72030740 -0.3282789 154    A
6 -0.08601891  0.3623756 110    A
7  0.30793594  0.2230840 154    A


scaling_fun <- function(alpha, beta, ind, site) {
    return(beta + ind^alpha + site*(ind^alpha))
}

# results in singularity in backsolve error

nlme(dep ~ trait_scaling_fun(alpha, beta, ind, site),
     data = scale_df,
     fixed = list(alpha + beta + site ~ 1), random = alpha ~ 1|spp,
     start = list(fixed = c(0.7, 0, 1)))


##############################
# simpler function converges #
##############################

scaling_fun <- function(alpha, beta, ind) {
    return(beta + ind^alpha)
}

nlme(dep ~ scaling_fun(alpha, beta, ind),
     data = scale_df,
     fixed = list(alpha + beta ~ 1), random = alpha ~ 1|spp,
     start = list(fixed = c(0.7, 0)))

Roland · Answer 1 · 2016-07-25T15:52:21.313

Your model does not really make sense since site is a factor variable (and not a parameter). I suspect you actually want to stratify alpha by site:

library(nlme)

scaling_fun <- function(alpha, beta, ind) {
  return(beta + ind^alpha)
}

nlme(dep ~ scaling_fun(alpha, beta, ind),
     data = scale_df,
     fixed = list(alpha ~ site, beta ~ 1), random = alpha ~ 1|spp,
     start = list(fixed = c(0.487, rep(0, 19), -0.3)))
#Nonlinear mixed-effects model fit by maximum likelihood
#  Model: dep ~ scaling_fun(alpha, beta, ind) 
#  Data: scale_df 
#  Log-likelihood: -716.4634
#  Fixed: list(alpha ~ site, beta ~ 1) 
#alpha.(Intercept)       alpha.siteB       alpha.siteC       alpha.siteD       alpha.siteE 
#       0.57671912       -0.61258632       -0.59244337       -0.25793558       -0.24572998 
#      alpha.siteF       alpha.siteG       alpha.siteH       alpha.siteI       alpha.siteJ 
#      -0.23615274       -0.31015393        0.17970575        0.01286117       -0.12539377 
#      alpha.siteK       alpha.siteL       alpha.siteM       alpha.siteN       alpha.siteO 
#       3.72445972       -0.08560994        0.13636185        0.31877456       -0.25952204 
#      alpha.siteQ       alpha.siteR       alpha.siteS       alpha.siteT       alpha.siteU 
#       0.15663989        0.66511079        0.10785082       -0.21547379       -0.23656126 
#             beta 
#      -0.30280707 
#
#Random effects:
# Formula: alpha ~ 1 | spp
#        alpha.(Intercept)  Residual
#StdDev:         0.6426563 0.4345844
#
#Number of Observations: 1031
#Number of Groups: 279

However, I also suspect that site should be a random effect.

Thanks Roland - this makes good sense. I'm interested in the following mathematical formulation: D = (β_1 + β_site + β_spp) • I^(α_1 + α_site + α_spp), where β_1 is the grand mean β term, β_site is the per-site deviation from that mean, and β_spp is a random term (same for α's). To that end, I code site with deviance (sum to zero) contrasts [e.g. contrasts(site) = contr.sum(levels(site))]. Using your scaling_fun, I run nlme with: fixed = list(alpha + beta ~ site), random = alpha + beta ~ 1|spp. I then interpret alpha.(Intercept) as α_1 and beta.(Intercept) as β_1. Do you agree? — Alexander Shenkin, Jul 28 '16 at 10:36
btw, I'm actually interested in the values of the per-site deviations, and hence, I've coded them as fixed effects. — Alexander Shenkin, Jul 28 '16 at 10:40
You know that you can extract random effects too and not only the fixed coefficients? — Roland, Jul 28 '16 at 11:03
Yes, but if i'm actually interested in the site deviations, fitting them as random effects constrains that fit to some extent (I think), and furthermore, doesn't 'count' those levels as degrees of freedom. So in a sense, it would seem to suggest that you have more statistical power that you actually do. My sense is that random effects should really be used for factors that you need to control for, but for whose actual levels you don't care about (aside from perhaps variance and other descriptions of the distribution of those levels...). Also, nlme doesn't like crossed random effects... — Alexander Shenkin, Jul 29 '16 at 07:57

Factors in nlme: singularity in backsolve error

1 Answers1