EXAMEN PARCIAL

Alumno: Giancarlo Javier Kuway Chero

DATA

resultados lab 58.7 A 62.7 B 55.9 C 60.7 D 61.4 A 64.5 B 56.1 C 60.3 D 60.9 A 63.1 B 57.3 C 60.9 D 59.1 A 59.2 B 55.2 C 61.4 D 58.2 A 60.3 B 58.1 C 62.3 D

PAQUETES UTILIZADOS

# Cargar paquetes
library(agricolae)
library(tidyverse)
library(apaTables)
library(gridExtra)
library(tseries) # jarque.bera
library(nortest) # Lilliefors (Kolmogorov-Smirnov)
library(car) # Prueba de leven
library(lmtest) # Prueba de durbin Watson

LECTURA DE DATOS

rm(list=ls())
library(readxl)
datos <- read_excel("a.xlsx")
View(datos)

PREGUNTA 1

VERIFICACIÓN DE SUPUESTOS

datos$lab<-as.factor(datos$lab)
modelo<-lm(resultados~lab,data=datos)

NORMALIDAD

METODO GRAFICO

e <- modelo$residuals
rstandard(modelo)
##          1          2          3          4          5          6          7 
## -0.7322060  0.5644088 -0.4728831 -0.3203401  1.3271234  1.9372951 -0.3203401 
##          8          9         10         11         12         13         14 
## -0.6254260  0.9457661  0.8694946  0.5949174 -0.1677972 -0.4271202 -2.1050923 
##         15         16         17         18         19         20 
## -1.0067833  0.2135601 -1.1135633 -1.2661062  1.2050891  0.9000032
residuals(modelo)
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## -0.96  0.74 -0.62 -0.42  1.74  2.54 -0.42 -0.82  1.24  1.14  0.78 -0.22 -0.56 
##    14    15    16    17    18    19    20 
## -2.76 -1.32  0.28 -1.46 -1.66  1.58  1.18
rstandard(modelo)
##          1          2          3          4          5          6          7 
## -0.7322060  0.5644088 -0.4728831 -0.3203401  1.3271234  1.9372951 -0.3203401 
##          8          9         10         11         12         13         14 
## -0.6254260  0.9457661  0.8694946  0.5949174 -0.1677972 -0.4271202 -2.1050923 
##         15         16         17         18         19         20 
## -1.0067833  0.2135601 -1.1135633 -1.2661062  1.2050891  0.9000032
e/summary(modelo)$sigma
##          1          2          3          4          5          6          7 
## -0.6549050  0.5048226 -0.4229595 -0.2865209  1.1870153  1.7327694 -0.2865209 
##          8          9         10         11         12         13         14 
## -0.5593980  0.8459189  0.7776997  0.5321103 -0.1500824 -0.3820279 -1.8828518 
##         15         16         17         18         19         20 
## -0.9004943  0.1910140 -0.9960013 -1.1324398  1.0778644  0.8049874
g1<-data.frame(e)|> # alt+124 y alt+62 
  ggplot(aes(x=e))+
  geom_histogram(bins = 6,# Numero de intervalos
                 colour = "black",
                 fill = "blue")+
  theme_bw()+
  labs(x = "Residuales", y = "Frecuencia")
g2<-data.frame(e)|>
  ggplot(aes(sample=e))+
  stat_qq(size = 2) +
  stat_qq_line(distribution = stats::qnorm)+
  labs(x = "Cuantil teórico",
       y = "Residuales")+
  theme_bw()

grid.arrange(g1,g2,ncol=2)

Gráficamente, el histograma no nos ayuda a apreciar si sigue o no una distribucion normal, mientas que el segundo grafico de linea con puntos, si nos demuestra que existe cierta cercania al comportamiento normal. Sin embargo, siempre sera mejor verificarlo con una prueba formal, como acontinuación.

METODO FORMAL

\[H0: Los\ residuos\ del\ modelo\ tienen\ distribución\ normal.\] \[Ha:Los\ residuos\ del\ modelo\ no\ tienen\ distribución\ normal\]

shapiro.test(modelo$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.97953, p-value = 0.9279
jarque.bera.test(modelo$residuals)
## 
##  Jarque Bera Test
## 
## data:  modelo$residuals
## X-squared = 0.36272, df = 2, p-value = 0.8341
lillie.test(modelo$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo$residuals
## D = 0.12257, p-value = 0.5996
ks.test(e, "pnorm", mean = mean(e), sd = sd(e))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  e
## D = 0.12257, p-value = 0.8895
## alternative hypothesis: two-sided

En cada una de las pruebas, los pvalores han sido mayor a 0.05 lo cual significa que no se rechaza la hipotesis nula. Eso quiere decir que los residuales tienen una distribucion normal.

HOMOGENEIDAD

METODO GRAFICO

y_est<-modelo$fitted.values
g3<-data.frame(y_est,e,lab=datos$lab)|>
  ggplot(aes(x=y_est,y=e))+
  geom_point(size=3.5)+
  geom_hline(yintercept = 0,col="red",linetype="dashed")+
  xlab("Valores ajustados")+ylab("Residuales") +
  ggtitle("Residuales vs Valores ajustados") +
  theme_bw()
g4<-data.frame(y_est,e,lab=datos$lab)|>
  ggplot(aes(x=lab,y=e))+
  geom_point(size=3.5)+
  geom_hline(yintercept = 0,col="red",linetype="dashed")+
  xlab("Dieta")+ylab("Residuales") +
  ggtitle("Residuales vs lab") +
  theme_bw()
grid.arrange(g3,g4,ncol=2)

Supuestamente los puntos no tienen forma de cono, pero sigue sin ser suficiente para determinar si cumple o no con el supuesto de homogeneidad, por lo cual se procedera con el metodo formal.

METODO FORMAL

\[H0:σ_{21}=σ_{22}=σ_{23}=σ_{24}=σ\ (homocedasticidad)\]

\[Ha: al\ menos\ un\ par\ σ_{2i}\ es\ diferente\ (heterocedasticidad)\]

bartlett.test(e~datos$lab)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  e by datos$lab
## Bartlett's K-squared = 3.8287, df = 3, p-value = 0.2806
leveneTest(y = e, group = datos$lab, center = "median")
## Levene's Test for Homogeneity of Variance (center = "median")
##       Df F value Pr(>F)
## group  3  1.0522 0.3967
##       16

En cada una de las pruebas, los pvalores han sido mayor a 0.05 lo cual significa que no se rechaza la hipotesis nula. Eso quiere decir que los residuales cumplen con el supuesto de homogeneidad.

INDEPENDENCIA

METODO GRAFICA

data.frame(e) |>
  ggplot(aes(x=1:nrow(datos),y=e))+
  geom_point(size = 3) +
  geom_line()+
  geom_hline(yintercept=0)+
  labs(x = "Orden", y = "Residual", 
       title = "Evaluación de independencia", 
       subtitle = "Modelo")+
  theme_minimal() 

plot(modelo$residuals, lag(modelo$residuals,1))

acf(e)

Supuestamente en el primer grafico, los puntos no siguen ningun patron mientras que en el tercer grafico ninguna bara supera los limites del grafico, por lo que deberian cumplir con el supuesto de independencia. Sin embargo, como se menciono anteriormente, mejor optar por la siguiente prueba formal.

METODO FORMAL

\[H0: Los\ errores\ del\ modelo\ son\ independientes\] \[Ha: Los\ errores\ del\ modelo\ no\ son\ independientes\]

dwtest(modelo,alternative =c("two.sided"))
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.3876, p-value = 0.285
## alternative hypothesis: true autocorrelation is not 0

En cada una de las pruebas, los pvalores han sido mayor a 0.05 lo cual significa que no se rechaza la hipotesis nula. Eso quiere decir que los residuales cumplen con el supuesto de independencia.

PREGUNTA 2

ANALISIS DE LA VARIANZA (ANOVA)

\[H0: u_1=u_2=u_3=u_4\]

\[H1: Al\ menos\ un\ u_i\ es\ diferente\ a\ los\ demàs\ (i=1,2,3,4)\]

anova(modelo)
## Analysis of Variance Table
## 
## Response: resultados
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## lab        3 85.926 28.6418  13.329 0.0001282 ***
## Residuals 16 34.380  2.1488                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como tiene Pvalor menor a 0.05 podemos decir que se rechaza la hipotesis nula, lo cual significa que existe diferencia en la media de los labs. Eso quiere decir que las medias de los analisis realizados en los diferentes laboratorios es significativamente diferente.

PREGUNTA 3

COMPARACIONES MULTIPLES

PRUEBA DE TUKEY

Tukey_man<-HSD.test(modelo,"lab",group=T,alpha=0.05,console=T)
## 
## Study: modelo ~ "lab"
## 
## HSD Test for resultados 
## 
## Mean Square Error:  2.14875 
## 
## lab,  means
## 
##   resultados       std r  Min  Max
## A      59.66 1.4081903 5 58.2 61.4
## B      61.96 2.1605555 5 59.2 64.5
## C      56.52 1.1627553 5 55.2 58.1
## D      61.12 0.7694154 5 60.3 62.3
## 
## Alpha: 0.05 ; DF Error: 16 
## Critical Value of Studentized Range: 4.046093 
## 
## Minimun Significant Difference: 2.652429 
## 
## Treatments with the same letter are not significantly different.
## 
##   resultados groups
## B      61.96      a
## D      61.12      a
## A      59.66      a
## C      56.52      b

Partiendo de este resultado, podemos decir que tanto el laboratorio A, B y D cuentan con medias en la medicion similares mientras que el laboratorio C difienre con todas las demas.

plot(Tukey_man,variation = "range")

Graficamente podemos observar que en efecto, el laboratorio C es el que difiere de los demas laboratorios.

Tukey_man<-HSD.test(modelo,"lab",group=F,alpha=0.05,console=T)
## 
## Study: modelo ~ "lab"
## 
## HSD Test for resultados 
## 
## Mean Square Error:  2.14875 
## 
## lab,  means
## 
##   resultados       std r  Min  Max
## A      59.66 1.4081903 5 58.2 61.4
## B      61.96 2.1605555 5 59.2 64.5
## C      56.52 1.1627553 5 55.2 58.1
## D      61.12 0.7694154 5 60.3 62.3
## 
## Alpha: 0.05 ; DF Error: 16 
## Critical Value of Studentized Range: 4.046093 
## 
## Comparison between treatments means
## 
##       difference pvalue signif.        LCL        UCL
## A - B      -2.30 0.1015         -4.9524292  0.3524292
## A - C       3.14 0.0178       *  0.4875708  5.7924292
## A - D      -1.46 0.4195         -4.1124292  1.1924292
## B - C       5.44 0.0001     ***  2.7875708  8.0924292
## B - D       0.84 0.8019         -1.8124292  3.4924292
## C - D      -4.60 0.0007     *** -7.2524292 -1.9475708

Individualmente tambien se realiza la comparación y verificamos que solamente es significativa la diferencia, solo cuando se incluye al laboratorio C en la comparacion.

PREGUNTA 4

PRUEBA DE DUNNETT

\[H_0: Di = D1\] \[H_1: Di \neq D1\]

library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
dunnet_man<-glht(modelo,linfct = mcp(lab="Dunnett"))
summary(dunnet_man)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: lm(formula = resultados ~ lab, data = datos)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)  
## B - A == 0   2.3000     0.9271   2.481   0.0621 .
## C - A == 0  -3.1400     0.9271  -3.387   0.0100 *
## D - A == 0   1.4600     0.9271   1.575   0.3006  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(dunnet_man,cex.axis=0.9,las=1)

Si el grupo control es el laboratorio A, entonces podemos decir lo siguiente:

  • Laboratorio A vs B

Su pvalor es cercano al 0.05 sin embargo esta no es menor a este, tiene un valor de 0.062 lo cual significa que no hay suficiente evidencia para decir que son significativamente diferentes.

  • Laboratorio A vs C

Su pvalor es mucho mas pequeño, es de 0.01 lo cual hace que si se rechaze la hipotesis nula y se diga que las medias de obteniadas en los laboratorios son significativamente distintos.

  • Laboratorio A vs D

Su pvalor es superior al 0.05, contando con un valor de 0.30 es el que mayor similitud tiene con el laboratorio control, por lo que tambien NO se rechaza la hipotesis nula por lo que no hay evidencia suficiente para decir que las medias son distintas.