## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144
)
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}

## ----personality-simulate-----------------------------------------------------
set.seed(20260715)
n_individual <- 80L
n_each <- 6L

individual_info <- data.frame(
  individual = factor(seq_len(n_individual)),
  sex = factor(
    rep(c("female", "male"), each = n_individual / 2),
    levels = c("female", "male")
  )
)

mean_by_sex <- c(female = 0.35, male = 0.70)
between_sd_by_sex <- c(female = 0.65, male = 0.40)
within_sd_by_sex <- c(female = 0.35, male = 0.60)

individual_effect <- stats::rnorm(
  n_individual,
  sd = between_sd_by_sex[individual_info$sex]
)

personality <- individual_info[rep(seq_len(n_individual), each = n_each), ]
personality$exploration_score <-
  mean_by_sex[personality$sex] +
  individual_effect[as.integer(personality$individual)] +
  stats::rnorm(
    nrow(personality),
    sd = within_sd_by_sex[personality$sex]
  )

## ----personality-fit----------------------------------------------------------
fit_personality <- drmTMB(
  bf(
    exploration_score ~ sex + (1 | individual),
    sigma ~ sex,
    sd(individual) ~ sex
  ),
  family = gaussian(),
  data = personality
)

check_drm(fit_personality)
round(coef(fit_personality, "mu"), 3)
round(coef(fit_personality, "sigma"), 3)
round(coef(fit_personality, "sd(individual)"), 3)

## ----personality-interval-targets---------------------------------------------
personality_targets <- profile_targets(fit_personality)
personality_targets[
  personality_targets$parm %in% c(
    "fixef:sigma:sexmale",
    "fixef:sd(individual):sexmale"
  ),
  c("parm", "estimate", "profile_ready", "profile_note")
]

## ----personality-component-intervals------------------------------------------
personality_ci <- confint(
  fit_personality,
  parm = c(
    "fixef:sigma:sexmale",
    "fixef:sd(individual):sexmale"
  ),
  method = "wald"
)
personality_ci

## ----personality-ratio-intervals----------------------------------------------
interval_ok <- personality_ci$conf.status == "wald" &
  stats::complete.cases(personality_ci[c("lower", "upper")])
data.frame(
  parm = personality_ci$parm[interval_ok],
  sd_ratio_lower = exp(personality_ci$lower[interval_ok]),
  sd_ratio_upper = exp(personality_ci$upper[interval_ok]),
  method = personality_ci$method[interval_ok],
  conf.status = personality_ci$conf.status[interval_ok]
)

## ----personality-data-figure, echo = FALSE, fig.width = 7, fig.height = 4.1, fig.cap = "Repeated exploration scores for females and males. Grey points are observations, blue ticks are individual means, and the vermillion line is the fitted sex-specific mean. The spread of individual means represents between-individual variation; scatter around each individual mean represents within-individual variation. Both panels use the same vertical scale.", fig.alt = "Two panels show repeated exploration scores for female and male individuals. Grey observations cluster around blue individual-mean ticks, while a vermillion horizontal line marks the fitted sex-specific mean. The female panel has wider spread among individual means, while the male panel has more scatter within individuals."----
individual_summary <- stats::aggregate(
  exploration_score ~ individual + sex,
  data = personality,
  FUN = mean
)
individual_summary <- individual_summary[
  order(individual_summary$sex, individual_summary$exploration_score),
]
individual_summary$display_id <- ave(
  individual_summary$exploration_score,
  individual_summary$sex,
  FUN = seq_along
)

plot_data <- merge(
  personality,
  individual_summary[c("individual", "sex", "display_id")],
  by = c("individual", "sex"),
  sort = FALSE
)

sex_grid <- data.frame(
  sex = factor(c("female", "male"), levels = levels(personality$sex))
)
sex_grid$fitted_mean <- predict(
  fit_personality,
  newdata = sex_grid,
  dpar = "mu"
)

ggplot2::ggplot(
  plot_data,
  ggplot2::aes(x = display_id, y = exploration_score)
) +
  ggplot2::geom_point(
    position = ggplot2::position_jitter(width = 0.08, height = 0, seed = 1),
    colour = "grey55",
    alpha = 0.42,
    size = 0.8
  ) +
  ggplot2::geom_segment(
    data = individual_summary,
    ggplot2::aes(
      x = display_id - 0.26,
      xend = display_id + 0.26,
      y = exploration_score,
      yend = exploration_score
    ),
    inherit.aes = FALSE,
    colour = "#0072B2",
    linewidth = 0.65
  ) +
  ggplot2::geom_hline(
    data = sex_grid,
    ggplot2::aes(yintercept = fitted_mean),
    inherit.aes = FALSE,
    colour = "#D55E00",
    linewidth = 0.8
  ) +
  ggplot2::geom_text(
    data = sex_grid,
    ggplot2::aes(
      x = Inf,
      y = fitted_mean,
      label = "fitted sex mean"
    ),
    inherit.aes = FALSE,
    colour = "#D55E00",
    hjust = 1.05,
    vjust = -0.6,
    size = 3
  ) +
  ggplot2::facet_wrap(
    ~sex,
    nrow = 1,
    labeller = ggplot2::as_labeller(c(female = "Female", male = "Male"))
  ) +
  ggplot2::labs(
    x = "Individuals, ordered within sex by their observed mean",
    y = "Exploration score"
  ) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    panel.grid.major.x = ggplot2::element_blank(),
    axis.text.x = ggplot2::element_blank(),
    axis.ticks.x = ggplot2::element_blank()
  )

## ----personality-components---------------------------------------------------
sex_grid$mean_score <- predict(
  fit_personality, newdata = sex_grid, dpar = "mu"
)
sex_grid$between_individual_sd <- predict(
  fit_personality, newdata = sex_grid, dpar = "sd(individual)"
)
sex_grid$within_individual_sd <- predict(
  fit_personality, newdata = sex_grid, dpar = "sigma"
)

sex_grid$repeatability <- with(
  sex_grid,
  between_individual_sd^2 /
    (between_individual_sd^2 + within_individual_sd^2)
)

repeatability_table <- sex_grid[c(
  "sex",
  "mean_score",
  "between_individual_sd",
  "within_individual_sd",
  "repeatability"
)]
repeatability_table[-1] <- lapply(
  repeatability_table[-1],
  round,
  digits = 3
)
repeatability_table

## ----personality-component-figure, echo = FALSE, fig.width = 7, fig.height = 2.9, fig.cap = "Model-implied values for females and males. Panels show the expected exploration score, between-individual SD from `sd(individual)`, and within-individual residual SD from `sigma`. Lines aid comparison and are not uncertainty intervals.", fig.alt = "Three small panels compare fitted female and male values. The expected score and within-individual standard deviation are higher for males, while the between-individual standard deviation is higher for females. No uncertainty intervals are shown."----
component_surface <- rbind(
  data.frame(
    sex = sex_grid$sex,
    component = "Expected score",
    estimate = sex_grid$mean_score
  ),
  data.frame(
    sex = sex_grid$sex,
    component = "Between-individual SD",
    estimate = sex_grid$between_individual_sd
  ),
  data.frame(
    sex = sex_grid$sex,
    component = "Within-individual SD",
    estimate = sex_grid$within_individual_sd
  )
)
component_surface$component <- factor(
  component_surface$component,
  levels = c(
    "Expected score",
    "Between-individual SD",
    "Within-individual SD"
  )
)

ggplot2::ggplot(
  component_surface,
  ggplot2::aes(x = sex, y = estimate, group = 1)
) +
  ggplot2::geom_line(colour = "grey65", linewidth = 0.55) +
  ggplot2::geom_point(colour = "#0072B2", size = 2.4) +
  ggplot2::geom_text(
    ggplot2::aes(label = sprintf("%.2f", estimate)),
    vjust = -0.8,
    colour = "grey20",
    size = 3.1
  ) +
  ggplot2::facet_wrap(~component, scales = "free_y", nrow = 1) +
  ggplot2::expand_limits(y = 0) +
  ggplot2::scale_y_continuous(
    expand = ggplot2::expansion(mult = c(0.02, 0.18))
  ) +
  ggplot2::scale_x_discrete(labels = c(female = "Female", male = "Male")) +
  ggplot2::labs(x = NULL, y = "Fitted value") +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    panel.grid.major.x = ggplot2::element_blank(),
    strip.text = ggplot2::element_text(face = "bold")
  )