Donde ~ nos indica vectores, y mayusculas nos indican matrices. Es decir que \( \widetilde{y} \) es el vector de respuesta, \( \widetilde{\beta} \) el vector de coeficientes del modelo, \( \widetilde{\varepsilon} \) el vector de errores aleatorios.
Entonces para este caso nos interesa \( \mathbf{X} \), i.e. la matriz de diseño.
Suponemos un diseño con 5 niveles en cada factor, y una sola réplica por tratamiento.Lo que corresponde a una matriz diseño de 25 renglones y 16 columnas.
clat <- matrix(nrow = 25, ncol = 1 + 5 * 3)
clat
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA NA NA NA NA NA NA NA
## [,14] [,15] [,16]
## [1,] NA NA NA
## [2,] NA NA NA
## [3,] NA NA NA
## [4,] NA NA NA
## [5,] NA NA NA
## [6,] NA NA NA
## [7,] NA NA NA
## [8,] NA NA NA
## [9,] NA NA NA
## [10,] NA NA NA
## [11,] NA NA NA
## [12,] NA NA NA
## [13,] NA NA NA
## [14,] NA NA NA
## [15,] NA NA NA
## [16,] NA NA NA
## [17,] NA NA NA
## [18,] NA NA NA
## [19,] NA NA NA
## [20,] NA NA NA
## [21,] NA NA NA
## [22,] NA NA NA
## [23,] NA NA NA
## [24,] NA NA NA
## [25,] NA NA NA
A la primer columna se le introducen 1, correspondientes a la media general
clat[, 1] <- 1
Las siguentes cinco columnas corresponden a los cinco niveles del primer factor de bloqueo, así a los primeros cinco renglones de la columna 2 se le introducen 1 y al resto ceros, en la columna 3 los renglones del 6 al 10 se introducen 1 y al resto ceros, y se continua así con las columnas siguientes hasta la 6.
clat[, 2] <- c(rep(1, 5), rep(0, 20))
clat[, 3] <- c(rep(0, 5), rep(1, 5), rep(0, 15))
clat[, 4] <- c(rep(0, 10), rep(1, 5), rep(0, 10))
clat[, 5] <- c(rep(0, 15), rep(1, 5), rep(0, 5))
clat[, 6] <- c(rep(0, 20), rep(1, 5), rep(0, 0))
Las siguientes cinco columnas, de la 7 a la 11, corresponden a el segundo factor de bloqueo, para este caso el primer renglon del primer nivel del primer factor se le agrega un 1 para la primer columna, correspondiente al primer nivel del segundo factor. La siguiente columna corresponde a el segundo nivel del segundo factor, y en este caso se agrega un 1 al segundo renglon de cada nivel del primer factor, asi sucesivamente hasta la columna 11.
clat[, 7] <- rep(c(1, rep(0, 4)), 5)
clat[, 8] <- rep(c(rep(0, 1), 1, rep(0, 3)), 5)
clat[, 9] <- rep(c(rep(0, 2), 1, rep(0, 2)), 5)
clat[, 10] <- rep(c(rep(0, 3), 1, rep(0, 1)), 5)
clat[, 11] <- rep(c(rep(0, 4), 1), 5)
Por último nos restan las columnas correspondientes al tercer factor,i.e. el factor de interes, en este caso consideramos heurísticamente mejor considerar cada conjunto de renglones correspondientes al primer factor de cada nivel como una matriz independiente, y utilizamos matrices identidad con renglones intercambiados para lograr que por renglon no se repitiera alguna de las cobinaciones lineales, y de esta manera obtener la matriz de diseño para cuadros latinos de 5x5x5.
clat[1:5, 12:16] <- diag(5)
m <- diag(5)
m[1:4, ] <- diag(5)[2:5, ]
m[5, ] <- diag(5)[1, ]
clat[6:10, 12:16] <- m
m <- diag(5)
m[1:3, ] <- diag(5)[3:5, ]
m[4:5, ] <- diag(5)[1:2, ]
clat[11:15, 12:16] <- m
m <- diag(5)
m[1:2, ] <- diag(5)[4:5, ]
m[3:5, ] <- diag(5)[1:3, ]
clat[16:20, 12:16] <- m
m <- diag(5)
m[1, ] <- diag(5)[5, ]
m[2:5, ] <- diag(5)[1:4, ]
clat[21:25, 12:16] <- m
Obtenemos así la matriz de diseño
clat
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 1 0 0 0 0 1 0 0 0 0 1 0
## [2,] 1 1 0 0 0 0 0 1 0 0 0 0 1
## [3,] 1 1 0 0 0 0 0 0 1 0 0 0 0
## [4,] 1 1 0 0 0 0 0 0 0 1 0 0 0
## [5,] 1 1 0 0 0 0 0 0 0 0 1 0 0
## [6,] 1 0 1 0 0 0 1 0 0 0 0 0 1
## [7,] 1 0 1 0 0 0 0 1 0 0 0 0 0
## [8,] 1 0 1 0 0 0 0 0 1 0 0 0 0
## [9,] 1 0 1 0 0 0 0 0 0 1 0 0 0
## [10,] 1 0 1 0 0 0 0 0 0 0 1 1 0
## [11,] 1 0 0 1 0 0 1 0 0 0 0 0 0
## [12,] 1 0 0 1 0 0 0 1 0 0 0 0 0
## [13,] 1 0 0 1 0 0 0 0 1 0 0 0 0
## [14,] 1 0 0 1 0 0 0 0 0 1 0 1 0
## [15,] 1 0 0 1 0 0 0 0 0 0 1 0 1
## [16,] 1 0 0 0 1 0 1 0 0 0 0 0 0
## [17,] 1 0 0 0 1 0 0 1 0 0 0 0 0
## [18,] 1 0 0 0 1 0 0 0 1 0 0 1 0
## [19,] 1 0 0 0 1 0 0 0 0 1 0 0 1
## [20,] 1 0 0 0 1 0 0 0 0 0 1 0 0
## [21,] 1 0 0 0 0 1 1 0 0 0 0 0 0
## [22,] 1 0 0 0 0 1 0 1 0 0 0 1 0
## [23,] 1 0 0 0 0 1 0 0 1 0 0 0 1
## [24,] 1 0 0 0 0 1 0 0 0 1 0 0 0
## [25,] 1 0 0 0 0 1 0 0 0 0 1 0 0
## [,14] [,15] [,16]
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 1 0 0
## [4,] 0 1 0
## [5,] 0 0 1
## [6,] 0 0 0
## [7,] 1 0 0
## [8,] 0 1 0
## [9,] 0 0 1
## [10,] 0 0 0
## [11,] 1 0 0
## [12,] 0 1 0
## [13,] 0 0 1
## [14,] 0 0 0
## [15,] 0 0 0
## [16,] 0 1 0
## [17,] 0 0 1
## [18,] 0 0 0
## [19,] 0 0 0
## [20,] 1 0 0
## [21,] 0 0 1
## [22,] 0 0 0
## [23,] 0 0 0
## [24,] 1 0 0
## [25,] 0 1 0
Para este caso obtamos por una función que realice la operación requerida.
#llamaremos a la funcion "mincuadrados"
mincuadrados<-function(X,y){ # requiere una matriz o un data.frame "X", que incluya las variables explicatorias; y un vector "y" con la variable explicada.
X<-as.matrix(X) #en caso de una data.frame este paso lo convierte en matriz
X_1<-matrix(nrow=nrow(X),ncol=ncol(X)+1) #creamos una nueva matriz con columnas p+1
X_1[,1]<-1 #a la primera columna de la nueva matriz le imputamos 1 para obtener el primer sumando del modelo (la media o intercepto)
X_1[,2:ncol(X_1)]<-X # al resto de la nueva matriz introducimos las variables explicatorias
t(solve(t(X_1)%*%X_1)%*%t(X_1)%*%y)} #ahora aplicamos la fórmula para obener el vector de coeficientes
Ahora falta verificar que efectivamente la funcion funcione, para esto simularemos algunos datos, 25 casos con 4 variables de interes, todas ellas normales:
set.seed(1984)
a <- rnorm(25, 20, 2)
B <- matrix(nrow = 25, ncol = 4)
set.seed(1)
B[, 1] <- rnorm(25, 5, 2) * a
set.seed(2)
B[, 2] <- rnorm(25, 9, 2) * a
set.seed(3)
B[, 3] <- rnorm(25, 8, 1) * a
set.seed(4)
B[, 4] <- rnorm(25, 7, 1) * a
Aplicamos ahora la función a los datos simulados
mincuadrados(B, a)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 8.372 -0.00481 -0.001084 0.03952 0.0425
lm(a ~ B)
##
## Call:
## lm(formula = a ~ B)
##
## Coefficients:
## (Intercept) B1 B2 B3 B4
## 8.37185 -0.00481 -0.00108 0.03952 0.04250
Aparentemente coincide, pero solo para estar seguros:
mincuadrados(B, a) == lm(a ~ B)$coefficients
## [,1] [,2] [,3] [,4] [,5]
## [1,] FALSE FALSE FALSE FALSE FALSE
Pues resulta que R dice que los vectores no son iguales, culparemos a alguna diferencia en decimales, tras verificar que efectivamente son iguales ambas estimaciones.
Cabe mencionar que la función “mincuadrados” debe funcionar para cualquier número de casos n y p variables explicatorias, pero eso ya lo verificaría alguien más.