Clase 3 Multivariado, factores y bloques
Javier Borbon
2024-02-29
#quedamos en taxonomia de suelos, recordemos que llamamos factores a las variables que controlo, ahora como las aleatorizo en campo, esto tiene que ver con lo el medio, si necesito bloquear, son bloques
#uno monta el ensayo para evaluar el efecto de los factores, lo que controlo
#factorail simple, compeltamente al azar
set.seed(2022)
rend1 = rnorm(48, 3, 0.3)
variedad = gl ( n = 3, k = 16, length = 48, #16 X 3 =48 GL es generador de niveles
labels = c( "V1", "V2", "V3")) #tres variedades, son mis tres factores
df = data.frame(variedad, rend1)
#GL no es grados de libertad
head(dt) #me da los primeros datos de rendimiento##
## 1 function (x, df, ncp, log = FALSE)
## 2 {
## 3 if (missing(ncp))
## 4 .Call(C_dt, x, df, log)
## 5 else .Call(C_dnt, x, df, ncp, log)
## 6 }library(collapsibleTree)
collapsibleTree(df, hierarchy = c("variedad", "rend1")) #pongo mi dataframe y le pongo la jerarquia#experimento balanceado, tengo el mismo numero de unidades en el trt, aca se decidio que para colocar 3 trt no hay problema, y los aleatorizo, pero si tengo problemas, seguramente debería bloquear
#estas 48 parcelas las puedo colocar como quiera#colocar mis parcelas con los diferentes genotipos
#expandgrid es para ahcer cuadricula
data = expand.grid(x = seq(0, 210, 30),
y = seq(0, 210, 40))
data$var = sample(1:3, 48, replace = T)
plot(data$x, data$y, pch = 8, col = data$var, cex =2)
grid(8,6)
#otra forma de mostrar la cuadricula
library(ggplot2)
#se le mete la data, coord x y y, si no tiene coord no importa, tiene fila 1, 2, 3...
#
ggplot(data, aes( x = data$x, y = data$y,
fill = as.factor(data$var)))+
geom_tile(color = "black")+
geom_text(aes(x = x, y = y,
label = data$var))
#vamos a hacer un ejemplo 2 donde no no requiera hacer aleatorizacion espacial, como en publicidad, mi variable respuesta es ver si la persona detecto uno de los puntos de la publicidad, dicotomia si y no
#ditocotmia si y no
#detección, detectó tal cosa
set.seed(123)
df = data.frame(detec = round(runif(n = 48, min = 0, max = 0.6)))
df$detec## [1] 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0
## [39] 0 0 0 0 0 0 0 0 0 0#depende,la ambiguedad no da peso.
#respuestas en 0 y 1#aca le voy a poner mas variables, para confiar mas en la respuesta, puse a las personas a diferentes segundos de publicidad
set.seed(123)
df2 = data.frame(detec = round(runif(n = 48, min = 0, max = 0.6)),
t_exp = gl(n = 3, k = 16, length = 48,
labels = c("30s", "50s", "70s")), persona = seq(1:48)) #30, 50 y 70 segundos de publicidad
df2$detec## [1] 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0
## [39] 0 0 0 0 0 0 0 0 0 0head(df2)## detec t_exp persona
## 1 0 30s 1
## 2 0 30s 2
## 3 0 30s 3
## 4 1 30s 4
## 5 1 30s 5
## 6 0 30s 6#anteriormente lo tengo todo ordenado, vamos a aleatorizar con sample
#lo ultimo, simplifica escritura de codigo, desarrollada por asesores de amazon Wickham
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondf3 = df2 %>%
sample_frac()
#Pipe, es un tubo, programacion tubular, con la salida de mi conjunto de datos, lo meto en otro analisis
#simple_frac es para organizarlos
#un conjunto de datos, de ahi saco ejemplo la media, y de esas medias, saco otra media...
#no todo lo de R esta bien, como tirar un Rcuadrado a cualquier modelo, siendo que R cuadrado no aplica en modelos no lineales
df3## detec t_exp persona
## 1 0 50s 23
## 2 0 50s 27
## 3 0 30s 7
## 4 0 70s 47
## 5 1 50s 32
## 6 0 70s 38
## 7 0 50s 25
## 8 0 70s 34
## 9 0 50s 29
## 10 1 30s 5
## 11 1 30s 8
## 12 0 30s 12
## 13 0 30s 13
## 14 0 50s 18
## 15 0 70s 33
## 16 0 70s 45
## 17 0 70s 42
## 18 0 30s 6
## 19 1 50s 21
## 20 0 30s 15
## 21 0 30s 9
## 22 0 70s 40
## 23 0 50s 26
## 24 1 30s 16
## 25 1 50s 20
## 26 1 50s 31
## 27 1 30s 11
## 28 0 70s 43
## 29 0 70s 46
## 30 0 70s 44
## 31 0 50s 17
## 32 0 70s 35
## 33 0 30s 2
## 34 1 30s 4
## 35 0 70s 36
## 36 0 70s 39
## 37 0 70s 48
## 38 0 30s 3
## 39 0 50s 28
## 40 0 70s 41
## 41 0 70s 37
## 42 0 30s 1
## 43 0 50s 30
## 44 1 50s 24
## 45 0 50s 22
## 46 0 30s 10
## 47 0 50s 19
## 48 0 30s 14#tengo grupos de personas con diferentes exposiciones a publicidad y respuestas no y si, 0 y 1#hay analisis de variables por separado: chirimoya, por dulzor, por aparte de la acidez, diferente al color, como si una variable no estuviera afectada por la otra.
#el sexo puede cambiar la percepcion en una encuesta
#FACTORIAL SIMPLE EN BLOQUES
#lo que miro, va a tener efecto de otra variable, que debo tenerla en cuenta para representar mas cercanamente la realidad
#Estrato vs bloque, es similar en publicidad, o factor de clasificacion, o factor experimental como la variedad
set.seed(2022)
rend2 = rnorm(48, 3, 0.3)
fertili = gl ( n = 3, k = 16, length = 48,
labels = c( "F0", "F10", "F20")) #dosis de fertilizante
df5 = data.frame(fertili, rend2)
# agregando variable de bloqueo
#estoy bloqueando en el invernadero, dado que el techo de plastico esta diferente: zona anterior y posterior del invernadero
#este ensayo se hizo para comparar la fertilizacion, no para ver el estado del plastico, no me interesa el plastico
#primera razon para bloquear
df5$bloque = gl(2,4,48,
c("anterior", "posterior"))
df5$x = data$x
df5$y = data$y
#segunda razon para bloquear, mañana y tarde
df5$turno = gl(2, 8, 48, c("mañana", "tarde"))
#hay 3 trt, tres factores, aunque se ve por bloque
collapsibleTree(df5, hierarchy = c("bloque", "turno", "fertili", "rend2"))#diseño, factorial simple en bloques,
#cuadrado latino, cuando se tiene doble razon para bloquear: bloque de turno y disposicion, mañana y tarde
#cuadrado de Jouden, tengo doble nivel de bloques?#Bloque dañado, hay un 8 en lugar de un 4
set.seed(2022)
rend2 = rnorm(48, 3, 0.3)
fertili = gl(n = 3, k = 16, length = 48,
labels = c("F0", "F10", "F20")) #dosis de fertilizante
df5 = data.frame(fertili, rend2)
# agregando variable de bloqueo
df5$bloque = gl(2, 8, 48, c("anterior", "posterior")) ##aca como tengo un 4 en lugar de un 4, me cambió que en el arbol no apareciera el bloque
# Suponiendo que data$x y data$y están definidos en algún lugar
df5$x = data$x
df5$y = data$y
# segunda razón para bloquear, mañana y tarde
df5$turno = gl(2, 8, 48, c("mañana", "tarde"))
# hay 3 trt, tres factores, aunque se ve por bloque
collapsibleTree(df5, hierarchy = c("bloque", "turno", "fertili", "rend2"))#cuando cada bloque tiene un resultado diferente, debo analizarlos por separado?
#si se analiza por separado, y se evalua el efecto del bloque si es positivo o negativo? una correlación?#FACTORIAL COMPLETAMENTE AL AZAR
set.seed(2022)
sineresis = rnorm(36, 3, 0.3) #lactosuero
fuente = gl(2, 18, 36,
c("cabra", "vaca"))
tiempo_exposicion = gl(2, 9, 36,
c("Temp12horas", "temp24horas"))
#bloqueo_laboratorio = gl(,), 36, C()) terminar con datos del DRIVE
df6 = data.frame(fuente, tiempo_exposicion, sineresis)
#Hay factores, 12 y 24 horas, completamente probados, significa que todos los factores del nivele anterior tienen los mismos factores
#si tengo una fuente -Invernadero adelante y atras, y tengo 12 y 24 horas en ambos lugares, es factorial completo
#factorial incompleto, cuando los niveles no son iguales, tengo invernadero y 12 y 24 horas, pero en algun lugar no pude medir 12 y 24
head(df6)## fuente tiempo_exposicion sineresis
## 1 cabra Temp12horas 3.270043
## 2 cabra Temp12horas 2.647996
## 3 cabra Temp12horas 2.730754
## 4 cabra Temp12horas 2.566650
## 5 cabra Temp12horas 2.900696
## 6 cabra Temp12horas 2.129811#sigue el arbolclorofila = rnorm(36, 8, 0.4)
ferti = gl(3, 12, 36, c("gallinaza", "pollinaza", "mixta")) #primer factor
dosis = gl(3, 4, 36, c("D0", "D5", "D10"))#segundo factor, dosis
df7 = data.frame(clorofila, ferti, dosis)
collapsibleTree(df7, c("ferti", "dosis", "clorofila"))#si hubiese dicho que tengo dosis de urea... si este fuera mi segundo factor, el factorial seria completo o incompleto? completo, xk son para los mismos
#Aca en este arbol, es ICOMPLETO, la dosis depende del factor anterior (EN EL MIXTO)
#si se trata del mismo producto de la dosis anterior, se vuelve incompleto###PARCELAS DIVIDIDAS, a menudo involucra maquinaria o drones. F1 3 genotipos f2 fertilizacion, dos tipos., aerea y terrestre (suelo) rendimiento de racimos
df8 = data.frame(F1=gl(3,24,72, c("G1", "G2", "G3")),
F2=gl(2,12,72, c("Fa", "Fs")),
rto = rnorm(72,3,0.3))
df8## F1 F2 rto
## 1 G1 Fa 3.130961
## 2 G1 Fa 3.019929
## 3 G1 Fa 3.391398
## 4 G1 Fa 2.937222
## 5 G1 Fa 3.305476
## 6 G1 Fa 3.409810
## 7 G1 Fa 3.443009
## 8 G1 Fa 3.266192
## 9 G1 Fa 2.695227
## 10 G1 Fa 3.561261
## 11 G1 Fa 3.322860
## 12 G1 Fa 2.677678
## 13 G1 Fs 2.341327
## 14 G1 Fs 3.160363
## 15 G1 Fs 3.403120
## 16 G1 Fs 3.415501
## 17 G1 Fs 3.824078
## 18 G1 Fs 2.986217
## 19 G1 Fs 3.222926
## 20 G1 Fs 3.078127
## 21 G1 Fs 3.128485
## 22 G1 Fs 2.889522
## 23 G1 Fs 3.866227
## 24 G1 Fs 2.817878
## 25 G2 Fa 2.436240
## 26 G2 Fa 3.215507
## 27 G2 Fa 3.075428
## 28 G2 Fa 3.140107
## 29 G2 Fa 2.934633
## 30 G2 Fa 2.634706
## 31 G2 Fa 2.773776
## 32 G2 Fa 3.057369
## 33 G2 Fa 3.314185
## 34 G2 Fa 3.326949
## 35 G2 Fa 3.208642
## 36 G2 Fa 2.842689
## 37 G2 Fs 3.079205
## 38 G2 Fs 3.372360
## 39 G2 Fs 3.192989
## 40 G2 Fs 2.700971
## 41 G2 Fs 2.399310
## 42 G2 Fs 2.857248
## 43 G2 Fs 3.013254
## 44 G2 Fs 2.585659
## 45 G2 Fs 3.074256
## 46 G2 Fs 3.067257
## 47 G2 Fs 3.258444
## 48 G2 Fs 2.721971
## 49 G3 Fa 2.903867
## 50 G3 Fa 3.201907
## 51 G3 Fa 2.620311
## 52 G3 Fa 2.844444
## 53 G3 Fa 3.336748
## 54 G3 Fa 2.905207
## 55 G3 Fa 2.472510
## 56 G3 Fa 3.371176
## 57 G3 Fa 2.500814
## 58 G3 Fa 3.556720
## 59 G3 Fa 3.317399
## 60 G3 Fa 3.011181
## 61 G3 Fs 3.001382
## 62 G3 Fs 2.416850
## 63 G3 Fs 2.529113
## 64 G3 Fs 2.918011
## 65 G3 Fs 3.140656
## 66 G3 Fs 3.309762
## 67 G3 Fs 2.574481
## 68 G3 Fs 2.973637
## 69 G3 Fs 3.473589
## 70 G3 Fs 3.073387
## 71 G3 Fs 3.449128
## 72 G3 Fs 3.274772#72 lotes
#uno de los factores me afecta el bloque, aca es el plan de vuelo, o la fertilizacion, o el sistema de riego, y dentro pongo el otro factor
#parcelas divididas, uso de los propios factores para bloquearcollapsibleTree(df8, c("F1", "F2", "rto"))# diseño factorial parcialmente anidado en 3 etapas
#Ejemplo de 2 especies de cerdo, dos organos, e implante y sin implante
df9= data.frame(F1=gl(2,12,24,c("Sp1","Sp2")),
F2=gl(2,6,24,c("H","R")), #higado y riñon
F3=gl(2,3,24,c("C","S")), #con implante, y sin implante
rto=rnorm(24,3,0.3))
collapsibleTree(df9,c("F1","F2","F3","rto"))#El higado no es ligado, varia por la especie, factorial incompleto
#parcialmente dado que el higado de cada especie me afecta el experimentoKriging, interpolacion con kriging
#cambio de propiedades de manera espacial
#el concepto de bloque es buscar areas lo mas homogeneas para la propiedad,
#el pH tiene un patron y materia organica es diferente
#superposicion de mapas? Nop
#bloques espaciales, absurdos, dado la alta variabilidad del suelo, toca manejar sustratos para manejar, bajar variabilidad#diametro ecuatorial y diametro polar
set.seed(2)
df11 = data.frame(
DE = sort.int(rnorm(36,4.5,0.6), partial = 20),
DP = sort.int(rnorm(36,5,0.65), partial = 20))
#view
plot(df11$DE, df11$DP, pch = 16)
mean(df11$DE)## [1] 4.622295mean(df11$DP)## [1] 4.891214\[H_0:\mu_{DE} = \mu_{DP}\]
plantear mi hipotesis \[H_a:\mu_{DE} \neq \mu_{DP}\]
#los bayecianos dicen, para probar que hay diferencia, pero usan estadisticamente para probar que no son iguales, deberia mas bien probar que son diferentes a partir de mostrar que es diferente
#el contraste de hipotesis es absurdo
#trato mis hipotesis conociendo temas biologicos, y propongo que una cosa es mejor porque es diferente
#El uso de x trt es mejor
#prueba t
t.test(x = df11$DE, y = df11$DP)##
## Welch Two Sample t-test
##
## data: df11$DE and df11$DP
## t = -1.5869, df = 68.52, p-value = 0.1171
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.60702501 0.06918774
## sample estimates:
## mean of x mean of y
## 4.622295 4.891214#me quedo con la Hipotesiis nula, pero es preocupante, porque no evaluo lo que me mostró la prueba# t.test(x = df2$de, y = df2$dl)
t.test(x = df11$DE, y = df11$DP)##
## Welch Two Sample t-test
##
## data: df11$DE and df11$DP
## t = -1.5869, df = 68.52, p-value = 0.1171
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.60702501 0.06918774
## sample estimates:
## mean of x mean of y
## 4.622295 4.891214boxplot(df11$DP, df11$DE)
points(1:2,c(mean(df11$DE),
mean(df11$DP)),
col = "red", pch = 16, cex = 3)
df11[df11$DE == max(df11$DE),][1] = 6.4#boxplot, comportamiento de los dos diametros
boxplot(df11$DE, df11$DP)
points(1:2, c (mean(df11$DE), mean(df11$DP)), col = "red", pch = 16, cex = 3)
#tienen ue correrse al tiempo, para poner los puntosdiam = c(df11$DE, df11$DP)
val = c(rep("DE", 36), rep("DP", 36))
df12 = data.frame(diam,val)
df12 %>%
ggplot(aes(x = df12$diam))+
geom_density(aes(fill = df12$val),
alpha = 0.4)+
theme_bw()## Warning: Use of `df12$val` is discouraged.
## ℹ Use `val` instead.## Warning: Use of `df12$diam` is discouraged.
## ℹ Use `diam` instead.
#tomar decisiones con mas datos además de solo el promedio, veamos como se comportan los datos