Transcript Modelos lineales generales y mixtos El modelo lineal mixto

Modelos lineales
generales y mixtos
El modelo lineal mixto
y  X   Zb  
 
 b1 
 
b 
 1
y  X  2   Z  2   ε
 
 
 
 
br 

 t 
Varianza de (Y)
Var (ε)  R
Var(b)  G   I
2
b r
Cov(b,ε)  0
Var (y)  V  ZGZ ' R
Estimación
1
n
n
ML : l (G, R)   log  V   log  r 'V 1r   1  log(2 / n) 
2
2
2


1
n
REML : lR (G, R)   log  V   log X 'V 1 X
2
2
n p
n p

log  r 'V 1r  
1  log(2 / (n  p))
2
2
r  y  X  X 'V X  X ´V 1 y
1

BLUEs-BLUPs

βˆ  X 'Vˆ 1 X



X ´Vˆ 1y ( EBLUE  empirical BLUE )
ˆ 'Vˆ 1 y  Xβˆ
bˆ  GZ

( EBLUP  empirical BLUP)
1
1

X
´
R
X
X
´
R
Z 
ˆ
ˆ
Covar (β  β, b  b)  
1
1
1 
Z
´
R
X
´
Z
´
R
Z

G




Var (bˆ )   Z ' Rˆ Z  Gˆ X  Cˆ
ˆ ´Vˆ X   X 'Vˆ X 
Cˆ    GZ
Var (βˆ )  X 'Vˆ 1 X
1
1
21


1
1
21

X ´Vˆ 1ZGˆ

Pruebas F en modelos mixtos
H0 : H β  0


βˆ  H  H X 'Vˆ 1 X
 
F
r(H )
v=N-Rango(X|Z)
(usualmente)


H
 H βˆ 

Fr ( H ),v
Theory and Computational Methods for LME Models
2.4 Hypothesis Tests and Confidence Intervals pag. 91
Análisis de diseños
experimentales clásicos
mediante modelos mixtos
El modelo clásico
Efectos de
tratamiento
Efectos de
bloque
Yij=  + i + j + ij
i=1,...,a; j=1,...,b
a:número de tratamientos,
b: número de bloques
DBCA
yij     i  bj   ij

Suposiciones sobre los errores
Var ( )   e2 Itr

i  1,..., t , j  1,..., r
¿Son los bi efectos fijos o aleatorios?


Fijos: inferencia restringida
Aleatorios: extensión de la inferencia
DBCA
ZGZ’+ R
DBCA
1t 0
0 1
t

Z


0 0
1t 0
0 1
t
ZGZ    b2 


0 0
0
0 


1t 
0  1t 0
0   0 1t


1t   0 0
G   b2 I r
 0
0 
11
t t
 0 11
0 
t t
  b2 




1t 
0
0
0
0 



11
t t
DBCA
ZGZ’ R
 Jt
0
Var (y )   b2 
0

0
0
Jt
0
0
0
0
Vb 0
0 V
b
Var (y )  
0

0 0
0
0
0
1 0
0 1
0 
  e2 

0


Jt 
0 0
0
0 


1
rt
0
0
0
0  Vb   b2 Jt   e2 It
0 
   

 
 V  



Vb 


2
b
2
b
2
e
2
b
2
b
b
0
2
b
2
b
2
e
 b2
 b2





 b2   e2 
Varianzas

Bloques fijos
var( yi ) 

 e2
r
2 e2
Var ( yi  y j  ) 
r
Bloques aleatorios
var( yi ) 
 e2   b2
r
Costo de la inferencia extendida
2( e2   b2 )
Var ( yi  y j ) 
r
Ejemplo
Efecto de la fertilización sobre la
resistencia de la fibra de algodón
Modelando la
heteroscedasticidad
Modelo lineal general
yij     i  ij
Ejemplo



Rendimientos de Maíz de 8 híbridos comerciales
Diseño completamente aleatorizado con tres repeticiones
Datos
D1
1496.7
1585.0
1583.0
D2
1540.0
1446.3
1490.0
N1
1790.0
1884.9
1712.0
N2
1581.4
1593.8
1633.5
M1
1617.5
1708.6
1767.1
M2
1688.7
1620.0
1590.0
S1
1304.0
1690.0
1463.5
S2
1350.0
1740.0
1698.0
¿Heteroscedasticidad?
yij     i  ij
1

2 
var(y )  



yij    i   ij g (ij , zij , δsij )
 12

 22

var(y ) 







2
 n 
1





1
¿Heteroscedasticidad?
yij    i   ij g (ij , zij , δsij )


var( y )  




CRITERIOS DE AGRUPAMIENTO







¿Heteroscedasticidad?
yij     i  ij
1

2 
var(y )  



modelo1=gls(Rinde14~Hib-1,data=datos)
1





1
modelo2=gls(Rinde14~Hib-1,data=datos,weight=varIdent(form=~1|Hib))
anova(modelo1,modelo2)
 12 I

 22 I

var(y ) 







2 
 s I 
varFunc()
var(i )   2 g 2 (i , zi , δ)
g (): función devarianza






varFixed(form)
varIdent(value, form,fixed)
varPower(value=0,form=~fitted,fixed)
varExp(value=0,form=~fitted,fixed)
varConstPower(const=1,power=0,form=~fitted,fixed)
varComb(lista de funciones de varianza)
varFixed(form)
var( i )   zi
2
g ( zi )  zi
Ejemplo: varFixed(~edad)
varIdent(value,form,fixed)
var( i )    s
2
2
i
g (δ)   s
i
Ejemplo: varIdent(form=~1|Estratos)
value=c(1,..,1)
fixed=NULL
varPower(value,form,fixed)
var( i )   zi
2
g ( zi ,  si )  zi
2 si
 si
Ejemplo: varPower(form~fitted(.)|Estratos)
value=value=c(0,..,0)
fixed=NULL
Ejemplo con datos de
diálisis
varExp(value,form,fixed)

var( i )   2 exp 2 si zi

g ( zi ,  si )  exp  si zi


Ejemplo: varExp(form~edad|Estratos)
value=c(0,..,0)
fixed=NULL
varConstPower(const,power,form,fixed)
var( i )  
2

1, si
 zi
g ( zi ,  si )  1, si  zi
 2,si

2
 2,si
Ejemplo: varConstPower(form~fitted(.)|Estratos)
const=1,power=0
fixed=list(const,power)
varComb(value,form,fixed)


var( i )   212si exp 2 2 si zi ....
g 2 ()  g12 () g 22 ()...g k2 ()
Ejemplo:
varComb(varIdent(form~1|Estrato1),varExp(form~edad|Estrato2))
Modelando la
correlación

cor ( ij ,  ij ' )  h jj ' d  pij , pij '  , ρ

Ejemplo de DBCA analizado
utilizando la correlación
espacial
Para ejemplificar la modelación de la
correlación espacial, utilizaremos el
archivo Variación Espacial
 El encabezamiento de este archivo se
presenta a continuación:

Ejemplo correlación espacial

El diseño es un diseño en bloque para
evaluar comparativamente el rendimiento
de un número grande de variedades de un
cultivo….

La inclusión de los bloques como efecto
aleatorio para dar cuenta de la variabilidad
entre parcelas no parece ayudar a distinguir
entre variedades.
Modelos lineales generales y mixtos
modeloRend_0_REML<-lme(Rend~1+Variedad
,random=list(Bloque=pdIdent(~1))
,method="REML"
,na.action=na.omit
,data=R.data)
Variable dependiente:Rend
Medidas resumen
AIC
516.77
BIC
550.46
logLik
-240.39
Sigma
28.73
Tabla de ANOVA (suma de cuadrados marginal)
(Intercept)
Variedad
numDF
1
15
denDF
45
45
F-value
16.75
0.38
p-value
0.0002
0.9783
Parámetros de los efectos aleatorios
Modelo de covarianzas de los efectos aleatorios: pdIdent
Formula: ~1|Bloque
Desvíos estándares relativos al residual y correlaciones
(Intercept)
(Intercept)
7.8E-05
Levelplot de los residuos en
función de la posición de la parcela
8
1.5
1.0
6
PosY
0.5
0.0
4
-0.5
-1.0
2
-1.5
2
4
6
PosX
8
Modelos de correlación

No serial



Simetría compuesta
(corCompSymm)
Serial

Desestructurada
(corSymm)

Autorregresivo de orden 1
(corAR1)

AR continuo de orden 1
(corCAR1)

Autorregresivo-Medias móviles
(corARMA)
Espacial (isotrópica)

Exponencial
(corExp)

Gaussiano
(corGaus)

Lineal
(corLin)

Rational quadratic
(corRatio)

Esférica
(corSpher)
Simetría compuesta
cor(ij , ij ' )  hjj ' ρ   (correlación intra clase); ρ  {}


La correlación entre cualquier par de errores
dentro de un estrato es la misma.
Ejemplo especificado en R (nlme)

correlation=corCompSymm(form = ~ 1 | Estrato )

Form=~1 implica que R va a calcular una correlación
0.3

común para
todas las
observaciones
dentro del estrato
1
0.3
0.3

Si hay 4 observaciones
por estrato yla correlación
0.3
1 correlación
0.3 resultante es:
intraclasees
0.3 la0.3
matriz de

1
0.3 0.3 0.3

0.3 0.3 0.3

1
Modelos de correlación
Serial
Desestructurada
(corSymm)
Autorregresivo
(corAR1)
AR
de orden 1
continuo de orden 1
Autorregresivo-Medias
móviles
(corCAR1)
(corARMA)
Desestructurada
cor(ij , ij ' )  hjj ' ρ   jj ' ; ρ  { jj '}
En este modelo todas las correlaciones
pueden ser diferentes
 Es el modelo más complejo
 Requiere la estimación de un número
grande de parámetros

Desestructurada

cor(ij , ij ' )  hjj ' ρ   jj ' ; ρ  { jj '}
Ejemplo especificado en R (nlme)

correlation=corSymm(form=~1|Estrato).

form=~1 implica que R toma la posición de las
observaciones en el archivo para calcular las
correlaciones
Si tengo 4 observaciones en el grupo i-ésimo, y sus
correlaciones son:


12=0.2, 13= 0.1, 14=-0.1, 23=0,24=0.2, 34= 0
0.2 0.1 0.1
 1
 0.2

1
0.0
0.2


 0.1 0.0 1
0.0 


1 
 0.1 0.2 0.0
Autoregresivo de orden 1
cor ( ij ,  ij ' )  h jj '  k jj ' ρ   



k jj '
; k jj '  0,1,....; ρ  { }
kjj’ representa la separación entre las
observaciones en una secuencia,
usualmente, temporal
 es la correlación entre cualquier par de
errores cuyas observaciones están
“separadas” una unidad (k=1).
La correlación entre errores decrece con k
Autoregresivo de orden 1

Ejemplo especificado en R (nlme)

correlation=corAR1(form = ~ 1 | Subject )

Form=~1 implica que R toma la posición de las
observaciones en el archivo para calcular la
separación k. De otra forma hay que especificar una
covariable que indica el orden de observación.

Si =0.8 y tengo 4 observaciones en el grupo o
estrato i-ésimo, su matriz de correlación sería:
0.8 0.64 0.512 
 1
 0.8

1
0.8
0.64


 0.64 0.8
1
0.8 


0.512
0.64
0.8
1


Auto-regresivo continuo
cor ( ,  )  h  s , ρ    ; s  0; ρ  {  0}
s jj '
ij
ij '
jj '
jj '
jj '

En este caso sjj’ representa la “separación” entre
observaciones en una escala, usualmente,
temporal.

sjj’ puede variar en una escala continua, a
diferencia con el modelo autoregresivo “clásico”
que lo hace en una escala entera.
Modelos de
correlación espacial
Modelos de correlación
espacial (isotrópicos)

Exponencial
(corExp)

Gaussiano
(corGaus)

Lineal
(corLin)

Rational quadratic
(corRatio)

Esférica
(corSpher)
Semivariograma

 
1
2

 s  p ,  p , ρ  Var  p   p
i
j
si var( pi )  1 i

 

i
j

 
Varianza de la diferencia/2
 s  p , p ,ρ  1 h s  p , p , ρ
i
j
i
j
Still
h(s, ρ) es la función
de correlación, ρ es el
parámetro de
correlación y “s” la
distancia
h(0, ρ) = 1
γ(0, ρ) = 0.
Nugget
1  nugget  h  s, ρ  , s  0
hnugget  s, ρ   
1, caso contrario
Distancia entre parcelas (s)
Semivariogramas de la
Correlación espacial
corExp
corExp(form = ~ x + y, metric=“man”, nugget = T )
si range=1 y nugget=0.2 la matriz de covarianzas podría
ser la siguiente
corGaus
corRatio
corLin
corSpher
Modelos alternativos para la
covarianza

Los modelos no son directamente
comparables mediante la prueba del
cociente de verosimilitudes

Los criterios de Akaike o de Schwarz
se pueden usar para selección de
modelo
Ejemplo de DBCA analizado
utilizando la correlación
espacial
Para ejemplificar la modelación de la
correlación espacial, utilizaremos el
archivo Variación Espacial
 El encabezamiento de este archivo se
presenta a continuación:

Ejemplo correlación espacial

El diseño es un diseño en bloque para
evaluar comparativamente el rendimiento
de un número grande de variedades de un
cultivo….

La inclusión de los bloques como efecto
aleatorio para dar cuenta de la variabilidad
entre parcelas no parece ayudar a distinguir
entre variedades.
Modelos lineales generales y mixtos
modeloRend_0_REML<-lme(Rend~1+Variedad
,random=list(Bloque=pdIdent(~1))
,method="REML"
,na.action=na.omit
,data=R.data)
Variable dependiente:Rend
Medidas resumen
AIC
516.77
BIC
550.46
logLik
-240.39
Sigma
28.73
Tabla de ANOVA (suma de cuadrados marginal)
(Intercept)
Variedad
numDF
1
15
denDF
45
45
F-value
16.75
0.38
p-value
0.0002
0.9783
Parámetros de los efectos aleatorios
Modelo de covarianzas de los efectos aleatorios: pdIdent
Formula: ~1|Bloque
Desvíos estándares relativos al residual y correlaciones
(Intercept)
(Intercept)
7.8E-05
Levelplot de los residuos en
función de la posición de la parcela
8
1.5
1.0
6
PosY
0.5
0.0
4
-0.5
-1.0
2
-1.5
2
4
6
PosX
8
La especificación en R

Modelos lineales generales y mixtos

modeloRend_9_REML<-gls(Rend~1+Variedad
,correlation=corGaus(form=~PosX+PosY
,metric="euclidean"
,nugget=FALSE)
,method="REML"
,na.action=na.omit
,data=R.data)






Variable dependiente:Rend
Medidas resumen
AIC
436.16
BIC
469.84
logLik
-200.08
Sigma
21.26
Tabla de ANOVA (suma de cuadrados marginal (tipo III))
numDF
(Intercept)
Variedad
F-value
1
15
p-value
63.25
17.90
<0.0001
<0.0001
Estructura de correlación
Modelo de correlacion: Gaussian spatial correlation
Formula: ~ PosX + PosY
Metrica: euclidean
Parámetros del modelo
range
Estimación
1.55
0.6
0.4
0.2
0.0
exp(-(x/1.55)^2)
0.8
1.0
La correlación entre parcelas
como función de la distancia es
0
1
2
x
3
4
Medidas ajustadas y errores estándares para Variedad
Variedad Medias E.E.
F
73.98
7.34
A
G
71.09
6.63
A
B
J
69.97
6.63
A
B
N
68.37
6.93
A
B
C
M
68.10
6.91
A
B
C
L
67.39
7.39
A
B
C
O
64.39
6.64
A
B
C
C
64.33
7.05
A
B
C
E
64.04
6.91
B
C
H
61.68
7.16
B
C
I
59.73
6.93
C
D
57.67
6.80
C
P
57.59
6.95
C
K
57.38
6.63
A
52.05
6.54
B
51.82
6.67
Letras distintas indican diferencias significativas(p<= 0.05
D
D
D
D
D
D
D
D
D
E
E
E
E
E
E
E

El mejor ajuste obtenido por la modelación
de la correlación espacial podría haberse
obtenido igualmente incluyendo como
efectos fijos la posición x e y de la parcela
actuando como regresoras.

No siempre, sin embargo, es posible
corregir tan fácilmente las respuesta
utilizando un modelo de efectos fijos.
Otro ejemplo de correlación
espacial

Los datos en el archivo Otra variación
espacial presenta un caso similar al
anterior pero donde el patrón de
variación espacial es menos regular.
Aquí la modelación de la correlación
es una clara ventaja.