Comandos_6
Sergio Camiz
2024-03-09
Aquí se definen parámetros generales y el espacio de trabajo
options(max.print = 100000, width=300) # dimensiones salida
getwd() # donde estamos ahora## [1] "/home/sergio/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_6/R"#setwd("~/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_6 /R") # donde vamos#Cruzado de variables cuantiativas
Se colgan las funciones: la primera dibuja un gráfico con caja y bigotes con promedio y desvíos.
boxplot.ext <- function(vec,range=1.5,nsd=2,at=1,...) { # boxplot
n <- length(vec) # vector length
m <- mean(vec) # mean
sd <- sqrt(sum(vec^2)/n-m^2) # standard deviation
sm <- (sum(abs(vec-m)))/n # average deviation
infd <- m-nsd*sd # sd extremes
supd <- m+nsd*sd
infm <- m-nsd*sm # sm extremes
supm <- m+nsd*sm
pm <- "\u00b1" # pm symbol
Encoding(pm) <- "UTF-8" # Unicode
namleg <- c("Mean", paste("Mean",pm, nsd,"std.deviations",sep=" "), # legend names
paste("Mean",pm, nsd, "avg.deviations",sep=" "))
boxplot(vec,...) # drws the boxplot
points(1,m,pch=20, col="red") # plots the mean
points(c(1,1),c(infd,supd),pch=20, col="blue") # plots std. devs.
points(c(1,1),c(infm,supm),pch=20,col="green") # plots avg. devs.
legend("topleft",pch=20,col=c("red","blue","green"), legend=namleg, cex=.8) # draws the legend
}La segunda cálcula estadísticas descriptivas por variables cualitativas
stat.con <- function (vec,boxplot=TRUE,quant=NULL,...) {
res <- list()
name <- deparse(substitute(vec)) # recoje el nombre
res$Name <- name # y lo coloca
if(is.null(quant)) quant <- c(.01,.05,.25,.50,.75,.95,.99) # cuantiles estándar
res$Quantiles <- quantile(vec,quant,type=1) # cuantiles
n <- length(vec) # número de unidades
med <- quantile(vec,.50) # mediana
m <- sum(vec)/n # promedio
ss <- (sum((vec-m)^2)) # desvíos cuadráticos
s2 <- ss/n # varianza
s <- sqrt(s2) # desvio estandar
cv <- s / m # coeficiente de variación
sm <- (sum(abs(vec-m)))/n # desvío promedio
stat <- round(cbind(n,med,m,ss,s2,s,cv,sm),3) # resultados
stat <- t(stat)
rownames(stat) <- c("número","mediana","promedio","desvíos_cuadr.","varianza","desvío_estándar",
"coef.var","desv.prom.")
colnames(stat) <- name # nombra stat
res$Statistics <- stat # estadísticas
########## boxplot con promedio y desvíos
if (boxplot) {
boxplot.ext(vec,main=name)
}
res
}##Relaciones entre un caracter cualitativo y uno cuantitativo.
Se trabaja con los datos de Iris. Se leen los datos y se calculan estadísticas descritivas.
iris <- read.csv("Iris.csv",header=TRUE,row.names=1) # data
iris$Species <- factor(iris$Species) # define species as factor
str(iris) # structure of iris## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...attach(iris) # se conecta iris
a=stat.con(iris[,2],) # ejemplo de juntar dos salidas
b=stat.con(iris[,3])
cbind(a$Statistics,b$Statistics)## iris[, 2] iris[, 3]
## número 150.000 150.000
## mediana 3.000 4.350
## promedio 3.057 3.758
## desvíos_cuadr. 28.307 464.325
## varianza 0.189 3.096
## desvío_estándar 0.434 1.759
## coef.var 0.142 0.468
## desv.prom. 0.337 1.563#estadísticas totales
total <- stat.con(Sepal.length) # estadísticas globales
total## $Name
## [1] "Sepal.length"
##
## $Quantiles
## 1% 5% 25% 50% 75% 95% 99%
## 4.4 4.6 5.1 5.8 6.4 7.3 7.7
##
## $Statistics
## Sepal.length
## número 150.000
## mediana 5.800
## promedio 5.843
## desvíos_cuadr. 102.168
## varianza 0.681
## desvío_estándar 0.825
## coef.var 0.141
## desv.prom. 0.688# mismas estadísticas pero compartidas por cada modalidad
n_unit <- table(Species); n_unit # números## Species
## setosa versicolor virginica
## 50 50 50nm <- length(n_unit); nm # número de modos ## [1] 3med <- tapply(Sepal.length,Species,median) ;median # medianas ## function (x, na.rm = FALSE, ...)
## UseMethod("median")
## <bytecode: 0x653342ee13a8>
## <environment: namespace:stats>mean <- tapply(Sepal.length,Species,mean) ;mean # promedios ## setosa versicolor virginica
## 5.006 5.936 6.588ss <- tapply(Sepal.length,Species,var)*(n_unit-1) # sumas de cuadrados
var <- ss/n_unit # varianzas
stdev <- sqrt(var); stdev # desvíos estándar ## setosa versicolor virginica
## 0.3489470 0.5109834 0.6294887cv <- stdev/mean; cv # coef. de variación ## setosa versicolor virginica
## 0.06970575 0.08608210 0.09555080qcla <- as.data.frame(split(Sepal.length,Species)) # divide el vector segun los modos
dev <- abs(sweep(qcla,2,mean,"-")) # cálcula los desvíos absolutos
desm <- apply(dev,2,mean) # desvíos promedios
brk <- round(rbind(n_unit,med,mean,ss,var,stdev,cv,desm),3); brk # tabla## setosa versicolor virginica
## n_unit 50.000 50.000 50.000
## med 5.000 5.900 6.500
## mean 5.006 5.936 6.588
## ss 6.088 13.055 19.813
## var 0.122 0.261 0.396
## stdev 0.349 0.511 0.629
## cv 0.070 0.086 0.096
## desm 0.271 0.421 0.503brtot <- cbind(brk,total) # incluye el total## Warning in cbind(brk, total): number of rows of result is not a multiple of vector length (arg 2)colnames(brtot) <- c(colnames(brk),"total")
brtot## setosa versicolor virginica total
## n_unit 50 50 50 "Sepal.length"
## med 5 5.9 6.5 numeric,7
## mean 5.006 5.936 6.588 numeric,8
## ss 6.088 13.055 19.813 "Sepal.length"
## var 0.122 0.261 0.396 numeric,7
## stdev 0.349 0.511 0.629 numeric,8
## cv 0.07 0.086 0.096 "Sepal.length"
## desm 0.271 0.421 0.503 numeric,7# razón empírica de correlación
sst <- total$Statistics[4] # toma la suma total
e <- (sst-sum(ss))/sst; e # calcula la razón de correlación ## [1] 0.6187045Se pueden construir funciones que dibujan gráficos por modos y estadísticas por modo también.
boxplot.mext <- function(quant,qual,range=1.5,nst=2,at=1:nm,...) {
n <- table(qual) # total
mean <- tapply(quant,qual,mean) # promedios
stdev <- sqrt(tapply(quant,qual,var)*(n-1)/n) # desvíos estándar
qcla <- as.data.frame(split(quant,qual)) # divide el vector segun los modos
dev <- abs(sweep(qcla,2,mean,"-")) # cálcula los desvíos absolutos
desm <- apply(dev,2,mean) # promedios
infd <- mean-nst*stdev # extremos de stdev
supd <- mean+nst*stdev
infm <- mean-nst*desm # y de desvio promedio
supm <- mean+nst*desm
pm <- "\u00b1" # pm symbol
Encoding(pm) <- "UTF-8" # Unicode
namleg <- c("Mean", paste("Mean",pm, nst,"std.deviations",sep=" "), # legend names
paste("Mean",pm, nst, "avg.deviations",sep=" "))
boxplot(quant~qual,...) # boxplot
points(1:nm,mean,pch=20, col="red") # dibuja el promedio
points(1:nm,supd,pch=20, col="blue") # dibuja los st.dev.
points(1:nm,infd,pch=20, col="blue")
points(1:nm,supm,pch=20, col="green") # dibuya los desv. promedios
points(1:nm,infm,pch=20, col="green")
legend("topleft",pch=20,col=c("red","blue","green"), legend=namleg, cex=.8) # draws the legend
}Breakdown <- function(quant,qual){
res <- list() # lista de salida
namev <- deparse(substitute(quant)) # nombre variable
nameq <- deparse(substitute(qual)) # nombre factor
name <- paste0("Boxplot of ",namev," broken down by ", nameq) # titulo
total <- stat.con(quant,boxplot=FALSE)$Statistics # estadísticas totales
sst <- total[4] # recoje sst
# estadísticas por clase
n_unit <- table(qual) # números
nm <- length(n_unit); nm # número de modos
med <- tapply(quant,qual,median) # medianas
mean <- tapply(quant,qual,mean) # promedis
ss <- tapply(quant,qual,var)*(n_unit-1) # sumas desvíos cuadráticos
var <- ss/n_unit # varianzas
stdev <- sqrt(var) # desvíos estándar
cv <- stdev/mean # coeficientes de variación
qcla <- as.data.frame(split(quant,qual)) # divide el vector
dev <- abs(sweep(qcla,2,mean,"-")) # cálcula los desvíos absolutos
desm <- apply(dev,2,mean) # promedios
brk <- rbind(n_unit,med,mean,ss,var,stdev,cv,desm);brk # construye tabla
brtot <- cbind(brk,total) # ajunta total
colnames(brtot) <- c(colnames(brk),"total") # nombres
res$breakdown <- brtot
# razón de correlación # tabla
e <- (sst-sum(ss))/sst # calcula razón de correlación
res$correlation_ratio <- e
# boxplot con promedio y 2 desvíos estándar
boxplot.mext(quant,qual, range=1.5, nst=2, at=1:nm,main=name,xlab=nameq,ylab=namev)
res
}Se pueden correr las estadísticas de las cuatro mediciones de Iris.
Breakdown(Sepal.length,Species)
## $breakdown
## setosa versicolor virginica total
## n_unit 50.00000000 50.0000000 50.0000000 150.000
## med 5.00000000 5.9000000 6.5000000 5.800
## mean 5.00600000 5.9360000 6.5880000 5.843
## ss 6.08820000 13.0552000 19.8128000 102.168
## var 0.12176400 0.2611040 0.3962560 0.681
## stdev 0.34894699 0.5109834 0.6294887 0.825
## cv 0.06970575 0.0860821 0.0955508 0.141
## desm 0.27072000 0.4214400 0.5025600 0.688
##
## $correlation_ratio
## [1] 0.6187045Breakdown(Sepal.width,Species)
## $breakdown
## setosa versicolor virginica total
## n_unit 50.0000000 50.0000000 50.0000000 150.000
## med 3.4000000 2.8000000 3.0000000 3.000
## mean 3.4280000 2.7700000 2.9740000 3.057
## ss 7.0408000 4.8250000 5.0962000 28.307
## var 0.1408160 0.0965000 0.1019240 0.189
## stdev 0.3752546 0.3106445 0.3192554 0.434
## cv 0.1094675 0.1121460 0.1073488 0.142
## desm 0.2873600 0.2548000 0.2421600 0.337
##
## $correlation_ratio
## [1] 0.4007843Breakdown(Petal.length,Species)
## $breakdown
## setosa versicolor virginica total
## n_unit 50.0000000 50.0000000 50.0000000 150.000
## med 1.5000000 4.3500000 5.5500000 4.350
## mean 1.4620000 4.2600000 5.5520000 3.758
## ss 1.4778000 10.8200000 14.9248000 464.325
## var 0.0295560 0.2164000 0.2984960 3.096
## stdev 0.1719186 0.4651881 0.5463479 1.759
## cv 0.1175914 0.1091991 0.0984056 0.468
## desm 0.1315200 0.3792000 0.4400000 1.563
##
## $correlation_ratio
## [1] 0.9413717Breakdown(Petal.width,Species)
## $breakdown
## setosa versicolor virginica total
## n_unit 50.0000000 50.0000000 50.0000000 150.000
## med 0.2000000 1.3000000 2.0000000 1.300
## mean 0.2460000 1.3260000 2.0260000 1.199
## ss 0.5442000 1.9162000 3.6962000 86.570
## var 0.0108840 0.0383240 0.0739240 0.577
## stdev 0.1043264 0.1957652 0.2718897 0.760
## cv 0.4240911 0.1476359 0.1342002 0.633
## desm 0.0825600 0.1571200 0.2280800 0.658
##
## $correlation_ratio
## [1] 0.928883##Dato cuantitativos
Con dos datos cuantitativo se calcula la covarianza
n <- dim(iris)[1] # total observaciones
cov(iris$Sepal.length,iris$Sepal.width) # covarianzas de las muestras## [1] -0.042434cov(iris$Petal.length,iris$Petal.width) # que estiman la población## [1] 1.295609cov(iris[,1:4]) # matriz de covarianza## Sepal.length Sepal.width Petal.length Petal.width
## Sepal.length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.width 0.5162707 -0.1216394 1.2956094 0.5810063cov(iris[,1:4])*(n-1)/n # covarianzas de la población ## Sepal.length Sepal.width Petal.length Petal.width
## Sepal.length 0.68112222 -0.04215111 1.2658200 0.5128289
## Sepal.width -0.04215111 0.18871289 -0.3274587 -0.1208284
## Petal.length 1.26582000 -0.32745867 3.0955027 1.2869720
## Petal.width 0.51282889 -0.12082844 1.2869720 0.5771329 # datos sin inferencia
corr <- cor(iris[,1:4]) # matriz de correlación
corr <- round(corr,3)
corr## Sepal.length Sepal.width Petal.length Petal.width
## Sepal.length 1.000 -0.118 0.872 0.818
## Sepal.width -0.118 1.000 -0.428 -0.366
## Petal.length 0.872 -0.428 1.000 0.963
## Petal.width 0.818 -0.366 0.963 1.000plot(iris[,1:4]) # gráfico total
plot(Sepal.length,Sepal.width,pch=20,sub=paste0("correlation = ",corr[1,2])) # gráficos individuales
plot(Petal.length,Petal.width,pch=20,sub=paste0("correlation = ",corr[3,4]))
Correlaciones para datos ordinales
Se trabaja con datos de rango de notas de alumnos
mrk <- read.csv("Marks_Ranks.csv",header=TRUE); mrk## M_English M_Maths R_English R_Maths
## 1 56 66 9.0 4
## 2 75 70 3.0 2
## 3 45 40 10.0 10
## 4 71 60 4.0 7
## 5 61 65 6.5 5
## 6 64 56 5.0 9
## 7 58 59 8.0 8
## 8 80 77 1.0 1
## 9 76 67 2.0 3
## 10 61 63 6.5 6cor(mrk[,3:4], method="spearman") # Spearman se emplea sobre rangos## R_English R_Maths
## R_English 1.0000000 0.6686961
## R_Maths 0.6686961 1.0000000cor(mrk[,3:4], method="pearson") # igual a Pearson## R_English R_Maths
## R_English 1.0000000 0.6686961
## R_Maths 0.6686961 1.0000000cor(mrk[,1:2], method="pearson") ## M_English M_Maths
## M_English 1.0000000 0.8005387
## M_Maths 0.8005387 1.0000000cor(mrk[,1:2], method="spearman") # pero también sobre valores## M_English M_Maths
## M_English 1.0000000 0.6686961
## M_Maths 0.6686961 1.0000000cor(iris[,1:4],method="spearman") # correlación más robusta## Sepal.length Sepal.width Petal.length Petal.width
## Sepal.length 1.0000000 -0.1667777 0.8818981 0.8342888
## Sepal.width -0.1667777 1.0000000 -0.3096351 -0.2890317
## Petal.length 0.8818981 -0.3096351 1.0000000 0.9376668
## Petal.width 0.8342888 -0.2890317 0.9376668 1.0000000También se emplea la correlación de Spearman para datos cuantitativos. Aquí un ejemplo de datos de SudAmerica.
sudam <- read.csv("SudAmerica.csv",row.names=1); sudam## Tajo_urbano IDH
## Argentina 86 0.833
## Bolivia 51 0.394
## Brasil 75 0.739
## Chile 86 0.863
## Colombia 70 0.758
## Ecuador 56 0.641
## Guyana 35 0.539
## Panama 54 0.731
## Paraguay 48 0.637
## Peru 70 0.600
## Suriname 48 0.749
## Uruguay 86 0.880
## Venezuela 84 0.824plot(sudam)
text(sudam,labels=rownames(sudam))
cor(sudam,method="pearson")## Tajo_urbano IDH
## Tajo_urbano 1.0000000 0.7441166
## IDH 0.7441166 1.0000000cor(sudam,method="spearman")## Tajo_urbano IDH
## Tajo_urbano 1.0000000 0.7978115
## IDH 0.7978115 1.0000000stat.con(sudam[,2])
## $Name
## [1] "sudam[, 2]"
##
## $Quantiles
## 1% 5% 25% 50% 75% 95% 99%
## 0.394 0.394 0.637 0.739 0.824 0.880 0.880
##
## $Statistics
## sudam[, 2]
## número 13.000
## mediana 0.739
## promedio 0.707
## desvíos_cuadr. 0.237
## varianza 0.018
## desvío_estándar 0.135
## coef.var 0.191
## desv.prom. 0.111