#problema Se desea evaluar 4 dosis de nitrógeno en el cultivo de arroz, se usó como respuesta un índice espectral conocido como ndvi. La variabilidad en el suelo obligó al uso de 8 bloques, donde los 4 tratamientos fueron aleatorizados. Realizar un análisis descriptivo e inferencial y concluir cual es la mejor dosis recomendada y el papel del bloqueo en este estudio.
El ndvi mide el verdor y la salud de la vegetación a partir de imágenes satelitales o aéreas. Se calcula a partir de la diferencia entre la reflectancia en las bandas de infrarrojo cercano (NIR) y rojo, dividida por su suma.
rm(list = ls())Se instalan las librerías, se siembra la semilla para que siempre se tengan los mismos datos, se crean las variables de las 4 dosis de nitrogeno 4 y los 8 bloques al azar.
library(ggplot2)
library(dplyr)
library(plotly)
set.seed(123)
dosis=c(sample(gl(4,1,4,c(0,10,20,30)),replace = F),
sample(gl(4,1,4,c(0,10,20,30)),replace = F),
sample(gl(4,1,4,c(0,10,20,30)),replace = F),
sample(gl(4,1,4,c(0,10,20,30)),replace = F))
bloque =gl(8,4,32,paste0("B_",1:8))
xy=expand.grid(x=seq(1,2),y=seq(1,2))
datos=data.frame(dosis,xy,bloque)Esta parte del código sirve para mostrar gráficamente como quedarían repartidos los 4 tratamientos en los 8 bloques
q=datos %>%
ggplot(aes(x=x,y = y,fill = dosis,label=dosis)) +
geom_tile()+
geom_text()+
facet_wrap(~bloque,nrow = 1)+
theme_bw()
ggplotly(q)Aquí se añadió la variable del ndvi dependiendo de la cantidad de nitrógeno
set.seed(123)
ndvi_0=runif(8,0.4,0.45)#n, min, max
ndvi_10=runif(8,0.45,0.55)
ndvi_20=runif(8,0.5,0.58)
ndvi_30=runif(8,0.58,0.65)
datos$ndvi[datos$dosis==0]=ndvi_0
datos$ndvi[datos$dosis==10]=ndvi_10
datos$ndvi[datos$dosis==20]=ndvi_20
datos$ndvi[datos$dosis==30]=ndvi_30
head(datos)# Graficar
q2=datos %>%
ggplot(aes(x = dosis,y = ndvi,fill = dosis))+
geom_point()+
facet_wrap(~bloque,nrow=1)+
theme(legend.position = 'none')+
theme_bw()
ggplotly(q2)q3=datos %>%
group_by(bloque) %>%
summarise(mean_ndvi=mean(ndvi)) %>%
ggplot(aes(x=bloque,y=mean_ndvi,colour=bloque))+
geom_point(size=8)+
scale_y_continuous(limits = c(0,1))+
theme_bw()
ggplotly(q3)Tabular estadísticas
datos %>%
group_by(dosis) %>%
summarise(mean_ndvi=mean(ndvi), # Promedio del NDVI por dosis
sd_ndvi=sd(ndvi),# Desviación estándar del NDVI por dosis
var_ndvi=var(ndvi), # Varianza del NDVI por dosis
cv_ndvi=sd_ndvi/mean_ndvi*100)# Coeficiente de variación(%)\[H_o:\mu_{D1}=\mu_{D2}=\mu_{D3}=\mu_{D4}\\H_a:H_o\text{ Es Falsa }\]
Significa que: Hipótesis nula (H0):Los promedios del NDVI (u otra variable de interés) son iguales entre todos los tratamientos de dosis. Hipótesis alternativa (Ha):Al menos uno de los tratamientos tiene un promedio significativamente diferente.
# Análisis de varianza
datos$dosis=as.factor(datos$dosis)
mod=aov(ndvi~bloque+dosis,data=datos)
s_mod=summary(mod)
s_mod## Df Sum Sq Mean Sq F value Pr(>F)
## bloque 7 0.00456 0.00065 1.267 0.313
## dosis 3 0.15403 0.05134 99.763 1.38e-12 ***
## Residuals 21 0.01081 0.00051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p valor de la dosis es mayor a 0,05 por lo tanto se rechaza la hipótesis nula (H0)
mod$residuals## 1 2 3 4 5 6
## -0.017431576 0.014622363 -0.004702967 0.007512180 -0.034486288 -0.002702304
## 7 8 9 10 11 12
## 0.017587634 0.019600958 -0.022220351 -0.009968411 0.036716599 -0.004527836
## 13 14 15 16 17 18
## -0.020379335 0.006987598 -0.007768531 0.021160268 0.016275563 -0.028804850
## 19 20 21 22 23 24
## 0.010172533 0.002356754 0.020767030 0.012093371 -0.018517736 -0.014342666
## 25 26 27 28 29 30
## 0.009061655 0.002263136 -0.042399149 0.031074357 0.006801885 -0.010010181
## 31 32
## -0.003665401 0.006873697Es un diseño factorial simple en bloques al azar
shapiro.test(mod$residuals)##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.98649, p-value = 0.9511Existe evidencia estadística con un nivel de significancia de 0.05 para no rechazar la hipótesis nula, por lo tanto los residuos del modelo presentan una distribución normal.
bartlett.test(mod$residuals, bloque)##
## Bartlett test of homogeneity of variances
##
## data: mod$residuals and bloque
## Bartlett's K-squared = 4.9709, df = 7, p-value = 0.6635Existe evidencia estadística con un nivel de significancia de 0.05 para no rechazar la hipótesis nula, por lo tanto las varianzas de los residuos son homogéneas entre los bloques.
datos %>%
group_by(dosis) %>%
summarise(mean_ndvi=mean(ndvi))Se observa una tendencia creciente en los valores medios de NDVI conforme aumenta la dosis del tratamiento. La dosis de 30 unidades presenta el mayor valor medio de NDVI (0.6220), siendo 0.1922 unidades superior al valor observado en la dosis 0 (0.4298).
Esta diferencia sugiere que la dosis 30 induce una mayor respuesta en el NDVI y, por tanto, podría considerarse como la más recomendable entre las evaluadas.
mean_total = mean(datos$ndvi)
mean_total## [1] 0.5270355medias = datos%>%
group_by(dosis) %>%
summarise(mean_ndvi = mean(ndvi))
mediasmedias = medias%>%
mutate(efecto=mean_ndvi - mean_total)
mediasLa dosis de 30 de nitrógeno es el mejor tratamiento en términos de efecto sobre el ndvi, ya que incrementó el ndvi aproximadamente 0.095 unidades por encima de promedio, indicando más vigor vegetal o cobertura verde en el cultivo de arroz.
# Análisis de varianza
datos$dosis=as.factor(datos$dosis)
mod=aov(ndvi~bloque+dosis,data=datos)
s_mod=summary(mod)
s_mod## Df Sum Sq Mean Sq F value Pr(>F)
## bloque 7 0.00456 0.00065 1.267 0.313
## dosis 3 0.15403 0.05134 99.763 1.38e-12 ***
## Residuals 21 0.01081 0.00051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p-valor es 0.313, que es mayor a 0.05 indica que el bloque no tiene un efecto significativo sobre el ndvi.