Skip to content

glmmTMB- fitting splines

9 messages · Ben Bolker, dani, John Maindonald +1 more

Original
dani
# 
Hello everyone,


I am working with a glmmTMB model with two random effects. Some of my covariates have non-parametric associations with my dependent variable so I would like to fit splines for them. I am not sure how my code should look like.


Could someone point me towards an example using glmmTMB with splines? I am not really sure how to interpret such a model.


Thanks!

Best regards,

Dani



<http://aka.ms/weboutlook>
Ben Bolker
# 
I don't know of an example offhand, but

https://stats.stackexchange.com/questions/301666/using-splines-in-r-lme4glmer-scale-issues

  gives an example of using splines::ns().  Basically, you can use ns()
as a drop-in term within a formula; unlike the magical s() function in
mgcv, you have to specify the number of knots/degrees of freedom
yourself (splines::ns fits regression splines, mgcv::s fits *penalized*
regression splines).

  Perhaps not known to everyone, mgcv can handle some forms of
zero-inflation (although I think it does ZIP but not ZINB), so
https://www.fromthebottomoftheheap.net/2017/05/04/compare-mgcv-with-glmmTMB/

might also be useful.

Here's an example. It is in principle possible to use
(ns(Days,5)|Subject) as the random effect (i.e. let curves vary among
individuals), but it didn't work in this case -- too complex for this
medium-size data set.

library(glmmTMB)
data(sleepstudy,package="lme4")

library(splines)
m1 <- glmmTMB(Reaction~ns(Days,5)+(1|Subject), data=sleepstudy)
sleepstudy$pred <- predict(m1)

library(ggplot2)
ggplot(sleepstudy,aes(x=Days))+geom_point(aes(y=Reaction))+geom_line(aes(y=pred,group=Subject))
On 2018-05-21 07:36 PM, dani wrote:
dani
# 
Hello Dr Bolker,


Thank you so much for your prompt and helpful answer! This is great!

Best regards,

Dani


Sent from Outlook<http://aka.ms/weboutlook>
dani
# 
Hi John,

Thank you so much! This is very helpful! I managed to run it  but I am not sure how to interpret the results as I get this:


# Conditional model:
#                                                    Estimate Std. Error z value Pr(>|z|)
# (Intercept)                             -8.40461    1.58077  -5.317 1.06e-07 ***
# splines::ns(newage, 2)1     -1.89262    0.57246  -3.306 0.000946 ***
# splines::ns(newage, 2)2      0.10296    0.47268   0.218 0.827575


I am not sure what to make of the two different spline results.

Best regards,

D
John Maindonald
# 
The spline coefficients multiply the two basis terms.  Most times, one wants to
work with predicted values and standard errors.  The predict method seems
not yet to have been implemented for glmmTMB models with a betabinomial
error family.  One can use fitted() to get just the fitted probabilities, and do a
complementary log-log transform (in this instance) to get predictions on the
scale of the linear predictor (NB, linear in the sense that it is a linear combination
of the basis functions, plus intercept).

The output suggests that the first basis term might be enough on its own.
Observe, however, to choose a simple case:
[1] 0.0000000 0.3570466 0.5662628 0.5290951 0.3440969

This is a very nonlinear function of x, quite different from the linear function
of x that one gets by typing splines::ns(x, 1)

Regression thin plate splines, as implemented in mgcv. have the advantage
that the initial basis terms change only very slightly as one moves to a higher
degree of freedom basis.


John Maindonald             email: john.maindonald at anu.edu.au<mailto:john.maindonald at anu.edu.au>
On 22/05/2018, at 14:41, dani <orchidn at live.com<mailto:orchidn at live.com>> wrote:
Hi John,
Thank you so much! This is very helpful! I managed to run it  but I am not sure how to interpret the results as I get this:


# Conditional model:
#                                                    Estimate Std. Error z value Pr(>|z|)
# (Intercept)                             -8.40461    1.58077  -5.317 1.06e-07 ***
# splines::ns(newage, 2)1     -1.89262    0.57246  -3.306 0.000946 ***
# splines::ns(newage, 2)2      0.10296    0.47268   0.218 0.827575


I am not sure what to make of the two different spline results.
Best regards,
D
Farrar, David
# 
I think I used functions from Hmisc, the last time I did regression splines.  
(I did not use penalized splines - I specified enough knots to give an appropriate range of shapes, and used the default knot placements.)
David 

-----Original Message-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of dani
Sent: Monday, May 21, 2018 10:41 PM
To: John Maindonald <john.maindonald at anu.edu.au>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] glmmTMB- fitting splines

Hi John,

Thank you so much! This is very helpful! I managed to run it  but I am not sure how to interpret the results as I get this:


# Conditional model:
#                                                    Estimate Std. Error z value Pr(>|z|)
# (Intercept)                             -8.40461    1.58077  -5.317 1.06e-07 ***
# splines::ns(newage, 2)1     -1.89262    0.57246  -3.306 0.000946 ***
# splines::ns(newage, 2)2      0.10296    0.47268   0.218 0.827575


I am not sure what to make of the two different spline results.

Best regards,

D
dani
# 
Hello again,


Thank you so much for your detailed explanation!


Best,

D


Sent from Outlook<http://aka.ms/weboutlook>
dani
# 
Hi David,


Thanks. I am not sure Hmisc works with random effects - I will check it out, though.


Best,

D


Sent from Outlook<http://aka.ms/weboutlook>