1 El `glsm`

Es de autoría propia (Llinás H, Villalba J, Borja J, Tilano J (2025)) que permite calcular de forma sencilla y precisa el log-verosímil del modelo saturado cuando la variable respuesta Y es multinomial (tiene más de dos niveles). Además, incluye conjuntos de datos de ejemplo útiles para ilustrar su aplicación. Se puede acceder a él así:

library(glsm)

2 Sobre el `glsm`

El paquete glsm, desarrollado por Llinás H, Villalba J, Borja J, Tilano J (2025), es una contribución propia que puede descargarse desde el repositorio oficial de CRAN. Su propósito principal es facilitar la estimación del log-verosimilitud del modelo saturado, herramienta clave en la evaluación comparativa de modelos logísticos.
Cuando la variable de respuesta $Y$ es multinomial (tiene más de dos niveles), la función glsm() calcula automáticamente el valor del log-verosímil del modelo saturado. Este valor permite comparar otros modelos y evaluar su bondad de ajuste relativa.
El paquete glsm no solo incluye funciones especializadas, sino también un conjunto de dato integrado que permiten realizar pruebas y ejemplos de manera práctica: hsbdemo.
Este paquete resulta especialmente útil en cursos de regresión logística, modelado estadístico y evaluación de modelos con base en la verosimilitud.

3 Exploración del `glsm`

Para explorar el contenido del paquete directamente desde R, se pueden usar los siguientes comandos:

help(package = "glsm")    # Muestra la documentación general del paquete
print(get("glsm", envir = asNamespace("glsm"))) # ver código fuente
ls("package:glsm")        # Lista todas las funciones y datasets disponibles

El comando help(package = "glsm") en R abre la documentación general del paquete glsm, mostrando:

Una descripción general del paquete.
Un listado de todas las funciones disponibles en el paquete.
Un listado de los datasets incluidos (si los hay).
Información sobre autores, dependencias y versiones.

La función print(get("glsm", envir = asNamespace("glsm"))) sirve para ver el código fuente de la función glsm() (o de cualquier otra función en R).

La función ls("package:glsm") devuelve un vector de nombres de objetos que puedes usar. Eso significa que el paquete glsm contiene, entre otros, una función llamada glsm y el dataset: hsbdemo.

## [1] "glsm"    "hsbdemo"

4 Salida del modelo `glsm()`

4.0.1 Descripción

La función glsm() ajusta un modelo logístico multinomial (con $R \geq 3$ categorías en la respuesta) y, además, calcula los log-verosímiles para cuatro configuraciones: nulo, completo, saturado y logístico. El objeto retornado (clase "glsm") contiene componentes para: estimación de parámetros, comparación entre modelos, matrices de información/varianzas, y resúmenes a nivel de poblaciones $j=1,\dots,J$.

Convención de índices:

$r$ denota categoría de la respuesta (respecto a la categoría de referencia ref).
$j$ denota la población (combinación única de covariables) en el modelo saturado.

Definiciones útiles:

Probabilidades predichas por el modelo logístico, para la población $j$:
\[ p_{rj}= \Pr(Y=r \mid \text{población }j). \]
Si $\text{ref}$ es la categoría de referencia, los odds y el logit para $r\neq \text{ref}$ son:

\[ \text{odds}_{rj} \;=\; \frac{p_{rj}}{p_{\text{ref},\,j}}, \qquad \text{logit}_{rj} \;=\; \log\!\left(\frac{p_{rj}}{p_{\text{ref},\,j}}\right). \]

Principales componentes del objeto `glsm`

Al ajustar un modelo logístico con la función glsm() en R, el objeto resultante contiene múltiples elementos que permiten realizar un análisis detallado del ajuste, la estimación de parámetros y las comparaciones entre modelos. A continuación, se describe cada uno de los principales componentes que se obtienen:

coefficients: Estimaciones de los coeficientes (interceptos y pendientes) del modelo logístico multinomial.
coef: Alias de coefficients (mismo contenido).
Std.Error: Errores estándar asociados a cada coeficiente.
ExpB: Exponenciación de los coeficientes $\exp(\beta)$, útil como odds ratios.
Wald: Estadísticos de Wald $(\beta/\text{SE})^2$ (asintóticamente $\chi^2$ con 1 grado de libertad).
DF: Grados de libertad asociados a cada prueba de Wald (1 por coeficiente).
P.value: $p$-valores de las pruebas de Wald.
Log_Lik_Complete: Log-verosímil del modelo completo (${\cal L}_{\text{complete}}$).
Log_Lik_Null: Log-verosímil del modelo nulo (${\cal L}_{\text{null}}$).
Log_Lik_Logit: Log-verosímil del modelo logístico ajustado (${\cal L}_{\text{logit}}$).
Log_Lik_Saturate: Log-verosímil del modelo saturado (${\cal L}_{\text{saturate}}$).
Populations: Número $J$ de poblaciones (combinaciones únicas de covariables) en el modelo saturado.
Dev_Null_vs_Logit: Deviance $2({\cal L}_{\text{logit}}-{\cal L}_{\text{null}})$.
Dev_Logit_vs_Complete: Deviance $-2\,{\cal L}_{\text{logit}}$.
Dev_Logit_vs_Saturate: Deviance $2({\cal L}_{\text{saturate}}-{\cal L}_{\text{logit}})$.
Df_Null_vs_Logit: Grados de libertad para la comparación nulo vs logístico.
Df_Logit_vs_Complete: Grados de libertad para la comparación logístico vs completo.
Df_Logit_vs_Saturate: Grados de libertad para la comparación logístico vs saturado.
P.v_Null_vs_Logit: $p$-valor para la comparación nulo vs logístico.
P.v_Logit_vs_Complete: $p$-valor para la comparación logístico vs completo.
P.v_Logit_vs_Saturate: $p$-valor para la comparación logístico vs saturado.
Logit_r: Matriz de log-odds $\text{logit}_{rj}$ respecto a la categoría de referencia.
p_hat_complete: Matriz de indicadores $u_{r}$ usada en el cálculo del modelo completo.
p_hat_null: Vector con proporciones globales por categoría (modelo nulo).
p_rj: Matriz $[J\times R]$ con probabilidades $\hat p_{rj}$ del modelo logístico.
odd: Matriz de odds $\text{odds}_{rj} = p_{rj}/p_{\text{ref},j}$.
OR: Odds ratios $\exp(\beta)$ de los coeficientes.
z_rj: Matriz de conteos $Z_{rj}$.
nj: Vector con tamaños $n_j$ por población.
p_rj_tilde: Matriz de proporciones empíricas $\tilde p_{rj}=Z_{rj}/n_j$ (modelo saturado).
v_rj: Matriz con varianzas Bernoulli $v_{rj}=\hat p_{rj}(1-\hat p_{rj})$.
m_rj: Matriz con esperanzas $m_{rj}=n_j\,\hat p_{rj}$.
V_rj: Matriz con varianzas $V_{rj}=n_j\,v_{rj}$.
V: Matriz de varianza–covarianza de $Z$.
S_p: Vector/matriz de scores bajo el saturado:
\[ S_{rj} \;\propto\; \frac{n_j(\tilde p_{rj}-\hat p_{rj})}{v_{rj}}. \]
I_p: Matriz de información de Fisher aproximada.
Zast_j: Valores estandarizados de $Z_{rj}$.
mcov: Matriz de varianza–covarianza de los coeficientes.
mcor: Matriz de correlaciones entre coeficientes.
Esm: Estimated Saturated Matrix: incluye covariables, $n_j$, $Z_{rj}$, $\tilde p_{rj}$ y log-verosímil ${\cal L}_{p,j}$.
Elm: Estimated Logit Matrix: incluye covariables, $n_j$, $Z_{rj}$, $p_{rj}$, $\text{logit}_{rj}$ y var_logit_rj.
call: Llamado original de la función.

Nota sobre la referencia `ref`

Si ref no se especifica, glsm() usa por defecto el primer nivel de la respuesta.
Es recomendable declarar explícitamente ref = "nombre_categoria" para interpretar de forma clara los odds ratios y log-odds.

5 Funciones aplicables al modelo `glsm()`

5.0.1 El modelo

El modelo se ajusta con:

library(glsm)

data("hsbdemo", package = "glsm")
modelo <- glsm(prog ~ ses + gender, data = hsbdemo, ref = "academic")

modelo      # imprime resumen breve (usa print.glsm)
str(modelo) # inspecciona todos los componentes del objeto

5.0.2 Salida de `print.glsm()`

En ella se muetra:

El Call.
El número de Populations en el saturado.
Una tabla con Coef(B), Std.Error, Exp(B).
Un bloque con los log-verosímiles (Complete, Null, Logit, Saturate).
Si hay singularidades, indica cuántos coeficientes quedaron en NA.
Si hubo datos faltantes, indica cuántas observaciones fueron omitidas.

5.0.3 Funciones

El objeto resultante (modelo) puede ser utilizado con funciones adicionales para análisis y visualización:

summary(modelo),
confint(modelo).

Estas funciones aprovechan la estructura de ese objeto para facilitar la interpretación y comunicación de los resultados.

summary(modelo)
confint(modelo)

5.0.4 `summary(modelo)`

Muestra un resumen detallado del modelo, incluyendo coeficientes, errores estándar, estadísticos, P-valores y las correspondientes pruebas de comparación del modelo logístico con el nulo, completo y saturado.

5.0.5 `confint(modelo)`

`confint()` para objetos `glsm`

confint(modelo)

Esta función calcula intervalos de confianza (IC) para los parámetros estimados del modelo. Devuelve una lista con dos componentes:

confint: IC en la escala del parámetro (coeficientes $\delta$ y todos los $\beta$).
OR: IC en la escala de odds ratio para las variables explicativas, obtenido como $\exp(\beta)$.

Argumentos

Esta función tiene los siguientes argumentos:

object : objeto de clase "glsm" (resultado de glsm()).
parm : (opcional) nombre del coeficiente para el cual se desea el IC (por ejemplo, "(Intercept)" o "ses:middle:academic" si ese es un nombre real en nuestro objeto).
- Debe ser cadena de texto y pertenecer a names(object$coef).
level: define el nivel de confianza (por defecto 0.95, es decir, 95%).
... : no se admite. La función genera error si se pasan argumentos adicionales.

Si no se indica parm, se devuelven los IC para todos los coeficientes.

Validaciones internas
- Si level no es un escalar en $(0,1)$, se produce error.
- Si parm se pasa y no es carácter, o no está en names(object$coef), se produce error.
- Si el error estándar de parm es NA, se produce error.

IC para un coeficiente específico (usando un nombre real del modelo)

Primero, se recomienda chequear nombres disponibles en el modelo:

names(modelo$coefficients)     # chequear
coef_name <- names(modelo$coefficients)[1]  # ejemplo

Luego, calcular el IC deseado.

confint(modelo, 
        parm = "(Intercept):general", # IC para el intercepto
        level = 0.90)                 # Grado de confianza

confint(modelo, 
        parm = "gendermale:general",  # IC para la pendiente
        level = 0.90)                 # Grado de confianza

Extraer directamente matrices del objeto retornado

Esto permite extraer directamente el tipo de intervalo deseado. Es decir,

confint(modelo)$confint # IC de coeficientes (escala beta)
confint(modelo)$OR      # IC de OR (exp(beta)); sin intercepto en el caso "todos"

6 Acceso directo a los componentes de `glsm()`

# Cargar paquete y ajustar un ejemplo reproducible
library(glsm)
data("hsbdemo", package = "glsm")
datos <- glsm::hsbdemo

# Ajuste (cambia 'ref' si deseas otra categoría de referencia)
modelo <- glsm(prog ~ ses + gender, data = datos, ref = "academic")

modelo 
summary(modelo)

# --- Estimación y pruebas ---
modelo$coefficients      # Estimaciones de coeficientes (interceptos y pendientes)
modelo$coef              # Alias de coefficients (mismo contenido)
modelo$Std.Error         # Errores estándar por coeficiente
modelo$ExpB              # exp(beta) = Odds Ratios
modelo$Wald              # Estadísticos de Wald (beta/SE)^2
modelo$DF                # Gl asignados a cada prueba de Wald (1 por coeficiente)
modelo$P.value           # p-valores de las pruebas de Wald

# --- Log-verosímiles de los cuatro modelos ---
modelo$Log_Lik_Complete  # Log-verosímil modelo completo
modelo$Log_Lik_Null      # Log-verosímil modelo nulo
modelo$Log_Lik_Logit     # Log-verosímil modelo logístico
modelo$Log_Lik_Saturate  # Log-verosímil modelo saturado

# --- Log-verosímiles de los cuatro modelos ---
modelo$Log_Lik_Complete  # Log-verosímil modelo completo
modelo$Log_Lik_Null      # Log-verosímil modelo nulo
modelo$Log_Lik_Logit     # Log-verosímil modelo logístico
modelo$Log_Lik_Saturate  # Log-verosímil modelo saturado

# --- Poblaciones (J) en el modelo saturado ---
modelo$Populations

# --- Deviances (comparaciones entre modelos) ---
modelo$Dev_Null_vs_Logit     # 2 (L_logit - L_null)
modelo$Dev_Logit_vs_Complete # -2 * L_logit
modelo$Dev_Logit_vs_Saturate # 2 (L_saturate - L_logit)

# --- Grados de libertad de los contrastes ---
modelo$Df_Null_vs_Logit
modelo$Df_Logit_vs_Complete
modelo$Df_Logit_vs_Saturate

# --- Grados de libertad de los contrastes ---
modelo$Df_Null_vs_Logit
modelo$Df_Logit_vs_Complete
modelo$Df_Logit_vs_Saturate

# --- p-valores de los contrastes ---
modelo$P.v_Null_vs_Logit
modelo$P.v_Logit_vs_Complete
modelo$P.v_Logit_vs_Saturate

# --- Salidas a nivel de categorías/poblaciones ---
modelo$Logit_r        # Matriz de log-odds (respecto a la categoría de referencia)
modelo$p_hat_complete # Matriz de indicadores u_r (usada en el cálculo "completo")
modelo$p_hat_null     # Proporciones globales por categoría (modelo nulo)
modelo$p_rj           # Matriz JxR de probabilidades estimadas (modelo logístico)
modelo$odd            # Matriz de odds p_rj / p_ref,j
modelo$OR             # exp(beta) por coeficiente (odds ratios)
modelo$z_rj           # Matriz de conteos Z_rj por población j y categoría r
modelo$nj             # Vector de tamaños n_j por población j
modelo$p_rj_tilde     # Proporciones empíricas (saturado): Z_rj / n_j
modelo$v_rj           # Varianzas Bernoulli v_rj = p_rj (1 - p_rj)
modelo$m_rj           # Esperanzas m_rj = n_j * p_rj
modelo$V_rj           # Var(Z_rj) = n_j * v_rj
modelo$V              # Matriz Var-Cov de Z (todas las categorías/poblaciones)
modelo$S_p            # Score bajo saturado: n_j (p~_rj - p_rj) / v_rj
modelo$I_p            # Información de Fisher aproximada (a partir de S_p)
modelo$Zast_j         # Z estandarizado por categoría (valores tipificados)

# --- Matrices de varianza-covarianza y correlación de coeficientes ---
modelo$mcov           # Var-Cov de coeficientes beta
modelo$mcor           # Matriz de correlaciones entre coeficientes

# --- Tablas estimadas (saturado y logístico) ---
modelo$Esm            # Estimated Saturated Matrix (covariables, n_j, Z_rj, p~_rj, Lp_tilde)
modelo$Elm            # Estimated Logit Matrix (covariables, n_j, Z_rj, p_rj, logit_rj, var_logit_rj)

# --- Llamado original ---
modelo$call

## Intervalos de confianza con `confint.glsm()` {.unnumbered .unlisted}

# IC (95%) para TODOS los coeficientes (escala beta) y sus OR (exp(beta))
ci_all <- confint(modelo)   # usa level = 0.95 por defecto

# Inspección rápida (usa el método print.confint.glsm)
ci_all

# Extraer matrices directamente
ci_all$confint    # IC en escala de coeficiente (beta): columnas "lower %", "upper %"
ci_all$OR         # IC en escala de OR (exp(beta)); en tu implementación, excluye el intercepto
ci_all$level      # nivel de confianza en porcentaje (ej. 95)

# IC (90%) para un coeficiente específico (reemplaza por un nombre válido de tu modelo)
names(modelo$coefficients)      # para ver los nombres disponibles
coef_name <- names(modelo$coefficients)[1]
ci_one <- confint(modelo, parm = coef_name, level = 0.90)
ci_one

# Ejemplo robusto: calcular IC para TODOS los coeficientes, fila a fila
# Devuelve una tabla "ordenada" con beta, IC_beta, OR y IC_OR (OR para intercepto se muestra como "-")
coef_names <- names(modelo$coefficients)

# IC beta (todos)
ci_beta <- as.data.frame(ci_all$confint)
ci_beta$term <- rownames(ci_all$confint)
names(ci_beta)[1:2] <- c("beta_low","beta_high")

# IC OR (puede no incluir intercepto, así que alineamos por nombre)
ci_or <- as.data.frame(ci_all$OR)
if (nrow(ci_or) > 0) {
  ci_or$term <- rownames(ci_all$OR)
  names(ci_or)[1:2] <- c("OR_low","OR_high")
}

# Construir tabla final
out_ci <- data.frame(
  term      = coef_names,
  beta_hat  = as.numeric(modelo$coefficients),
  se_beta   = as.numeric(modelo$Std.Error),
  OR_hat    = as.numeric(modelo$ExpB)
)

# Unir con IC beta
out_ci <- merge(out_ci, ci_beta[, c("term","beta_low","beta_high")],
                by = "term", all.x = TRUE, sort = FALSE)

# Unir con IC OR (si existe)
if (exists("ci_or") && nrow(ci_or) > 0) {
  out_ci <- merge(out_ci, ci_or[, c("term","OR_low","OR_high")],
                  by = "term", all.x = TRUE, sort = FALSE)
} else {
  out_ci$OR_low <- NA_real_
  out_ci$OR_high <- NA_real_
}

# Para el intercepto, OR no aplica en tu print.confint.glsm() cuando se pide 'parm';
# aquí lo dejamos numérico si existe (ExpB), y los IC de OR quedan NA si no vinieron en ci_all$OR.
out_ci <- out_ci[match(coef_names, out_ci$term), ]

# Mostrar la tabla
out_ci

# (Opcional) Presentación bonita con kable + kableExtra
# install.packages("kableExtra")  # si no lo tienes
library(kableExtra)

kable(out_ci,
      digits = 4,
      align = "lrrrrrrrr",
      col.names = c("Coeficiente",
                    "Beta", "SE(Beta)", "OR",
                    "IC Beta (Inf)", "IC Beta (Sup)",
                    "IC OR (Inf)",  "IC OR (Sup)")) |>
  kable_styling(full_width = FALSE) |>
  kable_classic_2()


# Ejemplo: construir una tabla solo para OR y sus IC, excluyendo el intercepto de forma explícita
is_intercept <- out_ci$term %in% c("(Intercept)")
or_tab <- subset(out_ci, !is_intercept, 
                 select = c(term, OR_hat, OR_low, OR_high))

kable(or_tab,
      digits = 4,
      align = "lrrr",
      col.names = c("Coeficiente", "OR", "IC OR (Inf)", "IC OR (Sup)")) |>
  kable_styling(full_width = FALSE) |>
  kable_classic_2()

# Cambiar el nivel de confianza (ej. 90%) y rehacer tablas rápidamente
ci90 <- confint(modelo, level = 0.90)

ci90$confint
ci90$OR
ci90$level

# Aplicar confint() a un conjunto de coeficientes y recolectar resultados (ejemplo)
coefs_to_check <- head(names(modelo$coefficients), 5)

res_list <- lapply(coefs_to_check, function(p) {
  ci_p <- confint(modelo, parm = p, level = 0.95)
  data.frame(
    term = p,
    beta_low = ci_p$confint[1,1],
    beta_high = ci_p$confint[1,2],
    OR_low = if (is.data.frame(ci_p$OR)) ci_p$OR[1,1] else NA_real_,
    OR_high = if (is.data.frame(ci_p$OR)) ci_p$OR[1,2] else NA_real_
  )
})

res_ci_subset <- do.call(rbind, res_list)
res_ci_subset

Bibliografía

Consultar el documento RPubs :: Regresión logística (bibliografía).

Post en las redes

2025, septiembre 2: Post de Storymodelers en Linkedin.

If you found any ERRORS or have SUGGESTIONS, please report them to my email. Thanks.

EL CASO MULTINOMIAL

Paquete `glsm` en R

Dr. rer. nat. Humberto LLinás Solano

1 El `glsm`

2 Sobre el `glsm`

3 Exploración del `glsm`

4 Salida del modelo `glsm()`

4.0.1 Descripción

Principales componentes del objeto `glsm`

Nota sobre la referencia `ref`

5 Funciones aplicables al modelo `glsm()`

5.0.1 El modelo

5.0.2 Salida de `print.glsm()`

5.0.3 Funciones

5.0.4 `summary(modelo)`

5.0.5 `confint(modelo)`

`confint()` para objetos `glsm`

Argumentos

IC para un coeficiente específico (usando un nombre real del modelo)

Extraer directamente matrices del objeto retornado

6 Acceso directo a los componentes de `glsm()`

Bibliografía

Post en las redes

EL CASO MULTINOMIAL

Paquete glsm en R

Dr. rer. nat. Humberto LLinás Solano

1 El glsm

2 Sobre el glsm

3 Exploración del glsm

4 Salida del modelo glsm()

4.0.1 Descripción

Principales componentes del objeto glsm

Nota sobre la referencia ref

5 Funciones aplicables al modelo glsm()

5.0.1 El modelo

5.0.2 Salida de print.glsm()

5.0.3 Funciones

5.0.4 summary(modelo)

5.0.5 confint(modelo)

confint() para objetos glsm

Argumentos

IC para un coeficiente específico (usando un nombre real del modelo)

Extraer directamente matrices del objeto retornado

6 Acceso directo a los componentes de glsm()

Bibliografía

Post en las redes

Paquete `glsm` en R

1 El `glsm`

2 Sobre el `glsm`

3 Exploración del `glsm`

4 Salida del modelo `glsm()`

Principales componentes del objeto `glsm`

Nota sobre la referencia `ref`

5 Funciones aplicables al modelo `glsm()`

5.0.2 Salida de `print.glsm()`

5.0.4 `summary(modelo)`

5.0.5 `confint(modelo)`

`confint()` para objetos `glsm`

6 Acceso directo a los componentes de `glsm()`