Análisis longitudinal de crecimiento de semillas con modelo mixto¶

Curso: Herramientas Estadísticas y Computacionales para el Análisis de Datos en Investigación Bioquímica
Sesión 3: Caso de estudio
Docente: Andrés García Medina
Contacto: andgarm.n@gmail.com
Sitio: https://sites.google.com/view/andresgm/home

Propósito del caso¶

Estudiaremos un experimento sintético de crecimiento radicular medido durante 15 días bajo 10 tratamientos. La pregunta central es:

¿Cambian los tratamientos la trayectoria o tasa de crecimiento de las semillas a lo largo del tiempo?

Se deben reconocer dos características del diseño:

  1. existen mediciones repetidas a través del tiempo;
  2. varias observaciones pertenecen a una misma unidad identificada como semilla.

Estas características motivan el paso desde la regresión lineal ordinaria hacia un modelo lineal mixto.

Caracteristica de los datos¶

  • Este notebook carga el archivo datos_crecimiento_semillas.xlsx, donde cada hoja representa una semilla distinta.
  • Cada hoja contiene 15 días de medición y 10 tratamientos.
  • El análisis transforma los datos de formato ancho a formato largo y ajusta un modelo mixto de efectos lineales.

Descripción del modelo¶

El modelo principal es:

$$ \text{crecimiento}_{ijt} = \beta_0 + \beta_1 \text{día}_t + \beta_2 \text{tratamiento}_i + \beta_3(\text{día}_t \times \text{tratamiento}_i) + u_j + \varepsilon_{ijt}. $$

Los subíndices identifican:

  • $i$: tratamiento o grupo experimental.
  • $j$: semilla individual.
  • $t$: momento o día de medición.
  • $\text{crecimiento}_{ijt}$: crecimiento observado de la semilla $j$, bajo el tratamiento $i$, en el día $t$.

Interpretación de los parámetros¶

  • $\beta_0$ — Intercepto: representa el crecimiento promedio esperado del grupo de referencia cuando $\text{día}=0$.

  • $\beta_1$ — Efecto del tiempo: representa el cambio promedio en el crecimiento por cada día adicional en el grupo de referencia. Corresponde a la pendiente de crecimiento del grupo control.

  • $\beta_2$ — Efecto del tratamiento: representa la diferencia promedio entre el tratamiento y el grupo de referencia cuando $\text{día}=0$.

  • $\beta_3$ — Interacción día × tratamiento: representa cuánto cambia la pendiente de crecimiento debido al tratamiento. Permite evaluar si las semillas sometidas al tratamiento crecen a una velocidad diferente de las semillas del grupo de referencia.

  • $u_j$ — Efecto aleatorio de la semilla: representa la desviación particular de la semilla $j$ respecto al comportamiento promedio. Permite considerar que cada semilla puede tener características individuales y que las mediciones repetidas sobre una misma semilla no son independientes.

  • $\varepsilon_{ijt}$ — Error residual: representa la variabilidad del crecimiento que no es explicada por el tiempo, el tratamiento, su interacción ni las diferencias individuales entre semillas.

Usualmente se supone que:

$$ u_j \sim N(0,\sigma_u^2), \qquad \varepsilon_{ijt} \sim N(0,\sigma^2), $$

donde:

  • $\sigma_u^2$ representa la variabilidad entre semillas.
  • $\sigma^2$ representa la variabilidad residual de las observaciones.

Interpretación de las pendientes¶

Uno de los principales objetivos del modelo es determinar si el tratamiento modifica la velocidad de crecimiento de las semillas a través del tiempo.

Partimos del modelo:

$$ \text{crecimiento}_{ijt} = \beta_0 + \beta_1 \text{día}_t + \beta_2 \text{tratamiento}_i + \beta_3(\text{día}_t \times \text{tratamiento}_i) + u_j + \varepsilon_{ijt}. $$

Para interpretar las pendientes, suponemos que la variable tratamiento fue codificada como:

$$ \text{tratamiento}= \begin{cases} 0, & \text{grupo control},\\ 1, & \text{grupo con tratamiento}. \end{cases} $$

Grupo control¶

Para el grupo control, $\text{tratamiento}=0$. Sustituyendo en el modelo:

$$ \text{crecimiento}_{jt} = \beta_0 + \beta_1\text{día}_t + \beta_2(0) + \beta_3(\text{día}_t\times0) + u_j + \varepsilon_{jt}. $$

Por tanto:

$$ \text{crecimiento}_{jt} = \beta_0+\beta_1\text{día}_t+u_j+\varepsilon_{jt}. $$

La pendiente respecto al tiempo es:

$$ \text{pendiente}_{\text{control}}=\beta_1 $$

Esto significa que $\beta_1$ representa el cambio esperado en el crecimiento por cada día adicional para las semillas del grupo control.

Por ejemplo, si:

$$ \beta_1=1.2, $$

el modelo estima que las semillas del grupo control aumentan su crecimiento, en promedio, en 1.2 unidades por día.

Grupo con tratamiento¶

Para las semillas sometidas al tratamiento:

$$ \text{tratamiento}=1. $$

Sustituyendo:

$$ \text{crecimiento}_{jt} = \beta_0 + \beta_1\text{día}_t + \beta_2 + \beta_3\text{día}_t + u_j + \varepsilon_{jt}. $$

Agrupando términos:

$$ \text{crecimiento}_{jt} = (\beta_0+\beta_2) + (\beta_1+\beta_3)\text{día}_t + u_j + \varepsilon_{jt}. $$

Por tanto, la pendiente del grupo con tratamiento es:

$$ \text{pendiente}_{\text{tratamiento}} = \beta_1+\beta_3 $$

El parámetro $\beta_3$ indica entonces cuánto cambia la velocidad de crecimiento debido al tratamiento:

$$ \beta_3= \text{pendiente}_{\text{tratamiento}} - \text{pendiente}_{\text{control}} $$

¿Cómo interpretar $\beta_3$?¶

El signo de $\beta_3$ permite determinar la dirección del efecto del tratamiento sobre la tasa de crecimiento:

  • Si $\beta_3>0$, el grupo con tratamiento crece más rápidamente que el grupo control.
  • Si $\beta_3<0$, el grupo con tratamiento crece más lentamente que el grupo control.
  • Si $\beta_3=0$, ambos grupos tienen la misma pendiente de crecimiento.

Hipótesis estadística de interés¶

Para determinar si la diferencia entre las pendientes puede considerarse estadísticamente distinta de cero se contrasta:

$$ H_0:\beta_3=0 $$

frente a:

$$ H_1:\beta_3\neq0. $$

La hipótesis nula establece que el tratamiento no modifica la pendiente de crecimiento.

Si se rechaza $H_0$, existe evidencia estadística de que la evolución temporal del crecimiento es diferente entre los grupos.

En consecuencia, la interacción día × tratamiento responde directamente a la pregunta:

¿El tratamiento modifica la velocidad con la que crecen las semillas a través del tiempo?

A. Entender los datos antes de modelar¶

Estructura codificada en el archivo¶

Elemento Papel en el caso
Respuesta crecimiento_mm: longitud radicular en milímetros
Tiempo dia: variable cuantitativa, de 1 a 15
Tratamiento 10 niveles; T01_Control funciona como referencia
Agrupación semilla: 10 niveles según las hojas del Excel
Observaciones $10\times15\times10=1500$ registros en formato largo

Los 1500 registros son mediciones, no necesariamente 1500 unidades biológicas independientes.

Unidad experimental, unidad de observación y réplica¶

En este estudio se analiza el crecimiento de semillas sometidas a diferentes tratamientos y observadas repetidamente a través del tiempo.

  • Unidad experimental: corresponde a cada semilla individual, suponiendo que el tratamiento se asigna de manera independiente a cada semilla. Cada semilla constituye una réplica biológica independiente dentro de su tratamiento.

  • Unidad de observación: corresponde a una medición del crecimiento de una semilla en un día determinado. Por tanto, una misma semilla genera varias unidades de observación a lo largo del experimento.

  • Medición repetida: corresponde a observar nuevamente el crecimiento de la misma semilla en otro día. Estas observaciones aportan información sobre su trayectoria de crecimiento, pero no constituyen nuevas réplicas biológicas.

  • Réplica biológica: corresponde a una semilla diferente sometida independientemente al mismo tratamiento. Las distintas semillas permiten estimar la variabilidad biológica existente dentro de cada grupo experimental.

  • Réplica técnica: sería la repetición del procedimiento de medición sobre la misma semilla en el mismo momento, por ejemplo medir dos veces su longitud el mismo día. Estas repeticiones permiten evaluar variabilidad del proceso de medición, pero no aumentan el número de unidades biológicas independientes.

Estructura longitudinal de los datos¶

Si la semilla $j$ es observada en $T$ momentos, produce una secuencia de mediciones:

$$ y_{j1},y_{j2},\ldots,y_{jT}, $$

donde $y_{jt}$ representa el crecimiento observado de la semilla $j$ en el día $t$.

Por ejemplo, si una misma semilla se mide durante cinco días:

$$ y_{j1},\;y_{j2},\;y_{j3},\;y_{j4},\;y_{j5}, $$

tenemos cinco observaciones, pero solamente una unidad experimental.

Las mediciones realizadas sobre una misma semilla suelen ser más parecidas entre sí que las mediciones provenientes de semillas diferentes. En consecuencia:

$$ \operatorname{Corr}(y_{jt},y_{jt'}) \neq 0, \qquad t\neq t'. $$

Por tanto, las observaciones longitudinales de una misma semilla no deben considerarse independientes.

Esta dependencia se incorpora en el modelo mediante el efecto aleatorio de la semilla:

$$ u_j, $$

que permite que cada semilla tenga su propio nivel de crecimiento alrededor del comportamiento promedio de su grupo.

B. Fundamentos de inferencia estadística¶

Conceptos básicos¶

  • Población: conjunto de unidades sobre el que se desea concluir.
  • Muestra: unidades efectivamente observadas.
  • Parámetro: cantidad poblacional desconocida, como una pendiente $\beta_1$ o una varianza $\sigma^2$.
  • Estadístico o estimador: cantidad calculada con la muestra, como $\widehat\beta_1$.
  • Distribución muestral: variación que tendría un estimador si el experimento se repitiera bajo condiciones comparables.
  • Error estándar (EE): desviación estándar de la distribución muestral de un estimador.

Una estimación puntual aislada no expresa su precisión. Por eso se acompaña con un error estándar o un intervalo de confianza. Un intervalo aproximado al 95 % suele escribirse como (asumiendo normalidad)

$$ \widehat\theta \pm 1.96\,\operatorname{EE}(\widehat\theta), $$

cuando la aproximación normal es razonable.

Hipótesis, valor $p$ y errores de decisión¶

Una prueba contrasta dos afirmaciones:

$$ H_0:\ \text{no existe el efecto o la diferencia especificada}, \qquad H_1:\ \text{sí existe}. $$

El valor $p$ es la probabilidad, suponiendo que $H_0$ y el modelo son correctos, de obtener un resultado al menos tan incompatible con $H_0$ como el observado. No es la probabilidad de que $H_0$ sea verdadera.

  • Nivel de significancia $\alpha$: umbral de decisión elegido antes de observar los datos; comúnmente 0.05.
  • Error tipo I: rechazar $H_0$ cuando es verdadera; su probabilidad se controla con $\alpha$.
  • Error tipo II: no rechazar $H_0$ cuando existe un efecto.
  • Potencia: probabilidad de detectar un efecto de magnitud relevante cuando realmente existe.

C. Regresión lineal¶

Modelo lineal simple¶

Para una respuesta cuantitativa $Y$ y un predictor cuantitativo $x$:

$$ Y_i=\beta_0+\beta_1x_i+\varepsilon_i, \qquad \varepsilon_i\sim N(0,\sigma^2). $$
  • $\beta_0$ es la respuesta media esperada cuando $x=0$.
  • $\beta_1$ es el cambio medio esperado en $Y$ por cada unidad adicional de $x$.
  • $\widehat Y_i=\widehat\beta_0+\widehat\beta_1x_i$ es el valor ajustado.
  • $e_i=Y_i-\widehat Y_i$ es el residuo.

En mínimos cuadrados ordinarios se eligen los coeficientes que minimizan

$$ \operatorname{SSE}=\sum_{i=1}^{n}(Y_i-\widehat Y_i)^2. $$

En notación matricial, $\mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon$ y, si $\mathbf X$ tiene rango completo,

$$ \widehat{\boldsymbol\beta} =(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y. $$

Regresión múltiple y variables categóricas¶

Con varios predictores,

$$ Y_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\varepsilon_i. $$

Un tratamiento con $K$ niveles se representa mediante $K-1$ variables indicadoras. Si T01_Control es la referencia,

$$ I_{ik}= \begin{cases} 1,&\text{si la observación }i\text{ pertenece al tratamiento }k,\\ 0,&\text{en otro caso.} \end{cases} $$

D. Modelo lineal mixto¶

Un modelo lineal mixto se representa como:

$$ \mathbf y=\mathbf X\boldsymbol\beta+\mathbf Z\mathbf b+\boldsymbol\varepsilon, $$

donde

$$ \mathbf b\sim N(\mathbf 0,\boldsymbol\Sigma_b), \qquad \boldsymbol\varepsilon\sim N(\mathbf 0,\boldsymbol\Sigma_\varepsilon), \qquad \mathbf b\perp\boldsymbol\varepsilon. $$
  • $\mathbf X\boldsymbol\beta$ contiene los efectos fijos, que describen la trayectoria media poblacional y las comparaciones de interés.
  • $\mathbf Z\mathbf b$ contiene los efectos aleatorios, que describen cómo se apartan las unidades particulares de esa trayectoria.
  • $\boldsymbol\varepsilon$ representa la variación no explicada dentro de la unidad.

¿Por qué no basta con mínimos cuadrados ordinarios?¶

En una regresión lineal convencional se supone:

$$ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, $$

con

$$ \boldsymbol\varepsilon\sim N(\mathbf 0,\sigma^2\mathbf I). $$

La matriz $\mathbf I$ implica que las observaciones son independientes y tienen una varianza común:

$$ \operatorname{Var}(\mathbf y)=\sigma^2\mathbf I. $$

Bajo estas condiciones, los parámetros pueden estimarse mediante mínimos cuadrados ordinarios (OLS) o método QR:

$$ \hat{\boldsymbol\beta}_{OLS} = (\mathbf X^\top\mathbf X)^{-1} \mathbf X^\top\mathbf y $$

Esta solución surge de minimizar la suma de cuadrados:

$$ S(\boldsymbol\beta) = (\mathbf y-\mathbf X\boldsymbol\beta)^\top (\mathbf y-\mathbf X\boldsymbol\beta). $$

¿Qué cambia en el experimento de las semillas?¶

En nuestro caso, una misma semilla se mide repetidamente a través del tiempo.

Por ejemplo, para una semilla $j$ tenemos:

$$ y_{j1},y_{j2},y_{j3},\ldots,y_{jT}. $$

Estas observaciones no son independientes, porque todas provienen de la misma semilla.

El modelo incorpora esta estructura mediante:

$$ \mathbf y = \mathbf X\boldsymbol\beta + \mathbf Z\mathbf b + \boldsymbol\varepsilon, $$

donde:

$$ \mathbf b\sim N(\mathbf 0,\boldsymbol\Sigma_b), \qquad \boldsymbol\varepsilon\sim N(\mathbf 0,\boldsymbol\Sigma_\varepsilon). $$

Los efectos aleatorios $\mathbf b$ hacen que las mediciones de una misma semilla estén relacionadas.

Al integrar los efectos aleatorios:

$$ E(\mathbf y)=\mathbf X\boldsymbol\beta, $$

pero ahora:

$$ \operatorname{Var}(\mathbf y) = \mathbf V = \mathbf Z\boldsymbol\Sigma_b\mathbf Z^\top + \boldsymbol\Sigma_\varepsilon $$

Por tanto, en general:

$$ \mathbf V\neq\sigma^2\mathbf I. $$

Esta es la diferencia fundamental respecto a la regresión lineal ordinaria.

Caso particular: intercepto aleatorio por semilla¶

En nuestro modelo:

$$ \text{crecimiento}_{ijt} = \beta_0 +\beta_1\text{día}_t +\beta_2\text{tratamiento}_i +\beta_3(\text{día}_t\times\text{tratamiento}_i) +u_j +\varepsilon_{ijt}, $$

suponemos:

$$ u_j\sim N(0,\sigma_u^2) $$

y

$$ \varepsilon_{ijt}\sim N(0,\sigma^2). $$

Dos observaciones de semillas diferentes tienen covarianza cero:

$$ \operatorname{Cov}(Y_{jt},Y_{ks})=0, \qquad j\neq k. $$

Sin embargo, dos mediciones realizadas sobre la misma semilla comparten el mismo $u_j$:

$$ Y_{jt}=\cdots+u_j+\varepsilon_{jt}, $$$$ Y_{js}=\cdots+u_j+\varepsilon_{js}. $$

Por tanto:

$$ \operatorname{Cov}(Y_{jt},Y_{js}) = \sigma_u^2 \qquad t\neq s. $$

Esta correlación aparece simplemente porque ambas mediciones contienen el mismo efecto aleatorio $u_j$.

La correlación intraclase (ICC) es:

$$ \rho= \frac{\sigma_u^2} {\sigma_u^2+\sigma^2} $$

y cuantifica qué tan semejantes son las observaciones realizadas sobre una misma semilla.

¿Qué ocurriría si ignoramos esta dependencia?¶

Podríamos aplicar OLS tratando cada combinación semilla × día como si fuera una observación independiente.

El problema es que estaríamos suponiendo incorrectamente:

$$ \operatorname{Cov}(Y_{jt},Y_{js})=0. $$

El modelo ignoraría que ambas observaciones pertenecen a la misma semilla.

En consecuencia, aunque bajo ciertas condiciones los coeficientes OLS de los efectos fijos pueden seguir siendo insesgados, sus errores estándar, intervalos de confianza y pruebas de hipótesis pueden ser incorrectos, además de que OLS no estima los componentes de variabilidad entre y dentro de semillas.

Por ello, el problema no es simplemente estimar $\beta_0,\beta_1,\beta_2,\beta_3$. También debemos estimar:

$$ \sigma_u^2 \qquad\text{y}\qquad \sigma^2. $$

Es decir, debemos determinar simultáneamente:

  1. la trayectoria promedio de crecimiento;
  2. el efecto del tratamiento;
  3. la diferencia entre semillas;
  4. la variabilidad dentro de cada semilla.

De OLS a mínimos cuadrados generalizados¶

Si conociéramos la matriz de covarianza $\mathbf V$, los efectos fijos podrían estimarse mediante mínimos cuadrados generalizados (GLS):

$$ \hat{\boldsymbol\beta} = (\mathbf X^\top\mathbf V^{-1}\mathbf X)^{-1} \mathbf X^\top\mathbf V^{-1}\mathbf y $$

Obsérvese que esta expresión es similar a OLS, pero incorpora $\mathbf V^{-1}$ para considerar la dependencia entre las observaciones.

El problema es que:

$$ \mathbf V\text{ es desconocida} $$

porque depende de parámetros desconocidos como:

$$ \sigma_u^2 \quad\text{y}\quad \sigma^2. $$

¿Por qué aparecen métodos numéricos?¶

Para resolver el problema se utilizan normalmente máxima verosimilitud (ML) o máxima verosimilitud restringida (REML).

Como:

$$ \mathbf y \sim N(\mathbf X\boldsymbol\beta,\mathbf V), $$

la log-verosimilitud contiene términos de la forma:

$$ \ell = -\frac{1}{2} \left[ \log|\mathbf V| + (\mathbf y-\mathbf X\boldsymbol\beta)^\top \mathbf V^{-1} (\mathbf y-\mathbf X\boldsymbol\beta) \right] +C. $$

La dificultad está en que los parámetros de varianza aparecen dentro de:

$$ \mathbf V^{-1} \qquad\text{y}\qquad \log|\mathbf V|. $$

Por ello, a diferencia del OLS convencional, en general no existe una fórmula cerrada única que permita obtener directamente todos los parámetros del modelo mixto.

El procedimiento se vuelve iterativo:

  1. Se proponen valores iniciales para los componentes de varianza.
  2. Con ellos se construye $\mathbf V$.
  3. Se estiman los efectos fijos $\boldsymbol\beta$.
  4. Se evalúa la verosimilitud ML o REML.
  5. El algoritmo modifica los componentes de varianza.
  6. Se repite el proceso hasta alcanzar convergencia.

De manera esquemática:

$$ (\sigma_u^2,\sigma^2)^{(0)} \rightarrow \mathbf V^{(0)} \rightarrow \hat{\boldsymbol\beta}^{(0)} \rightarrow \ell^{(0)} $$$$ \downarrow $$$$ (\sigma_u^2,\sigma^2)^{(1)} \rightarrow \mathbf V^{(1)} \rightarrow \hat{\boldsymbol\beta}^{(1)} \rightarrow \ell^{(1)} $$$$ \downarrow $$$$ \cdots \rightarrow \text{convergencia}. $$

Para esta optimización pueden utilizarse algoritmos numéricos como Newton--Raphson, Fisher scoring, EM, BFGS o variantes relacionadas, dependiendo del software utilizado.

En suma¶

El paso de regresión lineal a modelo lineal mixto puede resumirse como:

$$ \underbrace{\sigma^2\mathbf I}_{\text{OLS}} \quad\longrightarrow\quad \underbrace{ \mathbf Z\boldsymbol\Sigma_b\mathbf Z^\top+ \boldsymbol\Sigma_\varepsilon }_{\text{modelo mixto}} $$

El modelo mixto debe estimar no solamente los coeficientes:

$$ \beta_0,\beta_1,\beta_2,\beta_3, $$

sino también la estructura de dependencia y los componentes de varianza.

E. Implementación¶

In [7]:
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from scipy.stats import chi2
from statsmodels.stats.multicomp import pairwise_tukeyhsd

¿Qué función cumple cada biblioteca?¶

Biblioteca Uso en el caso
pathlib Maneja la ruta del archivo Excel.
numpy Realiza operaciones numéricas
pandas Lee, transforma, agrupa y resume los datos.
matplotlib Genera las gráficas.
statsmodels Ajusta el modelo lineal mixto y el contraste de Tukey.
scipy.stats.chi2 Convierte el estadístico de razón de verosimilitud en un valor $p$.

Las importaciones no modifican los datos: únicamente ponen las funciones a disposición del notebook.

1. Cargar el archivo Excel¶

El archivo debe estar en la misma carpeta que este notebook.

Diccionario de datos y comprobaciones esperadas¶

Cada hoja Semilla_01, ..., Semilla_10 contiene:

  • una columna Dia con 15 mediciones;
  • diez columnas de tratamiento, de T01_Control a T10_PosCtrl;
  • columnas auxiliares posteriores que el notebook no incorpora porque lee únicamente A:K y las primeras 15 filas.

Antes de analizar datos reales se comprobarían, al menos: nombres y unidades, ausencia de duplicados, valores faltantes, rango válido de días, número de observaciones por combinación y unicidad del identificador de la unidad experimental.

In [8]:
from google.colab import files
uploaded = files.upload()

#buscar "datos_crecimiento_semillas.xlsx"
Upload widget is only available when the cell has been executed in the current browser session. Please rerun this cell to enable. 
Saving datos_crecimiento_semillas.xlsx to datos_crecimiento_semillas (1).xlsx
In [9]:
EXCEL_PATH = Path("datos_crecimiento_semillas.xlsx")
xlsx = pd.ExcelFile(EXCEL_PATH)
xlsx.sheet_names
Out[9]:
['Semilla_01',
 'Semilla_02',
 'Semilla_03',
 'Semilla_04',
 'Semilla_05',
 'Semilla_06',
 'Semilla_07',
 'Semilla_08',
 'Semilla_09',
 'Semilla_10']

2. Transformar de formato ancho a formato largo¶

El formato original tiene una hoja por semilla.
Cada hoja tiene los días en filas y los tratamientos en columnas. Para el modelo mixto se necesita una tabla en formato largo:

| semilla | dia | tratamiento | crecimiento_mm | |---|---:|---|---:|

¿Por qué transformar de ancho a largo?¶

En formato ancho, cada tratamiento ocupa una columna. En formato largo, cada fila representa una medición y una columna indica el tratamiento. El formato largo es apropiado para fórmulas estadísticas porque separa claramente:

Variable Tipo
semilla identificador de agrupación
dia predictor cuantitativo
tratamiento predictor categórico
crecimiento_mm respuesta cuantitativa
In [10]:
# Lista donde guardaremos los datos transformados de cada semilla
long_parts = []

# Cada hoja del Excel corresponde a una semilla
for sheet in xlsx.sheet_names:

    # Leer las primeras 15 filas y las columnas A:K de la hoja
    tmp = pd.read_excel(
        EXCEL_PATH,
        sheet_name=sheet,
        usecols="A:K",
        nrows=15
    )

    # Guardar el nombre de la hoja como identificador de la semilla
    tmp["semilla"] = sheet

    # Transformar los datos de formato ancho a formato largo:
    # cada tratamiento pasa a ser una fila
    tmp_long = tmp.melt(
        id_vars=["semilla", "Dia"],     # columnas que se conservan
        var_name="tratamiento",         # nombre de la nueva columna
        value_name="crecimiento_mm"     # valores observados
    )

    # Guardar los datos de esta semilla
    long_parts.append(tmp_long)
In [11]:
# Unir todas las semillas en una sola tabla
data = pd.concat(long_parts, ignore_index=True)

# Cambiar "Dia" a "dia" para mantener nombres consistentes
data = data.rename(columns={"Dia": "dia"})
In [12]:
# Ordenar tratamientos como categorías
treatment_order = list(pd.read_excel(EXCEL_PATH, sheet_name=xlsx.sheet_names[0], usecols="B:K", nrows=0).columns)
treatment_order
Out[12]:
['T01_Control',
 'T02_5mgL',
 'T03_10mgL',
 'T04_25mgL',
 'T05_50mgL',
 'T06_75mgL',
 'T07_100mgL',
 'T08_150mgL',
 'T09_200mgL',
 'T10_PosCtrl']
In [13]:
data["tratamiento"] = pd.Categorical(data["tratamiento"], categories=treatment_order, ordered=True)
data.head()
Out[13]:
semilla dia tratamiento crecimiento_mm
0 Semilla_01 1 T01_Control 2.58
1 Semilla_01 2 T01_Control 3.16
2 Semilla_01 3 T01_Control 3.84
3 Semilla_01 4 T01_Control 4.40
4 Semilla_01 5 T01_Control 4.87
In [14]:
print(data.shape)
data.groupby(["semilla", "tratamiento"], observed=False).size().head(12)

(1500, 4)
Out[14]:
0
semilla tratamiento
Semilla_01 T01_Control 15
T02_5mgL 15
T03_10mgL 15
T04_25mgL 15
T05_50mgL 15
T06_75mgL 15
T07_100mgL 15
T08_150mgL 15
T09_200mgL 15
T10_PosCtrl 15
Semilla_02 T01_Control 15
T02_5mgL 15

3. Exploración gráfica¶

Primero se grafican las trayectorias promedio de crecimiento por tratamiento.
Las bandas representan aproximadamente ±1 error estándar.

Desviación estándar, error estándar y banda gráfica¶

Para cada día y tratamiento, el código calcula

$$ \bar y=\frac{1}{n}\sum_{j=1}^{n}y_j, \qquad s=\sqrt{\frac{1}{n-1}\sum_{j=1}^{n}(y_j-\bar y)^2}, \qquad \operatorname{EE}(\bar y)=\frac{s}{\sqrt n}. $$
  • $s$ describe la dispersión de observaciones individuales.
  • $\operatorname{EE}(\bar y)$ describe la precisión de la media estimada.
  • La banda de la figura es $\bar y\pm1\operatorname{EE}$; no es un intervalo de confianza al 95 %.

Las curvas permiten examinar dirección, forma, separación y posibles interacciones. No sustituyen la inferencia formal y, al resumir por día, tampoco muestran toda la dependencia longitudinal.

In [15]:
summary = (
    data
    .groupby(["dia", "tratamiento"], observed=False)
    .agg(
        media=("crecimiento_mm", "mean"),
        sd=("crecimiento_mm", "std"),
        n=("crecimiento_mm", "size")
    )
    .reset_index()
)

summary["se"] = summary["sd"] / np.sqrt(summary["n"])
summary
Out[15]:
dia tratamiento media sd n se
0 1 T01_Control 2.494 0.526987 10 0.166648
1 1 T02_5mgL 2.496 0.416205 10 0.131616
2 1 T03_10mgL 2.524 0.396238 10 0.125301
3 1 T04_25mgL 2.704 0.560044 10 0.177101
4 1 T05_50mgL 2.810 0.399944 10 0.126474
... ... ... ... ... ... ...
145 15 T06_75mgL 16.863 1.118432 10 0.353679
146 15 T07_100mgL 18.674 1.155049 10 0.365259
147 15 T08_150mgL 17.405 1.307459 10 0.413455
148 15 T09_200mgL 7.909 1.330066 10 0.420604
149 15 T10_PosCtrl 15.168 0.745234 10 0.235664

150 rows × 6 columns

In [16]:
plt.figure(figsize=(11, 6))

for trt in treatment_order:
    g = summary[summary["tratamiento"] == trt]
    x = g["dia"].to_numpy()
    y = g["media"].to_numpy()
    se = g["se"].to_numpy()

    plt.plot(x, y, marker="o", linewidth=1.8, label=trt)
    plt.fill_between(x, y - se, y + se, alpha=0.15)

plt.xlabel("Día")
plt.ylabel("Crecimiento / longitud radicular (mm)")
plt.title("Crecimiento promedio por tratamiento")
plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.show()

Pausa de interpretación¶

Observe la figura y responda antes de continuar:

  1. ¿Las trayectorias son aproximadamente lineales entre los días 1 y 15?
  2. ¿Qué tratamiento parece tener la pendiente mayor?
  3. ¿Existe una dosis a partir de la cual el crecimiento deja de aumentar?
  4. ¿Las diferencias parecen constantes o se amplían con el tiempo?
  5. ¿Las bandas representan variación entre semillas o incertidumbre de la media?

La gráfica sugiere que T07_100mgL presenta la trayectoria más pronunciada y que T09_200mgL queda por debajo del control al avanzar el tiempo. El modelo cuantificará (estadísticamente) esas diferencias.

4. Modelo mixto de efectos lineales¶

Se ajusta un modelo con:

  • efecto fijo de dia,
  • efecto fijo de tratamiento,
  • interacción dia × tratamiento,
  • intercepto aleatorio por semilla.

La interacción es la parte más importante si se desea saber si los tratamientos cambian la velocidad o trayectoria de crecimiento.

Supuestos específicos del modelo ajustado¶

  1. La relación media entre crecimiento y día es lineal dentro de cada tratamiento.
  2. Los interceptos aleatorios siguen aproximadamente una distribución normal.
  3. Los residuos tienen media cero y varianza constante.
  4. Una vez considerado el intercepto aleatorio, los residuos se tratan como independientes.
  5. Las semillas o grupos son independientes entre sí.
  6. El identificador de agrupación corresponde a la unidad experimental pertinente.

La inferencia debe presentarse como válida bajo esta especificación y acompañarse, en una aplicación real, de diagnósticos de residuos y análisis de sensibilidad.

In [17]:
# Definir una función para ajustar un modelo lineal mixto.
# `formula` especifica la relación estadística que se desea modelar.
# `df` es la tabla de datos que se utilizará.
def fit_mixedlm(formula, df):

    # Crear la estructura del modelo lineal mixto.
    model = smf.mixedlm(
        formula,              # Fórmula del modelo, por ejemplo:
                              # "crecimiento_mm ~ dia * C(tratamiento)".
        data=df,              # DataFrame que contiene las variables de la fórmula.
        groups=df["semilla"], # Indica que las observaciones se agrupan por semilla.
                              # Las mediciones de una misma semilla están relacionadas.
        re_formula="1"        # Especifica un intercepto aleatorio por semilla.
                              # Cada semilla puede iniciar con un nivel propio de crecimiento.
    )

    # Estimar los parámetros del modelo.
    result = model.fit(
        reml=False,           # Usar máxima verosimilitud (ML), no REML.
                              # ML permite comparar modelos con diferente número de
                              # efectos fijos mediante pruebas de razón de verosimilitud.
        method="powell",      # Usar el algoritmo numérico Powell para la optimización.
                              # Suele ser estable en conjuntos de datos pequeños.
        maxiter=500,          # Permitir hasta 500 iteraciones al algoritmo.
        disp=False            # No mostrar en pantalla el avance de la optimización.
    )

    # Devolver el modelo ya ajustado, con coeficientes, errores estándar y diagnósticos.
    return result
In [18]:
# Ajustar el modelo principal utilizando la función definida anteriormente.
modelo_full = fit_mixedlm(
    # Variable respuesta: crecimiento_mm.
    # Predictores: dia, tratamiento y su interacción.
    #
    # `C(tratamiento)` indica que tratamiento debe analizarse como variable categórica.
    # El símbolo `*` equivale a incluir:
    #   dia + C(tratamiento) + dia:C(tratamiento)
    "crecimiento_mm ~ dia * C(tratamiento)",

    # Tabla de datos en formato largo, con una fila por medición.
    data
)
In [19]:
print(modelo_full.summary())

                    Mixed Linear Model Regression Results
=============================================================================
Model:                   MixedLM      Dependent Variable:      crecimiento_mm
No. Observations:        1500         Method:                  ML            
No. Groups:              10           Scale:                   0.1772        
Min. group size:         150          Log-Likelihood:          -858.4480     
Max. group size:         150          Converged:               Yes           
Mean group size:         150.0                                               
-----------------------------------------------------------------------------
                                  Coef.  Std.Err.    z    P>|z| [0.025 0.975]
-----------------------------------------------------------------------------
Intercept                          1.859    0.191   9.723 0.000  1.484  2.234
C(tratamiento)[T.T02_5mgL]        -0.082    0.102  -0.800 0.424 -0.282  0.119
C(tratamiento)[T.T03_10mgL]       -0.056    0.102  -0.544 0.587 -0.256  0.145
C(tratamiento)[T.T04_25mgL]       -0.027    0.102  -0.267 0.789 -0.228  0.173
C(tratamiento)[T.T05_50mgL]       -0.062    0.102  -0.605 0.545 -0.262  0.139
C(tratamiento)[T.T06_75mgL]       -0.020    0.102  -0.199 0.842 -0.221  0.180
C(tratamiento)[T.T07_100mgL]       0.129    0.102   1.262 0.207 -0.071  0.330
C(tratamiento)[T.T08_150mgL]      -0.035    0.102  -0.338 0.735 -0.235  0.166
C(tratamiento)[T.T09_200mgL]      -0.242    0.102  -2.367 0.018 -0.443 -0.042
C(tratamiento)[T.T10_PosCtrl]      0.051    0.102   0.495 0.621 -0.150  0.251
dia                                0.549    0.008  69.036 0.000  0.534  0.565
dia:C(tratamiento)[T.T02_5mgL]     0.092    0.011   8.172 0.000  0.070  0.114
dia:C(tratamiento)[T.T03_10mgL]    0.177    0.011  15.699 0.000  0.155  0.199
dia:C(tratamiento)[T.T04_25mgL]    0.281    0.011  24.975 0.000  0.259  0.303
dia:C(tratamiento)[T.T05_50mgL]    0.373    0.011  33.178 0.000  0.351  0.395
dia:C(tratamiento)[T.T06_75mgL]    0.439    0.011  39.050 0.000  0.417  0.461
dia:C(tratamiento)[T.T07_100mgL]   0.563    0.011  50.067 0.000  0.541  0.585
dia:C(tratamiento)[T.T08_150mgL]   0.488    0.011  43.380 0.000  0.466  0.510
dia:C(tratamiento)[T.T09_200mgL]  -0.123    0.011 -10.920 0.000 -0.145 -0.101
dia:C(tratamiento)[T.T10_PosCtrl]  0.337    0.011  29.921 0.000  0.315  0.359
Group Var                          0.313    0.335                            
=============================================================================

Interpretación del modelo ajustado¶

El modelo usa T01_Control como grupo de referencia. Por ello, los coeficientes de dia:C(tratamiento) no son pendientes completas: representan la diferencia en la pendiente de cada tratamiento respecto al control.

La ecuación media estimada para el control es:

$$ \widehat{\text{crecimiento}}_{\text{Control}} = 1.859 + 0.549\,\text{día}. $$

Así, para una semilla con efecto aleatorio igual a cero, el crecimiento esperado en el control aumenta aproximadamente:

$$ 0.549\ \text{mm por día}. $$

El intercepto $1.859$ corresponde al día 0.

Ejemplo: tratamiento T07_100mgL¶

Para T07_100mgL, el modelo estima:

$$ \widehat{\text{crecimiento}}_{\text{T07}} = 1.859 + 0.129 + (0.549+0.563)\,\text{día}. $$

Por tanto:

$$ \widehat{\text{crecimiento}}_{\text{T07}} = 1.988+1.112\,\text{día}. $$

La pendiente total de este tratamiento es:

$$ 1.112\ \text{mm por día}. $$

Esto significa que T07_100mgL crece, en promedio, aproximadamente $1.112$ mm por día, mientras que el control crece $0.549$ mm por día.

La diferencia entre ambas tasas es:

$$ 1.112-0.549=0.563\ \text{mm/día}. $$

Como el valor $p$ de la interacción es menor que $0.001$, existe evidencia muy fuerte de que T07_100mgL incrementa la tasa de crecimiento respecto al control.

Es importante notar que el coeficiente de intercepto de T07_100mgL ($0.129$, $p=0.207$) no es significativo. Es decir, el modelo no detecta una diferencia clara respecto al control en el día 0; la diferencia aparece y aumenta conforme transcurren los días.

Pendientes estimadas de todos los tratamientos¶

Tratamiento Pendiente estimada (mm/día) Diferencia respecto al control
T01_Control $0.549$ $0.000$
T02_5mgL $0.549+0.092=0.641$ $+0.092$
T03_10mgL $0.549+0.177=0.726$ $+0.177$
T04_25mgL $0.549+0.281=0.830$ $+0.281$
T05_50mgL $0.549+0.373=0.922$ $+0.373$
T06_75mgL $0.549+0.439=0.988$ $+0.439$
T07_100mgL $0.549+0.563=1.112$ $+0.563$
T08_150mgL $0.549+0.488=1.037$ $+0.488$
T09_200mgL $0.549-0.123=0.426$ $-0.123$
T10_PosCtrl $0.549+0.337=0.886$ $+0.337$

Todas las diferencias de pendiente mostradas en la tabla son estadísticamente significativas ($p<0.001$).

Caso contrastante: T09_200mgL¶

Para T09_200mgL, la ecuación estimada es:

$$ \widehat{\text{crecimiento}}_{\text{T09}} = (1.859-0.242) + (0.549-0.123)\,\text{día}. $$

Por tanto:

$$ \widehat{\text{crecimiento}}_{\text{T09}} = 1.617+0.426\,\text{día}. $$

Este tratamiento presenta dos diferencias respecto al control:

  1. En el día 0, su valor esperado es menor por aproximadamente $0.242$ mm ($p=0.018$).
  2. Su crecimiento diario es menor por aproximadamente $0.123$ mm/día ($p<0.001$).

En consecuencia, la diferencia esperada entre T09_200mgL y el control en el día $t$ es:

$$ \Delta_{\text{T09}}(t) = -0.242-0.123t. $$

La desventaja estimada de este tratamiento aumenta conforme avanza el experimento.

Significado de las columnas inferenciales¶

Columna Interpretación
Coef. estimación puntual del parámetro
Std.Err. error estándar de la estimación
z cociente aproximado Coef./Std.Err.
`P> z ` valor $p$ bilateral para $H_0:\beta=0$
[0.025, 0.975] intervalo de confianza aproximado al 95 %
Group Var varianza estimada de los interceptos entre semillas
Scale varianza residual estimada

Con Group Var = 0.313 y Scale = 0.1772, la ICC aproximada basada en los valores redondeados es

$$ \frac{0.313}{0.313+0.1772}\approx0.64. $$

Esto indica agrupación sustancial según el modelo: ignorar la identidad de semilla produciría errores estándar inapropiados. La estimación descansa en solo 10 grupos, por lo que conviene ser prudente con aproximaciones asintóticas y componentes de varianza.

Interpretación de la variabilidad entre semillas¶

El modelo estima una varianza del intercepto aleatorio de:

$$ \widehat{\sigma_u^2}=0.313, $$

mientras que la varianza residual es:

$$ \widehat{\sigma^2}=0.177. $$

Esto indica que existe variabilidad apreciable entre las semillas o unidades agrupadas. Dos mediciones de la misma semilla tienden a parecerse entre sí, por lo que no es apropiado tratarlas como observaciones completamente independientes.

La proporción aproximada de variabilidad atribuible a diferencias entre semillas es:

$$ \text{ICC} = \frac{0.313}{0.313+0.177} \approx 0.64. $$

Es decir, aproximadamente el $64\%$ de la variabilidad marginal estimada se asocia con diferencias persistentes entre semillas, y el resto con variación residual dentro de cada semilla.

5. Pruebas de razón de verosimilitud¶

Se comparan modelos anidados:

  1. Modelo completo: dia * tratamiento
  2. Modelo sin interacción: dia + tratamiento
  3. Modelo solo con día: dia

Interpretación:

  • Si la interacción dia × tratamiento es significativa, los tratamientos generan trayectorias de crecimiento diferentes.
  • Si el efecto de tratamiento es significativo, hay diferencias promedio entre tratamientos.

¿Qué significa que los modelos estén anidados?¶

Un modelo reducido está anidado en el completo si puede obtenerse fijando ciertos parámetros del completo en valores específicos —aquí, cero—:

Completo:        dia + tratamiento + dia:tratamiento
Sin interacción: dia + tratamiento
Solo día:        dia

Para que la comparación sea coherente deben usarse los mismos datos, la misma respuesta, la misma estructura aleatoria y el mismo método ML.

Prueba de razón de verosimilitud (LRT)¶

La prueba de razón de verosimilitud permite comparar dos modelos anidados: un modelo reducido, con menos parámetros, y un modelo completo, que incorpora términos adicionales.

Se calcula:

$$ D=2\left( \ell_{\text{completo}} - \ell_{\text{reducido}} \right), $$

donde $\ell$ es la log-verosimilitud de cada modelo. Si los términos adicionales mejoran el ajuste, el modelo completo tendrá una log-verosimilitud mayor y (D) será grande.

Bajo la hipótesis nula y condiciones regulares:

$$ D\ \dot\sim\ \chi^2_{\Delta q}, $$

donde $\Delta q$ es el número de parámetros adicionales del modelo completo respecto al reducido.

Primer contraste: ¿las pendientes difieren entre tratamientos?¶

El modelo completo contiene:

  • el efecto de dia;
  • las diferencias entre tratamientos;
  • las interacciones dia × tratamiento.

El modelo reducido elimina las interacciones. Por tanto, supone que todos los tratamientos tienen la misma pendiente temporal.

La hipótesis es:

$$ H_0: \delta_2=\delta_3=\cdots=\delta_{10}=0, $$

donde cada $\delta_k$ representa la diferencia entre la pendiente del tratamiento $k$ y la pendiente del control.

Bajo $H_0$, todos los tratamientos cambian a la misma velocidad:

$$ \text{pendiente}_{1} = \text{pendiente}_{2} = \cdots = \text{pendiente}_{10}. $$

Si el valor $p$ de la LRT es pequeño, se rechaza $H_0$. Existe entonces evidencia de que al menos un tratamiento presenta una trayectoria temporal o tasa de crecimiento diferente respecto al control.

Segundo contraste: ¿persisten diferencias entre tratamientos?¶

Si el primer contraste indica que las interacciones no son necesarias, se ajusta un modelo sin dia × tratamiento. En ese modelo todos los grupos comparten una pendiente, pero todavía pueden diferir en su nivel medio de crecimiento.

Entonces se compara:

  • un modelo que incluye C(tratamiento), y
  • un modelo reducido que elimina C(tratamiento).

La hipótesis es:

$$ H_0: \gamma_2=\gamma_3=\cdots=\gamma_{10}=0. $$

Esto equivale a afirmar que, una vez asumida una pendiente común, no existen diferencias sistemáticas entre los tratamientos.

Si se rechaza $H_0$, al menos uno de los tratamientos presenta un crecimiento medio distinto del control, aunque todos evolucionen con la misma pendiente a través de los días.

In [20]:
modelo_sin_interaccion = fit_mixedlm(
    "crecimiento_mm ~ dia + C(tratamiento)",
    data
)

modelo_solo_dia = fit_mixedlm(
    "crecimiento_mm ~ dia",
    data
)

def lr_test(full_model, reduced_model, label):
    lr_stat = 2 * (full_model.llf - reduced_model.llf)
    df_diff = int(full_model.df_modelwc - reduced_model.df_modelwc)
    p_value = chi2.sf(lr_stat, df_diff)
    return {
        "comparacion": label,
        "LR_stat": lr_stat,
        "df": df_diff,
        "p_value": p_value
    }

lr_results = pd.DataFrame([
    lr_test(modelo_full, modelo_sin_interaccion, "Interacción día × tratamiento"),
    lr_test(modelo_sin_interaccion, modelo_solo_dia, "Efecto global de tratamiento")
])

lr_results
Out[20]:
comparacion LR_stat df p_value
0 Interacción día × tratamiento 2597.147907 9 0.0
1 Efecto global de tratamiento 2085.174342 9 0.0

6. Identificar el tratamiento con mayor efecto¶

Una manera práctica de resumir el efecto es comparar el crecimiento esperado entre el día 1 y el día 15 usando el modelo ajustado.

$$ \Delta_k=E(Y\mid d=15,k)-E(Y\mid d=1,k). $$

La expresión representa el cambio esperado en la respuesta Y para el tratamiento k, desde el día 1 hasta el día 15.

In [21]:
# Crear una tabla con todas las combinaciones de tratamiento y día
# que se usarán para obtener predicciones del modelo.
pred_grid = pd.DataFrame({

    # Repetir cada tratamiento dos veces:
    # una fila para el día 1 y otra para el día 15.
    "tratamiento": np.repeat(treatment_order, 2),

    # Repetir la secuencia [1, 15] para cada tratamiento.
    "dia": np.tile([1, 15], len(treatment_order))
})
In [22]:
# Convertir "tratamiento" en una variable categórica ordenada.
# Esto asegura que los nombres y el orden de los tratamientos
# coincidan con los usados al ajustar el modelo.
pred_grid["tratamiento"] = pd.Categorical(
    pred_grid["tratamiento"],   # Columna que se convertirá a categoría.
    categories=treatment_order, # Orden deseado de las categorías.
    ordered=True                # Indica que las categorías tienen un orden definido.
)
In [23]:
# Calcular la predicción media del crecimiento para cada combinación
# tratamiento–día. `predict()` usa los efectos fijos del modelo:
# representa la trayectoria de una semilla típica (efecto aleatorio = 0).
pred_grid["prediccion_mm"] = modelo_full.predict(pred_grid)
In [24]:
# Transformar la tabla de formato largo a formato ancho:
# cada tratamiento tendrá una sola fila y columnas separadas
# para las predicciones del día 1 y del día 15.
pred_wide = (
    pred_grid

    # Usar tratamiento como filas, día como columnas y predicción como valores.
    .pivot(index="tratamiento", columns="dia", values="prediccion_mm")

    # Convertir el índice "tratamiento" nuevamente en una columna normal.
    .reset_index()

    # Dar nombres más descriptivos a las columnas de los días.
    .rename(columns={1: "pred_dia_1", 15: "pred_dia_15"})
)
In [25]:
# Eliminar el nombre técnico "dia" que pandas conserva sobre las columnas.
pred_wide.columns.name = None


# Calcular el crecimiento esperado acumulado entre el día 1 y el día 15.
# Un valor positivo indica crecimiento durante el periodo.
pred_wide["incremento_predicho_15_1"] = (
    pred_wide["pred_dia_15"] - pred_wide["pred_dia_1"]
)
In [26]:
# Entre el día 1 y el día 15 transcurren 14 intervalos diarios.
# Dividir el incremento acumulado entre 14 produce la pendiente
# o tasa promedio de crecimiento predicha en mm por día.
pred_wide["pendiente_mm_dia"] = (
    pred_wide["incremento_predicho_15_1"] / 14
)
In [27]:
# Extraer el crecimiento acumulado predicho para el tratamiento control.
#
# `.loc[...]` selecciona la fila cuyo tratamiento es "T01_Control"
# y únicamente la columna "incremento_predicho_15_1".
# `.iloc[0]` obtiene el primer y único valor seleccionado.
control_incremento = pred_wide.loc[
    pred_wide["tratamiento"] == "T01_Control",
    "incremento_predicho_15_1"
].iloc[0]
In [28]:
# Calcular cuánto mayor o menor es el crecimiento acumulado de cada
# tratamiento respecto al control.
#
# Valor positivo: mayor crecimiento que el control.
# Valor negativo: menor crecimiento que el control.
# Para el propio control, el resultado será cero.
pred_wide["diferencia_vs_control"] = (
    pred_wide["incremento_predicho_15_1"] - control_incremento
)
In [29]:
# Ordenar los tratamientos desde el mayor hasta el menor crecimiento
# acumulado predicho entre los días 1 y 15.
ranking_modelo = (
    pred_wide
    .sort_values("incremento_predicho_15_1", ascending=False)

    # Crear un índice nuevo consecutivo: 0, 1, 2, ...
    .reset_index(drop=True)
)

# Mostrar la tabla final con el ranking de tratamientos.
ranking_modelo
Out[29]:
tratamiento pred_dia_1 pred_dia_15 incremento_predicho_15_1 pendiente_mm_dia diferencia_vs_control
0 T07_100mgL 3.100375 18.674825 15.57445 1.112461 7.88575
1 T08_150mgL 2.861508 17.382758 14.52125 1.037232 6.83255
2 T06_75mgL 2.827000 16.666200 13.83920 0.988514 6.15050
3 T05_50mgL 2.719400 15.633800 12.91440 0.922457 5.22570
4 T10_PosCtrl 2.795258 15.196608 12.40135 0.885811 4.71265
5 T04_25mgL 2.661667 14.284067 11.62240 0.830171 3.93370
6 T03_10mgL 2.529058 12.690408 10.16135 0.725811 2.47265
7 T02_5mgL 2.418142 11.393992 8.97585 0.641132 1.28715
8 T01_Control 2.408050 10.096750 7.68870 0.549193 0.00000
9 T09_200mgL 2.043092 8.011842 5.96875 0.426339 -1.71995
In [30]:
plt.figure(figsize=(9, 5))
plt.bar(
    ranking_modelo["tratamiento"].astype(str),
    ranking_modelo["incremento_predicho_15_1"]
)
plt.xticks(rotation=45, ha="right")
plt.ylabel("Incremento predicho día 15 - día 1 (mm)")
plt.title("Ranking de tratamientos según incremento predicho")
plt.tight_layout()
plt.show()

7. Análisis complementario¶

Con 10 tratamientos existen

$$ \binom{10}{2}=45 $$

comparaciones por pares. Si cada una se probara a nivel 0.05 sin ajuste, aumentaría la probabilidad de al menos un falso positivo. Tukey HSD controla la tasa de error familiar (FWER) para el conjunto de comparaciones.

En la salida:

  • meandiff es la diferencia group2 - group1;
  • lower y upper forman el intervalo simultáneo;
  • p-adj es el valor $p$ ajustado;
  • reject=True indica rechazo de igualdad al nivel familiar 0.05.

Si el intervalo incluye cero, no se detecta una diferencia compatible con el control familiar establecido.

In [31]:
# Crear una lista vacía para guardar una pendiente por cada
# combinación de semilla y tratamiento.
slopes = []


# Agrupar la tabla por semilla y tratamiento.
#
# En cada iteración:
# - `seed` contiene el nombre o identificador de la semilla;
# - `trt` contiene el nombre del tratamiento;
# - `g` contiene únicamente las 15 observaciones de esa combinación.
#
# `observed=True` evita incluir combinaciones categóricas que no aparecen
# realmente en los datos.
for (seed, trt), g in data.groupby(
    ["semilla", "tratamiento"],
    observed=True
):

    # Ajustar una recta a los datos de esta semilla y tratamiento:
    #
    # crecimiento_mm = intercepto + pendiente × dia
    #
    # `np.polyfit(..., 1)` ajusta un polinomio de grado 1, es decir,
    # una regresión lineal simple. Devuelve:
    # [pendiente, intercepto]
    #
    # `[0]` extrae solamente la pendiente, en mm por día.
    slope = np.polyfit(
        g["dia"],
        g["crecimiento_mm"],
        1
    )[0]

    # Guardar los resultados como un diccionario:
    # identificador de semilla, tratamiento y pendiente estimada.
    slopes.append({
        "semilla": seed,
        "tratamiento": trt,
        "pendiente_mm_dia": slope
    })


# Convertir la lista de diccionarios en una tabla de pandas.
# Cada fila representa una combinación semilla × tratamiento.
slopes_df = pd.DataFrame(slopes)


# Mostrar la tabla final de pendientes individuales.
slopes_df
Out[31]:
semilla tratamiento pendiente_mm_dia
0 Semilla_01 T01_Control 0.593643
1 Semilla_01 T02_5mgL 0.679500
2 Semilla_01 T03_10mgL 0.807179
3 Semilla_01 T04_25mgL 0.920179
4 Semilla_01 T05_50mgL 0.988929
... ... ... ...
95 Semilla_10 T06_75mgL 1.110536
96 Semilla_10 T07_100mgL 1.215036
97 Semilla_10 T08_150mgL 1.142036
98 Semilla_10 T09_200mgL 0.518464
99 Semilla_10 T10_PosCtrl 0.963786

100 rows × 3 columns

In [32]:
# Agrupar las pendientes individuales por tratamiento.
# Cada grupo reúne las pendientes estimadas de todas las semillas
# que recibieron el mismo tratamiento.
ranking_pendientes = (
    slopes_df
    .groupby("tratamiento", observed=True)

    # Calcular estadísticas descriptivas de las pendientes por tratamiento.
    .agg(
        # Promedio de las pendientes: tasa media de crecimiento en mm/día.
        pendiente_media=("pendiente_mm_dia", "mean"),

        # Desviación estándar: variabilidad de las pendientes entre semillas.
        sd=("pendiente_mm_dia", "std"),

        # Número de pendientes disponibles para ese tratamiento.
        n=("pendiente_mm_dia", "size")
    )

    # Convertir "tratamiento", que quedó como índice, en una columna normal.
    .reset_index()

    # Ordenar desde el tratamiento con mayor pendiente media
    # hasta el tratamiento con menor pendiente media.
    .sort_values("pendiente_media", ascending=False)

    # Crear un índice consecutivo: 0, 1, 2, ...
    .reset_index(drop=True)
)

# Mostrar la tabla de ranking.
ranking_pendientes
Out[32]:
tratamiento pendiente_media sd n
0 T07_100mgL 1.112461 0.073065 10
1 T08_150mgL 1.037232 0.072740 10
2 T06_75mgL 0.988514 0.063488 10
3 T05_50mgL 0.922457 0.049439 10
4 T10_PosCtrl 0.885811 0.055004 10
5 T04_25mgL 0.830171 0.053016 10
6 T03_10mgL 0.725811 0.060639 10
7 T02_5mgL 0.641132 0.060084 10
8 T01_Control 0.549193 0.050869 10
9 T09_200mgL 0.426339 0.058352 10
In [33]:
# Aplicar la prueba post hoc de Tukey HSD para comparar, por pares,
# las pendientes medias de crecimiento entre todos los tratamientos.
tukey = pairwise_tukeyhsd(

    # Variable numérica que se comparará:
    # contiene una pendiente estimada (mm/día) por cada semilla y tratamiento.
    endog=slopes_df["pendiente_mm_dia"],

    # Variable categórica que identifica a qué tratamiento pertenece
    # cada pendiente de la columna anterior.
    groups=slopes_df["tratamiento"],

    # Nivel de significancia familiar para todas las comparaciones simultáneas.
    # Un valor de 0.05 corresponde a un 95% de confianza.
    alpha=0.05
)
In [34]:
# Mostrar la tabla de resultados de Tukey HSD:
# diferencias de medias, intervalos de confianza, valores p ajustados
# y decisión para cada par de tratamientos.
print(tukey.summary())

     Multiple Comparison of Means - Tukey HSD, FWER=0.05      
==============================================================
   group1      group2   meandiff p-adj   lower   upper  reject
--------------------------------------------------------------
T01_Control    T02_5mgL   0.0919 0.0308  0.0046  0.1793   True
T01_Control   T03_10mgL   0.1766    0.0  0.0893  0.2639   True
T01_Control   T04_25mgL    0.281    0.0  0.1937  0.3683   True
T01_Control   T05_50mgL   0.3733    0.0  0.2859  0.4606   True
T01_Control   T06_75mgL   0.4393    0.0   0.352  0.5266   True
T01_Control  T07_100mgL   0.5633    0.0  0.4759  0.6506   True
T01_Control  T08_150mgL    0.488    0.0  0.4007  0.5754   True
T01_Control  T09_200mgL  -0.1229 0.0006 -0.2102 -0.0355   True
T01_Control T10_PosCtrl   0.3366    0.0  0.2493  0.4239   True
   T02_5mgL   T03_10mgL   0.0847 0.0652 -0.0026   0.172  False
   T02_5mgL   T04_25mgL    0.189    0.0  0.1017  0.2764   True
   T02_5mgL   T05_50mgL   0.2813    0.0   0.194  0.3686   True
   T02_5mgL   T06_75mgL   0.3474    0.0  0.2601  0.4347   True
   T02_5mgL  T07_100mgL   0.4713    0.0   0.384  0.5587   True
   T02_5mgL  T08_150mgL   0.3961    0.0  0.3088  0.4834   True
   T02_5mgL  T09_200mgL  -0.2148    0.0 -0.3021 -0.1275   True
   T02_5mgL T10_PosCtrl   0.2447    0.0  0.1574   0.332   True
  T03_10mgL   T04_25mgL   0.1044 0.0073   0.017  0.1917   True
  T03_10mgL   T05_50mgL   0.1966    0.0  0.1093   0.284   True
  T03_10mgL   T06_75mgL   0.2627    0.0  0.1754    0.35   True
  T03_10mgL  T07_100mgL   0.3866    0.0  0.2993   0.474   True
  T03_10mgL  T08_150mgL   0.3114    0.0  0.2241  0.3987   True
  T03_10mgL  T09_200mgL  -0.2995    0.0 -0.3868 -0.2121   True
  T03_10mgL T10_PosCtrl     0.16    0.0  0.0727  0.2473   True
  T04_25mgL   T05_50mgL   0.0923 0.0296   0.005  0.1796   True
  T04_25mgL   T06_75mgL   0.1583    0.0   0.071  0.2457   True
  T04_25mgL  T07_100mgL   0.2823    0.0   0.195  0.3696   True
  T04_25mgL  T08_150mgL   0.2071    0.0  0.1197  0.2944   True
  T04_25mgL  T09_200mgL  -0.4038    0.0 -0.4912 -0.3165   True
  T04_25mgL T10_PosCtrl   0.0556  0.555 -0.0317   0.143  False
  T05_50mgL   T06_75mgL   0.0661 0.3075 -0.0213  0.1534  False
  T05_50mgL  T07_100mgL     0.19    0.0  0.1027  0.2773   True
  T05_50mgL  T08_150mgL   0.1148 0.0019  0.0275  0.2021   True
  T05_50mgL  T09_200mgL  -0.4961    0.0 -0.5834 -0.4088   True
  T05_50mgL T10_PosCtrl  -0.0366 0.9355  -0.124  0.0507  False
  T06_75mgL  T07_100mgL   0.1239 0.0005  0.0366  0.2113   True
  T06_75mgL  T08_150mgL   0.0487 0.7272 -0.0386   0.136  False
  T06_75mgL  T09_200mgL  -0.5622    0.0 -0.6495 -0.4749   True
  T06_75mgL T10_PosCtrl  -0.1027 0.0089   -0.19 -0.0154   True
 T07_100mgL  T08_150mgL  -0.0752  0.154 -0.1626  0.0121  False
 T07_100mgL  T09_200mgL  -0.6861    0.0 -0.7734 -0.5988   True
 T07_100mgL T10_PosCtrl  -0.2266    0.0  -0.314 -0.1393   True
 T08_150mgL  T09_200mgL  -0.6109    0.0 -0.6982 -0.5236   True
 T08_150mgL T10_PosCtrl  -0.1514    0.0 -0.2387 -0.0641   True
 T09_200mgL T10_PosCtrl   0.4595    0.0  0.3721  0.5468   True
--------------------------------------------------------------

Ejemplos de lectura de Tukey¶

  • T01_Control vs. T07_100mgL: la pendiente media de T07_100mgL es mayor por $0.5633$ mm/día. El intervalo simultáneo no incluye cero y reject=True; se detecta una diferencia bajo el ajuste de Tukey.

  • T07_100mgL vs. T08_150mgL: la diferencia de pendientes medias es $-0.0752$ mm/día. El intervalo simultáneo incluye cero y reject=False; no hay evidencia suficiente de una diferencia entre estos tratamientos.

  • T04_25mgL vs. T10_PosCtrl: la diferencia de pendientes medias es $0.0556$ mm/día. El intervalo simultáneo incluye cero y reject=False; no hay evidencia suficiente de que sus tasas medias de crecimiento difieran.

F. Alcances y limitaciones del caso¶

  • Diez grupos pueden ser pocos para caracterizar con gran precisión una distribución de efectos aleatorios.
  • Los valores $p$ se basan en aproximaciones asintóticas.
  • Las bandas exploratorias son $\pm1$ EE, no intervalos de confianza al 95 %.
  • El ranking no incluye barras de incertidumbre.
  • Tukey se aplica a pendientes estimadas en un procedimiento complementario de dos etapas.
  • Revisar gráficos de diagnóstico del modelo

Referencias y materiales de apoyo¶

  • Pinheiro, J. C. y Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.
  • Wood, S. N. (2017). Generalized Additive Models: An Introduction with R. CRC Press.
  • Tukey, John W. “Comparing Individual Means in the Analysis of Variance.” Biometrics, vol. 5, no. 2, 1949, pp. 99-114. https://www.jstor.org/stable/3001913