Primero cargo los datos en cuestión:

require(agricolae)
## Loading required package: agricolae
## Warning: package 'agricolae' was built under R version 4.0.5
require(faraway)
## Loading required package: faraway
## Warning: package 'faraway' was built under R version 4.0.5
require(ggplot2)
## Loading required package: ggplot2
data("rats")

Según esta base, el experimento consistió en un ensayo de 48 observaciones, en las que se sometieron un grupo de ratones a diferentes tipos de agentes tóxicos con la finalidad de comparar el efecto de diferentes toxinas y tratamientos sobre el tiempo de supervivencia de los ratones.

Las variables en cuestión son:

time = tiempo de supervivencia en términos de horas

poison = el tipo de veneno, un factor con los niveles I, II y III

treat = el tratamiento, un factor con los niveles A, B, C y D

Exploratorio Bivariado (tiempo ~ veneno)

Definimos para un solo par de variables un diagrama de cajas donde:

y=tiempo de vida

x= tipo de veneno

attach(rats)
dats=data.frame(poison,treat)
tapply(time,dats,mean,na.rm=TRUE)
##       treat
## poison      A     B      C      D
##    I   0.4125 0.880 0.5675 0.6100
##    II  0.3200 0.815 0.3750 0.6675
##    III 0.2100 0.335 0.2350 0.3250

Según este resumen, se evidencian las medias en el tiempo de vida de los ratones según el tipo de veneno empleado, en este caso se destaca que el veneno III es el que más incide en el poco tiempo de vida de los ratones inoculados.Además,los ratones bajo el efecto del veneno II y con el tratamiento B, exhiben mayor tiempo de vida tras inoculación.

Para comparar este resumen gráficamente, se hace un gráfico de cajas con los tres tipos de venenos empleados contra el tiempo de vida de los ratones de ensayo:

ggplot(rats,aes(x="veneno",y=time,fill=poison))+geom_boxplot()+theme_bw()

Según se observa, el veneno I es aquel con menor mortalidad en los ratones, a medida que el factor crece de II a III, la mortalidad aumenta, haciendo que el tiempo de vida medio de los ratones decrezca considerablemente.

##Modelo de Diseño - ANOVA de un solo factor

Se considera realizar una regresión simple entre el tiempo y el tipo de veneno.

modelo1=lm(time~poison,data = rats)#Este es mi modelo lineal
summary(modelo1)
## 
## Call:
## lm(formula = time ~ poison, data = rats)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31438 -0.15922 -0.03125  0.08594  0.69563 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61750    0.05234  11.799 2.28e-15 ***
## poisonII    -0.07313    0.07401  -0.988    0.328    
## poisonIII   -0.34125    0.07401  -4.611 3.32e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2093 on 45 degrees of freedom
## Multiple R-squared:  0.3438, Adjusted R-squared:  0.3146 
## F-statistic: 11.79 on 2 and 45 DF,  p-value: 7.656e-05

Según este modelo, existen sólo dos factores de la variable poison que si son significativos a la hora de relacionarlos con la variable de respuesta, al obtener p-value de 0.328, se añade que este factor no afecta fuertemente en sí el comportamiento del experimento. Hya que añadir que el p-value de toda la variable poison es bajo y por ende se relaciona fuertemente sobre el tiempo de vida de los ratones.

Tabla post-ANOVA de un solo factor

anova(modelo1) #Si en el ANOVA hay asterisco, se hace un análisis PostANOVA
## Analysis of Variance Table
## 
## Response: time
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## poison     2 1.0330 0.51651  11.786 7.656e-05 ***
## Residuals 45 1.9721 0.04382                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test1=LSD.test(modelo1,"poison")
test1
## $statistics
##      MSerror Df     Mean       CV  t.value       LSD
##   0.04382375 45 0.479375 43.66962 2.014103 0.1490704
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none poison   3  0.05
## 
## $means
##         time        std  r       LCL       UCL  Min  Max    Q25   Q50   Q75
## I   0.617500 0.20942779 16 0.5120913 0.7229087 0.31 1.10 0.4500 0.625 0.730
## II  0.544375 0.28936641 16 0.4389663 0.6497837 0.23 1.24 0.3575 0.420 0.635
## III 0.276250 0.06227627 16 0.1708413 0.3816587 0.18 0.38 0.2275 0.270 0.315
## 
## $comparison
## NULL
## 
## $groups
##         time groups
## I   0.617500      a
## II  0.544375      a
## III 0.276250      b
## 
## attr(,"class")
## [1] "group"

Según el análisis ANOVA con el test LSD de Fisher, si hay diferencias en el tiempo de vida medio de los ratones con los venenos I y III comparado al veneno de tipo II, ya que el veneno II al arrojar un p-value tan alto (0.328), se rechaza la hipótesis de que su media sea significativa.

Exploratorio Bivariado (tiempo ~ tratamiento)

Definimos para un solo par de variables un diagrama de cajas donde:

y=tiempo de vida

x= factor de tratamiento

attach(rats)
## The following objects are masked from rats (pos = 3):
## 
##     poison, time, treat
tapply(time,treat,mean,na.rm=TRUE)
##         A         B         C         D 
## 0.3141667 0.6766667 0.3925000 0.5341667

Este resumen indica que el tratamiento B aumenta la esperanza de supervivencia de los ratones de una forma moderadamente mayor que en resto de tratamientos; aunque el tratamiento D también es efectivo en este propósito. Los tratamientos A y C son poco efectivos contra las dosis de inoculación de venenos en los ensayos.

ggplot(rats,aes(x="tratamiento",y=time,fill=treat))+geom_boxplot()+theme_bw()

Según se observa, con la gráfica es mejor visualizar cómo los tratamientos B y D son los mejores en aumentar la esperanza de supervivencia de los ratones tras inoculación

Modelo de Diseño - ANOVA de un solo factor

Se considera realizar una regresión simple entre el tiempo y el tipo de tratamiento:

modelo2=lm(time~treat,data = rats)#Este es mi modelo lineal
summary(modelo2)
## 
## Call:
## lm(formula = time ~ treat, data = rats)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38667 -0.15292 -0.01417  0.12833  0.56333 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.31417    0.06282   5.001 9.62e-06 ***
## treatB       0.36250    0.08885   4.080 0.000186 ***
## treatC       0.07833    0.08885   0.882 0.382739    
## treatD       0.22000    0.08885   2.476 0.017196 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2176 on 44 degrees of freedom
## Multiple R-squared:  0.3065, Adjusted R-squared:  0.2593 
## F-statistic: 6.484 on 3 and 44 DF,  p-value: 0.0009921

A contrario de lo que se pensaba, el tratamiento C no posee un valor significativo sobre la regresión lineal, ya que su p-value es muy alto (0.382739), por lo que se descarta que este tratamiento haga diferencias en las medias de los tiempos de supervivencia. El modelo anterior al arrojar varios valores significativos, dan cabida a un análisis Post-ANOVA ya que indica que si hay diferencias significativas en los resultados.

Tabla post-ANOVA de un solo factor

anova(modelo2) #Si en el ANOVA hay asterisco, se hace un análisis PostANOVA
## Analysis of Variance Table
## 
## Response: time
##           Df  Sum Sq  Mean Sq F value    Pr(>F)    
## treat      3 0.92121 0.307069  6.4836 0.0009921 ***
## Residuals 44 2.08388 0.047361                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test2=LSD.test(modelo2,"treat")
test2
## $statistics
##     MSerror Df     Mean       CV  t.value       LSD
##   0.0473608 44 0.479375 45.39773 2.015368 0.1790557
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none  treat   4  0.05
## 
## $means
##        time       std  r       LCL       UCL  Min  Max    Q25   Q50    Q75
## A 0.3141667 0.1022882 12 0.1875552 0.4407782 0.18 0.46 0.2275 0.300 0.4075
## B 0.6766667 0.3208323 12 0.5500552 0.8032782 0.29 1.24 0.3775 0.665 0.8900
## C 0.3925000 0.1670125 12 0.2658885 0.5191115 0.22 0.76 0.2475 0.375 0.4425
## D 0.5341667 0.2194397 12 0.4075552 0.6607782 0.30 1.02 0.3525 0.505 0.6725
## 
## $comparison
## NULL
## 
## $groups
##        time groups
## B 0.6766667      a
## D 0.5341667     ab
## C 0.3925000     bc
## A 0.3141667      c
## 
## attr(,"class")
## [1] "group"

Un coeficiente de variación (CV) del 45.39% indica fuerte heterogeneidad de los datos, lo que da probabilidad a que las mediciones no sean del todo centradas a la media, son datos muy dispersos.Por ejemplo, para el tratamiento B, el tiempo máximo de vida se prolonga hasta 1.24 horas, con el tratamiento D, el tiempo se prolonga hasta 1.02 horas.

Para concluir todo y resumir las anteriores afirmaciones, se realiza un modelo de regresión lineal donde se incluyen los factores que afectan el tiempo de supervivencia de los ratones según el tipo de veneno y el tipo de tratamiento.

LinearModel <- lm(time ~ treat +poison, data=rats)
summary(LinearModel)
## 
## Call:
## lm(formula = time ~ treat + poison, data = rats)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25167 -0.09625 -0.01490  0.06177  0.49833 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.45229    0.05592   8.088 4.22e-10 ***
## treatB       0.36250    0.06458   5.614 1.43e-06 ***
## treatC       0.07833    0.06458   1.213  0.23189    
## treatD       0.22000    0.06458   3.407  0.00146 ** 
## poisonII    -0.07312    0.05592  -1.308  0.19813    
## poisonIII   -0.34125    0.05592  -6.102 2.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1582 on 42 degrees of freedom
## Multiple R-squared:  0.6503, Adjusted R-squared:  0.6087 
## F-statistic: 15.62 on 5 and 42 DF,  p-value: 1.123e-08

Tal como se dijo antes, el tratamiento C y la inoculación con el veneno II no son variables significaticas sobre el experimento, dado su p-value tan alto. Si existen diferencias significativas en los ensayos entre todos los grupos al existir solapamientos de los intervalos de confianza.

kruskal.test(time ~ treat, data=rats)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  time by treat
## Kruskal-Wallis chi-squared = 13.118, df = 3, p-value = 0.004387