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

[Ben Bolker]

   I agree that this ought to be more visible: I note that it's included 
in the print method [print(m1)] but should be included in the summary 
output as well.

   As far as making family_params() more visible, maybe we could link to 
it from ?glmmTMB and ?sigma.glmmTMB ...


On 10/29/24 10:18, Henrik Singmann wrote:
> Thanks Mollie, exactly what I was looking for!
> 
> Am Di., 29. Okt. 2024 um 13:13 Uhr schrieb Mollie Brooks <
> mollieebrooks using 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 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
>>
>>
> 

-- 
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
* E-mail is sent at my convenience; I don't expect replies outside of 
working hours.