Para ilustrar este caso se tiene la siguiente tabla de frecuencias:

Tabla_1<-as.table(rbind(c(6,3,7,0),c(1,6,0,8)))
dimnames(Tabla_1) <- list(Mediana=c(">89","<89"),
                     Variedad=c("Variedad 1", "Variedad 2", "Variedad 3", "Variedad 4"))
print(Tabla_1)
##        Variedad
## Mediana Variedad 1 Variedad 2 Variedad 3 Variedad 4
##     >89          6          3          7          0
##     <89          1          6          0          8

Para la elaboración del test aleatorizado es necesario asignarle valores aleatorios a algunas celdas, y el valor del resto dependera de dichas asignaciones, por está razón se le asignaron nombres a cada una de las casillas.

##          Variedad
## Mediana   Variedad 1 Variedad 2 Variedad 3 Variedad 4 Total.C
##   >89     a          b          c          d          16     
##   <89     e          f          g          h          15     
##   Total.F 7          9          7          8          31

Los valores que tomen, las celdas c, d, e, f, g, y h; van a depender de a, b y c. Por lo cual se obtendran valores de a, b y c mediante una distribución uniforme que va desde 0 hasta la suma total de la marginal columna; de la siguiente forma:

por lo cual el resto de valores se obtendrá así:

##          Variedad
## Mediana   Variedad 1 Variedad 2 Variedad 3 Variedad 4         Total.C
##   >89     a          b          c          16-(a+b+c)         16     
##   <89     7-a        9-b        7-c        31-(a+b+c+d+e+f+g) 15     
##   Total.F 7          9          7          8                  31

Ya que es posible que al crear las tablas algunas columnas resulten con valores negativos y no será posible realizar la prueba, se crea una función que va a permitir evaluar cuando esto sucede.

Positiva <- function(tabla){
  tabla2<-tabla
  for (i in 1:2) {
    for (j in 1:4) {
      if (tabla[i,j]>=0){
        tabla2[i,j]<-0
      }
      else{
        tabla2[i,j]<-1
      }
    }
    resul<-sum(tabla2) # esta suma será diferente de cero cuando haya almenos un valor negativo
  }
  return(resul)
}

Ahora se realiza la estimación de los valores de la tabla.

a<-c()
b<-c()
c<-c()

Nsim<-10000
chisq_sim<-c() ## 
for (i in 1:Nsim) { # esto nos va a generar Nsim tablas
  a[i]<-round(runif(1,0,7),0) #Estimación de a
  b[i]<-round(runif(1,0,9),0) #EStimación de b
  c[i]<-round(runif(1,0,7),0) #Estimación de c
  d<-16-(a[i]+b[i]+c[i]) # Marginal fila menos a+b+c
  e<-7-a[i] 
  f<-9-b[i] 
  g<-7-c[i]
  h<-31-(a[i]+b[i]+c[i]+d+e+f+g)
  tabla<-matrix(c(a[i],b[i],c[i],d,e,f,g,h),ncol=4,nrow = 2,TRUE) #Ceación de la tabla
  if(Positiva(tabla)==0){ #Se valida que la tabla no tenga valores <0
    chisq_sim[i]<-chisq.test(tabla)$statistic #se almacena el estadístico para las 1000 tablas
  }
}

Hallamos el valor del estadístico para la tabla_1, sobre el cual se hallará la probabilidad de obtener un valor mayor al obtenido, en pocas palabras el valor-p:

\(Valor-p=P(\chi^2_{aleatorizado}>\chi^2_{calculado})\)

chisq_obs<-chisq.test(Tabla_1)$statistic
chisq_obs
## X-squared 
##  19.55952
Exito<-0 # Se crea el vector que identificará que se cumpla la condición
chisq_final<-na.omit(chisq_sim) #Eliminamos los NA´s estos corresponden a los casos en que no se pudo aplicar la prueba

valor_p <- mean(chisq_final>chisq_obs)
valor_p
## [1] 0.04413377
hist(chisq_final,freq=FALSE,col = "darkviolet", xlab = "Valores X-cuadrado permutación aleatorizada", main = "Distribución de referencia")
abline(v=chisq_obs, col=c("turquoise"), lty=c(8,3),lwd=2)
legend("topright",paste("Valor-P = ",round(valor_p,4)))

length(chisq_final) # número total de pruebas realizadas
## [1] 7296
## [1] "72.96 %"

Aquí se muestra la gráfica de frecuencias para los valores de a, b y c, considerando los casos en que se pudo realizar la prueba y los que no.

# Con NA´s
par(mfrow=c(2,3)) 
plot(table(a),col=4)
plot(table(b),col=2)
plot(table(c),col=3)
### Sin NA´s
plot(table(a[chisq_sim!="NA"]),ylab="a",col=4)
plot(table(b[chisq_sim!="NA"]),ylab="b",col=2)
plot(table(c[chisq_sim!="NA"]),ylab="C",col=3)