marzo 21, 2017

ANOVA combinado (2) — mediciones temporales repetidas


“El ANOVA combinado [de medidas repetidas] se lleva a cabo considerando el tiempo de observación como un factor adicional .... página 172


library(ggplot2)
# Ejemplo en página 176
pd <- position_dodge(0.1) # mover barras de errores .05 hacia la izquierda y la derecha
plot.interaction <- ggplot(data = medias, aes(x        = TRAT , 
                                              y        = LOPROP, 
                                              group    = TIEMPO, 
                                              color    = TIEMPO, 
                                              linetype = TIEMPO, 
                                              shape    = TIEMPO)) +
  geom_errorbar(aes(ymin = LOPROP-se, ymax = LOPROP+se), width = .1, position=pd) +
  geom_line(position = pd) + geom_point(position = pd) + ggtitle("Gráfico con errores típicos") +
  xlab("TRATAMIENTOS") + ylab("Log10(No. Propágalos)"); plot.interaction


Ambas colecciones de códigos cargan las tablas 13.2 y 13.2c directamente desde el escritorio. Fue necesario crear dos archivos porque la organización de variables para ejecutar ANOVAs separados y combinados es diferente.



Códigos para R


#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA                              
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010       
# CAPÍTULO 13: ANÁLISIS DE LA VARIANZA COMBINADO 
# (páginas 167 - 172) 
#************************************************************

#### Importar Tabla 13.2
datos <- read.csv("C:/Users/Administrator/Desktop/tabla 13.2.csv", header=TRUE)
# Las transformaciones logarítmicas fueron hechas en Excel
# '=LOG10(PROPT1+1)' y '=LOG10(PROPT2+1)' 
attach(datos)

#### Etiquetar variable 'TRAT'
TRAT   <- factor(TRAT, levels = c(0,1,2), labels = c("Sin solarizar", "Solarizado 1", "Solarizado 2"))
REP    <- factor(REP) 
MUES   <- factor(MUES)
datID  <- data.frame(TRAT, REP, MUES, PROPT1, PROPT2, LOPROPT1, LOPROPT2)
attach(datID); detach(datos); str(datID)
#### ANOVA de tiempos individuales
# ANOVA para LOPROP (T=1)
MODELOT1 <- aov(LOPROPT1 ~ TRAT * REP + Error(TRAT + REP), datID)
summary(MODELOT1) # opción 1
## 
## Error: TRAT
##      Df Sum Sq Mean Sq
## TRAT  2  15.44   7.721
## 
## Error: REP
##     Df  Sum Sq Mean Sq
## REP  3 0.07266 0.02422
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## TRAT:REP   6  0.503 0.08385     1.5  0.194
## Residuals 60  3.354 0.05590
MODELOT1 <- aov(LOPROPT1 ~ TRAT * REP + TRAT*REP, datID)
summary(MODELOT1) # opción 2
##             Df Sum Sq Mean Sq F value Pr(>F)    
## TRAT         2 15.442   7.721 138.126 <2e-16 ***
## REP          3  0.073   0.024   0.433  0.730    
## TRAT:REP     6  0.503   0.084   1.500  0.194    
## Residuals   60  3.354   0.056                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA para LOPROP (T=2)
MODELOT2 <- aov(LOPROPT2 ~ TRAT * REP + Error(TRAT + REP), datID) 
summary(MODELOT2) # opción 1
## Error: TRAT
##      Df Sum Sq Mean Sq
## TRAT  2  9.755   4.878
## 
## Error: REP
##     Df  Sum Sq Mean Sq
## REP  3 0.05174 0.01725
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## TRAT:REP   6  0.464 0.07732   1.061  0.396
## Residuals 60  4.373 0.07288
MODELOT1 <- aov(LOPROPT2 ~ TRAT * REP + TRAT*REP, datID)
summary(MODELOT1) # opción 2
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## TRAT         2  9.755   4.878  66.931 5.24e-16 ***
## REP          3  0.052   0.017   0.237    0.870    
## TRAT:REP     6  0.464   0.077   1.061    0.396    
## Residuals   60  4.373   0.073                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
detach(datID)
rm(list = ls()) # remueve (casi) todo en el ambiente de trabajo; útil para cargar otra tabla
#### Importar Tabla 13.2c con variable 'TIEMPO" y 'LOPROP' en columnas individuales 
datosC <- read.csv("C:/Users/Administrator/Desktop/tabla 13.2c.csv", header=TRUE)
attach(datosC)
#### Etiquetar variable 'TRAT' 
TRAT   <- factor(TRAT, levels = c(0,1,2), labels = c("Sin solarizar", "Solarizado 1", "Solarizado 2"))
TIEMPO <- factor(TIEMPO, levels = c(1,2)  , labels = c("Tiempo 1", "Tiempo 2"))
REP    <- factor(REP)
MUES   <- factor(MUES)
datIDc <- data.frame(TRAT, TIEMPO, REP, MUES, PROP, LOPROP)
attach(datIDc)
# Mostrando las primeras y las últimas 5 líneas de la tabla 13.2c
head(datIDc, 5)
##            TRAT   TIEMPO REP MUES PROP   LOPROP
## 1 Sin solarizar Tiempo 1   1    1 18.4 1.287802
## 2 Sin solarizar Tiempo 1   1    2 19.8 1.318063
## 3 Sin solarizar Tiempo 1   1    3 23.6 1.390935
## 4 Sin solarizar Tiempo 1   1    4 60.2 1.786751
## 5 Sin solarizar Tiempo 1   1    5 34.5 1.550228
tail(datIDc, 5)
##             TRAT   TIEMPO REP MUES PROP    LOPROP
## 140 Solarizado 2 Tiempo 2   4    2  6.4 0.8692317
## 141 Solarizado 2 Tiempo 2   4    3  1.3 0.3617278
## 142 Solarizado 2 Tiempo 2   4    4  2.1 0.4913617
## 143 Solarizado 2 Tiempo 2   4    5  0.7 0.2304489
## 144 Solarizado 2 Tiempo 2   4    6  0.0 0.0000000
#### Estadísticas descriptivas por grupo
source("C:/Users/Administrator/Desktop/summarySE.R") # Función summarySE
medias <- summarySE(datIDc, measurevar="LOPROP", groupvars=c("TRAT", "TIEMPO"))
attach(medias); medias
##            TRAT   TIEMPO  N    LOPROP        sd         se         ci
## 1 Sin solarizar Tiempo 1 24 1.2957658 0.2613633 0.05335056 0.11036404
## 2 Sin solarizar Tiempo 2 24 1.2651218 0.2503973 0.05111214 0.10573351
## 3  Solarizado 1 Tiempo 1 24 0.4502163 0.2626273 0.05360857 0.11089777
## 4  Solarizado 1 Tiempo 2 24 0.5118877 0.2552353 0.05209969 0.10777641
## 5  Solarizado 2 Tiempo 1 24 0.2180824 0.1832201 0.03739964 0.07736704
## 6  Solarizado 2 Tiempo 2 24 0.4593394 0.2910055 0.05940126 0.12288086
#### Código de gráfico de interacción
library(ggplot2)
# Ejemplo anterior invertido 
pd <- position_dodge(0.1) # mover barras de errores .05 hacia la izquierda y la derecha
plot.interaction <- ggplot(data = medias, aes(x        = TIEMPO , 
                                              y        = LOPROP, 
                                              group    = TRAT, 
                                              color    = TRAT, 
                                              linetype = TRAT, 
                                              shape    = TRAT)) +
  geom_errorbar(aes(ymin = LOPROP-se, ymax = LOPROP+se), width = .1, position=pd) +
  geom_line(position = pd) + geom_point(position = pd) + ggtitle("Gráfico con errores típicos") +
  xlab("TIEMPOS") + ylab("Log10(No. Propágalos)"); plot.interaction
detach(medias)
#### ANOVA combinado
# Modelo TRAT REP TRAT*REP(E) TIEMPO TRAT*TIEMPO TRAT*REP*TIEMPO(E)
MODELO <- aov(LOPROP ~ TRAT  * REP * TIEMPO + TIEMPO*TRAT, datIDc) 
summary(MODELO)
##                  Df Sum Sq Mean Sq F value Pr(>F)    
## TRAT              2 24.738  12.369 192.109 <2e-16 ***
## REP               3  0.036   0.012   0.184 0.9069    
## TIEMPO            1  0.297   0.297   4.606 0.0339 *  
## TRAT:REP          6  0.772   0.129   1.999 0.0710 .  
## TRAT:TIEMPO       2  0.459   0.229   3.563 0.0314 *  
## REP:TIEMPO        3  0.089   0.030   0.460 0.7110    
## TRAT:REP:TIEMPO   6  0.195   0.032   0.504 0.8042    
## Residuals       120  7.726   0.064                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo TRAT TRAT*REP(E) TIEMPO TRAT*TIEMPO
MODELO <- aov(LOPROP ~ TRAT + REP + TRAT*REP + TRAT*TIEMPO, datIDc) 
summary(MODELO)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## TRAT          2 24.738  12.369 199.206 <2e-16 ***
## REP           3  0.036   0.012   0.191 0.9023    
## TIEMPO        1  0.297   0.297   4.776 0.0307 *  
## TRAT:REP      6  0.772   0.129   2.073 0.0607 .  
## TRAT:TIEMPO   2  0.459   0.229   3.695 0.0275 *  
## Residuals   129  8.010   0.062                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### Contrastes pendientes...



Códigos para SAS


/***********************************************************
 EXPERIMENTACIÓN EN AGRICULTURA                             
 Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010      
 CAPÍTULO 13: ANÁLISIS DE LA VARIANZA COMBINADO
 (páginas 172 - 176)              
************************************************************/

* Creando una libraria en donde se encuentran datos de SAS;
LIBNAME lib1 'C:\Users\Administrator\Desktop';
RUN;

* Importando tabla 13.2.csv;
PROC IMPORT OUT= datos
            DATAFILE="C:\Users\Administrator\Desktop\tabla 13.2.csv"
            DBMS=csv
            REPLACE;
     GETNAMES=YES;
     DATAROW=2;
RUN;
ODS HTML;

 * 'Imprimiento' datos in la pantalla del Output;
PROC PRINT DATA = datos (OBS=20);
RUN;

* Enlistando variables;
PROC CONTENTS DATA = datos;
RUN;

* Etiquetar TRATAMIENTOS;
PROC FORMAT;
       VALUE TRATCODIGO
       0 = 'Sin solarizar'
       1 = 'Solarizado 1'
       2 = 'Solarizado 2';
RUN;

DATA datosID;
       SET datos;
       FORMAT TRAT TRATCODIGO. REP MUES PROPT1 LOPROP1 LOPROP2;
RUN;

* Datos con ID ;
PROC PRINT DATA = datosID;
RUN;

* Ver resultados en página 173;
* Opción 1: ANOVA LOPROP (T=1);
PROC GLM DATA = datosID;
       CLASS TRAT REP MUES;
       MODEL LOPROPT1 = TRAT REP TRAT*REP / SS3;
       TEST  H =TRAT REP E = TRAT*REP;
RUN;

* Opción 2: ANOVA LOPROP (T=1);
PROC MIXED DATA = datosID METHOD = TYPE3;
       CLASS TRAT REP MUES;
       MODEL LOPROPT1 = TRAT REP /DDFM=SATTERTH;
       RANDOM REP*TRAT;
RUN;

* Ver primeros resultados en página 174;
* Opción 1: ANOVA LOPROP (T=2);
PROC GLM DATA = datosID;
       CLASS TRAT REP MUES;
       MODEL LOPROPT2 = TRAT REP TRAT*REP / SS3;
       TEST  H =TRAT REP E = TRAT*REP;
RUN;

* Ver segundos resultados en página 174;
* Opción 2: ANOVA LOPROP (T=2);
PROC MIXED DATA = datosID METHOD = TYPE3;
       CLASS TRAT REP MUES;
       MODEL LOPROPT2 = TRAT REP /DDFM=SATTERTH;
       RANDOM REP*TRAT;
RUN;

* Variable TIEMPO agregada en una sola columna;
* Importando tabla 13.2c.csv;
PROC IMPORT OUT= datos2
            DATAFILE="C:\Users\Administrator\Desktop\tabla 13.2c.csv"
            DBMS=csv
                     REPLACE;
     GETNAMES=YES;
     DATAROW=2;
RUN;

* Etiquetar TRATAMIENTOS;
PROC FORMAT;
       VALUE TRATCODIGO
       0 = 'Sin solarizar'
       1 = 'Solarizado 1'
       2 = 'Solarizado 2';
RUN;

DATA datosIDc;
       SET datos2;
       FORMAT TRAT TRATCODIGO. REP MUES TIEMPO PROP LOPROP;
RUN;

* Enlistando variables;
PROC CONTENTS DATA = datosIDc;
RUN;

* Datos con ID ;
PROC PRINT DATA = datosIDc;
RUN;

*Ordenando datos por TRAT;
PROC SORT DATA = datosIDc; BY TRAT; RUN;
PROC UNIVARIATE DATA = datosIDc;
       VAR LOPROP;
       BY TRAT TIEMPO;
RUN;

* Modelo TRAT REP TRAT*REP(E) TIEMPO TRAT*TIEMPO TRAT*REP*TIEMPO(E)
  Opción 1: ANOVA Combinado LOPROP (página 174);
PROC GLM DATA = datosIDc;
       CLASS TRAT REP TIEMPO MUES;
       MODEL LOPROP = TRAT REP TRAT*REP TIEMPO TRAT*TIEMPO TRAT*REP*TIEMPO / SS3;
       TEST  H =TRAT REP TIEMPO E = TRAT*REP*TIEMPO;
RUN;

* Opción 2: ANOVA Combinado LOPROP (página 174);
PROC MIXED DATA = datosIDc METHOD = TYPE3;
       CLASS TRAT REP TIEMPO MUES;
       MODEL LOPROP = TRAT REP TRAT*REP TIEMPO TRAT*TIEMPO / DDFM=SATTERTH;
       RANDOM REP*TRAT*TIEMPO;
RUN;

* Modelo TRAT TRAT*REP(E) TIEMPO TRAT*TIEMPO
  Opción 1: ANOVA Combinado LOPROP (página 175);
PROC GLM DATA = datosIDc;
       CLASS TRAT REP TIEMPO;
       MODEL LOPROP = TRAT REP TRAT*REP TIEMPO TRAT*TIEMPO / SS3;
       TEST  H =TRAT TIEMPO E = TRAT*TIEMPO;
RUN;

* Opción 2: ANOVA Combinado LOPROP (página 175);
PROC MIXED DATA = datosIDc METHOD = TYPE3;
       CLASS TRAT REP TIEMPO;
       MODEL LOPROP = TRAT REP TRAT*REP TIEMPO / DDFM=SATTERTH;
       RANDOM TRAT*TIEMPO;
RUN;

*Contrastes ortogonales;
PROC GLM DATA = datosIDc;
CLASS TRAT TIEMPO;
MODEL LOPROP = TRAT TIEMPO TRAT*TIEMPO / SOLUTION;
*Tratamaientos e interacciones        S0 S1 S2            S0T1 S0T2 S1T1 S1T2 S2T1 S2T2;         
CONTRAST 'T1 : S0 vs. (S1 + S2)' TRAT  2 -1 -1 TRAT*TIEMPO   2   0   -1    0   -1   0 ;
CONTRAST 'T1 : S1 vs. S2'        TRAT  0  1 -1 TRAT*TIEMPO   0   0    1    0   -1   0 ;
CONTRAST 'T2 : S0 vs. (S1 + S2)' TRAT  2 -1 -1 TRAT*TIEMPO   0   2    0   -1    0  -1 ;
CONTRAST 'T2 : S1 vs. S2'        TRAT  0  1 -1 TRAT*TIEMPO   0   0    0    1    0  -1 ;
*Tiempos (Tiempo 1 = T1, Tiempo 2 = T2)  T1 T2;
CONTRAST 'S0 : T1 vs. T2'        TIEMPO   1 -1 TRAT*TIEMPO   1  -1    0    0    0   0 ;
CONTRAST 'S1 : T1 vs. T2'        TIEMPO   1 -1 TRAT*TIEMPO   0   0    1   -1    0   0 ;
CONTRAST 'S2 : T1 vs. T2'        TIEMPO   1 -1 TRAT*TIEMPO   0   0    0    0    1  -1 ;
RUN;
QUIT;

ODS HTML CLOSE;


-------------------------------------------------------


Si desean repetir o ajustar códigos a datos de sus experimentos, bajen los archivos y cambien directorios “C:/Users/Administrator/Desktop/tabla 13.2.csv” al de su computador, dejando únicamente el nombre de la tabla a usar.


Recomendamos, además, vean esta guía sobre cómo crear contrastes ortogonales en SAS. 




No hay comentarios:

Publicar un comentario