[R-sig-ME] glmmTMB output for t-family

[Mollie Brooks]

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 using 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]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models