Ejercicio cuadrado de youden
Lady Diana Arias
3/8/2021
library(collapsibleTree)
b<-8 #columnas
k<-7 #filas
row1 <- factor(rep(1:k, times <- b))
col1 <- factor(c(rep(1,k),rep(2,k),rep(3,k),rep(4,k),rep(5,k),rep(6,k),rep(7,k),rep(8,k)))
treat<-factor(c("D","H","C","B","E","A","G",
"F","E","G","A","H","B","C",
"B","C","H","D","G","F","E",
"A","G","E","F","C","D","H",
"C","B","D","H","A","E","F",
"E","F","A","G","D","C","B",
"G","A","F","E","B","H","D",
"H","D","B","C","F","G","A"))
y<-c(16.6,16.9,17.4,17.4,15.8,18.2,15.7,
15.9,16.4,15.8,19.0,17.6,17.8,18.9,
17.1,16.8,19.2,16.6,15.8,17.8,18.4,
17.7,15.9,16.3,16.0,17.6,17.8,18.1,
17.4,17.0,16.8,19.2,20.3,18.4,15.9,
16.5,16.0,16.9,15.9,17.1,17.5,17.4,
15.8,16.9,15.9,16.5,17.6,19.4,17.1,
18.6,17.4,17.4,19.2,16.8,15.7,17.4)
dat2<-data.frame(row1, col1, treat, y)library(collapsibleTree)
collapsibleTreeSummary(dat2, hierarchy = c('col1', 'row1', 'treat','y'))N<-xtabs(~treat + col1, data <- dat2);N## col1
## treat 1 2 3 4 5 6 7 8
## A 1 1 0 1 1 1 1 1
## B 1 1 1 0 1 1 1 1
## C 1 1 1 1 1 1 0 1
## D 1 0 1 1 1 1 1 1
## E 1 1 1 1 1 1 1 0
## F 0 1 1 1 1 1 1 1
## G 1 1 1 1 0 1 1 1
## H 1 1 1 1 1 0 1 1v<-nlevels(treat);
n<-b*k;
T1<-aggregate(dat2$y, by<-list(treat<-dat2$treat),FUN<-sum);
Tt<-matrix(T1$x,v,1);# totales de Trat
T2<-aggregate(dat2$y, by<-list(row1<-dat2$row1),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
T3<-aggregate(dat2$y, by<-list(col1<-dat2$col1),FUN<-sum);
Tc<-matrix(T3$x,b,1)#totales Columna
i<-1:k
j<-1:b
if((b==k)&&(N[i,j]==1)){
print("dado el diseño cuadrado latino");
#si es un DCL uno corre el código antes considerado.
}else{
x<-0;z<-0;
if(v>=2 && k<v)
{
A<-N%*%t(N);
for(i in 2:v){
x1<-A[1,1]-A[i,i];
x<-x+abs(x1);
}
sum1<-sum(x);sum1
for(i in 1:v){
for(j in 1:v){
if(i!=j){
x2<-A[1,2]-A[i,j];
z<-z+abs(x2)}
}}
sum2<-sum(z);
if(sum1==0 && sum2==0){
r<-A[1,1];
lambda<-A[1,2];
cat("dado el diseño cuadrado de Youden \n con r y lambda respectivamente")
print(r);
print(lambda);
#se define el vector columna de k entradas todas en 1
ek<-c(matrix(rep(1,k),k,1));
Dk<-diag(k)-(ek%*%t(ek))/v;
Q<-Tt-(N%*%Tc)/k;
alpha_hat<-(k/(lambda*v))*Q;
beta_hat<-(Tc/k)-(t(N)%*%Q)/(lambda*v);
gamma_hat<-Tr/v;
print("Las estimaciones son");
print("mu_hat");
print("0");
print("alpha_hat");
print(alpha_hat);
print("beta_hat")
print(beta_hat)
print("gamma_hat")
print(gamma_hat);
#probando hipótesis H_alpha, H_beta y H_gamma
myfit <- lm(y~treat+row1+col1, dat2);
anova(myfit);
}}}## dado el diseño cuadrado de Youden
## con r y lambda respectivamente[1] 7
## [1] 6
## [1] "Las estimaciones son"
## [1] "mu_hat"
## [1] "0"
## [1] "alpha_hat"
##
## treat [,1]
## A 0.8750000
## B 0.1416667
## C 0.5895833
## D -0.1520833
## E -0.2895833
## F -0.9666667
## G -1.3604167
## H 1.1625000
## [1] "beta_hat"
##
## col1 [,1]
## 1 16.71905
## 2 17.32113
## 3 17.51071
## 4 17.07738
## 5 17.66280
## 6 16.92321
## 7 17.11280
## 8 17.45863
## [1] "gamma_hat"
## [,1]
## [1,] 16.9500
## [2,] 16.6625
## [3,] 16.9625
## [4,] 17.4750
## [5,] 17.3250
## [6,] 17.8250
## [7,] 17.3625## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 7 38.494 5.4992 10.5175 4.866e-07 ***
## row1 6 7.299 1.2164 2.3265 0.05386 .
## col1 7 4.927 0.7039 1.3462 0.25867
## Residuals 35 18.300 0.5229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Comparación de los tratamientos cuando el diseño es de tipo youden (modelo diferente)
El efecto de los tratamientos: tienen mejor respuesta los tratamientos A, B, H y c Segun el p-Valor de los tratamientos, se rechaza la H_o ya que el valor es <0.05
#Analisis descriptivo
med_trat=tapply(dat2$y, dat2$treat, mean)
barplot(med_trat, ylim = c(0,20), main = 'Distribucion de la respuesta por tratamiento')
abline(h=mean(dat2$y))
Al ver el grafico podemos decir que 4 de los tratamientos tienen efecto negativo y 4 de los tratamientos tienen efecto positivo
med_trat=tapply(dat2$y, dat2$row1, mean)
barplot(med_trat, ylim = c(0,20), main = 'Distribucion de la respuesta por bloque fila')
abline(h=mean(dat2$y))
med_trat=tapply(dat2$y, dat2$col1, mean)
barplot(med_trat, ylim = c(0,20), main = 'Distribucion de la respuesta por bloque columna')
abline(h=mean(dat2$y))