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

[Henrik Singmann]

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. Henrik Singmann
Associate Professor, Experimental Psychology
University College London (UCL), UK
http://singmann.org

	[[alternative HTML version deleted]]