Es una expresión matemática que indica cómo una variable de respuesta aleatoria, con una distribución de probabilidad dada, se puede relacionar con una o más variables predictoras aleatorias o no, consideradas en el diseño experimental.

Es posible expresar el modelo estadístico de la forma:

\[Y_i = f(x_{1i},...,x_{pi}) + \varepsilon_i \] donde \(f(x_{1i},...,x_{pi})\) es una función matemática de depende de \(p\) variables predictoras y \(\epsilon_i\) representa una variable aleatoria.

• Modelos lineales en los parámetros: La variable de respuesta es descripta a través de una combinación lineal de variables predictoras como por ejemplo:

• Modelos no lineales:

En el análisis de regresión las variables predictoras pueden ser cuantitativas o cualitativas(dummy en este último caso). El objetivo es ajustar funciones y/o predecir la variable de respuesta

En el análisis de la varianza las variables predictoras se denominan factores y se las trata como cualitaticas. El objetivo se centra en la comparación de medias.

Supongamos que tenemos \(n\) observaciones independientes. Para cada individuo \(i\) se observa:

1. Una variable de respuesta \(Y_i\).

2. \(p\) covariables \(x_i =(x_{i1},...,x_{ip})^t\). Donde \(x_i\) es un vector columna de dimensión \(p<n\).

suponiendo que \(\varepsilon_i \overset{i.i.d.}{\sim} N(0, \sigma^2)\)

donde \(\beta\) es un vector columna con \(p\) parámetros y i.i.d en \(\varepsilon_i\) significa independientes e identicamente distribuidos

\(\varepsilon_i \overset{i.i.d.}{\sim} N(0, \sigma^2)\)

Implica:

• Varianza Constante

• Observaciones Independientes

• Distribución Normal

Cuando los supuestos no se cumplen, las estimaciones de los parámetros siguen siendo insesgadas, pero los errores estándar de ellos no, por los que los valores p y los Intervalos de Confianza no son confiables.

Los supuestos anteriores implican que:

1. \(E(Y_i|\mathbf{x}_i)=\mathbf{x}_i^t\mathbf{\beta}\)

2. \(Var(Y_i|\mathbf{x}_i)=\sigma^2\)

3. \(Y_i|\mathbf{x}_i\overset{indep.}{\sim}N(\mathbf{x}_i^t\mathbf{\beta}, \sigma^2)\)

Esto implica que para dos individuos \(i\) y \(j\), con \(i\neq j\) se tiene que \(Cov(Y_i,y_j)=0\) ya que asumimos errores independientes.

Es factible plantear el modelo anterior de forma matricial. Para ello se define la siguiente matriz del modelo de dimensión \(nxp\):

\[\mathbf{X}_{n\times p} = \begin{pmatrix} x_{1,1} & x_{1,2} & \cdots & x_{1,p} \\ x_{2,1} & x_{2,2} & \cdots & x_{2,p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n,1} & x_{n,2} & \cdots & x_{n,p} \end{pmatrix} \] Donde cada fila en \(\mathbf{X}\) representa al vector de covariables de cada uno de los \(n\) individuos.

Es decir: \[\mathbf{X} = \begin{pmatrix} \mathbf{x}_{1}^t \\ \mathbf{x}_{2}^t\\ \vdots \\ \mathbf{x}_{n}^t \end{pmatrix}\]

Por lo que el modelo (2.1) se puede escribir como: \[ \begin{equation} \mathbf{Y} = \mathbf{X}\mathbf{\beta} + \mathbf{\varepsilon} \tag{2.2} \end{equation} \]

donde:

• \(\mathbf{Y}=(Y_1,\ldots,Y_n)^t\)

• \(\mathbf{\varepsilon}=(\varepsilon_1,\ldots, \varepsilon_n)^t\)

Se asume tambien que: \(\mathbf{\varepsilon}\sim N(\mathbf{0}, \sigma^2 \mathbf{I}_n)\) donde \(\mathbf{I}_n\) es la matriz identidad de orden \(n\).

\[ \begin{equation} \mathbf{Y} = \mathbf{X}\mathbf{\beta} + \mathbf{\varepsilon} \tag{2.2} \end{equation} \] Los supuestos para el modelo (2.2) implican que:

1. \(E(\mathbf{Y}|\mathbf{X})=\mathbf{X}\mathbf{\beta}\)
2. \(Var(\mathbf{Y}|\mathbf{X})=\sigma^2\mathbf{I}_n\)
3. \(\mathbf{Y}|\mathbf{X}\sim N(\mathbf{X}\mathbf{\beta}, \sigma^2\mathbf{I}_n)\)

La función objetivo a ser minimizada es: \[\sum(y_i-\hat{y_i})^2=\sum e_i^2=\sum(Y-X^t \beta)^2 \] el estimador mínimo cuadrático obtenido se define como: \[\hat{\beta}=(X^tX)^{-1}X^tY\]

Si el diseño es desbalanceado las estimaciónes dejan de ser insesgadas y las sumas de cuadrados pierden ortogonalidad.

1. Modelos lineales generales con estimación por Mínimos Cuadrados: Pueden modelar la falta de independencia, pero no corresponde su aplicación ante falta de normalidad, heterocedasticidad o desbalanceo importante.

2. Modelos lineales generales con estimación por Mínimos Cuadrados generalizados: Posibilitan modelar heterocedasticidad, pero no corresponde su aplicación ante falta de normalidad o desbalanceo importante.

3. Modelos lineales generales con estimación por máxima verosimilitud: Posibilitan modelar heterocedasticidad y soportan desbalanceo. No corresponde su aplicación ante falta de normalidad.

4. Modelos lineales generalizados con estimación por máxima verosimilitud: Posibilitan modelar heterocedasticidad, soportan desbalanceo y consideran varias distribuciones.

• Es uno de los supuestos más importantes en los modelos lineales generales
• Consiste en suponer que todos los tratamientos tienen la misma variabilidad (\(\sigma^2\)) o alternativamente, que los errores tienen una variabilidad constante \(\sigma^2\). \[Var(\varepsilon_i)=\sigma^2 \]
• La violación al supuesto de igualdad de varianzas provoca:
  1. Estimaciones erróneas de los Errores Estandar de los estimadores
  2. P-valores e Intervalos de Confianza que no son válidos
• La solución clásica: transformar. Ello implica pérdida de información y trabajar con datos que no son los originales.

1. Biológicas
2. Presencia de outliers
3. Proceso generador alejado de la distribución normal:

• Distribución Poisson
• Distribución Binomial
• Otras distribuciones

• Si no existe información a priori sobre la naturaleza de la heterocedasticidad, es posible posible estimar el modelo sin considerarla para luego hacer un análisis de los residuos que se generan. -Gráficos de residuos vs esperados o predichos, esperando encontrar una distribución al azar y con variabilidad constante

• Pruebas Analíticas: Test de Levene

• Los residuos son fundamentales en el proceso de validación de modelos
• El Residuo (\(e_{ij}\)) es la diferencia entre el valor observado en \(y\) y el valor esperado según el modelo: \[e_{ij}=y_{ij}-\hat{y_{ij}} \]
• Cuando se modelan estructuras de varianza es necesario tener una medida relativa que considere a las mismas. Por lo que suele utilizarse el residuo estandarizado: \[e_{ij}=\frac{y_{ij}-\hat{y_{ij}}}{\sqrt{\sigma^2_{ij}}} \]

• Por Estimación Mínimo Cuadrática Ordinaria se puede incorporar una función de varianza de la forma: \[Var(\varepsilon_i)=\sigma^2 f(\mu_i,X,\gamma)\] Donde:
• \(\mu_i\): esperanza de la variable de respuesta
• \(X\): una covariable para la varianza
• \(\gamma\): parámetro estimado a partír de la función de varianza propuesta.

El estimador por mínimos cuadrados ponderados se puede obtener a partír de: \[\hat{\beta}=(X^t\hat{V}^{-1}X)^{-1}X^t \hat{V}^{-1} Y\] donde \(\hat{V}\) es una matriz de varianza-covarianza obtenida a partír de la función de varianza propuesta.

Eisenhart (1947) introduce modelos de Análisis de la Varianza para efectos fijos y aleatorios (Modelos I y II). Los criterios para definir si un efecto es fijo o aleatorio son:

1. Si al repetir el experimento los mismos niveles se van a estudiar nuevamente, el efecto es fijo (mecanismo de muestreo aleatorio).
2. Si las inferencias estarán limitadas a los niveles de los factores realmente empleados en el experimento, el factor es fijo. Si las conclusiones se aplicarán a poblaciones más generales, es aleatorio (espacio de inferencia).
3. Si los niveles estudiados forman una estructura jerárquica (multinivel), como alumnos en un colegio, lagos en una región, etc, entons los efectos son aleatorios.
4. Otras consideraciones: contexto, espacio de inferencia (sujeto-específico o promedio-poblacional), estrategia de modelación, etc.

Es posible plantear un modelo que contenga efectos fijos, efectos aleatorios y una función de varianza para los residuales de la forma:

\[Y=fijos+aleatorios+\varepsilon\] donde los efectos aleatorios \(\sim N(0,G )\) y \(\varepsilon \sim(N(0,R)\).

Los modelos Mixtos permiten un ambiente flexible para modelar correlaciones y varianzas heterogéneas modelando así:

• Modelos heterocedásticos.
• Modelos jerárquicos con efectos aleatorios
• Modelos de regresión específico de sujetos
• Modelos de covarianza y correlación (temporal, espacial)

La estimación se basa en verosimilitud, por lo que las pruebas son válidas en diseños complejos con datos desbalanceados.

Es posible, además, realizar estimaciones máximo verosimiles restringidas para poder estimar conponentes de varianza y realizar inferencia sore los efectos aleatorios a través de predictores encogidos hacia la media.

Primer Formulación de un Modelo Mixto

\[\mathbf{y} = \mathbf{X\beta+Zu+e},\] \[\begin{bmatrix} \mathbf{u}\\ \mathbf{e} \end{bmatrix}\sim N\left(\begin{bmatrix} \mathbf{0}\\ \mathbf{0} \end{bmatrix}, \begin{bmatrix} \mathbf{G} & \mathbf{0}\\ \mathbf{0} & \mathbf{R} \end{bmatrix}\right)\]

Segunda Formulación de un Modelo Mixto

\[\mathbf{y}|\mathbf{u} \overset{i.d.}{\sim} N(\mathbf{X\beta+Zu}, \mathbf{R})\] \[\mathbf{u} \overset{i.d.}{\sim} N(\mathbf{0}, \mathbf{G}).\]

\[E(\mathbf{y})=\mathbf{X\beta}, \quad Var(\mathbf{y})=\mathbf{V}=\mathbf{ZGZ^t+R} \quad \text{(dist. marginal)}.\]

\[E(\mathbf{y}|\mathbf{u})=\mathbf{X\beta+Zu}, \quad Var(\mathbf{y}|\mathbf{u})=\mathbf{R} \quad \text{(dist. condicional)}.\]

Si escribimos la matriz de varianzas y correlación serial \(R\) como: \[\mathbf{R} = \begin{bmatrix} \mathbf{R}_1 & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{R}_2 & \cdots & \mathbf{0} \\ \cdots & \cdots & \ddots & \cdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{R}_n \end{bmatrix}\] con: \[\mathbf{R}_i=Diag^{1/2}(\sigma_{ij}) \mathbf{Corr}_i Diag^{1/2}(\sigma_{ij}).\] formulamos el modelo mixto en R con la función lme de la libreria nlme:

library(nlme)
Modelo <-lme(Y~1+factor(Trat)+X,
             random=list(U=pdIdent(~1)), # efectos aleatorios
             weight=varComb(varIdent(form=~1|Trat), # varianza residual
             correlation=corCompSymm(form=~1|Trat), # estructura de correlación
method='REML', data=Datos)

• Sin Estructura(pdSymm):Varianzas distintas para cada componente aleatorio y correlaciones tambien distintas entre ellos

\[G = \begin{pmatrix} \sigma^2_{b1} & \sigma_{b12} & \sigma_{b13} \\ \sigma_{b21} & \sigma^2_{b2} & \sigma_{b23} \\ \sigma_{b31} & \sigma_{b32} & \sigma^2_{b3} \end{pmatrix} \]

• Diagonal(pdDiag): Varianzas distintas para cada componente aleatorio, sin correlación entre ellos \[G = \begin{pmatrix} \sigma^2_{b1} & 0 & 0 \\ 0 & \sigma^2_{b2} & 0 \\ 0 & 0 & \sigma^2_{b3} \end{pmatrix} \]

• Identidad(pdIdent): Varianza única para todos los componentes aleatorios, sin correlación entre ellos

\[G = \begin{pmatrix} \sigma^2_{b} & 0 & 0 \\ 0 & \sigma^2_{b} & 0 \\ 0 & 0 & \sigma^2_{b} \end{pmatrix} \]

• Componente de Simetría(pdCompSymm): Varianza única para todos los componentes aleatorios y única correlación entre ellos.

\[G = \begin{pmatrix} \sigma^2_{b1} & \sigma_{b2} & \sigma_{b2} \\ \sigma_{b2} & \sigma^2_{b1} & \sigma_{b2} \\ \sigma_{b2} & \sigma_{b2} & \sigma^2_{b1} \end{pmatrix} \]

-Fija(varFixed): La varianza como función lineal de alguna covariable \[Var(\varepsilon_i) \sim N(0,\sigma^2 X_i)\] -Identidad(varIdent): Una varianza distinta para cada grupo \[Var(\varepsilon_i) \sim N(0,\sigma^2_i)\]

-Exponencial(varExp): La varianza como función exponencial de alguna covariable \[Var(\varepsilon_i) \sim N(0,e^{2 \gamma \sigma X_i } )\] -Potencial(varPower): La varianza como función de potencia de alguna covariable \[Var(\varepsilon_i) \sim N(0,\sigma^2 \lvert{X_i}\rvert^{2 \gamma})\]

• Requiere calculos Complejos
• Para la estimación es fundamental asumir una distribución de probabilidades de la variable
• Si se cumplen los supuestos de independencia, normalidad y homocedasticidad y el diseño es balanceado, las estimaciones coinciden con las obtenidas por mínimos cuadrados.

• Se arrojan 5 monedas equilibradas y se cuentra la cantidad de caras.
• La distribución de probabilidad en cada ensayo es binomial: \[P(X=x) \sim Bi(n,P) \] donde \(n\)=5 se corresponde a la cantidad de ensayos y \(P\)=0,5 es la probabilidad de obtener cara en cada tirada.

• Si en una experiencia dada se arrojan 5 monedas equilibradas 10 veces y se cuenta la cantidad de caras: 5 3 2 2 4 1 1 3 3 4

• Para cada tirada se puede calcular, utilizando la distribución binomial, la probabilidad de obtener ese resultado: 0.03125, 0.03125, 0.03125, 0.03125, 0.15625, 0.15625, 0.15625, 0.31250, 0.31250, 0.15625

• La función de distribución conjunta de los valores observados, asumiendo independencia, se calcula como el producto de las probabilidades individuales:

\[f(X| n,p,f)=P_{(5)}P_{(3)}...P_{(4)}=5.55x10^{-8}\] - La función de verosimilitud se obtiene de la misma forma que la distribución conjunta, pero mientras que en esta \(X\) es variable y los parámetros son fijos, en la función de verosimilitud \(X\) es fijo y los parámetros son variables:

\[L=f(p|X,f)=P_{(5)}P_{(3)}...P_{(4)}=5.55x10^{-8}\] -En general: \[L(\theta|datos,modelo)\]

• Si ahora en vez de presuponer que la moneda es equilibrada, estamos dispuestos a aceptar la posibilidad de que no lo sea….. es decir….

• Ante la misma experiencia realizada:5 3 2 2 4 1 1 3 3 4 , cuál es la probabilidad de los valores observados. Podríamos calcularlos asumiendo la distribución de probabilidad antes dada, pero…. con que valor de \(P\)?

• Como no lo sabemos, podemos obtener la función de verosimilitud a partír de ensayo y error obteniendo algo como:

  

• El valor mas próbable de obtener cada de la moneda a partír de los datos recabados es \(0.74\), este es el valor estimado por máxima verosimilitud.
• Al ser un producto de probabilidades, \(L(\theta|X )\) toma valores entre 0 y 1.
• Para facilitar los cálculos se aplica logarinto natural a los efectos de trabajar con sumas en vez de con productos.
• Para hallas el máximo se deriva la función con respecto a \(\theta\) y se la iguala a 0: \[\frac{dlogL(\theta|X)}{d(\theta_j)=0}\]
• \(logL(\theta|X)\) toma siempre valores negativos dado que \(0<L<1\)

• Elegiremos el valor de \(\theta\) para el cual \(L(\theta)\) nos de el máximo
• Si el tamaño de la muestra es grande, las estimaciones por ML proveen estimadores insesgados (la esperanza del estimador coincide con el parámetro) y consistentes (su varianza tiende a cero cuando \(n\) tiende a infinito).
• Sin embargo, las estimaciones de las varianzas están sesgadas, por lo lo que en la práctica se aplica una corrección a los grados de libertad o se obtienen las estimaciones a través de máxima verosimilitud restringida (REML).

Cuando los supuesto usuales del modelo clásico se cumplen la verosimilitud resulta maximizada por:

\[ log L=-(n/2) log 2 \pi - (n/2) log \sigma ^2 - (1/2) \sum_{j=1}^{n} (y_j - X^T_j \beta)^2/ \sigma^2 \]

• La estimación REML es un camino alternativo para estimar los parámetros de varianza-covarianza.
• REML es preferido a la estimación ML porque produce estimadores insesgados debido a que toma en cuenta la pérdida en grados de libertar que resulta de estimar los efectos fijos.
• Debido a que en la verosilimitud las estimaciones de los efectos fijos y los aleatorios están relacionadas, las estimaciones obtenidas por ML y REML son diferentes
• La estimación REML se utiliza a los efectos de presentar el modelo final
• La estimación ML se utiliza a los efectos de comparar modelos cuando se plantean cambios en la componente fija del modelo.

Criterios de Información

• Akaike: \(AIC= -2logL(\theta)+2p\)

• Bayesiano de Schwartz:\(BIC= -2logL(\theta)+p\;ln(n)\)

• Representa un resumen de la información de un modelo, teniendo en cuenta \(L(\theta)\)(cuanto mayor, mejor) y el número de parámetros a estimar del modelo ($p)(cuanto mayor, peor)

• El valor individual no es interpretable, solo sirve con fines comparativos: cuanto menor, mejor el modelo

• BIC penaliza más que AIC los modelos con más parámetros a estimar

• Algunos autores sugieren que AIC mide la capacidad predictiva mientras que BIC lo hace con la bondad de ajuste, por lo que, al proveer estimaciones de métricas diferentes no deberían considerarse como criterios en competencia

Prueba t F y de Wald para Efectos Fijos

La prueba t para evaluar el efecto de un único efecto aleatorio puede realizarse a partír del estadístico t: \[t=\frac{\hat{\beta}} {se(\hat{\beta})} \]

• Pero en el contexto de los Modelos Lineales Mixtos la distribución de este estadístico no tiene una forma exacta, el mismo planteo se presenta para las pruebas F. Para salvar ese inconveniente se proponen métodos aproximados que consideran los eventuales efectos aleatorios y residuales correlacionados.

• También es posible plantear la prueba de Wald: \[X^2=(\frac{\hat{\beta}} {se(\hat{\beta})})^2 \] la cual asintóticamente tiene una distribución \(\chi^2_{gl}\) con \(gl\) correspondiente a la varianza del error.

Prueba de Razón de Verosimilitudes (RV) para Efectos Fijos

Consiste en una prueba asintótica basada en el principio de verosimilitud:

\[G= -2 log(\frac{L(\theta_1)}{L(\theta_2)} =2logL(\theta_2)-2logL(\theta_1) \sim \chi^2 {(p_1-p_2)}\]

• Es un método para comparar modelos anidados por máxima verosimilitud

• Un modelo estadístico se dice anidado dentro de otro si este representa un caso especial del otro modelo

• \(L(\theta)_2\) es el logaritmo de la verosimilitud para el modelo mas general que tiene asociados \(p_2\) parámetros estimados

• \(L(\theta)_1\) es el logaritmo de la verosimilitud para el modelo con restricciones que tiene asociados \(p_1\) estimados

Prueba RV para Efectos Aleatorios

• En el caso de probar hipótesis sobre parámetros de varianzas y covarianzas es posible también usar la prueba de Razón de Verosimilitudes

• En este caso podemos hacerlo con ML o REML, ya que estamos comparando modelos con la misma parte fija.

• Debemos tener en cuenta que al probar si una o más varianzas son cero estamos en la frontera del espacio de parámetros, y por lo tanto, la distribución asintótica del estadístico de la prueba chi-cuadrado.

• (Molenberghs and Verbeke 2007) recomienda usar como distribución asintótica una mezcla de distribuciones chi-cuadrado. Por ejemplo, para probar \(\sigma^2_b=0\),el cociente de verosimilitud se prueba con una mezcla de \(\chi^2_0\) y \(\chi^2_1\).

A efectos de verificar el cumplimiento de supuestos se proponen diversos residuales

Condicionales

• Residual Condicional Crudo \[\hat{\varepsilon_i}= y_i - X_i \hat{\beta} - Z_i \hat{u_i} \]
• Residual Condicional Estandarizado de PEARSON \[\hat{\varepsilon_i}= \frac{ y_i - X_i \hat{\beta} - Z_i \hat{u_i}} {\hat{\sigma_i}} \]

Marginales

• Residual Marginal Crudo \[\hat{\varepsilon_i}= y_i - X_i \hat{\beta} \]

1- Comenzar con un modelo donde la parte fija contenga todas las variables explicatorias y sus interacciones, se asume normalidad, homocedasticidad e independencia de los errores. Utilizar REML como método de estimación.

2- Seleccionar la estructura aleatoria óptima. Para ello utilizar AIC o BIC y la prueba RV (esta última solo para el caso de modelos anidados)

3- Analizar gráficamente el supuesto de homocedasticidad. Si se detecta heterocedasticidad elegir una estructura de varianzas apropiada. Para ello utilizar gráficos, AIC o BIC o RV

4- Chequear la normalidad

5- Si no se logra homocedasticidad y/o normalidad en los residuales considerar modelos lineales generalizados o transformación de la variable

6- Estudiar el componente fijo del modelo. Para ello utilizar las pruebas F o la prueba RV, si se recurre a esta última utilizar estimación ML ya que las estimaciones REML no son válidad para comparar modelos con efectos fijos diferentes

7- Comenzar analizando las interacciones de mayor orden y sólo en caso de no resultar significativas, analizar las interacciones de menor orden y/o efectos principales

8- Efectuar las comparaciones adecuadas según la significación de los efectos fijos

9- Presentar el modelo final utilizado utilizando estimación REML

Consideraciones

• No comparar AIC con BIC

• No comparar los Criterios de Información obtenidos por REML con los obtenidos con ML

• Es posible comparar modelos con distinta estructura aleatorio pero la misma estructura fija mediante Criterios de Información, pero la prueba RV sólo se puede aplicar para modelos con estructura aleatoria anidada

• Es posible comparar modelos con distinta estructura fija pero la misma aleatoria mediante AIC y BIC, pero la estimación debe efectuarse por ML

El modelo mixto para este diseño puede plantearse como:

\[y_{ijk}=\mu+\tau_i+\delta_k+\mu+\tau_i*\delta_k+b_j+p_{ik}+\varepsilon_{ijk}\] donde:

• \(\mu\) representa la media general de la respuesta
• \(\tau_i\) representa el efecto del i-ésimo nivel asociado a las parcelas principales (Factor Agua)
• \(\delta_k\) representa al efecto dek k-ésimo nivel asociado a las subparcelas (Factor Variedad)
• \(b_j\),\(p_{ik}\) y \(\varepsilon_{ijk}\) corresponden a efectos aleatorios de los bloques, de las parcelas dentro de los bloques y de los errores experimentales.

El modelo mixto para este diseño puede plantearse como:

\[y_{ijk}=\mu+\tau_i+\delta_k+\mu+\tau_i*\delta_k+b_j+p_{ik}+\varepsilon_{ijk}\]

Los supuestos sobre los componentes aleatorios son:

• \(b_j\sim N(0,\sigma^2_b)\)

• \(p_jk\sim N(0,\sigma^2_p)\)

• \(\varepsilon_{ijk} \sim N(0,\sigma^2_{\varepsilon})\)

• Los tres componentes aleatorios son independientes

MODELO_1<-lme(Rendimiento~1+Agua+Variedad+Agua:Variedad
    ,random=list(Bloque=pdIdent(~1)
                 ,Agua=pdIdent(~1))
    ,method="REML"
    ,control=lmeControl(niterEM=150
                        ,msMaxIter=200)
    ,na.action=na.omit
    ,data=TRIGO
    ,keep.data=FALSE)


summary(MODELO_1)

plot(MODELO_1)

Codigo en R para poder reproducir el ejemplo

Littell (2006) propone los siguientes criterios al momento de proponer un modelo:

1. Visualizar la forma en la cual se obtuvieron las observariones
2. Enumerar los criterios de bloqueo, si los hay.
3. Enumerar los factores de tratamiento y todas las posibles interacciones.
4. Identificar la unidad experimental de cada cada efecto principal y interacción.
5. Determinar los efectos fijos del modelo del paso 3 y del paso 2.
6. Determine los efectos aleatorios del modelo del paso 4 y 2. En general, habrá un efecto aleatorio por tamaño de unidad experimental.
7. Si hay unidades de muestreo dentro de la unidad experimental más pequeña (por ejemplo, si la clase es la unidad experimental, pero las medidas se toman en estudiantes individuales dentro de la clase), entonces el error residual corresponde al error de la unidad de muestreo. De lo contrario, el error residual se corresponde al error de unidad experimental más pequeño.

Pinheiro J. and Bate D 2000 Mixed-Effects Models in S and S-PLUS. Springer

Di Rienzo, J. and Machiavelli, R and Casanoves, F. 2017 Modelos Lineales Mixtos-Aplicaciones en Infostat. Grupo Infostat. Cordoba. Argentina

Molenberghs, G., and G. Verbeke. 2007. “Likelihood Ratio, Score, and Wald Tests in a Constrained Parameter Space.” The American Statistician 61 (1): 567–79.

Littel,R and Milliken, G and Strup, W and Wolfinges, R. 2006 Sas System for Mixed Models. Sas Institute. Cary, NC. USA. Faraway, J 2017 Extending the Linear Model with R. Generalized Linear, Mixed Effects and Nonparametric Regression Models. CRC Press. New York

Zuur, A and Ieno, E and Walker N et al 2009 Mixed Effects Models and Extension in Ecology with R. Springer-Verlag