---
title: "Práctica 4"
output: html_document
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# MODELOS CON DOS FACTORES

## Ejercicio 1
Nos interesa determinar si el tiempo de coagulación (en minutos) del plasma sanguíneo presenta diferencias para 3 tratamientos y 2 concentraciones de adrenalina mezclada con el plasma. Para cada combinación de tratamiento y concentración de adrenalina, se tomaron 3 observaciones independientes. Los resultados se encuentran en el archivo coagulacion.txt.


Se trata de un diseño de dos factores clasificación cruzada, ya que se han realizado observaciones para cada combinación de un tratamiento con una concentración de adrenalina. 

Vamos a importar los datos y a representar un gráfico de interacción que nos ayude a decidir si consideramos un modelo con o sin presencia de interacción entre ambos factores:

```{r gráfico de coagulacion, message=FALSE}
coagulacion<-read.table("./datos/coagulacion.txt", head=T)
colnames(coagulacion)

#convertimos las variables tratamiento y concentracion en factores

coagulacion$tratamiento<-as.factor(coagulacion$tratamiento)
coagulacion$concentracion<-as.factor(coagulacion$concentracion)

#dibujamos el gráfico de interacción
attach(coagulacion)
interaction.plot(tratamiento,concentracion,tiempo)
interaction.plot(concentracion, tratamiento,tiempo)
detach(coagulacion)
```
Aparentemente, no hay interacción entre ambos factores, no obstante, para quedarnos más seguros ajustaremos un modelo que considere dicha interacción como posible y contrastaremos la aditividad de dicho modelo:

$$Y_{ijk}=\mu+\tau_i+\delta_j+(\tau\delta)_{ij}+\varepsilon_{ijk}$$

Ajustemos dicho modelo lineal a estos datos:

```{r modelo coagulacion, message=FALSE}
ml1<-lm(tiempo~tratamiento+concentracion+tratamiento:concentracion,data=coagulacion)
```


Veamos que este modelo verifica las condiciones teóricas para proceder a su análisis. En primer lugar dibujamos los gráficos de residuos:

```{r gráficos residuos coagulacion, message=FALSE}
layout(matrix(1:4,2))
plot(ml1)
```

Contrastemos la normalidad de los datos:

```{r normalidad coagulacion, message=FALSE}
shapiro.test(rstandard(ml1))
```

Contrastemos también la homocedasticidad:

```{r homocedasticidad coagulacion, message=FALSE}
library(lmtest)
bptest(ml1)
```


Por último, la no correlación de errores:


```{r autocorrelación coagulacion, message=FALSE}
library(lmtest)
dwtest(ml1)
```


***Comentarios a los diagnósticos:***




Procedamos al análisis, realizando el contraste

$$H_0: (\tau\delta)^*_{ij}=0 \quad\mbox{ para todo } i,j$$
Calculamos la tabla de análisis de la varianza:

```{r anova coagulacion, message=FALSE}
anova(ml1)
```
```{r Tukey, message=FALSE}
TukeyHSD(aov(ml1))
```
***Comentarios al análisis:***




Puesto que la presencia de interacción es significativa, separaremos los datos de cada concentración de adrenalina y veremos, en cada uno de los nuevos conjuntos de datos, cómo se comporta el tiempo de coagulación frente al tratamiento mediante un diseño completamente aleatorizado.

$$Y_{ijk}=\mu+\tau_i+\varepsilon_{ijk}$$


Separamos en primer lugar los datos de cada concentración de adrenalina:

```{r filtrado de datos, message=FALSE}
attach(coagulacion)
c1<-coagulacion[concentracion=="1",]
c2<-coagulacion[concentracion=="2",]
detach(coagulacion)
```

Ajustamos un diseño completamente aleatorizado de un factor para cada uno de los nuevos conjuntos de datos y procedemos a su análisis.


#### Concentración 1
```{r ajuste concentración 1, message=TRUE}
ml2<-lm(tiempo~tratamiento,data=c1)
anova(ml2)
```


Calculemos las medias para cada tratamiento:

```{r medias concentración 1, message=TRUE}
tapply(c1$tiempo, c1$tratamiento, mean)
```


Realicemos las comparaciones múltiples de Tukey

```{r Tukey concentración 1, message=TRUE}
TukeyHSD(aov(ml2))
```

#### Concentración 2
```{r ajuste concentración 2, message=TRUE}
ml3<-lm(tiempo~tratamiento,data=c2)
anova(ml3)
```


Calculemos las medias para cada tratamiento:

```{r medias concentración 2, message=TRUE}
tapply(c2$tiempo, c2$tratamiento, mean)
```


Realicemos las comparaciones múltiples de Tukey

```{r Tukey concentración 2, message=TRUE}
TukeyHSD(aov(ml3))
```



***Comentarios al análisis:***







## Ejercicio 2
En una fábrica se produce un determinado tipo de fibra textil. Se piensa que tanto el telar como el proveedor de la seda puede tener influencia en la resistencia de la fibra. Se diseña un experimento en el que se elaboran dos muestras fibra en cada uno de los cuatro telares y utilizando seda de los tres proveedores de la fábrica.

Los datos de resistencia obtenidos se encuentran en el archivo fibra.dat. Analizar la influencia que el telar y el proveedor pueden tener sobre sobre la resistencia de la fibra producida.



Vamos a importar los datos y a representarlos en un gráfico de interacción:


```{r gráfico de fibra, message=FALSE}
fibra<-read.table("./datos/fibra.dat",head=T)
fibra$telar<-as.factor(fibra$telar)
colnames(fibra)

attach(fibra)
interaction.plot(telar,proveedor,resistencia)
interaction.plot(proveedor,telar,resistencia)
detach(fibra)
```

Como se utilizado material de todos los proveedores en todos los telares, tenemos un modelo de clasificación cruzada. De nuevo consideraremos la posibilidad de que exista interacción: 

$$Y_{ijk}=\mu+\tau_i+\delta_j+(\tau\delta)_{ij}+\varepsilon_{ijk}$$

Ajustemos un modelo lineal a estos datos:

```{r modelo fibra, message=FALSE}
ml4<-lm(resistencia~telar+proveedor+telar:proveedor,data=fibra)
```


Veamos que este modelo verifica las condiciones teóricas para proceder a su análisis. En primer lugar dibujamos los gráficos de residuos:

```{r gráficos residuos fibra, message=FALSE}
layout(matrix(1:4,2))
plot(ml4)
```

Contrastemos la normalidad de los datos:

```{r normalidad fibra, message=FALSE}
shapiro.test(rstandard(ml4))
```

Contrastemos también la homocedasticidad:

```{r homocedasticidad fibra, message=FALSE}
library(lmtest)
bptest(ml4)
```


Por último, la no correlación de errores:


```{r autocorrelación fibra, message=FALSE}
library(lmtest)
dwtest(ml4)
```



***Comentarios a los diagnósticos:***



Procedamos al análisis, realizando el contraste

$$H_0: (\tau\delta)^*_{ij}=0 \quad\mbox{ para todo } i,j$$

Calculamos la tabla de análisis de la varianza:

```{r anova fibra, message=FALSE}
anova(ml4)
```

***Comentarios al análisis:***




Puesto que el término de interacción es claramente no significativo, ajustamos un modelo sin interacción:

$$Y_{ijk}=\mu+\tau_i+\delta_j+\varepsilon_{ijk}$$

```{r modelo fibra 2, message=FALSE}
ml5<-lm(resistencia~telar+proveedor,data=fibra)
```


Veamos que este modelo verifica las condiciones teóricas para proceder a su análisis. En primer lugar dibujamos los gráficos de residuos:

```{r gráficos residuos fibra 2, message=FALSE}
layout(matrix(1:4,2))
plot(ml5)
```

Contrastemos la normalidad de los datos:

```{r normalidad fibra 2, message=FALSE}
shapiro.test(rstandard(ml5))
```

Contrastemos también la homocedasticidad:

```{r homocedasticidad fibra 2, message=FALSE}
library(lmtest)
bptest(ml5)
```


Por último, la no correlación de errores:


```{r autocorrelación fibra 2, message=FALSE}
library(lmtest)
dwtest(ml5)
```



***Comentarios a los diagnósticos:***


Procedamos al análisis, realizando los contrastes

$$H_0: \tau_1=\tau_2=\tau_3=\tau_4\quad,\quad H'_0:\delta_1=\delta_2=\delta_3$$

Calculamos la tabla de análisis de la varianza:

```{r anova fibra 2, message=FALSE}
anova(ml5)
```

***Comentarios al análisis:***

Calculamos la resistencia media de la fibra por proveedor y por telar:

```{r medias fibra, message=TRUE}
tapply(fibra$resistencia, fibra$telar, mean)
tapply(fibra$resistencia, fibra$proveedor, mean)
```


Realizamos las comparaciones múltiples de Tukey:

```{r Tukey fibra, message=TRUE}
TukeyHSD(aov(ml5))
```

***Comentarios al análisis:***



## Ejercicio 3
Se estudia el grosor del remate de la superficie de piezas metálicas elaboradas en 4 máquinas. En este experimento cada máquina es manejada por 3 operadores diferentes, entre los cuales también estamos interesados en detectar diferencias; de cada uno de los operarios se eligen dos piezas, registrándose el grosor del remate de cada una de ellas.

A causa de la localización de las máquinas, se utilizaron distintos operarios para cada una de ellas.Los datos aparecen en el archivo remate.dat.


Como los obreros solo trabajan con una máquina, el diseño es anidado o jerarquizado. El factor principal es la máquina y el anidado el obrero.

$$Y_{ijk}=\mu+\tau_i+(\tau\delta)_{ij}+\varepsilon_{ijk}$$
Vamos a importar los datos:

```{r remate, message=FALSE}
remate<-read.table("./datos/remate.dat", head=T, sep="\t")

colnames(remate)

remate$maquina<-as.factor(remate$maquina)
remate$obrero<-as.factor(remate$obrero)
```

Ajustemos el modelo de dos factores anidado:

```{r ajuste aguacate, message=FALSE}
ml6<-lm(grosor~maquina+maquina:obrero,data=remate)
```




Veamos que este modelo verifica las condiciones teóricas para proceder a su análisis. En primer lugar dibujamos los gráficos de residuos:

```{r gráficos residuos remate, message=FALSE}
layout(matrix(1:4,2))
plot(ml6)
```

Contrastemos la normalidad de los datos:

```{r normalidad aguacate, message=FALSE}
shapiro.test(rstandard(ml6))
```

Contrastemos también la homocedasticidad:

```{r homocedasticidad aguacate, message=FALSE}
library(lmtest)
bptest(ml6)
```


Por último, la no correlación de errores:


```{r autocorrelación aguacate, message=FALSE}
library(lmtest)
dwtest(ml6)
```



***Comentarios a los diagnósticos:***


Procedamos al análisis, realizando el contraste

$$H_0:(\tau\delta)_{i1}=(\tau\delta)_{i2}=(\tau\delta)_{i3}\quad,\quad i=1,2,3,4$$
mediante la tabla de análisis de la varianza:

```{r anova remate, message=FALSE}
anova(ml6)
```



Calculemos los estimadores del grosor promedio de las tres máquinas:
```{r medias remate, message=FALSE}
tapply(remate$grosor, remate$maquina, mean)
```

***Comentarios al análisis:***


Procedamos a realizar las comparaciones múltiples por el método de Tukey:

```{r Tukey remate, message=FALSE}
TukeyHSD(aov(ml6))
```

***Comentarios al análisis:***


Separamos los datos por máquina para ver diferencias entre los obreros:

```{r filtrado remate por máquinas}
attach(remate)
m1<-remate[maquina=="1",]
m2<-remate[maquina=="2",]
m3<-remate[maquina=="3",]
m4<-remate[maquina=="4",]
detach(remate)
```


Ajustamos modelos de un factor, diseño completamente aleatorizado para cada máquina:

#### Máquina 1
```{r remate máquina 1, message=FALSE}
ml7<-lm(grosor~obrero,data=m1)
anova(ml7)
```

#### Máquina 2
```{r remate máquina 2, message=FALSE}
ml8<-lm(grosor~obrero,data=m2)
anova(ml8)
TukeyHSD(aov(ml8))
```

#### Máquina 3
```{r remate máquina 3, message=FALSE}
ml9<-lm(grosor~obrero,data=m3)
anova(ml9)
TukeyHSD(aov(ml9))
```

#### Máquina 4
```{r remate máquina 4, message=FALSE}
ml10<-lm(grosor~obrero,data=m4)
anova(ml10)
```





## Ejercicio 4
La Confederación Hidrográfica del Guadiana realiza un estudio para determinar la concentración media de mercurio en dicho río a su paso por Extremadura. Para ello se seleccionan tres estaciones pertenecientes a la Confederación en las tres comarcas que cruza el río (Siberia, Vegas Altas y Vegas Bajas).
Cada uno de estos observatorios toma 6 mediciones (en ppm) de la sustancia en cuestión. Los datos aparecen en el archivo conc.dat.
Analizar las diferencias entre comarcas y entre estaciones.