packages

Instalar paquetes

Ver el video que demuestra como instalar paquetes enlace

cargar packages

Leyendo datos

1. Caso1 (y ~ x)

Para leer la base de datos directamente desde el archivo en Drive * Copiar el siguiente código, tal y como está * y luego pegarlo en la consola de R (ejecutarlo)

gs4_deauth()
ss= "https://docs.google.com/spreadsheets/d/1-be0wPZgFEKIGsQZyrVZgWq-P1XitcAMzMJqdGaTe14/edit?usp=sharing"
hoja= "Hoja5"
rango = "A1:B142"
caso1 <- read_sheet(ss,
                    sheet= hoja,
                    range= rango,
                    col_names= TRUE
                    )
caso1 <- caso1 %>%
  mutate_if(is.character, as.factor)

Paso 2. Explorar datos

hist(caso1$y, breaks = 10, probability = TRUE, col=NA)
lines(density(caso1$y), col="black", lwd=2)
lines(density(rnorm(100000, mean = mean(caso1$y), sd=sd(caso1$y))), col="red")

test de normalidad

shapiro.test(caso1$y)
## 
##  Shapiro-Wilk normality test
## 
## data:  caso1$y
## W = 0.9923, p-value = 0.644

Modelo

variable respuesta ? Variable predictora ?

ggplot( aes(y=y, x=x), data=caso1) +
  geom_hline(yintercept = mean(caso1$y), lty=2, lwd=1,col="blue") +
  geom_boxplot(alpha=0.2) +
  geom_point(
  position = position_jitter(width = 0.2),
             size=1) +
 # geom_curve(aes(x=x, y= aov.1$fitted.values, xend= x, yend= y), 
  #           curvature = 0.3) +
  stat_summary(geom = "pointrange", col= "brown2", linewidth= 0.6,
               fun.data = "mean_cl_normal" )

Paso 3. Análisis aov(y~x) o t.test(y~x) (si son dos grupos)

aov.1 <- aov(y ~ x, data=caso1)
summary(aov.1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## x             4   2703   675.6   8.106 6.83e-06 ***
## Residuals   136  11336    83.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Paso 4. Comparaciones no planeadas

plot(TukeyHSD(aov.1))

TukeyHSD(aov.1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y ~ x, data = caso1)
## 
## $x
##           diff         lwr       upr     p adj
## b-a -2.8793167 -10.1659098  4.407276 0.8101586
## c-a -0.4530643  -7.5284068  6.622278 0.9997775
## d-a  6.9723000  -1.0097629 14.954363 0.1175355
## e-a  7.8035167   0.7639989 14.843034 0.0217436
## c-b  2.4262524  -3.8540011  8.706506 0.8225342
## d-b  9.8516167   2.5650235 17.138210 0.0024906
## e-b 10.6828333   4.4429678 16.922699 0.0000533
## d-c  7.4253643   0.3500218 14.500707 0.0345836
## e-c  8.2565810   2.2647570 14.248405 0.0019335
## e-d  0.8312167  -6.2083011  7.870734 0.9975187

Paso 4. Conclusión ?

Por lo menos un grupo es diferente del resto (P= 6.83^{-6}). Las comparaciones no planeadas nos muestran las diferencias entre todos los grupos. Los p adj son los valores de P ajustados al número de comparaciones no planeadas posibles.

===============================================================

Aproximado como un lm()

lm.1 <- lm(y ~ x, data=caso1)
summary(lm.1)
## 
## Call:
## lm(formula = y ~ x, data = caso1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.4714  -6.0108   0.5879   5.4339  21.8349 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 153.2982     2.0415  75.091  < 2e-16 ***
## xb           -2.8793     2.6356  -1.092  0.27655    
## xc           -0.4531     2.5592  -0.177  0.85974    
## xd            6.9723     2.8871   2.415  0.01707 *  
## xe            7.8035     2.5462   3.065  0.00263 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.13 on 136 degrees of freedom
## Multiple R-squared:  0.1925, Adjusted R-squared:  0.1688 
## F-statistic: 8.106 on 4 and 136 DF,  p-value: 6.829e-06

Este método compara cada grupo contra el primero, de la forma en que se harîa comparando con la distribución de t. Esto es lo que se haría si se tuviera un diseño en el que el primer grupo es el grupo control. No es lo aconsejable si se quieren comparar ´todos contra todos´, por eso vemos que la comparación del a-d es más exagerada aquí (P= 0.017) que con las comparaciones de TukeyHSD.

El R-squared sí lo podemos calcular desde el lm(). Vemos que es R-squared= 0.19

ggplot( aes(y=y, x=x), data=caso1) +
  geom_hline(yintercept = mean(caso1$y[caso1$x=="a"]), lty=2, lwd=0.7,col="blue") +
#  geom_boxplot(alpha=0.2) +
 # geom_point(
#  position = position_jitter(width = 0.2),
 #            size=1) +
  geom_segment(aes(x=x[1], y= coef(lm.1)[c(1)] , xend= x[x=="b"][1], yend= mean(y[x=="b"])) )+
  geom_segment(aes(x=x[1], y= coef(lm.1)[c(1)] , xend= x[x=="c"][1], yend= mean(y[x=="c"])) )+
  geom_segment(aes(x=x[1], y= coef(lm.1)[c(1)] , xend= x[x=="d"][1], yend= mean(y[x=="d"])) )+
  geom_segment(aes(x=x[1], y= coef(lm.1)[c(1)] , xend= x[x=="e"][1], yend= mean(y[x=="e"])) )+
#, 
 #            curvature = 0.3) +
  stat_summary(geom = "pointrange", col= "brown2", linewidth= 0.6,
               fun.data = "mean_cl_normal" )