#crear cojuto de datos, 1. fijamos semilla
#acostumbramos a medir muchas cosas, y seguimos haciendo el analisis variable por variable
#Comencemos a usar analisis multivariante, usemos la prueba J y ver relacion entre variables
set.seed(2024)
misdatos = data.frame(clorofila = sort(rnorm(72,250,10)),
nitrogeno = sort(rnorm(72,2.5,0.12)),
met = gl(2,36,72, c("Met1", "Met2"))) #una variable que separa, metodo 1 y metodo 2
#tengo mi data.frame con 3 variables, 2 respuestas, 2 var cuantitativas y una factor, es la que yo controlo como metodo, voy a medir clorofila en 2 equipos diferentes.
#1. voy a hacer un diagrama de dispersion para ver la relacion en mis variables, a ver si hay alguna relacion
plot(misdatos$nitrogeno, misdatos$clorofila)
cor(misdatos$nitrogeno,misdatos$clorofila)
## [1] 0.9755675
#como ordene mis datos, hay relacion, lo hice intencionalmente
#2. como hago para meter la otra variable, meter los metodos?
plot(misdatos$nitrogeno,misdatos$clorofila,
col =misdatos$met, pch = 19)
#aca coloreé los datos que fueron tomados con diferentes metodos
#si hago cajas de bigotes, se separan las variables
#si tuviera que decir sobre los factores,
#aca sigo con la parte descriptiva de mi investigación
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, union
#
#
#libreria de resumenes estadisticos
#pipe, agrupo por metodo
#tengo la data completa de 72 columna, ahora particioneme los datos segun el metodo, los parte en 2
#Aca se intentó hacer una tuberia pero creo que se fracasó en el intento, lo colocó en #
#misdatos %>%
# group_by(met) %>%
#summarise( media = mean(clorofila), desv = sd(clorofila), min = min (clorofila),
#max = max (clorofila), cv = sd (clorofila)/mean (clorofila)*100, media_nitro = mean(nitrogeno),
#desv_nitro = sd(nitrogeno), min_nitro = min (nitrogeno), max_nitro = max (nitrogeno),
#cv = sd (nitrogeno)/mean (nitrogeno)*100, )
#aca nombre y distingui los de nitrogeno
#trasnponer tabla, para visualizar
View(t(misdatos))
View(misdatos)
#esta tabla nos muestra que no tenemos una alta correlacion, así en la imagen así lo veamos
#si hay baja o no hay correlacion, puedo analizarlas de manera diferente, en este analisis se pierde el analisis en separado.
#la opinion por separado, a menudo no cuenta
#a menudo, si no hay interaccion (ambos ), podemos hacer analisis de los datos marginales
#la omision de la correlacion, me afecta
#cuando hay estrecha relacion entre variables, debo tener en cuenta, debo medir todo en un conjunto y no variable por variable
#en encuestas, tengo 10 preguntas y las termino analizando una por una, el profe SE DEPRIME jajaja, no es lo que se debe hacer, pretender sacar conclusiones por partes, peligroso.pensamos que siempre se hace así.
#Indicadores de gestion, evaluados por areas: como esta el desempeño segun lo que dicen los indicadores 1 a 1
#Quiebra del banco: los indicadores individuales muestran que el banco esta bien, pero si miro los indicadores en conjunto, se veía que el banco estaba mal
##union datos estadisticos sociales (segregacion de poblacion en segmentos para publicidad), a quienes se les aplica la publicidad generada por IA
#solemos poner la prueba t student por cada datos
#corramos prueba t
t.test(misdatos$clorofila, misdatos$nitrogeno)
##
## Welch Two Sample t-test
##
## data: misdatos$clorofila and misdatos$nitrogeno
## t = 197.02, df = 71.016, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 244.2415 249.2357
## sample estimates:
## mean of x mean of y
## 249.226933 2.488317
#Prueba t es para cada respuesta, pues MAL, debe ser para ambas variables, no para una sola
#Aca me dio una prueba Welch, otras no las tienen
#T al cuadrado, es una prueba para todas las respuestas
#T tengo varias modalidades, a una muestra o 2 muestras, este segundo tiene otro camino, las 2 muestras son independientes o pareadas?
## las tomé de dos metodos diferentes, es posible que no esten pareadas, si fuera datos en pH,midiendo en 2 profundidades, aca se relaciona el pH de arriba con la de abajo. Similar a tomar datos de un fruto.
#MUESTRAS RELACIONADAS SON MUESTRAS RELACIONADAS, si no lo se, correr a prueba en R no es suficiente.
#si son independientes, me pregunto si tiene varianzas iguales o desiguales
##con simulacion puedo hacer pruebas
#Miremos la prueba que corremos, miremos sus argumentos, en HELP t.test
#Aca se intentó pasar por un pipe pero no estoy seguró el por qué, lo pongo en # ya que no me deja seguir
#df2.misdatos %>%
# filter(met)
#df2
t.test(misdatos$clorofila ~ misdatos$met)
##
## Welch Two Sample t-test
##
## data: misdatos$clorofila by misdatos$met
## t = -11.19, df = 59.987, p-value = 2.556e-16
## alternative hypothesis: true difference in means between group Met1 and group Met2 is not equal to 0
## 95 percent confidence interval:
## -19.92332 -13.88078
## sample estimates:
## mean in group Met1 mean in group Met2
## 240.7759 257.6780
#la virgulilla es funcion
#Este es un valor bastante bajo, un cero, rechazo la H0, pero aca hay varios argumentos que no he mirado, corri la prueba por defecto.
#ahora hagamoslo para nitrogeno
t.test(misdatos$nitrogeno ~ misdatos$met)
##
## Welch Two Sample t-test
##
## data: misdatos$nitrogeno by misdatos$met
## t = -11.857, df = 67.144, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Met1 and group Met2 is not equal to 0
## 95 percent confidence interval:
## -0.2152943 -0.1532564
## sample estimates:
## mean in group Met1 mean in group Met2
## 2.396179 2.580455
#en ambas variables las medi en t test y me dio diferencia, en ambos nos dio que rechazar, pero el 0.05 del valor p
#correciones, esta bonferroni
#en ambos concluí pensando en el valor p de 0.05%, pero al hacer multivariable, este valor es diferentes, ya no es 0.05
#¿POR QUE CAMBIA EL VALOR DE ALPHA?
#Enfoque bayesiano no usa el alpha?
#en la interaccion de dos vriables, un p valor puede bajar hasta 0.025, pero esto lo debo definir yo
#No hago muchas hipotesis, hago una sola
#prueba que no debemos hacer
t.test(misdatos$clorofila ~ misdatos$met)
##
## Welch Two Sample t-test
##
## data: misdatos$clorofila by misdatos$met
## t = -11.19, df = 59.987, p-value = 2.556e-16
## alternative hypothesis: true difference in means between group Met1 and group Met2 is not equal to 0
## 95 percent confidence interval:
## -19.92332 -13.88078
## sample estimates:
## mean in group Met1 mean in group Met2
## 240.7759 257.6780
#en ambos me dieron valores diferentes
#misdatos2 = data.frame(clorofila = sort(rnorm(36,250,10)),
#sort(rnorm(36,254.8, 10))
#nitrogeno = sort(rnorm(36,5,0.12)),
#sort(rnorm(36,5,0.3))
#met = gl(2,36,72, c("Met1", "Met2")))
#version multivariante t student para correr ensayos de varias variables
#PRUEBA T CUADRADO de Hotelling
#comparar trt para todas las respuestas y no uno x uno
#dos pacquetes hacer esta prueba
#ICSNP ES UNA de estas librerias
#la otra es Hotelling
library(ICSNP)
## Loading required package: mvtnorm
## Loading required package: ICS
library(Hotelling)
## Loading required package: corpcor
##
## Attaching package: 'Hotelling'
## The following object is masked from 'package:dplyr':
##
## summarise
set.seed(123456)
clorof<-misdatos$clorofila #var1
nitr<-misdatos$nitrogeno #var2
#Z<-rbind(clorof,nitr)
g = factor(misdatos$met) #grupo
#g<-factor(rep(c(1,2),c(25,25)))
datos_prueba1 = data.frame(clorof, nitr,g)
resp = as.matrix(datos_prueba1[ , 1:2])
HotellingsT2(resp ~ datos_prueba1$g ) #Mi valor de interes, correr la prueba T para 2 variables
##
## Hotelling's two sample T2-test
##
## data: resp by datos_prueba1$g
## T.2 = 69.386, df1 = 2, df2 = 69, p-value < 2.2e-16
## alternative hypothesis: true location difference is not equal to c(0,0)
#HotellingsT2(Z~g,mu=rep(-0.5,4))
#rm(.Random.seed)
#La prueba de T^2 muestra en T2 la magnitud de la diferencia, df1 y df2 dice los Grados de libertad de los datos
#2 variables y 72 observaciones, df1: 2 y df2 : 69, df2 es el tamaño de la muestra - # de grupos
#p valor muy pequeño, se rechaza H0 de igualdad, se acepta que hay diferencia
cov(resp)
## clorof nitr
## clorof 112.911063 1.17716946
## nitr 1.177169 0.01289516
#Coeficiente de variacion, es una varianza para medir cuanto varian mis datos
#Coef de variacion, algunos piensan que por encima de 20% no era bueno, PRECONCEPCION
#
##########################coeficiente de variacion
#datos de diametros de duraznos
#library(faux)
diametro = faux:: rnorm_multi(n=100,
vars = 2,
mu = c(4,5),
sd = c(0.1, 1.1),
r = 0.8,
varnames = c("DE", "DP"))
loc=gl(2,50,100, c("Pamplona", "Paipa")) #localidad
diametro=data.frame(loc,diametro)
#el código genera datos simulados con una estructura multivariada, donde las variables están correlacionadas, y agrega una variable de localidad a estos datos.
datos.boyaca=data.frame(loc,diametro)
datos.boyaca
## loc loc.1 DE DP
## 1 Pamplona Pamplona 4.117674 5.913395
## 2 Pamplona Pamplona 3.912286 4.701132
## 3 Pamplona Pamplona 4.043515 4.604249
## 4 Pamplona Pamplona 3.972292 5.098769
## 5 Pamplona Pamplona 4.129166 7.481222
## 6 Pamplona Pamplona 4.008283 5.922179
## 7 Pamplona Pamplona 4.181177 6.438111
## 8 Pamplona Pamplona 4.200868 7.752884
## 9 Pamplona Pamplona 4.151085 6.280862
## 10 Pamplona Pamplona 3.890653 4.536703
## 11 Pamplona Pamplona 3.971795 3.900493
## 12 Pamplona Pamplona 3.895757 3.775749
## 13 Pamplona Pamplona 4.000671 4.938321
## 14 Pamplona Pamplona 4.123650 6.289720
## 15 Pamplona Pamplona 4.100562 6.157360
## 16 Pamplona Pamplona 4.028450 5.061628
## 17 Pamplona Pamplona 3.997381 4.187348
## 18 Pamplona Pamplona 4.014193 6.027984
## 19 Pamplona Pamplona 4.158245 6.833238
## 20 Pamplona Pamplona 4.055715 5.614863
## 21 Pamplona Pamplona 3.868853 4.175788
## 22 Pamplona Pamplona 4.196822 6.375196
## 23 Pamplona Pamplona 4.138714 5.032448
## 24 Pamplona Pamplona 4.136998 5.199608
## 25 Pamplona Pamplona 4.075792 5.506029
## 26 Pamplona Pamplona 3.922874 4.533029
## 27 Pamplona Pamplona 4.007502 5.017794
## 28 Pamplona Pamplona 4.037040 5.776785
## 29 Pamplona Pamplona 4.136324 6.064770
## 30 Pamplona Pamplona 3.929835 4.318951
## 31 Pamplona Pamplona 3.928712 4.058744
## 32 Pamplona Pamplona 4.066552 5.072066
## 33 Pamplona Pamplona 3.868702 3.852645
## 34 Pamplona Pamplona 3.749501 1.978452
## 35 Pamplona Pamplona 4.002119 3.750397
## 36 Pamplona Pamplona 3.878131 4.056000
## 37 Pamplona Pamplona 4.137576 6.715163
## 38 Pamplona Pamplona 3.972146 6.124561
## 39 Pamplona Pamplona 4.073187 6.149157
## 40 Pamplona Pamplona 3.958558 3.769006
## 41 Pamplona Pamplona 3.797784 3.830055
## 42 Pamplona Pamplona 4.070421 6.065122
## 43 Pamplona Pamplona 4.034908 5.186589
## 44 Pamplona Pamplona 3.865250 4.018716
## 45 Pamplona Pamplona 3.943283 5.179180
## 46 Pamplona Pamplona 4.007402 4.444383
## 47 Pamplona Pamplona 3.876216 3.940869
## 48 Pamplona Pamplona 4.109942 4.866208
## 49 Pamplona Pamplona 4.108056 6.193009
## 50 Pamplona Pamplona 3.880050 3.668855
## 51 Paipa Paipa 3.842753 3.056645
## 52 Paipa Paipa 3.902983 4.463342
## 53 Paipa Paipa 4.055114 5.351471
## 54 Paipa Paipa 4.112326 6.607070
## 55 Paipa Paipa 4.132998 6.690249
## 56 Paipa Paipa 3.896789 4.632019
## 57 Paipa Paipa 3.956537 3.811680
## 58 Paipa Paipa 3.788411 3.366880
## 59 Paipa Paipa 3.928771 4.725697
## 60 Paipa Paipa 3.942328 4.869380
## 61 Paipa Paipa 3.968376 4.284632
## 62 Paipa Paipa 4.019128 5.344818
## 63 Paipa Paipa 4.000783 5.135752
## 64 Paipa Paipa 4.079425 4.068791
## 65 Paipa Paipa 4.068589 5.658970
## 66 Paipa Paipa 4.058749 4.723749
## 67 Paipa Paipa 4.099181 4.788555
## 68 Paipa Paipa 3.967270 5.027995
## 69 Paipa Paipa 4.063480 4.459420
## 70 Paipa Paipa 3.964569 4.186575
## 71 Paipa Paipa 4.198733 6.425874
## 72 Paipa Paipa 3.935190 5.628976
## 73 Paipa Paipa 3.796933 3.018715
## 74 Paipa Paipa 4.054878 5.480732
## 75 Paipa Paipa 3.813146 3.420010
## 76 Paipa Paipa 4.096357 6.046652
## 77 Paipa Paipa 3.932061 4.205450
## 78 Paipa Paipa 4.029591 5.051763
## 79 Paipa Paipa 3.757213 3.288687
## 80 Paipa Paipa 3.781881 2.469852
## 81 Paipa Paipa 4.131198 6.270957
## 82 Paipa Paipa 3.946258 4.224593
## 83 Paipa Paipa 3.806032 3.273747
## 84 Paipa Paipa 4.024329 5.571441
## 85 Paipa Paipa 3.886956 3.830599
## 86 Paipa Paipa 3.933238 5.056905
## 87 Paipa Paipa 4.152181 5.926708
## 88 Paipa Paipa 4.071935 5.473227
## 89 Paipa Paipa 4.081165 5.572221
## 90 Paipa Paipa 3.937384 4.723809
## 91 Paipa Paipa 4.057294 4.446380
## 92 Paipa Paipa 4.127851 6.384150
## 93 Paipa Paipa 3.978546 5.626345
## 94 Paipa Paipa 3.961247 4.621711
## 95 Paipa Paipa 4.079765 5.795708
## 96 Paipa Paipa 4.054590 5.951569
## 97 Paipa Paipa 4.008227 5.407729
## 98 Paipa Paipa 4.126052 6.742016
## 99 Paipa Paipa 4.018886 5.065383
## 100 Paipa Paipa 4.033645 5.139880
#Aca no pude hacer Knit así que colocaré la siguiente parte del codigo en ##
#datos.boyaca %>%
# group_by(loc) %>%
#summarise(media_DE = mean(DE),
# sd_DE = sd(DE),
# cv_DE = sd_DE / media_DE * 100,
# media_DP = mean(DP),
# sd_DP = sd(DP),
# cv_DP = sd_DP / media_DP * 100)
#a novel coefficient of variation
#aca el DE y el DP crecen relacionados, le debo usar el Coef. de variacion diferente dado que estan relacionados.