Cargando librerias que seran utiles en el analisis de la informacion:
library(dplyr)Asignamos los valores (signos \(+\) y \(-\)) que pueden tomar las variables A(Gap) B(\(C_2F_6\) Flow) y C(Power). Como son dos replicas se pondran los valores dos veces.
ABC_values <- c(-1,-1,-1,
1,-1,-1,
-1,1,-1,
1,1,-1,
-1,-1,1,
1,-1,1,
-1,1,1,
1,1,1,
-1,-1,-1,
1,-1,-1,
-1,1,-1,
1,1,-1,
-1,-1,1,
1,-1,1,
-1,1,1,
1,1,1)El vector lo convertiremos en una matriz para luego hacer el data frame.
ABC_values <- matrix(ABC_values, ncol = 3, byrow = TRUE)
ABC_values [,1] [,2] [,3]
[1,] -1 -1 -1
[2,] 1 -1 -1
[3,] -1 1 -1
[4,] 1 1 -1
[5,] -1 -1 1
[6,] 1 -1 1
[7,] -1 1 1
[8,] 1 1 1
[9,] -1 -1 -1
[10,] 1 -1 -1
[11,] -1 1 -1
[12,] 1 1 -1
[13,] -1 -1 1
[14,] 1 -1 1
[15,] -1 1 1
[16,] 1 1 1Asignamos los valores al vector de las replicas del experimento.
replicates <- c(247,470,429,435,837,551,775,660,400,446,405,445,850,670,865,530)Creamos el data frame con los valores anteriores.
plasma <- data.frame(cbind(A = ABC_values[,1], B = ABC_values[,2], C = ABC_values[,3], replicate = replicates))
head(plasma)Asignaremos los totales de las replicas en una matriz.
totals <- c(647,916,834,880,1687,1221,1640,1190)
totals <- matrix(totals)
rownames(totals) <- c("(1)","a","b","ab","c","ac","bc","abc")
totals [,1]
(1) 647
a 916
b 834
ab 880
c 1687
ac 1221
bc 1640
abc 1190Calculando el estimado de los efectos principales: A, B y C. Donde \(n\) es la cantidad de replicas que se hicieron. En este caso son 2. La parte que se encuentra en los corchetes \([\hspace{0.1cm}]\) se define como el contraste. \[A = \frac{1}{4n}[a + ab + ac + abc - (1) - b - c - bc]\]
rep_num <- 2
A_contrast <- totals["a",] + totals["ab",] + totals["ac",] + totals["abc",] - totals["(1)",] - totals["b",] - totals["c",] - totals["bc",]
A_contrast a
-601 A_est_effect <- (1/(4*rep_num)) * A_contrast
A_est_effect a
-75.125 \[B = \frac{1}{4n}[b + ab + bc + abc - (1) - a - c - ac]\]
B_contrast <- totals["b",] + totals["ab",] + totals["bc",] + totals["abc",] - totals["(1)",] - totals["a",] - totals["c",] - totals["ac",]
B_contrast b
73 B_est_effect <- (1/(4*rep_num)) * B_contrast
B_est_effect b
9.125 \[C = \frac{1}{4n}[c + ac + bc + abc - (1) - a - c - ab]\]
C_contrast <- totals["c",] + totals["ac",] + totals["bc",] + totals["abc",] - totals["(1)",] - totals["a",] - totals["b",] - totals["ab",]
C_contrast c
2461 C_est_effect <- (1/(4*rep_num)) * C_contrast
C_est_effect c
307.625 Ahora se obtendran los contrastes y el estimado del efecto de las interacciones de dos factores: AB, AC y BC.
\[AB = \frac{1}{4n}[abc + ab + c + (1) - bc - b - ac - a]\]
AB_contrast <- totals["abc",] + totals["ab",] + totals["c",] + totals["(1)",] - totals["bc",] - totals["b",] - totals["ac",] - totals["a",]
AB_contrast abc
-207 AB_est_effect <- (1/(4*rep_num)) * AB_contrast
AB_est_effect abc
-25.875 \[AC = \frac{1}{4n}[(1) + b + ac + abc - a - ab - c - bc]\]
AC_contrast <- totals["(1)",] + totals["b",] + totals["ac",] + totals["abc",] - totals["a",] - totals["ab",] - totals["c",] - totals["bc",]
AC_contrast (1)
-1231 AC_est_effect <- (1/(4*rep_num)) * AC_contrast
AC_est_effect (1)
-153.875 \[BC = \frac{1}{4n}[(1) + a + bc + abc - b - ab - c - ac]\]
BC_contrast <- totals["(1)",] + totals["a",] + totals["bc",] + totals["abc",] - totals["b",] - totals["ab",] - totals["c",] - totals["ac",]
BC_contrast (1)
-229 BC_est_effect <- (1/(4*rep_num)) * BC_contrast
BC_est_effect (1)
-28.625 Y por ultimo calculamos el contraste y el estimado del efecto de la interaccion de tres factores: ABC. \[ABC = \frac{1}{4n}[abc + c + b + a - bc - ac - ab - (1)]\]
ABC_contrast <- totals["abc",] + totals["c",] + totals["b",] + totals["a",] - totals["bc",] - totals["ac",] - totals["ab",] - totals["(1)",]
ABC_contrastabc
239 ABC_est_effect <- (1/(4*rep_num)) * ABC_contrast
ABC_est_effect abc
29.875 Ahora podemos calcular la Suma de Cuadrados (SS) de cada una de estas interacciones: \[SS = \frac{(contrast)^2}{8n}\]
SS_A <- (A_contrast^2)/(8*rep_num)
SS_B <- (B_contrast^2)/(8*rep_num)
SS_C <- (C_contrast^2)/(8*rep_num)
SS_AB <- (AB_contrast^2)/(8*rep_num)
SS_AC <- (AC_contrast^2)/(8*rep_num)
SS_BC <- (BC_contrast^2)/(8*rep_num)
SS_ABC <- (ABC_contrast^2)/(8*rep_num)
c(SS_A,SS_B,SS_C,SS_AB,SS_AC,SS_BC,SS_ABC) a b c abc (1) (1) abc
22575.0625 333.0625 378532.5625 2678.0625 94710.0625 3277.5625 3570.0625 Ahora se hara el analisis de la varianza, se obtendran los P-values y se concluira cuales factores tienen los efectos mas estadisticamente significativos. Se usara la funcion anova() Analysis of Variance. Esta función recibe una regresion, por lo que usaremos todos los efectos: A, B, C, AB, AC, BC y ABC.
fit_for_AOV <- lm(replicate ~ A + B + C + (A*B) + (A*C) + (B*C) + (A*B*C), data = plasma)
anova(fit_for_AOV)Analysis of Variance Table
Response: replicate
Df Sum Sq Mean Sq F value Pr(>F)
A 1 22575 22575 5.6446 0.044837 *
B 1 333 333 0.0833 0.780240
C 1 378533 378533 94.6465 1.042e-05 ***
A:B 1 2678 2678 0.6696 0.436881
A:C 1 94710 94710 23.6808 0.001246 **
B:C 1 3278 3278 0.8195 0.391772
A:B:C 1 3570 3570 0.8926 0.372419
Residuals 8 31996 3999
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1En la columna Sum Sq podemos ver que tenemos los mismos resultados de los calculos de la suma de cuadrados que se calcularon anteriormente. Tambien se puede ver que el error es de \(31,996\). Luego de este analisis se puede concluir que las variables mas significativas para la regresion son: A, C y AC. Ahora ya podemos hacer un modelo de regresion lineal con esas variables.
fit_plasma <- lm(replicate ~ A + C + (A*C), data = plasma)
summary(fit_plasma)
Call:
lm(formula = replicate ~ A + C + (A * C), data = plasma)
Residuals:
Min 1Q Median 3Q Max
-123.25 -23.44 11.75 33.62 67.25
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 563.44 14.76 38.162 6.75e-14 ***
A -37.56 14.76 -2.544 0.025743 *
C 153.81 14.76 10.418 2.30e-07 ***
A:C -76.94 14.76 -5.211 0.000218 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 59.06 on 12 degrees of freedom
Multiple R-squared: 0.9222, Adjusted R-squared: 0.9027
F-statistic: 47.39 on 3 and 12 DF, p-value: 6.304e-07Por lo que podemos ver que nuestras variables si son estadisticamente significativas. La regresion con estas variables tiene los siguientes coeficientes:
fit_plasma$coefficients(Intercept) A C A:C
563.4375 -37.5625 153.8125 -76.9375 Y como resultado tenemos: \[563.438 - 37.563x_1 + 153.813x_3 - 76.938x_1x_3\]
Otro metodo para juzgar la significancia de los efectos es por medio del calculo de el error estandar de los efectos y compararlo con las magnitudes de los efectos y sus errores estandar.
La varianza estimada esta dada por el error de los minimos cuadrados en el analisis de la varianza. La varianza de cada efecto estimado esta dada por: \[Var(Effect) = Var\bigg(\frac{Contrast}{n2^{k-1}}\bigg) = \frac{1}{(n2^{k-1})^2}Var(Contrast)\]
Cada contraste es una combinacion lineal de \(2^k\) tratamientos totales y cada total consiste en \(n\) observaciones. Entonces: \[Var(Contrast) = n2^k\sigma^2\]
Por lo que en la ecuacion anterior se sustituye en la ecuacion de \(Var(Effect)\) por lo que la varianza de un efecto es: \[Var(Effect)=\frac{1}{(n2^{k-1})^2}n2^k\sigma^2 = \frac{1}{n2^{k-2}}\sigma^2\]
El error estandar estimado se encuentra reemplazando \(\sigma^2\) por su estimado \(s^2\) y despejando \(s\): \[se(Effect) =\sqrt{\frac{1}{n2^{k-2}}s^2}\]
El valor del estimado del error de minimos cuadrados se puede obtener del analisis de varianza en la columna de Mean Sq, el ultimo valor.
MSE <- anova(fit_for_AOV)$`Mean Sq`[8]
MSE[1] 3999.438se_effect <- sqrt(1/(rep_num*(2^(3-2)))*(MSE))
se_effect[1] 31.62055Y los limites de los errores en el estimado de los efectos con dos desviaciones estandar son:
limits <- c(A_est_effect - (2*se_effect), A_est_effect + (2*se_effect),
B_est_effect - (2*se_effect), B_est_effect + (2*se_effect),
C_est_effect - (2*se_effect), C_est_effect + (2*se_effect),
AB_est_effect - (2*se_effect), AB_est_effect + (2*se_effect),
AC_est_effect - (2*se_effect), AC_est_effect + (2*se_effect),
BC_est_effect - (2*se_effect), BC_est_effect + (2*se_effect),
ABC_est_effect - (2*se_effect), ABC_est_effect + (2*se_effect))
limits <- matrix(limits, ncol = 2, byrow = TRUE)
rownames(limits) <- c("A","B","C","AB","AC","BC","ABC")
colnames(limits) <- c("low", "high")
limits low high
A -138.36611 -11.88389
B -54.11611 72.36611
C 244.38389 370.86611
AB -89.11611 37.36611
AC -217.11611 -90.63389
BC -91.86611 34.61611
ABC -33.36611 93.11611