1. Lee el fichero dieta.csv en el que se recogen las siguientes variables y resuelve los apartados propuestos:
dieta.df <- read.table("dieta.csv", sep = ";", dec = ".", header = T)

head(dieta.df)
##   id edad tipoDiet     peso0     peso1     peso2
## 1  1    3        1 122.01756 108.64895  90.94509
## 2  2    3        1 140.86780  89.24381  88.86856
## 3  3    2        3 110.52651 115.39116 100.41721
## 4  4    2        3 120.74844 106.17004 105.34358
## 5  5    1        1  82.73085  96.49512  91.90629
## 6  6    2        1 118.06524 104.84634  94.96339
summary(dieta.df)
##        id             edad      tipoDiet     peso0            peso1       
##  Min.   : 1.00   Min.   :1   Min.   :1   Min.   : 82.73   Min.   : 85.86  
##  1st Qu.: 6.75   1st Qu.:1   1st Qu.:1   1st Qu.:110.06   1st Qu.:100.22  
##  Median :12.50   Median :2   Median :2   Median :116.26   Median :106.05  
##  Mean   :12.50   Mean   :2   Mean   :2   Mean   :115.91   Mean   :107.54  
##  3rd Qu.:18.25   3rd Qu.:3   3rd Qu.:3   3rd Qu.:122.41   3rd Qu.:114.56  
##  Max.   :24.00   Max.   :3   Max.   :3   Max.   :146.63   Max.   :135.73  
##      peso2       
##  Min.   : 86.48  
##  1st Qu.: 94.20  
##  Median : 99.60  
##  Mean   :100.61  
##  3rd Qu.:108.36  
##  Max.   :118.83
  1. Convierte en factor las variables tipoDiet y edad.
dieta.df[,2:3] <- factor(c(dieta.df$edad, dieta.df$tipoDiet))
str(dieta.df)
## 'data.frame':    24 obs. of  6 variables:
##  $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ edad    : chr  "3" "3" "2" "2" ...
##  $ tipoDiet: chr  "1" "1" "3" "3" ...
##  $ peso0   : num  122 140.9 110.5 120.7 82.7 ...
##  $ peso1   : num  108.6 89.2 115.4 106.2 96.5 ...
##  $ peso2   : num  90.9 88.9 100.4 105.3 91.9 ...
  1. Comprueba si hay diferencias en la media de los pesos antes de comenzar la dieta (Variable peso0) según el grupo de edad al que pertenece y encuentra en qué grupos están las diferencias.

Primero compruebo la normalidad en cada grupo de edad, para ello uso el test de Shapiro-Wilk (hay menos de 50 datos).

grupos_edad <- unique(dieta.df$edad)

for(x in grupos_edad){
  print(shapiro.test(dieta.df$peso0[dieta.df$edad == x]))
}
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso0[dieta.df$edad == x]
## W = 0.92652, p-value = 0.4849
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso0[dieta.df$edad == x]
## W = 0.83071, p-value = 0.06038
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso0[dieta.df$edad == x]
## W = 0.853, p-value = 0.1022

En todos los casos el p-valor es superior a 0.05, por lo que no se rechaza la hipótesis nula: los pesos iniciales siguen una distribución normal en cada grupo de edad.

Ahora se comprueba si existe homocedasticidad mediante el test de Barlett (es muy robusto si existe normalidad).

bartlett.test(peso0 ~ edad, data = dieta.df)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  peso0 by edad
## Bartlett's K-squared = 0.43715, df = 2, p-value = 0.8037

De nuevo obtenemos un p-valor superior a 0.05: hay homocedasticidad.

Llevo a cabo un ANOVA de una vía.

fitDieta <- aov(peso0 ~ edad, data = dieta.df)
summary(fitDieta)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## edad         2   2544  1272.0   13.46 0.000173 ***
## Residuals   21   1985    94.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-valor inferior a 0.001 indica un efecto muy significativo de la edad sobre el peso inicial.

Para saber en qué grupo o grupos difiere el peso inicial se hace una prueba de contrastes post-hoc. Se puede hacer tanto el contraste de Bonferroni como el de Holm (indicados para menos de 6 niveles).

pairwise.t.test(dieta.df$peso0, dieta.df$edad, p.adj = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dieta.df$peso0 and dieta.df$edad 
## 
##   1       2      
## 2 0.33973 -      
## 3 0.00015 0.00751
## 
## P value adjustment method: bonferroni
pairwise.t.test(dieta.df$peso0, dieta.df$edad, p.adj = "holm")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dieta.df$peso0 and dieta.df$edad 
## 
##   1       2      
## 2 0.11324 -      
## 3 0.00015 0.00501
## 
## P value adjustment method: holm

Como podemos ver, al comparar el grupo 3 con el 1 y con el 2 obtenemos p-valores inferiores a 0.05, por lo que hay una diferencia significativa en el peso0 del grupo 3 con respecto a los otros dos, mientras que no se aprecia una diferencia significativa entre los grupos 1 y 2.

  1. Comprueba si hay diferencias significativas en el peso final (peso2) según la interacción del grupo de edad al que pertenece el sujeto y el tipo de dieta al que se somete. Representa el gráfico de interacción con la edad en el eje x y el tipoDiet en tres líneas diferentes.

En primer lugar comprobamos la normalidad de todos los subconjuntos de peso2 construidos en base a los diferentes grupos de edad y tipos de dieta.

for(x in grupos_edad){
  print(shapiro.test(dieta.df$peso2[dieta.df$edad == x]))
}
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$edad == x]
## W = 0.91386, p-value = 0.382
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$edad == x]
## W = 0.9739, p-value = 0.9267
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$edad == x]
## W = 0.91981, p-value = 0.4284
tipos_dieta <- unique(dieta.df$tipoDiet)
for(y in grupos_edad){
  print(shapiro.test(dieta.df$peso2[dieta.df$tipoDiet == y]))
}
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$tipoDiet == y]
## W = 0.96944, p-value = 0.8936
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$tipoDiet == y]
## W = 0.91767, p-value = 0.4113
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  dieta.df$peso2[dieta.df$tipoDiet == y]
## W = 0.90802, p-value = 0.3403

En todos los casos el p-valor es superior a 0.05, por lo que no se rechaza la hipótesis nula: los pesos finales siguen una distribución normal tanto dentro de cada grupo de edad como dentro de cada tipo de dieta.

Se comprueba la homocedasticidad mediante el test de Barlett.

bartlett.test(peso2 ~ edad, data = dieta.df)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  peso2 by edad
## Bartlett's K-squared = 3.4329, df = 2, p-value = 0.1797
bartlett.test(peso2 ~ tipoDiet, data = dieta.df)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  peso2 by tipoDiet
## Bartlett's K-squared = 1.5668, df = 2, p-value = 0.4568

Ambos valores son superiores a 0.05, por lo que se considera que hay homocedasticidad.

Ahora se lleva a cabo un ANOVA de dos vías para evaluar tanto el efecto de la edad como el del tipo de dieta sobre el peso final además del efecto de la interacción entre ambas.

Para ello primero debemos ver si los datos son balanceados, puesto que de no serlo el orden SÍ importa.

table(dieta.df$edad)
## 
## 1 2 3 
## 8 8 8
table(dieta.df$tipoDiet)
## 
## 1 2 3 
## 8 8 8

Los datos son balanceados. Podemos continuar de formal convencional.

fitDiet2 <- aov(peso2 ~ edad * tipoDiet, data = dieta.df)
summary(fitDiet2)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## edad           2   84.2    42.1   1.320    0.296    
## tipoDiet       2 1391.1   695.5  21.823 3.62e-05 ***
## edad:tipoDiet  4  160.0    40.0   1.255    0.331    
## Residuals     15  478.1    31.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La única variable que afecta de forma significativa (el Pr es tan cercano a 0 que podría decirse que es el factor determinante de las diferencias de peso) es el tipo de dieta.

Podemos visualizarlo mediante el gráfico de interacción:

interaction.plot(dieta.df$edad, dieta.df$tipoDiet, dieta.df$peso2,
                 col = c("blue", "green", "red"),
                 lty = c(1, 3, 5),
                 ylab = "media del peso final",
                 xlab = "edad", trace.label = "dieta")

  1. Comprueba si la media de los pesos difieren dependiendo del momento de la dieta que se esté. Reestructura el data frame con la función melt() del paquete reshape2. En caso de que difieran, encuentra en qué grupos están las diferencias.

Primero reestructuro el dataframe y compruebo que esté como quiero.

library(reshape2)
## Warning: package 'reshape2' was built under R version 3.3.3
pesoRe <- melt( dieta.df, id = colnames(dieta.df)[1:3],
                measure = colnames(dieta.df)[4:6],
                variable.name = "medida",
                value.name = "peso" )

head(pesoRe)
##   id edad tipoDiet medida      peso
## 1  1    3        1  peso0 122.01756
## 2  2    3        1  peso0 140.86780
## 3  3    2        3  peso0 110.52651
## 4  4    2        3  peso0 120.74844
## 5  5    1        1  peso0  82.73085
## 6  6    2        1  peso0 118.06524

Ahora realizo un ANOVA para medidas repetidas (en este caso las diferentes medidas tomadas para la variable peso).

library(ez)
## Warning: package 'ez' was built under R version 3.3.3
options(contrasts = c("contr.sum", "contr.poly"))
ezANOVA(data = pesoRe, dv = peso,
         wid = id, within = medida,
         type = 3)
## Warning: Converting "id" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2 medida   2  46 9.007468 0.0004999883     * 0.2143998
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 medida 0.9348576 0.4766505      
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 2 medida 0.9388416 0.0006806327         * 1.019405 0.0004999883         *

Observando el p-valor inferior a 0.01 concluimos que existen diferencias significativas según el momento de la dieta en el que se esté.

Además, el test de Mauchly da un p-valor superior a 0.05 que al no ser significativo nos lleva a rechazar la hipótesis nula de ausencia de esfericidad, por lo que es correcto llevar a cabo el ANOVA.

Para terminar se hace la comparación por grupos para saber entre qué pesos hay una diferencia significativa.

pairwise.t.test(pesoRe$peso, pesoRe$medida, p.adj = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pesoRe$peso and pesoRe$medida 
## 
##       peso0   peso1  
## peso1 0.06168 -      
## peso2 0.00015 0.16121
## 
## P value adjustment method: bonferroni
pairwise.t.test(pesoRe$peso, pesoRe$medida, p.adj = "holm")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pesoRe$peso and pesoRe$medida 
## 
##       peso0   peso1  
## peso1 0.04112 -      
## peso2 0.00015 0.05374
## 
## P value adjustment method: holm

Interesante resultado: tanto haciendo el ajuste de Bonferroni como el de Holm el peso final es significativamente diferente del peso inicial y las diferencias entre el peso 1 y el final no son significativas (aunque según Holm estamos cerca del p-valor 0.05 que se toma como corte). Sin embargo, mientras que con el ajuste de Bonferroni se considera la dierencia entre el peso inicial y el 1 como no significativa, según el ajuste de Holm la diferencia sí es significativa, esto se debe a que ambos ajustes de corrección dan lugar a que aparezcan p-valores muy próximos a 0.05.

Así, dependiendo del método de ajuste aplicado, las conclusiones podrían ser diferentes.

  1. Una empresa quiere saber si existe relación entre el salario de un trabajador y las ausencias del mismo al trabajo. Para el estudio se dividió el salario en distintas categorías y se eligió aleatoriamente un grupo de trabajadores para determinar el número de días que habían faltado en los últimos tres años. ¿Se puede construir un modelo que relacione la categoría del salario y las ausencias en el trabajo? (Trabaja con el fichero “william.csv”)
william <- read.table(file = "william.csv", sep = ";", dec = ".", header = TRUE)

head(william)
##   salario ausencias
## 1      11        18
## 2      10        17
## 3       8        29
## 4       5        36
## 5       9        11
## 6       7        28
  1. Comprueba la normalidad de la variable explicativa (salario) y calcula la correlación entre salario y ausencias.
shapiro.test(william$salario)
## 
##  Shapiro-Wilk normality test
## 
## data:  william$salario
## W = 0.93541, p-value = 0.3281

El p-valor es superior a 0.05, por lo que no se rechaza la hipótesis nula: los valores de la variable salario siguen una distribución normal.

cor.test(william$salario, william$ausencias)
## 
##  Pearson's product-moment correlation
## 
## data:  william$salario and william$ausencias
## t = -7.4737, df = 13, p-value = 4.672e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9668476 -0.7211085
## sample estimates:
##        cor 
## -0.9006674

Existe una correlación significativa y negativa entre el salario y las ausencias.

  1. Realiza el ajuste del modelo.
model <- lm (ausencias ~ salario, data = william)
summary(model)
## 
## Call:
## lm(formula = ausencias ~ salario, data = william)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.516 -3.053  1.428  2.961  5.475 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  47.6002     3.0789  15.460 9.50e-10 ***
## salario      -3.0094     0.4027  -7.474 4.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.294 on 13 degrees of freedom
## Multiple R-squared:  0.8112, Adjusted R-squared:  0.7967 
## F-statistic: 55.86 on 1 and 13 DF,  p-value: 4.672e-06

R^2 ajustado es 0.8: parece que la recta modelo se ajusta bastante bien a los datos (explica un 80% de la variabilidad), aunque hay que hacer un diagnóstico del modelo para poder explicar con precisión su bondad.

Visualización:

plot(william$salario, william$ausencias,
     xlab = "Salario",
     ylab = "Ausencias",
     main = "Salario vs Ausencias")

abline(model, col = "red")

  1. Estudia la bondad de ajuste con la función anova() y explica el significado de los coeficientes obtenidos en la recta de regresión.
anova(model)
## Analysis of Variance Table
## 
## Response: ausencias
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## salario    1 1030.01 1030.01  55.857 4.672e-06 ***
## Residuals 13  239.72   18.44                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-valor menor que 0.01 (y muy próximo a 0) indica un rechazo de la hipótesis nula, es decir, el modelo de regresión lineal permite explicar la correlación entre ausencias y salario.

SSM de 1030.01 es muy grande: se ha hecho una gran mejora a la hora de predecir la variable dependiente (las ausencias).

El estadístico F también muestra que MSm es 56 veces más grande que MSr.

  1. Realiza el diagnóstico del modelo.

Normalidad de los residuos. Mediante el test de Shapiro-Wilk.

william$fitted.model <- fitted(model)
william$residuals.model <- residuals(model)
william$rstudent.model <- rstudent(model)

shapiro.test(william$rstudent.model)
## 
##  Shapiro-Wilk normality test
## 
## data:  william$rstudent.model
## W = 0.91538, p-value = 0.1637

El p-valor es mayor que 0.05, por lo que los datos siguen una distribución normal.

Homogeneidad de varianzas. Test de Breusch-Pagan.

library(lmtest)
## Warning: package 'lmtest' was built under R version 3.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.3.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 0.52916, df = 1, p-value = 0.467

Como el p-valor es superior a 0.05 se rechaza la hipótesis nula y se concluye que hay homogeneidad de las varianzas.

Incorrelación de los residuos. Test de Durbin-Watson.

dwtest(ausencias ~ salario, alternative = "two.sided", data = william)
## 
##  Durbin-Watson test
## 
## data:  ausencias ~ salario
## DW = 2.452, p-value = 0.3763
## alternative hypothesis: true autocorrelation is not 0

Como la hipótesis nula es que no exista correlación entre los residuos y se rechaza puesto que el p-valor es mayor que 0.05, concluimos que sí hay correlación entre los residuos.