glmmTMB output for t-family

family_params(m1)

On 29 Oct 2024, at 11.48, Henrik Singmann <singmann at gmail.com> wrote:

Hi all,

It seems to me as if the estimate for the df of the t-family in glmmTMB are
somewhat hidden from the user. In particular, the only way I have found to
get to it is via exp(model$fit$par["psi"])

So my question is just if I am missing something obvious here.

To give some background, I am trying to fit some data that looks pretty
leptokurtic so t-distribution  seems appropriate which is luckily available
in glmmTMB. To get going with it, I first wanted to see whether I can
actually recover the df, which seems to work reasonably well for large N,
see code below.

Best,
Henrik

#### Simulation code (somewhat long) follows ###
library("mvtnorm")
library("extraDistr")
library("glmmTMB")

# Sample size
nsubjects <- 100
replicates_cell <- 20

# Fixed effects
ifixed <- 0.5 # intercept
sfixed1 <- 0.5 # slope

# Subject random effects
is <- 0.5 # sd intercept
ss <- 0.5 # sd slope

# Correlations
rs.is <- 0.5 # intercept and slope

# Residual
sigma <- 1.5
phi <- 3.4

set.seed(56771234)
cov.is <- rs.is*is*ss
cov.subjects <- matrix(c(is^2, cov.is, cov.is, ss^2),
                      nrow=2, byrow=TRUE)
# Random sample from bivariate normal, for each subject
re.subjects <- data.frame(cbind(1:nsubjects,
                               rmvnorm(nsubjects, mean=c(rep(0, 2)),
                                       sigma=cov.subjects)))
colnames(re.subjects) <- c("Subject", "IntSubject", "SlopeSubject")
# Observations
df <- expand.grid(unique(re.subjects$Subject), seq(replicates_cell),
                 c("A", "B"))
colnames(df) <- c("Subject", "Item", "Condition")
# Numerical coding (treatment contrasts)
df$Condition.num <- ifelse(df$Condition=="A", 0, 1)
# Fixed effects and sigma
df$Intercept <- ifixed
df$Slope <- sfixed1
df$Error <- rlst(nrow(df), phi, 0, sigma)
# Merge random effects
df <- merge(df, re.subjects)
# Response variable
df$Y <- with(df, (Intercept + IntSubject) + (Slope + SlopeSubject) *
Condition.num + Error)

m1 <- glmmTMB(Y ~ Condition.num + (1+Condition.num|Subject), df,
             family = t_family())
summary(m1)
sigma(m1) ## 1.558537
exp(m1$fit$par["psi"]) ## 3.633625

#### Simulation code end ###

-- 
Dr. Henrik Singmann
Associate Professor, Experimental Psychology
University College London (UCL), UK
http://singmann.org

	[[alternative HTML version deleted]]

Hi Henrik,

Thanks for the reproducible example and clear question. The function
family_params(m1) will give you

family_params(m1)

Student-t df
    3.633625

cheers,
Mollie

On 29 Oct 2024, at 11.48, Henrik Singmann <singmann at gmail.com> wrote:

Hi all,

It seems to me as if the estimate for the df of the t-family in glmmTMB

are

somewhat hidden from the user. In particular, the only way I have found

to

get to it is via exp(model$fit$par["psi"])

So my question is just if I am missing something obvious here.

To give some background, I am trying to fit some data that looks pretty
leptokurtic so t-distribution  seems appropriate which is luckily

available

in glmmTMB. To get going with it, I first wanted to see whether I can
actually recover the df, which seems to work reasonably well for large N,
see code below.

Best,
Henrik

#### Simulation code (somewhat long) follows ###
library("mvtnorm")
library("extraDistr")
library("glmmTMB")

# Sample size
nsubjects <- 100
replicates_cell <- 20

# Fixed effects
ifixed <- 0.5 # intercept
sfixed1 <- 0.5 # slope

# Subject random effects
is <- 0.5 # sd intercept
ss <- 0.5 # sd slope

# Correlations
rs.is <- 0.5 # intercept and slope

# Residual
sigma <- 1.5
phi <- 3.4

set.seed(56771234)
cov.is <- rs.is*is*ss
cov.subjects <- matrix(c(is^2, cov.is, cov.is, ss^2),
                      nrow=2, byrow=TRUE)
# Random sample from bivariate normal, for each subject
re.subjects <- data.frame(cbind(1:nsubjects,
                               rmvnorm(nsubjects, mean=c(rep(0, 2)),
                                       sigma=cov.subjects)))
colnames(re.subjects) <- c("Subject", "IntSubject", "SlopeSubject")
# Observations
df <- expand.grid(unique(re.subjects$Subject), seq(replicates_cell),
                 c("A", "B"))
colnames(df) <- c("Subject", "Item", "Condition")
# Numerical coding (treatment contrasts)
df$Condition.num <- ifelse(df$Condition=="A", 0, 1)
# Fixed effects and sigma
df$Intercept <- ifixed
df$Slope <- sfixed1
df$Error <- rlst(nrow(df), phi, 0, sigma)
# Merge random effects
df <- merge(df, re.subjects)
# Response variable
df$Y <- with(df, (Intercept + IntSubject) + (Slope + SlopeSubject) *
Condition.num + Error)

m1 <- glmmTMB(Y ~ Condition.num + (1+Condition.num|Subject), df,
             family = t_family())
summary(m1)
sigma(m1) ## 1.558537
exp(m1$fit$par["psi"]) ## 3.633625

#### Simulation code end ###

--
Dr. Henrik Singmann
Associate Professor, Experimental Psychology
University College London (UCL), UK
http://singmann.org

      [[alternative HTML version deleted]]

Thanks Mollie, exactly what I was looking for!

Am Di., 29. Okt. 2024 um 13:13 Uhr schrieb Mollie Brooks <
mollieebrooks at gmail.com>:

Hi Henrik,

Thanks for the reproducible example and clear question. The function
family_params(m1) will give you

family_params(m1)

Student-t df
     3.633625

cheers,
Mollie

On 29 Oct 2024, at 11.48, Henrik Singmann <singmann at gmail.com> wrote:

Hi all,

It seems to me as if the estimate for the df of the t-family in glmmTMB

are

somewhat hidden from the user. In particular, the only way I have found

to

get to it is via exp(model$fit$par["psi"])

So my question is just if I am missing something obvious here.

To give some background, I am trying to fit some data that looks pretty
leptokurtic so t-distribution  seems appropriate which is luckily

available

in glmmTMB. To get going with it, I first wanted to see whether I can
actually recover the df, which seems to work reasonably well for large N,
see code below.

Best,
Henrik

#### Simulation code (somewhat long) follows ###
library("mvtnorm")
library("extraDistr")
library("glmmTMB")

# Sample size
nsubjects <- 100
replicates_cell <- 20

# Fixed effects
ifixed <- 0.5 # intercept
sfixed1 <- 0.5 # slope

# Subject random effects
is <- 0.5 # sd intercept
ss <- 0.5 # sd slope

# Correlations
rs.is <- 0.5 # intercept and slope

# Residual
sigma <- 1.5
phi <- 3.4

set.seed(56771234)
cov.is <- rs.is*is*ss
cov.subjects <- matrix(c(is^2, cov.is, cov.is, ss^2),
                       nrow=2, byrow=TRUE)
# Random sample from bivariate normal, for each subject
re.subjects <- data.frame(cbind(1:nsubjects,
                                rmvnorm(nsubjects, mean=c(rep(0, 2)),
                                        sigma=cov.subjects)))
colnames(re.subjects) <- c("Subject", "IntSubject", "SlopeSubject")
# Observations
df <- expand.grid(unique(re.subjects$Subject), seq(replicates_cell),
                  c("A", "B"))
colnames(df) <- c("Subject", "Item", "Condition")
# Numerical coding (treatment contrasts)
df$Condition.num <- ifelse(df$Condition=="A", 0, 1)
# Fixed effects and sigma
df$Intercept <- ifixed
df$Slope <- sfixed1
df$Error <- rlst(nrow(df), phi, 0, sigma)
# Merge random effects
df <- merge(df, re.subjects)
# Response variable
df$Y <- with(df, (Intercept + IntSubject) + (Slope + SlopeSubject) *
Condition.num + Error)

m1 <- glmmTMB(Y ~ Condition.num + (1+Condition.num|Subject), df,
              family = t_family())
summary(m1)
sigma(m1) ## 1.558537
exp(m1$fit$par["psi"]) ## 3.633625

#### Simulation code end ###

--
Dr. Henrik Singmann
Associate Professor, Experimental Psychology
University College London (UCL), UK
http://singmann.org

       [[alternative HTML version deleted]]