library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Asignamos los valores:
ABC <- 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)
Convertimos el vector en matriz:
ABC <- matrix(ABC, ncol = 3, byrow = TRUE)
ABC
## [,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 1
Asignamos estos valores al vector de replica y creamos el dataframe:
replica <- c(247,470,429,435,837,551,775,660,400,446,405,445,850,670,865,530)
plasma <- data.frame(cbind(A = ABC[,1], B = ABC[,2], C = ABC[,3], replicate = replica))
head(plasma)
Asignamos totales de las replicas:
totales <- c(647,916,834,880,1687,1221,1640,1190)
totales <- matrix(totales)
rownames(totales) <- c("(1)","a","b","ab","c","ac","bc","abc")
totales
## [,1]
## (1) 647
## a 916
## b 834
## ab 880
## c 1687
## ac 1221
## bc 1640
## abc 1190
\(A = \frac{1}{4n}[a + ab + ac + abc - (1) - b - c - bc]\)
rep_num <- 2
Contraste_A <- totales["a",] + totales["ab",] + totales["ac",] + totales["abc",] - totales["(1)",] - totales["b",] - totales["c",] - totales["bc",]
Contraste_A
## a
## -601
A_est_effect <- (1/(4*rep_num)) * Contraste_A
A_est_effect
## a
## -75.125
\(B = \frac{1}{4n}[b + ab + bc + abc - (1) - a - c - ac]\)
Contraste_B <- totales["b",] + totales["ab",] + totales["bc",] + totales["abc",] - totales["(1)",] - totales["a",] - totales["c",] - totales["ac",]
Contraste_B
## b
## 73
B_est_effect <- (1/(4*rep_num)) * Contraste_B
B_est_effect
## b
## 9.125
\(C = \frac{1}{4n}[c + ac + bc + abc - (1) - a - c - ab]\)
Contraste_C <- totales["c",] + totales["ac",] + totales["bc",] + totales["abc",] - totales["(1)",] - totales["a",] - totales["b",] - totales["ab",]
Contraste_C
## c
## 2461
C_est_effect <- (1/(4*rep_num)) * Contraste_C
C_est_effect
## c
## 307.625
\(AB = \frac{1}{4n}[abc + ab + c + (1) - bc - b - ac - a]\)
Contraste_AB <- totales["abc",] + totales["ab",] + totales["c",] + totales["(1)",] - totales["bc",] - totales["b",] - totales["ac",] - totales["a",]
Contraste_AB
## abc
## -207
AB_est_effect <- (1/(4*rep_num)) * Contraste_AB
AB_est_effect
## abc
## -25.875
\(AC = \frac{1}{4n}[(1) + b + ac + abc - a - ab - c - bc]\)
Contraste_AC <- totales["(1)",] + totales["b",] + totales["ac",] + totales["abc",] - totales["a",] - totales["ab",] - totales["c",] - totales["bc",]
Contraste_AC
## (1)
## -1231
AC_est_effect <- (1/(4*rep_num)) * Contraste_AC
AC_est_effect
## (1)
## -153.875
\(BC = \frac{1}{4n}[(1) + a + bc + abc - b - ab - c - ac]\)
Contraste_BC <- totales["(1)",] + totales["a",] + totales["bc",] + totales["abc",] - totales["b",] - totales["ab",] - totales["c",] - totales["ac",]
Contraste_BC
## (1)
## -229
BC_est_effect <- (1/(4*rep_num)) * Contraste_BC
BC_est_effect
## (1)
## -28.625
\(ABC = \frac{1}{4n}[abc + c + b + a - bc - ac - ab - (1)]\)
Contraste_ABC <- totales["abc",] + totales["c",] + totales["b",] + totales["a",] - totales["bc",] - totales["ac",] - totales["ab",] - totales["(1)",]
Contraste_ABC
## abc
## 239
ABC_est_effect <- (1/(4*rep_num)) * Contraste_ABC
ABC_est_effect
## abc
## 29.875
Calculamos la suma de cuadrados por iteracción:
\(SS = \frac{(contrast)^2}{8n}\)
SS_A <- (Contraste_A^2)/(8*rep_num)
SS_B <- (Contraste_B^2)/(8*rep_num)
SS_C <- (Contraste_C^2)/(8*rep_num)
SS_AB <- (Contraste_AB^2)/(8*rep_num)
SS_AC <- (Contraste_AC^2)/(8*rep_num)
SS_BC <- (Contraste_BC^2)/(8*rep_num)
SS_ABC <- (Contraste_ABC^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)
## 22575.0625 333.0625 378532.5625 2678.0625 94710.0625 3277.5625
## abc
## 3570.0625
Analisis de Varianza:
fit_for_AOV <- lm(replicate ~ A + B + C + (A*B) + (A*C) + (B*C) + (A*B*C), data = plasma)
anova(fit_for_AOV)
Variables significativas A, C y AC, se procede hacer modelo de regresión lineal:
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-07
Coeficientes:
fit_plasma$coefficients
## (Intercept) A C A:C
## 563.4375 -37.5625 153.8125 -76.9375
Resultado:
\(563.43 - 37.56x_1 + 153.81x_3 - 76.93x_1x_3\)
\(se(Effect) =\sqrt{\frac{1}{n2^{k-2}}s^2}\)
MSE <- anova(fit_for_AOV)$`Mean Sq`[8]
MSE
## [1] 3999.438
se_effect <- sqrt(1/(rep_num*(2^(3-2)))*(MSE))
se_effect
## [1] 31.62055
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