1 El `lsm`

Es de autoría propia (Villalba JL, Llinás HJ y Fabregas OJ (2025)) que permite calcular de forma sencilla y precisa el log-verosímil del modelo saturado cuando la variable respuesta Y es dicotómica (0 o 1). Además, incluye conjuntos de datos de ejemplo útiles para ilustrar su aplicación. Se puede acceder a él así:

library(lsm)

2 Sobre el `lsm`

El paquete lsm, desarrollado por Villalba JL, Llinás HJ y Fabregas OJ (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 dicotómica (es decir, toma valores de 0 o 1), la función lsm() 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 lsm no solo incluye funciones especializadas, sino también conjuntos de datos integrados que permiten realizar pruebas y ejemplos de manera práctica: chdage, icu, lowbwt, pros, survey y uis.
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 `lsm`

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

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

El comando help(package = "lsm") en R abre la documentación general del paquete lsm, 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("lsm", envir = asNamespace("lsm"))) sirve para ver el código fuente de la función lsm() (o de cualquier otra función en R).

La función ls("package:lsm") devuelve un vector de nombres de objetos que puedes usar. Eso significa que el paquete lsm contiene, entre otros, una función llamada lsm y los datasets: chdage, icu, lowbwt, pros, survey y uis.

## [1] "chdage" "icu"    "lowbwt" "lsm"    "pros"   "survey" "uis"

4 Salida del modelo `lsm()`

Al ajustar un modelo logístico con la función lsm() 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: Vector con las estimaciones de los coeficientes (intercepto y pendientes).
coef: Vector con las estimaciones de los coeficientes (intercepto y pendientes).
Std.Error: Vector con los errores estándar de los coeficientes (intercepto y pendientes).
ExpB: Vector con el exponencial de los coeficientes (intercepto y pendientes).
Wald: Valor de la estadística de Wald (con distribución Chi-cuadrado).
DF: Grados de libertad para la distribución Chi-cuadrado.
P.value: P-valor calculado con la distribución Chi-cuadrado.
Log_Lik_Complete: Estimación del log-verosímil en el modelo completo.
Log_Lik_Null: Estimación del log-verosímil en el modelo nulo.
Log_Lik_Logit: Estimación del log-verosímil en el modelo logístico.
Log_Lik_Saturate: Estimación del log-verosímil en el modelo saturado.
Populations: Número de poblaciones en el modelo saturado.
Dev_Null_vs_Logit: Valor de la estadística de prueba (Hipótesis: modelo nulo vs modelo logístico).
Dev_Logit_vs_Complete: Valor de la estadística de prueba (Hipótesis: modelo logístico vs modelo completo).
Dev_Logit_vs_Saturate: Valor de la estadística de prueba (Hipótesis: modelo logístico vs modelo saturado).
Df_Null_vs_Logit: Grados de libertad para la estadística de prueba (Hipótesis: modelo nulo vs modelo logístico).
Df_Logit_vs_Complete: Grados de libertad para la estadística de prueba (Hipótesis: modelo logístico vs modelo completo).
Df_Logit_vs_Saturate: Grados de libertad para la estadística de prueba (Hipótesis: modelo logístico vs modelo saturado).
P.v_Null_vs_Logit: P-valor para la prueba de hipótesis: modelo nulo vs modelo logístico.
P.v_Logit_vs_Complete: P-valor para la prueba de hipótesis: modelo logístico vs modelo completo.
P.v_Logit_vs_Saturate: P-valor para la prueba de hipótesis: modelo logístico vs modelo saturado.
Logit: Vector con los log-odds.
p_hat_complete: Vector con las probabilidades estimadas de que la variable respuesta sea 1, dado la \(j\)-ésima población (estimadas con el modelo completo y sin el modelo logístico).
p_hat_null: Vector con las probabilidades estimadas de que la variable respuesta sea 1, dado la \(j\)-ésima población (estimadas con el modelo nulo y sin el modelo logístico).
p_j: Vector con las probabilidades estimadas de que la variable respuesta sea 1, dado la \(j\)-ésima población (estimadas con el modelo logístico).
odd: Vector con los valores de los odds en cada \(j\)-ésima población.
OR: Vector con los valores del odds ratio para cada coeficiente de las variables.
z_j: Vector con los valores de cada \(Z_j\) (la suma de observaciones en la \(j\)-ésima población).
n_j: Vector con los valores de \(n_j\) (el número de observaciones en cada \(j\)-ésima población).
p_j_tilde: Vector con la estimación de cada \(p_j\) (probabilidad de éxito en la \(j\)-ésima población) en el modelo saturado (sin estimar los parámetros logísticos).
v_j: Vector con la varianza de las variables Bernoulli en la \(j\)-ésima población.
m_j: Vector con los valores esperados de \(Z_j\) en la \(j\)-ésima población.
V_j: Vector con las varianzas de \(Z_j\) en la \(j\)-ésima población.
V: Matriz de varianza y covarianza de \(Z\), el vector que contiene todos los \(Z_j\).
S_p: Vector de puntuaciones (score vector) en el modelo saturado.
I_p: Matriz de información en el modelo saturado.
Zast_j: Vector con los valores estandarizados de \(Z_j\).
mcov: Matriz de varianza y covarianza para las estimaciones de los coeficientes.
mcor: Matriz de correlaciones para las estimaciones de los coeficientes.
Esm: Data frame con estimaciones en el modelo saturado. Contiene para cada población \(j\): el valor de las variables explicativas, \(n_j\), \(Z_j\), \(p_j\) y el log-verosímil \(L_{j}\).
Elm: Data frame con estimaciones en el modelo logístico. Para cada población \(j\), contiene: \(n_j\), \(Z_j\), \(g_j\) (el Logit de \(p_j\)), \(p_j\), el log-verosímil \(L_j\) y var.logit (la varianza del logit).

5 Funciones aplicables al modelo `lsm()`

5.0.1 Funciones

Una vez que se ha ajustado el modelo con

modelo <- lsm(y~x, data)

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

summary(modelo),
confint(modelo),
predict(modelo),
plot(modelo).

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

5.0.2 `summary(modelo)`

Muestra un resumen detallado del modelo, incluyendo coeficientes, errores estándar, estadísticos y P-valores.

5.0.3 `confint(modelo)`

Esta función calcula intervalos de confianza (IC) para los parámetros estimados del modelo. El resultado incluye:

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

Esta función tiene los siguientes argumentos:

parm: permite indicar un coeficiente específico (por ejemplo, "(Intercept)" o "AGE").
level: define el nivel de confianza (por defecto 0.95, es decir, 95%).

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

Ejemplo de sintaxis:

confint(modelo, 
        parm = "(Intercept)", # IC para el intercepto
        #parm = "AGE",        # IC para la pendiente
        level = 0.90)         # Grado de confianza

La función confint.lsm() devuelve una lista con dos componentes:

confint: IC para los coeficientes \(\beta\).
OR: IC para las odds ratio asociadas (\(\exp(\beta)\)).

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

confint(modelo)$confint # IC de β
confint(modelo)$OR      # IC de la OR

5.0.4 `confint(modelo)` para nuevos valores (`newdata`)

También es posible obtener IC para medidas derivadas del modelo (logit, probabilidad y odds) en puntos concretos de la variable explicativa. Por ejemplo,

# IC para el LOGIT en AGE = 50
confint(modelo, target = "link", newdata = data.frame(AGE = 50))

# IC para la PROBABILIDAD en AGE = 35 y 50
confint(modelo, target = "response", newdata = data.frame(AGE = c(35, 50)))

# IC para los ODDS en varias edades
ages <- data.frame(AGE = seq(30, 65, by = 5))
confint(modelo, target = "odd", newdata = ages, level = 0.90)  # IC al 90%

Se resalta que, con predict(modelo), también se obtienen estos mismos intervalos.

5.0.5 `predict(modelo)`

Calcula valores predichos (por ejemplo, probabilidades), en diferentes escalas, para nuevas observaciones o para las observaciones originales y sus correspondientes intervalos de confianzas. También, se puede cambiar el grado de confianza con level (por defecto, 0.95).

predict(modelo, type = "link")     # logit
predict(modelo, type = "response") # Probabilidades estimadas pj
predict(modelo, type = "odd")       # odds
predict(modelo, type = "OR")        # odds ratio para cada covariable (no depende de `newdata`).

5.0.6 `predict(modelo)` para nuevos valores (`newdata`)

Para un un solo valor de una variable, digamos X:

predict(modelo, newdata = data.frame(X = 35), type = "link")  
predict(modelo, newdata = data.frame(X = 35), type = "response")
predict(modelo, newdata = data.frame(X = 35), type = "odd")

Para varios valores, digamos de X:

x <- data.frame(X = c(30, 40, 50, 60))
predict(modelo, newdata = x, type = "response")

Comentarios:

Cuando se usa newdata, la salida corresponde a esas filas nuevas.
"OR" ignora newdata (es propiedad del coeficiente, no del punto de predicción).

5.0.7 `plot(modelo)`

Genera gráficos relacionados con el ajuste del modelo logístico:

plot(modelo, type = "scatter")     # Dispersión y variable respuesta
plot(modelo, type = "probability") # Probabilidades estimadas pj
plot(modelo, type = "Logit")       # gj (logit de pj)
plot(modelo, type = "odds")        # Odds = pj /(1-pj)

5.0.8 Personalización de gráficos

La función plot() permite incluir argumentos adicionales para personalizar la visualización, como título, colores, tamaño de los puntos o etiquetas de ejes. Por ejemplo:

plot(
  modelo,
  type  = "probability",
  color = "blue",
  size  = 2.5,
  title = "Title",
  xlab  = "Variable X",
  ylab  = "Variable Y"
     ) + 

facet_wrap(. ~ "Otro título") +    
theme_bw(base_size = 12)

En este ejemplo:

color cambia el color de los puntos.
size controla el tamaño de los puntos.
title define el título del gráfico.
xlab y ylab modifican las etiquetas de los ejes.
Después del plot(), se pueden añadir capas extra de ggplot2 (las que se muestran en el código, es solo un ejemplo).

Esta flexibilidad permite adaptar las gráficas generadas por plot() a necesidades específicas de análisis, publicación o presentación, aprovechando todo el potencial de ggplot2.

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

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

# Ajuste del modelo (ejemplo simple)
modelo <- lsm(CHD~AGE, data=datos)

modelo 
summary(modelo)

# --- Estimación y pruebas ---
modelo$coefficients      # Estimaciones de coeficientes (intercepto y pendientes)
modelo$coef              # Alias de coefficients
modelo$Std.Error         # Errores estándar
modelo$ExpB              # exp(beta) = Odds Ratios
modelo$Wald              # Estadísticos de Wald (beta/SE)^2
modelo$DF                # Grados de libertad (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

# --- Número de poblaciones (J) ---
modelo$Populations

# --- Deviances (comparaciones entre modelos) ---
modelo$Dev_Null_vs_Logit     # 2 (L_logit - L_null)
modelo$Dev_Logit_vs_Complete # 2 (L_complete - 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

# --- 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 poblaciones ---
modelo$Logit            # Vector de log-odds g_j = log(p_j / (1 - p_j))
modelo$p_hat_complete   # Probabilidades estimadas en el modelo completo
modelo$p_hat_null       # Probabilidades estimadas en el modelo nulo
modelo$p_j              # Probabilidades estimadas en el modelo logístico
modelo$odd              # Odds p_j / (1 - p_j)
modelo$OR               # exp(beta) (odds ratios por coeficiente)
modelo$z_j              # Conteos Z_j por población j
modelo$n_j              # Tamaños n_j por población j
modelo$p_j_tilde        # Proporciones empíricas (saturado): Z_j / n_j
modelo$v_j              # Varianzas Bernoulli v_j = p_j (1 - p_j)
modelo$m_j              # Esperanzas m_j = n_j * p_j
modelo$V_j              # Var(Z_j) = n_j * v_j
modelo$V                # Matriz Var-Cov de Z (todas las poblaciones)
modelo$S_p              # Score bajo saturado
modelo$I_p              # Matriz de información bajo saturado
modelo$Zast_j           # Valores tipificados de Z_j

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

# --- Tablas estimadas (saturado y logístico) ---
modelo$Esm              # Estimated Saturated Matrix (covariables, n_j, Z_j, p~_j, L_j)
modelo$Elm              # Estimated Logit Matrix (n_j, Z_j, g_j, p_j, L_j, var_logit_j)

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

# Inspección rápida (usa el método print.confint.lsm)
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)); según tu implementación, puede omitir el intercepto
ci_all$level      # nivel de confianza en porcentaje (p. ej., 95)

# IC (90%) para un coeficiente específico
names(modelo$coefficients)      # ver nombres disponibles
coef_name <- names(modelo$coefficients)[1]
ci_one <- confint(modelo, parm = coef_name, level = 0.90)
ci_one

## Tabla "ordenada" con beta, SE, OR y sus IC {.unnumbered .unlisted}
# Construye una tabla unificada como en glsm

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; 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")
}

# Tabla base con estimaciones puntuales
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_
}

# Mantener el orden original de coef_names
out_ci <- out_ci[match(coef_names, out_ci$term), ]

# Mostrar tabla cruda
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()

# Tabla solo para OR (excluyendo intercepto explícitamente)
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 rápidamente
ci90 <- confint(modelo, level = 0.90)
ci90$confint
ci90$OR
ci90$level

# Aplicar confint() a un subconjunto de coeficientes y recolectar resultados
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

#REVISAR: La función confint.lsm() no acepta target ni newdata. Solo para coeficientes. Debe mejorarse.

## IC para cantidades derivadas vía `confint.lsm()` con newdata {.unnumbered .unlisted}
# También puedes pedir IC para logit, probabilidad y odds en puntos específicos

# IC para el LOGIT en age = 50
confint(modelo, target = "link", newdata = data.frame(age = 50))

# IC para la PROBABILIDAD en age = 35 y 50
confint(modelo, target = "response", newdata = data.frame(age = c(35, 50)))

# IC para los ODDS en varias edades (al 90%)
ages <- data.frame(age = seq(30, 65, by = 5))
confint(modelo, target = "odd", newdata = ages, level = 0.90)

## Predicción y bandas de confianza con `predict.lsm()` {.unnumbered .unlisted}
# Predicciones en distintas escalas (sobre los datos originales)
predict(modelo, type = "link")      # logit g_j
predict(modelo, type = "response")  # probabilidad p_j
predict(modelo, type = "odd")       # odds p_j / (1 - p_j)
predict(modelo, type = "OR")        # OR por coeficiente (no depende de newdata)

# Predicciones para nuevos valores
predict(modelo, newdata = data.frame(AGE = 35), type = "link")
predict(modelo, newdata = data.frame(AGE = 35), type = "response")
predict(modelo, newdata = data.frame(AGE = 35), type = "odd")

# REVISAR. NO FUNCIONA. DEBE MEJORARSE

# Varios puntos de predicción
nd <- data.frame(AGE = c(30, 40, 50, 60))
predict(modelo, newdata = nd, type = "response")

# Cambiar nivel de confianza en predict (si la implementación lo soporta)
predict(modelo, newdata = nd, type = "response", level = 0.90)

plot(modelo, type = "scatter")     # Dispersión y variable respuesta
plot(modelo, type = "probability") # Probabilidades estimadas pj
plot(modelo, type = "Logit")       # gj (logit de pj)
plot(modelo, type = "odds")        # Odds = pj /(1-pj)

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 BINARIO

Paquete `lsm` en R

Dr. rer. nat. Humberto LLinás Solano

1 El `lsm`

2 Sobre el `lsm`

3 Exploración del `lsm`

4 Salida del modelo `lsm()`

5 Funciones aplicables al modelo `lsm()`

5.0.1 Funciones

5.0.2 `summary(modelo)`

5.0.3 `confint(modelo)`

5.0.4 `confint(modelo)` para nuevos valores (`newdata`)

5.0.5 `predict(modelo)`

5.0.6 `predict(modelo)` para nuevos valores (`newdata`)

5.0.7 `plot(modelo)`

5.0.8 Personalización de gráficos

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

Bibliografía

Post en las redes

EL CASO BINARIO

Paquete lsm en R

Dr. rer. nat. Humberto LLinás Solano

1 El lsm

2 Sobre el lsm

3 Exploración del lsm

4 Salida del modelo lsm()

5 Funciones aplicables al modelo lsm()

5.0.1 Funciones

5.0.2 summary(modelo)

5.0.3 confint(modelo)

5.0.4 confint(modelo) para nuevos valores (newdata)

5.0.5 predict(modelo)

5.0.6 predict(modelo) para nuevos valores (newdata)

5.0.7 plot(modelo)

5.0.8 Personalización de gráficos

6 Acceso directo a los componentes de lsm()

Bibliografía

Post en las redes

Paquete `lsm` en R

1 El `lsm`

2 Sobre el `lsm`

3 Exploración del `lsm`

4 Salida del modelo `lsm()`

5 Funciones aplicables al modelo `lsm()`

5.0.2 `summary(modelo)`

5.0.3 `confint(modelo)`

5.0.4 `confint(modelo)` para nuevos valores (`newdata`)

5.0.5 `predict(modelo)`

5.0.6 `predict(modelo)` para nuevos valores (`newdata`)

5.0.7 `plot(modelo)`

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