-
-
Save geofis/14b9284af187ad630950 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#SCRIPT CONSOLIDADO DE LA ASIGNATURA MÉTODOS DE INVESTIGACIÓN GEOGRÁFICA II, BASADO SOBRE TODO EN TRIOLA (2009) | |
#EL OBJETIVO DE ESTAS SENTENCIAS ES MOSTRAR CÓMO REALIZAR ALGUNOS EJERCICIOS DE EJEMPLO DE TRIOLA (2009) EN R. ¡NO SE EXPLICAN CONCEPTOS, NI SE INTERPRETAN RESULTADOS! SÓLO MOSTRAR LA SOLUCIÓN DE LOS MISMOS EN R | |
#SE RECOMIENDAN ESTAS OTRAS FUENTES: | |
#GUERRERO, 2010: http://geografiafisica.org/consigna/visitante/de_gd/introduccion_al_analisis_de_datos_con_R_cuaderno_de_ejercicios_andres_guerrero_BUENO.pdf | |
#PARADIS, 2003: https://cran.r-project.org/doc/contrib/rdebuts_es.pdf | |
#TANTO EN GUERRERO COMO EN PARADIS HAY MUCHOS TRUCOS SOBRE CÓMO MANIPULAR LOS DATOS EN R. LOS PAQUETES dplyr Y tidyr OFRECEN MUCHAS HERRAMIENTAS PARA MANIPULACIÓN ROBUSTA DE DATOS. SE RECOMIENDA APRENDER A USARLOS DESPUÉS (¡IMPORTANTE!) DE APRENDER LA MANIPULACIÓN BÁSICA SIN ELLOS | |
#ESTE SCRIPT SE DIVIDE EN TRES PARTES: | |
#PRIMERA: INTRODUCCIÓN, USANDO DATOS FICTICIOS GENERADOS EN R | |
#SEGUNDA: EJERCICIOS EXTRAÍDOS DE TRIOLA (2009) | |
#TERCERA: MANIPULACIÓN Y ANÁLISIS DE DATOS PROPIOS | |
#IMPORTANTE: CONSULTAR Y EJECUTAR, LÉASE BIEN, CON CONCIENCIA, LOS EJERCICIOS, TANTO DE TRIOLA (2009) COMO DE GUERRERO (2010) Y PARADIS (2003): | |
#CARGA DE PAQUETES. IMPORTANTE: DEBE CARGAR UNO A UNO, PARA VERIFICAR POSIBLES ERRORES | |
library(moments) #PARA CALCULAR EL SESGO Y OTROS ESTADÍSTICOS | |
library(dplyr) #MANIPULACIÓN DE DATOS | |
library(tidyr) #MANIPULACIÓN DE DATOS | |
library(ggplo2) #GRÁFICOS ESTILIZADOS | |
library(psych) #ÚTIL FUNCIÓN describe de este paquete | |
library(plotGoogleMaps) #CARGA PAQUETE PARA REALIZAR MAPAS EN GOOGLEMAPS | |
library(raster) #CARGA PAQUETE PARA TRABAJAR CON RASTERS | |
library(sp) #CARGA PAQUETE PARA MANEJAR OBJETOS ESPACIALES | |
library(maptools) #CARGA HERRAMIENTAS PARA MAPAS | |
library(plotKML) #CARGA HERRAMIENTAS PARA CREAR ARCHIVOS KML | |
#SI DURANTE LA EJECUCIÓN DEL SCRIPT APARECIERAN MENSAJES DEL TIPO "NO SE PUDO ENCONTRAR LA FUNCIÓN...", ES PROBABLE QUE DEBA VERIFICAR SI EL PAQUETE CORRESPONDIENTE SE HA CARGADO | |
#PRIMERA PARTE: INTRODUCCIÓN, USANDO DATOS FICTICIOS GENERADOS EN R | |
#OPERACIONES SENCILLAS SIN DATOS EXTERNOS | |
2+2 #REALIZA Y MUESTRA UNA SUMA | |
s<-2+2 #REALIZA UNA SUMA Y LA ASIGNA AL VECTOR s | |
c(1,2,3,4,5) #MUESTRA UN VECTOR | |
v<-c(1,2,3,4,5) #CREA UN VECTOR, PERO LO ASIGNA AL OBJETO v | |
ls() #MUESTRA LOS OBJETOS CREADOS. EN PRINCIPIO, SI NO HAY SESIONES PREVIAS RESTAURADAS, SÓLO DEBERÍAN APARECER c Y v | |
s+v #MUESTRA LA SUMA LOS OBJETOS s Y v | |
1:10 #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 1 AL 10, DE 1 EN 1 | |
1:100 #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 1 AL 100, DE 1 EN 1 | |
seq(1,500,by=5) #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 1 AL 500, DE 5 EN 5 | |
rep(1,100) #MUESTRA LA REPETICIÓN DEL NÚMERO "1" CIEN VECES | |
rep("mg",100) #MUESTRA LA REPETICIÓN DE LA CADENA "mg" CIEN VECES | |
rep(s,5) #MUESTRA LA REPETICIÓN DEL VECTOR s CINCO VECES | |
seq(1,100,by=5) #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 1 AL 100, DE 5 EN 5 | |
seq(1,15,by=3) #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 1 AL 15, DE 3 EN 3: 1,4,7,10,13. NO SE INCLUYE EL 15 | |
seq(0,15,by=3) #MUESTRA UNA SECUENCIA DE NÚMEROS DEL 0 AL 15, DE 3 EN 3: 0,3,6,9,12,15 | |
history(Inf) #MUESTRA UNA VENTANA CON EL HISTÓRICO DE SENTENCIAS | |
#SEGUNDA PARTE: EJERCICIOS EXTRAÍDOS DE TRIOLA (2009) | |
#CAPÍTULO 2. GRÁFICAS Y DATOS | |
#EDADES ACTRICES/ACTORES, PÁGINA 41 | |
m <- c(22,30,35,40,43,41,26,26,37,26,33,39,35,33,80,25,28,29,29,29,34,31,42,33,63,24,38,27,34,74,29,35,32,38,54,31,27,33,33,35,26,25,24,38,37,50,35,28,31,29,25,29,42,38,45,27,41,46,25,41,61,49,27,30,41,35,36,21,39,28,35,28,60,32,41,34) | |
h <- c(44,38,41,49,44,40,51,46,41,34,38,35,62,42,32,40,62,32,42,47,43,36,42,36,52,40,52,31,42,76,54,47,41,43,51,47,48,39,52,29,34,56,35,37,49,53,37,43,34,41,30,57,56,45,38,52,39,39,42,38,36,32,41,49,41,45,60,62,45,37,57,44,42,30,43,60) | |
table(m)#TABLAS DE FRECUENCIA. NO TIENE MUCHO SENTIDO SI NO SE AGRUPAN POR RANGOS | |
table(h)#ÍDEM ANTERIOR | |
#TABLA 2.2 DE LA PÁGINA 43 | |
cut(m,seq(21,81,by=10),include.lowest = TRUE, right=F) #INCLUYE EL VALOR MÁS BAJO DE TODOS, Y EL INTERVALO ES CERRADO POR LA IZQUIERDA, NO POR LA DERECHA. DE ESTOS RANGOS SE PUEDE HACER UNA TABLA DE FRECUENCIAS | |
table(cut(m,seq(21,81,by=10),include.lowest = TRUE,right=F)) | |
#HISTOGRAMAS, A PARTIR DE PÁGINA 51 | |
#PÁGINA 52, FIGURA 2.2 | |
par(mfrow=c(1,1)) | |
hist(m,breaks=seq(20.5,80.5,by=10),xaxt='n', xlab='Edades de las mejores actrices', ylab='Frecuencia', col='grey') #PÁGINA 52, FIGURA 2.2 | |
axis(side=1, at=seq(20.5,80.5,by=10), labels=seq(20.5,80.5,by=10)) | |
#PÁGINA 52, FIGURA 2.3. SE TRATA DE UN HISTOGRAMA DE FRECUENCIAS RELATIVAS. NO TIENE MUCHO SENTIDO, PORQUE O SON DE DENSIDADES O SON DE FRECUENCIAS, PERO DE PORCENTAJES NO SON COMUNES | |
hist(m, breaks=seq(20.5,80.5,by=10), freq=F, xaxt='n', yaxt='n', xlab='Edades de las mejores actrices', ylab='Frecuencia relativa', col='grey') | |
axis(side=1, at=seq(20.5,80.5,by=10), labels=seq(20.5,80.5,by=10)) | |
axis(side=2, at=seq(0,0.04,by=0.01), labels=paste(seq(0,40,by=10),"%",sep='')) | |
#PÁGINA 57, FIGURA 2.4. POLÍGONO DE FRECUENCIAS | |
plot(table(cut(m,seq(10.5,90.5,by=10))), type='o', xaxt='n', xlab='Edades de las mejores actrices', ylab='Frecuencia') | |
axis(side=1, at=seq(1,8,by=1), labels=c("",seq(25.5,75.5,by=10),"")) | |
#PÁGINA 57, FIGURA 2.5. POLÍGONO DE FRECUENCIAS RELATIVAS | |
plot(table(cut(m,seq(10.5,90.5,by=10)))/length(m)*100, type='o', xaxt='n', yaxt='n', xlab='Edades de las mejores actrices', ylab='Frecuencia relativa') | |
axis(side=1, at=seq(1,8,by=1), labels=c("",seq(25.5,75.5,by=10),"")) | |
axis(side=2, at=seq(0,40,by=10), labels=paste(seq(0,40,by=10),"%",sep="")) | |
lines(table(cut(h,seq(10.5,90.5,by=10)))/length(h)*100, type="o", lty=1, col="grey") | |
text(6,15,"actores",col='grey') | |
text(3.7,10,"actrices") | |
#GRÁFICO DE TALLO Y HOJAS | |
stem(m) #PÁGINA 59 | |
stem(h) #NO INCLUIDO EN EL LIBRO | |
#CAPÍTULO 3. ESTADÍSTICOS DESCRIPTIVOS, EXPLORAR Y COMPARAR DATOS | |
#MEDIDAS DE TENDENCIA CENTRAL. IMPRESCINDIBLE HABER CREADO LOS OBJETOS m Y h DEL CAPÍTULO ANTERIOR (VER ARRIBA) | |
pb <- c(5.4,1.1,0.42,0.73,0.48,1.1) #PÁGINA 79 | |
mean(pb) #MEDIA | |
median(pb) #MEDIANA | |
#MEJORES ACTRICES Y ACTORES. PÁGINA 81 | |
#MEJORES ACTRICES | |
mean(m) #MEDIA | |
median(m) #MEDIANA | |
#PARA LA MODA, CREAR FUNCIÓN O CARGAR PAQUETE RASTER U OTROS DISPONIBLES | |
moda <- function(x){y <- as.factor(x); freq <- summary(y); mode <- names(freq)[freq[names(freq)] == max(freq)]; as.numeric(mode)} #TOMADO DE Wei (2014), URL: http://stackoverflow.com/questions/2547402/standard-library-function-in-r-for-finding-the-mode | |
moda(m) #MODA | |
mean(range(m)) #MITAD DEL RANGO | |
#MEJORES ACTORES | |
mean(h) #MEDIA | |
median(h) #MEDIANA | |
moda(h) #MODA. ARRIBA FUE CREADA LA FUNCIÓN "moda", POR LO QUE AHORA SÓLO ES NECESARIO LLAMARLA | |
mean(range(h)) #MITAD DEL RANGO | |
#SESGO. PAQUETE moments YA CARGADO. VALOR DEL SESGO CERCANO 0 SIGNIFICA DISTRIBUCIÓN NO SESGADA. VALOR NEGATIVO, SESGADA A LA IZQUIERDA; VALOR POSITIVO, SESGADA A LA DERECHA | |
skewness(m) | |
hist(m) #COMPROBANDO GRÁFICAMENTE | |
skewness(h) | |
hist(h) | |
#MEDIDAS DE VARIACIÓN | |
#CREACIÓN DE TABLA Y GRÁFICOS | |
clientes <- data.frame(`1`=c(6,4,1),`2`=c(6,7,3),`3`=c(6,7,14),row.names=c('Banco 1: filas variables','Banco 2: una sola fila', 'Banco 3: múltiples filas'), check.names = F) #TOMADO DE TRIOLA (2009), PÁGINA 92, QUE CONSISTE EN UNA TABLA DE TIEMPOS DE ESPERA PARA LOS CLIENTES 1, 2 Y 3 EN 3 BANCOS CON DIFERENTES ESTILOS DE GESTIÓN DE FILAS | |
par(mfrow=c(1,3)) | |
barplot(t(clientes)[,1],ylim=c(0,15),yaxt='n') | |
axis(side=2,at=c(0,5,10,15),labels=c(0,5,10,15)) | |
abline(h=6,lwd=2) | |
barplot(t(clientes)[,2],ylim=c(0,15),yaxt='n') | |
axis(side=2,at=c(0,5,10,15),labels=c(0,5,10,15)) | |
abline(h=6,lwd=2) | |
barplot(t(clientes)[,3],ylim=c(0,15),yaxt='n') | |
axis(side=2,at=c(0,5,10,15),labels=c(0,5,10,15)) | |
abline(h=6,lwd=2) | |
#RANGO | |
diff(range(t(clientes)[,1])) | |
diff(range(t(clientes)[,2])) | |
diff(range(t(clientes)[,3])) | |
#DESVIACIÓN ESTÁNDAR | |
sd(t(clientes)[,1]) | |
sd(t(clientes)[,2]) | |
sd(t(clientes)[,3]) | |
#VARIANZA | |
var(t(clientes)[,1]) | |
var(t(clientes)[,2]) | |
var(t(clientes)[,3]) | |
#COEFICIENTE DE VARIACIÓN. SE PUEDE CALCULAR CON LA FUNCIÓN cv DEL PAQUETE raster, O SE PUEDE CREAR UNA SENCILLA FUNCIÓN DE R | |
(sd(m)/mean(m))*100 | |
(sd(h)/mean(h))*100 | |
#MEDIDAS DE POSICIÓN RELATIVA | |
#PUNTUACIONES z | |
mz <- round(scale(m),2) | |
hz <- round(scale(h),2) | |
data.frame(m,mz,h,hz) #VISTAS EN UNA TABLA JUNTO CON SU VALOR ORIGINAL | |
#CUARTILES | |
summary(m) #CUARTILES Y VALORES MÁXIMOS/MÍNIMOS PARA MEJORES ACTRICES | |
summary(h) #ÍDEM PARA ACTORES | |
#PERCENTILES | |
quantile(m,prob=0.2) #PÁGINAS 112 Y 113 | |
quantile(m,prob=seq(0.01,0.99,by=0.01)) #TODOS LOS PERCENTILES DEL 1% AL 99% | |
Hmisc::describe(m) #PERCENTILES HABITUALES | |
#RANGO INTERCUARTILAR (RIC) | |
paste("RIC =",as.vector(diff(summary(m)[c(2,5)]))) | |
paste("RIC =",quantile(m,0.75,names=F)-quantile(m,0.25,names=F)) | |
#RANGO SEMIINTERCUARTILAR | |
paste("RSI =",as.vector(diff(summary(m)[c(2,5)])/2)) | |
paste("RSI =",(quantile(m,0.75,names=F)-quantile(m,0.25,names=F))/2) #FORMA MÁS LARGA | |
#CUARTIL MEDIO | |
paste("CM =",as.vector(mean(summary(m)[c(2,5)]))) | |
paste("CM =",(quantile(m,0.75,names=F)+quantile(m,0.25,names=F))/2) #FORMA MÁS LARGA | |
#RANGO DE PERCENTILES 10-90 | |
paste("R10-90 =",quantile(m,0.9,names=F)-quantile(m,0.1,names=F)) | |
#DIAGRAMA DE CUADRO O DE CAJA. PÁGINA 124. EN EL LIBRO SE LE DENOMINA "GRÁFICO DE CUADRO MODIFICADO", PORQUE USA UMBRALES PARA CONSIDERAR VALORES EXTREMOS | |
exsalud <- read.csv('http://geografiafisica.org/sem_2015_02/geo301/triola_2009_apendice_conjunto_datos_1_examen_salud.csv') | |
exsalud$Sexo <- factor(exsalud$Sexo,levels=c('mujer','hombre')) | |
par(mfrow=c(1,1)) | |
boxplot(Pulso~Sexo,exsalud,horizontal=T) | |
#DIAGRAMAS DE CUADRO CON ACTRICES | |
par(mfrow=c(1,1)) | |
boxplot(m,h,names=c('actrices','actores'),horizontal=T) | |
#VARIOS ESTADÍSTICOS DESCRIPTIVOS, ELIGIENDO UNO A UNO LOS QUE INTERESAN. REQUIERE TENER CARGADOS LOS PAQUETES moments, dplyr Y tidyr. ES UNA BUENA SENTENCIA PARA RESUMIR ESTADÍSTICOS. ABAJO SE MUESTRA UNA FORMA CORTA, PERO CON OTROS ESTADÍSTICOS | |
clientes$banco <- row.names(clientes) | |
clientes <- clientes[,c(4,1:3)] | |
as.data.frame(clientes %>% gather(cliente,tiempo,-banco) %>% do(na.omit(.)) %>% group_by(banco) %>% summarise(cantidad=n(),min=min(tiempo),max=max(tiempo),media=mean(tiempo),mediana=median(tiempo),de=sd(tiempo),mitad_rango=mean(range(tiempo)),sesgo=skewness(tiempo),rango=diff(range(tiempo)),cv=(sd(tiempo)/mean(tiempo))*100)) | |
#VARIOS ESTADÍSTICOS DESCRIPTIVOS PREDETERMINADOS DE LA FUNCIÓN describe DE psych. REQUIERE TENER CARGADOS LOS PAQUETES psych, moments, dplyr Y tidyr. OTRA BUENA SENTENCIA PARA RESUMIR ESTADÍSTICOS | |
clientes %>% gather(cliente,tiempo,-banco) %>% group_by(banco) %>% do(describe(.$tiempo)) %>% select(-vars) #NOTAR QUE EL SESGO NO ES IGUAL ENTRE ESTA SENTENCIA Y LA ANTERIOR | |
#VARIOS ESTADÍSTICOS DESCRIPTIVOS PREDETERMINADOS DE LA FUNCIÓN describe DE psych Y summary. REQUIERE TENER CARGADOS LOS PAQUETES psych, moments, dplyr Y tidyr. OTRA BUENA SENTENCIA PARA RESUMIR ESTADÍSTICOS | |
clientes %>% gather(cliente,tiempo,-banco) %>% group_by(banco) %>% do(data.frame(psych::describe(.$tiempo),rbind(summary(.$tiempo)),check.names = F)) %>% dplyr::select(-vars) | |
#VARIOS ESTADÍSTICOS CON sapply | |
sapply(clientes[sapply(clientes,is.numeric)],function(x) {c(media=mean(x,na.rm=T),mediana=median(x,na.rm=T))}) #SE PUEDEN AGREGAR OTROS ESTADÍSTICOS | |
t(sapply(clientes[sapply(clientes,is.numeric)],function(x) {c(media=mean(x,na.rm=T),mediana=median(x,na.rm=T))})) #RESULTADO ANTERIOR TRANSPUESTO | |
#CAPÍTULO 6. DISTRIBUCIÓN DE PROBABILIDAD NORMAL | |
#USO DE LAS TABLAS DE PROBABILIDAD, TABLA A2 | |
#PÁGINA 251. ÁREA ACUMULADA DESDE LA IZQUIERDA DE LA PUNTUACIÓN z=1.58, EQUIVALENTE A "MENOR QUE 1.58" | |
#ÁREAS SOMBREADAS BASADAS EN: http://www.r-bloggers.com/creating-shaded-areas-in-r/ | |
#HAY CÓDIGOS QUE HACEN DE ESTOS GRÁFICOS DE ÁREAS SOMBREADAS MÁS "SEXY". VER: http://www.konstantinkashin.com/blog/2012/04/10/ggplot2-density-plots-and-histograms/ | |
pnorm(1.58) | |
par(mfrow=c(1,1)) | |
curve(dnorm(x,0,1),xlim=c(-3,3)) | |
polygon(c(-4,seq(-4,1.58,0.01),1.58),c(0,dnorm(seq(-4,1.58,0.01)),0),col='grey') | |
#PÁGINA 252. ÁREA A LA DERECHA A PARTIR DE UNA PUNTUACIÓN z | |
pnorm(-1.23,lower.tail = F) | |
par(mfrow=c(1,1)) | |
curve(dnorm(x,0,1),xlim=c(-3,3)) | |
polygon(c(-1.23,seq(-1.23,4,0.01),4),c(0,dnorm(seq(-1.23,4,0.01)),0),col='grey') | |
#PÁGINA 255. PUNTUACIÓN z A PARTIR DE UN ÁREA | |
qnorm(0.95) | |
par(mfrow=c(1,1)) | |
curve(dnorm(x,0,1),xlim=c(-3,3)) | |
polygon(c(-4, seq(-4,qnorm(0.95),0.01),qnorm(0.95)),c(0,dnorm(seq(-4,qnorm(0.95),0.01)),0),col='grey') | |
#ENTRE UN VALOR Y OTRO | |
pnorm(1.5)-pnorm(-2) #PÁGINA 252 | |
par(mfrow=c(1,1)) | |
curve(dnorm(x,0,1),xlim=c(-3,3)) | |
polygon(c(-2,seq(-2,1.5,0.01),1.5),c(0,dnorm(seq(-2,1.5,0.01)),0),col='grey') | |
#LA DISTRIBUCIÓN NORMAL COMO APROXIMACIÓN DE LA DISTRIBUCIÓN BINOMIAL. PÁGINA 291 Y SIGUIENTES (TEORÍA), PÁGINAS 293 Y SIGUIENTES (EJEMPLO) | |
#PRIMERO, VERIFICAR REQUISITOS (VER PÁGINA 291) | |
n=213 | |
x=122-0.5 #LA PROBABILIDAD DE QUE HAYA AL MENOS 122 PASAJEROS, ES EL ÁREA A LA DERECHA DE LA PUNTUACIÓN z CORRESPONDIENTE A ESTA CANTIDAD DE PASAJEROS MENOS 0.5 | |
p=0.5 | |
q=1-p | |
n*p>5 | |
n*q>5 | |
mu=n*p | |
mu | |
sigma=sqrt(n*p*q) | |
sigma | |
z=(x-mu)/sigma | |
z | |
pnorm(z,lower.tail = F) | |
#CAPÍTULO 7. ESTIMACIONES Y TAMAÑOS DE MUESTRA | |
#UNA NIÑA DE 9 AÑOS CONDUCE UN EXPERIMENTO (TRIOLA, 2009, PÁGINA 318). CONSULTA A 280 TERAPEUTAS DE CONTACTO PARA VER SI PUEDEN DESCUBRIR QUÉ MANO TIENE PUESTA SOBRE UNA MESA QUE SE OCULTA TRAS UNA MAMPARA. DE ESTOS, 123 ACIERTAN. CALCULAR EL INTERVALO DE CONFIANZA USANDO LA DISTRIBUCIÓN NORMAL COMO APROXIMACIÓN DE LA BINOMIAL PARA DISTINTOS NIVELES DE CONFIANZA | |
n=280 #MUESTRA | |
n | |
pm=123/280 #PROPORCIÓN DE ACIERTOS | |
pm | |
#SE FIJA EL NIVEL DE CONFIANZA EN 95%, LO CUAL FIJA EL NIVEL DE SIGNIFICANCIA EN 0.05. | |
#PARA OBTENER EL INTERVALO DE CONFIANZA, ES PRECISO PRIMERO CALCULAR EL ERROR: error=z*(sqrt(pm*qm/n)) | |
#DE LA FÓRMULA ANTERIOR SE DISPONE DE pm (PROPORCIÓN MUESTRAL DE ACIERTOS), qm (PROPORCIÓN MUESTRAL DE DESACIERTOS) Y n. p EQUIVALE A LA PROPORCIÓN DE ACIERTOS, Y qm=1-pm | |
qm=1-pm | |
#z SE OBTIENE A PARTIR DE DIVIDIR EL NIVEL DE SIGNIFICANCIA ENTRE 2 (0.05/2=0.025). PARA ESTE VALOR SE LOCALIZA, EN UNA TABLA PUNTUACIONES z Y ÁREAS DESDE LA IZQUIERDA (TABLA A-2 DE TRIOLA 2009), EL VALOR DE p (PROBABILIDAD) CORRESPIENTE A DICHO NIVEL DE SIGNIFICANCIA, QUE ES p=0.975, Y SE LOCALIZA LA PUNTUACIÓN z=1.96 | |
nc='95%' | |
z=1.96 | |
error=z*(sqrt(pm*qm/n)) | |
error | |
ui=pm-error | |
ui | |
us=pm+error | |
us | |
paste("intervalo de confianza para NC de ",nc,": ",round(ui,2)," < ",round(pm,2)," < ",round(us,2),sep="") #PÁGINA 327 | |
#SI SE HUBIESE TENIDO UNA MUESTRA MAYOR, POR EJEMPLO, DE 3000 PERSONAS, ENTONCES EL ERROR Y EL INTERVALO DE CONFIANZA SERÍA EL SIGUIENTE, PARA UN NIVEL DE CONFIANZA DE 95%: | |
n=3000 | |
error=z*(sqrt(pm*qm/n)) | |
error | |
ui=pm-error | |
ui | |
us=pm+error | |
us | |
paste("intervalo de confianza para NC de ",nc,": ",round(ui,2)," < ",round(pm,2)," < ",round(us,2),sep="") | |
#SI AUMENTAMOS EL NIVEL DE CONFIANZA A 99%, EL VALOR DE z CORRESPONDIENTE SERÍA z=2.575, PARA n=280 | |
nc='99%' | |
z=2.575 | |
n=280 | |
error=z*(sqrt(pm*qm/n)) | |
error | |
ui=pm-error | |
ui | |
us=pm+error | |
us | |
paste("intervalo de confianza para NC de ",nc,": ",round(ui,2)," < ",round(pm,2)," < ",round(us,2),sep="") | |
#Y PARA UNA MUESTRA DE n=3000 ELEMENTOS | |
nc='95%' | |
z=2.575 | |
n=3000 | |
error=z*(sqrt(pm*qm/n)) | |
error | |
ui=pm-error | |
ui | |
us=pm+error | |
us | |
paste("intervalo de confianza para NC de ",nc,": ",round(ui,2)," < ",round(pm,2)," < ",round(us,2),sep="") | |
#SI QUEREMOS DETERMINAR EL TAMAÑO DE LA MUESTRA NECESARIA FIJANDO UN MARGEN DE ERROR | |
#EN ESTE CASO, SE CREARÁ UNA FUNCIÓN | |
tam <- function(error,z,pm,qm){(z^2)*pm*qm/(error^2)} | |
#PARA UN ERROR DE 3% Y UN NC DE 95% | |
error=0.03 | |
z=1.96 | |
pm=123/280 | |
qm=1-pm | |
tam(error,z,pm,qm) | |
#PARA UN ERROR DE 3% Y UN NC DE 99% | |
error=0.05 | |
z=2.575 | |
tam(error,z,pm,qm) | |
#PÁGINA 337, EJERCICIO 49; SE REQUIERE REALIZAR TAMBIÉN EL EJERCICIO 44. CALCULAR EL ERROR PARA POBLACIÓN FINITA. SE EMPLEA UN FACTOR DE CORRECCIÓN. NOTA: SE REQUIERE HABER CREADO LA FUNCIÓN tam (VER ARRIBA) | |
nc='95%' | |
#EJERCICIO 44 | |
error=0.04 | |
z=qnorm(0.975) | |
pm=0.5 | |
qm=0.5 | |
n=tam(error,z,pm,qm) | |
N=1250 | |
errorpf=error*(sqrt((N-n)/(N-1))) | |
errorpf | |
#CAPÍTULO 8. PRUEBA DE HIPÓTESIS. ES IMPORTANTE ESTUDIAR ESTE CAPÍTULO ÍNTEGRO PARA COMPRENDER BIEN EL SIGNIFICADO DEL USO DEL ÁREA A LA DERECHA, A LA IZQUIERDA, DOS COLAS Y OTROS DETALLES | |
#PÁGINA 388, ASEVERACIONES SOBRE POROPORCIONES: NACIMIENTOS DE NIÑAS. EN ESTE CASO, ES NECESARIO EMPLEAR LA DISTRIBUCIÓN NORMAL COMO APROXIMACIÓN DE LA DISTRIBUCIÓN BINOMIAL. PÁGINA 291 Y SIGUIENTES (TEORÍA), PÁGINAS 293 Y SIGUIENTES (EJEMPLO) | |
n=100 | |
x=52-0.5 #LA PROBABILIDAD DE QUE NAZCAN AL MENOS 52 NIÑAS, ES EL ÁREA A LA DERECHA DE LA PUNTUACIÓN z CORRESPONDIENTE A 52 NIÑAS DE 100 MENOS 0.5 | |
p=0.5 | |
q=1-p | |
mu=n*p | |
n*p>5 | |
n*q>5 | |
mu | |
sigma=sqrt(n*p*q) | |
sigma | |
z=(x-mu)/sigma | |
z | |
pnorm(z,lower.tail = F) #PROBABILIDAD DE QUE NAZCAN AL MENOS 52 NIÑAS | |
#A LA INVERSA (VER HISTOGRAMA DE LA FIGURA 8-1, PÁGINA 389): PROPORCIÓN DE NIÑAS QUE PODRÍA CONSIDERARSE COMO ALGO NO ATRIBUIBLE AL AZAR, ALGO SIGNIFICATIVAMENTE GRANDE | |
z=qnorm(0.05,lower.tail = F) #ESTA ES LA PUNTUACIÓN z CORRESPONDIENTE A UNA PROBABILIDAD DE 0.05 | |
#DESPEJAR x | |
x=(z*sigma)+mu #EL RESULTADO DA 58.23, QUE EQUIVALE "POR LO MENOS 59 NIÑAS POR CADA 100 NACIMIENTOS" | |
#PÁGINA 392, OTRO EJEMPLO SOBRE PORPORCIONES. EL 61% DE 703 ENCUESTADOS ASEGURA QUE CONSIGUIÓ TRABAJO POR MEDIO DE REDES DE CONTACTOS. DETERMINAR SI ESTE VALOR ES SIGNIFICATIVAMENTE MAYOR QUE LO SE PODRÍA ESPERAR POR AZAR | |
p=0.5 #ESTA ES DE HECHO LA HIPÓTESIS NULA, QUE CONTIENE LA IGUALDAD Y ES LO ESPERADO POR AZAR. LA HIPÓTESIS ALTERNA SERÍA P>0.5 | |
n=703 | |
pm=0.61 #LA PROBABILIDAD DE QUE 61 O MÁS DE CADA 100 PERSONAS CONSIGA TRABAJO POR AZAR, ES EL ÁREA A LA DERECHA DE LA PUNTUACIÓN z CORRESPONDIENTE A ESTA PROPORCIÓN | |
x=n*pm | |
q=1-p | |
#VERIFICAR REQUISITOS (VER PÁGINA 291) | |
n*p>5 | |
n*q>5 | |
#CÁLCULOS | |
mu=n*p | |
mu | |
sigma=sqrt(n*p*q) | |
sigma | |
z=(x-mu)/sigma | |
z | |
pnorm(z,lower.tail = F) #PROBABILIDAD DE QUE 61 O MÁS DE CADA 100 PERSONAS CONSIGA TRABAJO POR AZAR, LO CUAL REVELA QUE LAS REDES DE CONTACTOS SÍ FUNCIONAN | |
#NÓTESE QUE EN ESTE CASO SE HA USADO UNA VARIANTE DE LA FÓRMULA ORIGINAL PARA EL ESTADÍSTICO DE PRUEBA DE PROPORCIONES (PÁGINA 391). LA FÓRMULA YA EMPLEADA ARRIBA EN EL CASO DEL EJEMPLO DE LA PÁGINA 388, ES LA MISMA QUE LA EMPLEADA AQUÍ. SE COMPRUEBA A CONTINUACIÓN: | |
p=0.5 #ESTA ES DE HECHO LA HIPÓTESIS NULA, QUE CONTIENE LA IGUALDAD Y ES LO ESPERADO POR AZAR. LA HIPÓTESIS ALTERNA SERÍA P>0.5 | |
n=703 | |
pm=0.61 #LA PROBABILIDAD DE QUE 61 O MÁS DE CADA 100 PERSONAS CONSIGA TRABAJO POR AZAR, ES EL ÁREA A LA DERECHA DE LA PUNTUACIÓN z CORRESPONDIENTE A ESTA PROPORCIÓN | |
q=1-p | |
z=(pm-p)/sqrt(p*q/n) | |
z #SE VERIFICA QUE ES EL MISMO VALOR | |
z>qnorm(0.01,lower.tail = F) #SE COMPRUEBA QUE z ES MAYOR QUE LA PUNTUACIÓN CRÍTICA PARA UN NIVEL DE SIGNIFICANCIA DE INCLUSO 0.01, QUE ES ca. 2.33 | |
#LOS FLUJOGRAMAS QUE APARECEN EN LAS PÁGINAS SIGUIENTES SON MUY ÚTILES. SE RECOMIENDA ESTUDIARLOS CON DETENIMIENTO Y EJECUTAR LOS EJEMPLOS | |
#CAPÍTULO 9. INFERENCIAS A PARTIR DE DOS MUESTRAS | |
#INFERENCIAS A PARTIR DE PROPORCIONES. PÁGINA 458 Y SIGUIENTES. ¿LA CIRUGIA ES MEJOR QUE EL ENTABLILLADO? LA HIPÓTESIS NULA EN ESTE CASO ES: LA PROPORCIÓN DE ÉXTIOS CON LA CIRUGÍA ES MAYOR QUE CON EL ENTABLILLADO | |
ex1=67 #TASA DE ÉXITOS DE CIRUGIA | |
tr1=73 #TOTAL DE SUJETOS TRATADOS POR CIRUGIA | |
ex2=60 #TASA DE ÉXITOS DE ENTABLILLADO | |
tr2=83 #TOTAL DE SUJETOS TRATADOS POR ENTABLILLADO | |
pma=(ex1+ex2)/(tr1+tr2) #PROPORCIÓN MUESTRAL DE ÉXITOS AGRUPADA | |
pma | |
qma=1-pma #PROPORCIÓN MUESTRAL DE FRACASOS AGRUPADA | |
qma | |
pr1=ex1/tr1 #PROPORCIÓN DE ÉXITOS DE CIRUGÍA | |
pr2=ex2/tr2 #PROPORCIÓN DE ÉXITOS DE ENTABLILLADO | |
z=(pr1-pr2)/sqrt((pma*qma/tr1)+(pma*qma/tr2)) | |
z | |
pnorm(z,lower.tail = F) #SE RECHAZA LA HIPÓTESIS NULA, POR TRATARSE DE UN VALOR DE p MENOR A 0.05 | |
#INFERENCIAS A PARTIR DE DOS MEDIAS DE MUESTRAS INDEPENDIENTES. PÁGINA 471 Y SIGUIENTES. | |
#ESTUDIAR LA VERIFICACIÓN DE REQUISITOS (PÁGINA 471) | |
#HIPÓTESIS NULA: LOS SOLICITANTES SIN ÉXITO PROVIENEN DE UNA POBLACIÓN CUYA MEDIA DE EDAD ES IGUAL A LA DE LOS SOLICITANTES CON ÉXITO | |
#HIPÓTESIS ALTERNA: LOS SOLICITANTES SIN ÉXITO PROVIENEN DE UNA POBLACIÓN CUYA MEDIA DE EDAD ES MAYOR QUE LA DE LOS SOLICITANTES CON ÉXITO | |
#CARGAR DATOS | |
edadessinexito <- c(34,37,37,38,41,42,43,44,44,45,45,45,46,48,49,53,53,54,54,55,56,57,60) | |
edadesconexito <- c(27,33,36,37,38,38,39,42,42,43,43,44,44,44,45,45,45,45,46,46,47,47,48,48,49,49,51,51,52,54) | |
#GRÁFICOS Y ESTADÍSTICOS DESCRIPTIVOS | |
par(mfrow=c(1,2)) | |
hist(edadessinexito) | |
hist(edadesconexito) | |
qqnorm(edadessinexito) | |
qqnorm(edadesconexito) | |
estd <- function(x) {c(n=length(x),media=mean(x,na.rm=T),mediana=median(x,na.rm=T),sd=sd(x),`p prueba ST`=shapiro.test(x)$p.value)} | |
data.frame(`Sin éxito`=estd(edadessinexito),check.names = F) | |
data.frame(`Con éxito`=estd(edadesconexito),check.names = F) | |
#APLICACIÓN DE LA PRUEBA | |
#PRIMERO: USANDO FÓRMULAS DE TRIOLA (2009), ESTADÍSTICO DE PRUEBA Y VALOR DE P RESULTANTE. LOS GRADOS DE LIBERTAD EQUIVALEN AL TAMAÑO DE LA MUESTRA MÁS PEQUEÑO MENOS 1 | |
gl=length(edadessinexito)-1 | |
t=(mean(edadessinexito)-mean(edadesconexito))/sqrt((var(edadessinexito)/length(edadessinexito))+(var(edadesconexito)/length(edadesconexito))) | |
t | |
p=pt(t,df=gl,lower.tail = F) | |
p | |
#SEGUNDO: FUNCIÓN INCORPORADA EN PAQUETE stats DE R. NÓTESE LA DIFERENCIA EN EL RESULTADO RESPECTO DE LA IMPLEMENTACIÓN ANTERIOR | |
t.test(edadessinexito,edadesconexito,alternative = 'greater') #NOTAR QUE EL VALOR DEL ESTADÍSTICO DE PRUEBA (t) ES IGUAL AL DE LA IMPLEMENTACIÓN MANUAL, PERO HAY DIFERENCIAS EN EL VALOR DE P. ESTO OCURRE PORQUE ESTA FUNCIÓN (t.test) DEL PAQUETE STATS CALCULA LOS GRADOS DE LIBERTAD MEDIANTE UNA FÓRMULA DIFERENTE, A LA QUE TRIOLA DENOMINA "Fórmula 9-1" | |
#TERCERO: CREANDO UNA FUNCIÓN PERSONALIZADA. ÉSTA UNE LAS DOS MODALIDADES ANTERIORES, Y SÓLO ES NECESARIO ESPECIFICIAR LAS MUESTRAS | |
#UNA FUNCIÓN PARA REALIZAR LA PRUEBA, TANTO MEDIANTE LA FÓRMULA COMO MEDIANTE LA FUNCIÓN INCORPORADA DE R, QUE INCLUYE CORRECCIONES | |
#PROGRAMACIÓN DE LA FUNCIÓN | |
pr.test.per <- function(m1,m2){ | |
estd <- function(x) {c(n=length(x),media=mean(x,na.rm=T),mediana=median(x,na.rm=T),sd=sd(x),`p prueba ST`=shapiro.test(x)$p.value)} | |
cat("\nEstadísticos\n") | |
print(data.frame(`Sin éxito`=estd(m1),`Con éxito`=estd(m2),check.names = F)) | |
gl=min(length(m1),length(m2))-1 | |
t=(mean(m1)-mean(m2))/sqrt((var(m1)/length(m1))+(var(m2)/length(m2))) | |
p=pt(t,df=gl,lower.tail = F) | |
cat("\nSegún fórmula original\n","Estadístico de prueba z = ",t,"\nvalor de p = ",p,sep="") | |
cat("\n\nPrueba T de Student incorporada") | |
t.test(m1,m2,alternative='greater') | |
} | |
#EJECUCIÓN | |
pr.test.per(edadessinexito,edadesconexito) | |
#TERCERA PARTE: MANIPULACIÓN Y ANÁLISIS DE DATOS PROPIOS (APLICACIÓN DE LO ANTERIOR) | |
#SE DISPONE DE MEDICIONES E IDENTIFICACIONES DE CLASTOS DEL RÍO OCOA A LA ALTURA DE LA LOCALIDAD "LOS PILONES". CADA MUESTRA SE COMPONE DE 100 ELEMENTOS, Y ÉSTAS PUEDEN SER DE TRES TIPOS: A) RODAL CON IDENTIFICACIÓN LITOLÓGICA; B) RODAL SIN IDENTIFICACIÓN LITOLÓGICA; C) TRANSECTO WOLMANN | |
#REALIZAR LOS SIGUIENTES EJERCICIOS: | |
#1) MAPA CON LA LOCALIZACIÓN DE UN PUNTO QUE REPRESENTE A CADA MUESTRA. EN EL CASO DE LOS RODALES, SÓLO EXISTE UNA COORDENADA, POR LO QUE SE UTILIZA ÉSTA. EN EL CASO DE LOS TRANSECTOS, BASTA CON COLOCAR LA COORDENADA INICIAL | |
#2) ANÁLISIS EXPLORATORIO DE DATOS (AED), QUE INCLUYE: | |
#2.0) TABLA DE LUGAR(=CÓDIGO)/TIPO DE MUESTRA/UNIDAD GEOMORFOLÓGICA | |
#2.1) TABLA DE FRECUENCIAS SEGÚN: LUGAR DE LA MUESTRA, TIPO DE MUESTRA (RODAL CON Y SIN IDENTIFICACIÓN, TRANSECTO), TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#2.2) MEDIAS DE CADA EJE | |
#2.2.1) MEDIAS DE CADA EJE DE LAS MUESTRAS LPRO_01 Y LPRO_02 SEGÚN TIPO DE ROCA | |
#2.2.2) CON DATOS COMBINADOS (POOLED), USANDO UN ÚNICO CRITERIO A LA VEZ, CALCULAR LAS MEDIAS DE CADA EJE SEGÚN: LUGAR DE MUESTRA, TIPO DE MUESTRA (RODAL CON Y SIN IDENTIFICACIÓN, TRANSECTO), TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#2.3) VALORES EXTREMOS | |
#2.3.1) VALORES EXTREMOS DE CADA EJE DE LAS MUESTRAS LPRO_01 Y LPRO_02 SEGÚN TIPO DE ROCA | |
#2.3.2) CON DATOS COMBINADOS (POOLED), USANDO UN ÚNICO CRITERIO A LA VEZ, CALCULAR LOS VALORES MÁXIMO Y MÍNIMO DE CADA EJE SEGÚN: LUGAR DE MUESTRA, TIPO DE MUESTRA, TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#3) ANÁLISIS EXPLORATORIO DE DATOS (AED) GRÁFICAMENTE, QUE INCLUYE: | |
#3.1) GRÁFICOS DE BARRA SEGÚN FRECUENCIAS DE 2.1 | |
#3.2) HISTOGRAMA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN LUGAR DE MUESTRA | |
#3.3) HISTOGRAMA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN UNIDAD GEOMORFOLÓGICA | |
#3.4) DIAGRAMAS DE CAJA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN LUGAR DE MUESTRA | |
#3.5) DIAGRAMAS DE CAJA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN UNIDAD GEOMORFOLÓGICA | |
#4) PARA LAS MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03, CALCULAR EL MARGEN DE ERROR Y EL INTERVALO DE CONFIANZA, PARA UN 95% DE NIVEL DE CONFIANZA DEL EJE a | |
#5) COMPARAR LAS PROPORCIONES DE MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03. HACER LO MISMO CON LAS CALIZAS DE LAS MUESTRAS LPRO_01 Y LPRO_02. ¿QUÉ CONCLUSIONES SE OBTIENEN DE ESTE RESULTADO? | |
#6) COMPARAR MEDIAS DEL EJE a DE MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03. HACER LO MISMO CON LAS CALIZAS DE LAS MUESTRAS LPRO_01 Y LPRO_02. ¿QUÉ CONCLUSIONES SE OBTIENEN DE ESTE RESULTADO? | |
#PASOS COMUNES A TODAS LAS RESPUESTAS: | |
#LECTURA DE LOS DATOS DESAGREGADOS: UNA TABLA DONDE CADA FILA CONSTITUYE UN CLASTO MEDIDO EN SUS TRES EJES, a, b Y c. EN ALGUNOS CASOS, LOS EJES a Y c NO ESTÁN PRESENTES. ALGUNOS MUESTREOS INCLUYERON IDENTIFICACIÓN DEL TIPO DE ROCA | |
d<-read.csv('http://geografiafisica.org/sem_2015_02/geo301/integrados_ok.csv') | |
head(d) | |
tail(d) | |
str(d) #LA ESTRUCTURA ¡DEBE SER ESTUDIADA DETENIDAMENTE! NO SE TRATA DE MANIPULAR DATOS CUYO CONTENIDO NO SE CONOCE | |
d.env<-read.csv('http://geografiafisica.org/sem_2015_02/geo301/integrados_env_ok.csv') #DESCARGA Y CONVIERTE A UNA TABLA (data.frame) EL ARCHIVO CONTENIENDO LA INFORMACIÓN AMBIENTAL DE LOS MUESTREOS | |
head(d.env) | |
tail(d.env) | |
str(d.env) | |
#RESPUESTAS | |
#1) MAPA CON LA LOCALIZACIÓN DE UN PUNTO QUE REPRESENTE A CADA MUESTRA. EN EL CASO DE LOS RODALES, SÓLO EXISTE UNA COORDENADA, POR LO QUE SE UTILIZA ÉSTA. EN EL CASO DE LOS TRANSECTOS, BASTA CON COLOCAR LA COORDENADA INICIAL | |
#VÍA googlemaps.com | |
d.env.sp <- d.env | |
coordinates(d.env.sp)<-~x_inicial+y_inicial #INDICA AL PROGRAMA QUE LAS COORDENADAS DEL OBJETO t.env.ok ESTÁN EN LOS CAMPOS x_inicial E y_inicial. ESTO CONVERTIRÁ A d.env.ok EN UN OBJETO ESPACIAL | |
proj4string(d.env.sp) <- CRS("+init=epsg:32619") #LE ASIGNA EL SISTEMA DE REFERENCIA DE COORDENADAS CORRESPONDIENTE AL OBJETO | |
plotGoogleMaps(d.env.sp,iconMarker='',zcol="codigo","cantometria.html") #GENERA UN ARCHIVO .html EN LA CARPETA DE TRABAJO, LA CUAL SE SABE ESCRIBIENDO getwd(). EL PROGRAMA LANZARÁ EL NAVEGADOR DE MANERA AUTOMÁTICA, Y CARGARÁ LOS PUNTOS. NOTA: SI R NO TIENE PERMISOS DE ESCRITURA EN EL DIRECTORIO DE TRABAJO, NO PODRÁ CREAR EL ARCHIVO Y DARÁ ERROR. EN ESE CASO, SE RECOMIENDA FIJAR UNA CARPETA DE TRABAJO EN UN DIRECTORIO CON TODOS LOS PRIVILEGIOS, USANDO EL COMANDO setwd(CARPETA). POR EJEMPLO, setwd("C:/"), O setwd("D:/"), O setwd("C:/Users/Public") | |
#KML DE GoogleEarth | |
d.env.sp.g<-reproject(d.env.sp) #REPROYECTA EL OBJETO d.env.ok A COORDENADAS GEOGRÁFICAS | |
kmlPoints(d.env.sp.g,kmlfile="d.env.sp.g.kml",name=d.env.sp.g$codigo) #ESTO GUARDARÁ UN ARCHIVO kml EN EL DIRECTORIO DE TRABAJO; EN DICHO DIRECTORIO R DEBE TENER PERMISOS DE ESCRITURA, PORQUE DE LO CONTRARIO DARA ERROR -SE PUEDE SABER CUÁL ES EL DIRECTORIO DE TRABAJO ESCRIBIENDO getwd(); TAMBIÉN SE PUEDE CREAR UNA CARPETA, MEDIANTE EL EXPLORADOR DE WINDOWS, Y LUEGO FIJARLA COMO DIRECTORIO DE TRABAJO EN R CON setwd("RUTA"), DONDE "RUTA" SERÁ LA LOCALIZACIÓN DE LA CARPETA CREADA, POR EJEMPLO "C:/PRACTICAS_CON_R")-. EL ARCHIVO html CREADO SE ABRIRÁ AUTOMÁTICAMENTE EN UNA VENTANA DEL NAVEGADOR PREDETERMINADO. EL NOMBRE DE CADA MARCADOR SERÁ OBTENIDO DE LA COLUMNA "codigo" | |
#HASTA AQUÍ LA REPRESENTACIÓN ESPACIAL DE LOS PUNTOS DE MUESTREO CON COORDENADAS VALIDADAS | |
#OTRAS OPCIONES SON FUNCIONES DE LOS PAQUETES leafletr Y ggmap | |
#2) ANÁLISIS EXPLORATORIO DE DATOS (AED), QUE INCLUYE: | |
#2.0) TABLA DE LUGAR(=CÓDIGO)/TIPO DE MUESTRA | |
d.env[,c('codigo','tipo_de_muestreo','unidad_geom')] | |
#2.1) TABLA DE FRECUENCIAS SEGÚN: LUGAR DE LA MUESTRA, TIPO DE MUESTRA (RODAL CON Y SIN IDENTIFICACIÓN, TRANSECTO), TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#FREQ POR LUGARES DE CADA MUESTRA | |
table(d$codigo) | |
#FREQ POR TIPOS DE ROCAS | |
table(d$tipo) #CON CÓDIGOS DE ROCAS NUMÉRICOS | |
table(d$tipot) #CON NOMBRES DE ROCAS | |
#FREQ POR TIPO DE MUESTRA | |
table(d.env$tipo_de_muestreo) | |
#FREQ POR UNIDAD GEOMORFOLÓGICA | |
table(d.env$unidad_geom) | |
#USANDO CÓDIGOS SE REALIZAN LAS TABLAS DE FRECUENCIAS ANTERIORES (Y MÁS) CON MAYOR EFICIENCIA | |
#TABLA DE FRECUENCIAS DE TODAS LAS COLUMNAS FACTORIALES DE LA TABLA DE CLASTOS (d) | |
sapply(d[sapply(d,is.factor)],table) | |
#TABLA DE FRECUENCIAS DE TODAS LAS COLUMNAS FACTORIALES DE LA TABLA DE INFORMACIÓN AMBIENTAL (d) | |
sapply(d.env[sapply(d.env,is.factor)],table) | |
#2.2) MEDIAS DE CADA EJE | |
#2.2.1) MEDIAS DE CADA EJE DE LAS MUESTRAS LPRO_01 Y LPRO_02 SEGÚN TIPO DE ROCA | |
#APLICANDO FUERZA BRUTA, CON tapply, CADA TIPO DE ROCA POR CADA EJE | |
tapply(d[d$codigo=='LPRO_01','a'],droplevels(d[d$codigo=='LPRO_01','tipot']), function(x) mean(x,na.rm=T)) | |
tapply(d[d$codigo=='LPRO_01','b'],droplevels(d[d$codigo=='LPRO_01','tipot']), function(x) mean(x,na.rm=T)) | |
tapply(d[d$codigo=='LPRO_01','c'],droplevels(d[d$codigo=='LPRO_01','tipot']), function(x) mean(x,na.rm=T)) | |
tapply(d[d$codigo=='LPRO_02','a'],droplevels(d[d$codigo=='LPRO_02','tipot']), function(x) mean(x,na.rm=T)) | |
tapply(d[d$codigo=='LPRO_02','b'],droplevels(d[d$codigo=='LPRO_02','tipot']), function(x) mean(x,na.rm=T)) | |
tapply(d[d$codigo=='LPRO_02','c'],droplevels(d[d$codigo=='LPRO_02','tipot']), function(x) mean(x,na.rm=T)) | |
#APLICANDO sapply, CON UNA SOLA LÍNEA DE CÓDIGO, PERO CON SINTÁXIS COMPLEJA | |
sapply(c('LPRO_01','LPRO_02'),function(x) sapply(levels(droplevels(d[d$codigo==x,'tipot'])), function(y) sapply(d[d$codigo==x&d$tipot==y,c('a','b','c')], function(z) mean(z,na.rm=T))),simplify=F) | |
#EN UNA ÚNICA TABLA MÁS SIMPLIFICADA | |
t(as.data.frame(sapply(c('LPRO_01','LPRO_02'),function(x) sapply(levels(droplevels(d[d$codigo==x,'tipot'])), function(y) sapply(d[d$codigo==x&d$tipot==y,c('a','b','c')], function(z) mean(z,na.rm=T))),simplify=F))) | |
#CON dplyr MUUUUUUCHO MÁS SENCILLO | |
d %>% filter(grepl("LPRO_01|LPRO_02",codigo)) %>% group_by(codigo,tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(mean(.,na.rm=T))) | |
#MÁS ESTADÍSTICOS CON dplyr | |
d %>% filter(grepl("LPRO_01|LPRO_02",codigo)) %>% group_by(codigo,tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(mean(.,na.rm=T),median(.,na.rm=T),length,sd)) | |
#TAMBIÉN SE PUEDEN ORDENAR, EN ESTE CASO, DESCENDENTEMENTE POR EL EJE a | |
d %>% filter(grepl("LPRO_01|LPRO_02",codigo)) %>% group_by(codigo,tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(mean(.,na.rm=T))) %>% arrange(desc(a)) | |
#2.2.2) CON DATOS COMBINADOS (POOLED), USANDO UN ÚNICO CRITERIO A LA VEZ, CALCULAR LAS MEDIAS DE CADA EJE SEGÚN: LUGAR DE MUESTRA, TIPO DE MUESTRA (RODAL CON Y SIN IDENTIFICACIÓN, TRANSECTO), TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#MEDIAS POR LUGARES DE MUESTRA | |
t(sapply(levels(d$codigo),function(y) sapply(d[d$codigo==y,c('a','b','c')],function(x) c(media=round(mean(x,na.rm=T),2))))) | |
#CON dplyr | |
d %>% group_by(codigo) %>% dplyr::select(a,b,c) %>% summarise_each(funs(round(mean(.,na.rm=T),2))) | |
#MEDIAS POR TIPO DE ROCA | |
sapply(levels(d$tipot),function(y) sapply(d[d$tipot==y,c('a','b','c')],function(x) c(media=round(mean(x,na.rm=T),2)))) | |
#CON dplyr | |
d %>% group_by(tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(round(mean(.,na.rm=T),2))) | |
#MEDIAS POR TIPO DE MUESTRA | |
duni <- merge(d,d.env,by='codigo') | |
sapply(levels(duni$tipo_de_muestreo),function(y) sapply(duni[duni$tipo_de_muestreo==y,c('a','b','c')],function(x) c(media=round(mean(x,na.rm=T),2)))) | |
#CON dplyr | |
inner_join(d,d.env,by='codigo') %>% group_by(tipo_de_muestreo) %>% dplyr::select(a,b,c) %>% summarise_each(funs(round(mean(.,na.rm=T),2))) | |
#MEDIAS POR UNIDAD GEOMORFOLÓGICA | |
sapply(levels(duni$unidad_geom),function(y) sapply(duni[duni$unidad_geom==y,c('a','b','c')],function(x) c(media=round(mean(x,na.rm=T),2)))) | |
#CON dplyr | |
inner_join(d,d.env,by='codigo') %>% group_by(unidad_geom) %>% dplyr::select(a,b,c) %>% summarise_each(funs(round(mean(.,na.rm=T),2))) | |
#2.3) VALORES EXTREMOS | |
#2.3.1) VALORES EXTREMOS DE CADA EJE DE LAS MUESTRAS LPRO_01 Y LPRO_02 SEGÚN TIPO DE ROCA | |
#SÓLO CON dplyr | |
d %>% filter(grepl("LPRO_01|LPRO_02",codigo)) %>% group_by(codigo,tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(min,max)) | |
#2.3.2) CON DATOS COMBINADOS (POOLED), USANDO UN ÚNICO CRITERIO A LA VEZ, CALCULAR LOS VALORES MÁXIMO Y MÍNIMO DE CADA EJE SEGÚN: LUGAR DE MUESTRA, TIPO DE MUESTRA, TIPO DE ROCA, UNIDAD GEOMORFOLÓGICA | |
#SÓLO CON dplyr | |
#MEDIAS POR LUGARES DE MUESTRA | |
d %>% group_by(codigo) %>% dplyr::select(a,b,c) %>% summarise_each(funs(min,max)) | |
#MEDIAS POR TIPO DE ROCA | |
d %>% group_by(tipot) %>% dplyr::select(a,b,c) %>% summarise_each(funs(min,max)) | |
#MEDIAS POR TIPO DE MUESTRA | |
inner_join(d,d.env,by='codigo') %>% group_by(tipo_de_muestreo) %>% dplyr::select(a,b,c) %>% summarise_each(funs(min,max)) | |
#MEDIAS POR UNIDAD GEOMORFOLÓGICA | |
inner_join(d,d.env,by='codigo') %>% group_by(unidad_geom) %>% dplyr::select(a,b,c) %>% summarise_each(funs(min,max)) | |
#3) ANÁLISIS EXPLORATORIO DE DATOS (AED) GRÁFICAMENTE, QUE INCLUYE: | |
#3.1) GRÁFICOS DE BARRA SEGÚN FRECUENCIAS DE 2.1 | |
#NOTA: TANTO CON LAS FUNCIONES PROPIAS DEL PAQUETE graphics DE R (PAQUETE DE GRÁFICOS POR DEFECTO), Y CON ggplot2 SE PUEDEN VIRTUALMENTE PERONALIZAR CUALQUIER PARTE DE UN GRÁFICO. EN ESTE SCRIPT SÓLO SE MUESTRAN LOS FORMATOS BÁSICOS | |
par(mfrow=c(1,4))#CONFIGURA LA VENTANA GRÁFICA PARA QUE HAYA UN PANEL CON UNA FILA Y CUATRO COLUMNAS | |
par(mar=c(8,2,2,2))#CONFIGURA LA VENTANA GRÁFICA PARA QUE EL MARGEN INFERIOR SEA MAYOR | |
#POR LUGARES DE CADA MUESTRA | |
barplot(table(d$codigo),las='2') | |
#FREQ POR TIPOS DE ROCAS | |
barplot(table(d$tipot),las='2') | |
#FREQ POR TIPO DE MUESTRA | |
text(barplot(table(d.env$tipo_de_muestreo),xaxt='n'),y=-0.6, gsub(" ","\n",names(table(d.env$tipo_de_muestreo))),srt=0,adj = 0.5, xpd = T, font = 1, cex = 1) | |
#FREQ POR UNIDAD GEOMORFOLÓGICA | |
text(barplot(table(d.env$unidad_geom),xaxt='n'),y=-0.6, gsub(" ","\n",names(table(d.env$unidad_geom))),srt=0,adj = 0.5, xpd = T, font = 1, cex = 1) | |
#3.2) HISTOGRAMA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN LUGAR DE MUESTRA | |
#ROCAS REPRESENTADAS POR MÁS DE 10 ELEMENTOS POR LUGAR DE MUESTRA: areniscas, calizas, margas, volcánicas | |
d %>% group_by(codigo,tipot) %>% mutate(n=n()) %>% filter(n>10&!tipot=="NA") %>% summarise(n=n()) %>% ungroup() %>% arrange(tipot) | |
#ASIGNANDO LOS VALORES OBTENIDOS EN EL EJE a phist | |
phistcodtipo <- d %>% group_by(codigo,tipot) %>% mutate(n=n()) %>% filter(n>10&!tipot=="NA") %>% dplyr::select(codigo,tipot,a) %>% gather(eje,valor,-(codigo:tipot)) | |
#TODOS LOS HISTOGRAMAS | |
dev.off() | |
ggplot(phistcodtipo,aes(x=valor)) + geom_histogram(binwidth=5) + xlim(30,80) + facet_wrap(codigo~tipot,scales='free') | |
#3.3) HISTOGRAMA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN UNIDAD GEOMORFOLÓGICA | |
#ÍDEM CASO ANTERIOR | |
inner_join(d,d.env) %>% group_by(unidad_geom,tipot) %>% mutate(n=n()) %>% filter(n>10&!tipot=="NA") %>% summarise(n=n()) %>% ungroup() %>% arrange(tipot) | |
phistugtipo <- inner_join(d,d.env) %>% group_by(unidad_geom,tipot) %>% mutate(n=n()) %>% filter(n>10&!tipot=="NA") %>% dplyr::select(unidad_geom,tipot,a) %>% gather(eje,valor,-(unidad_geom:tipot)) | |
ggplot(phistugtipo,aes(x=valor)) + geom_histogram(binwidth=5) + xlim(30,80) + facet_wrap(unidad_geom~tipot,scales='free') | |
#3.4) DIAGRAMAS DE CAJA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN LUGAR DE MUESTRA | |
ggplot(phistcodtipo,aes(tipot, fill=tipot, valor)) + geom_boxplot() + coord_cartesian(ylim=c(30,80)) + facet_grid(~codigo) | |
#3.5) DIAGRAMAS DE CAJA DEL EJE a DE CADA TIPO DE ROCA REPRESENTADA POR MÁS DE 10 ELEMENTOS SEGÚN UNIDAD GEOMORFOLÓGICA | |
ggplot(phistugtipo,aes(tipot, fill=tipot, valor)) + geom_boxplot() + coord_cartesian(ylim=c(30,80)) + facet_grid(~unidad_geom) | |
#4) PARA LAS MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03, CALCULAR EL MARGEN DE ERROR Y EL INTERVALO DE CONFIANZA, PARA UN 95% DE NIVEL DE CONFIANZA DEL EJE a | |
d %>% dplyr::select(codigo,tipot,a) %>% filter(grepl('LPRO_02|LPRO_03',codigo),grepl('margas/lutitas',tipot)) %>% group_by(codigo) %>% summarise(n=n()) | |
#NO HAY MÁS DE 30 ELEMENTOS EN UNA DE LAS MUESTRAS, PERO ESTÁ BASTANTE CERCA, POR LO QUE SE ELEGIRÍA LA DISTRIBUCIÓN t | |
#VER SELECCIÓN DE DISTRIBUCIÓN (SI t O NORMAL) EN DIAGRAMA DE FLUJO DE LA PÁGINA 354 DE TRIOLA (2009) | |
#VER EJEMPLOS DE APLICACIÓN EN LA PÁGINA 355 DE TRIOLA (2009) | |
#USANDO sapply. SE PUEDE HACER MANUALMENTE, UNO A UNO, USANDO subset PARA AISLAR LAS MARGAS Y CADA MUESTRA, ACUMULANDO MÁS PASOS Y GENERANDO OBJETOS INNECESARIOS | |
sapply(c('LPRO_02','LPRO_03'),function(x) t.test(d[d$codigo==x,'a']),simplify=F) | |
#5) COMPARAR LAS PROPORCIONES DE MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03. HACER LO MISMO CON LAS CALIZAS DE LAS MUESTRAS LPRO_01 Y LPRO_02. ¿QUÉ CONCLUSIONES SE OBTIENEN DE ESTE RESULTADO? | |
#ESTE ES UN EJEMPLO DE APLICACIÓN DEL CAPÍTULO 9 DE TRIOLA, PARECIDO AL MOSTRADO EN PÁGINAS 458 Y SIGUIENTES. | |
#PROPORCIONES DE MARGAS ENTRE LPRO_02 y LPRO_03 | |
d %>% dplyr::select(codigo,tipot,a) %>% filter(grepl('LPRO_02|LPRO_03',codigo),grepl('margas/lutitas',tipot)) %>% group_by(codigo) %>% summarise(n=n()) #CANTIDAD DE MARGAS SEGÚN MUESTRAS | |
d %>% dplyr::select(codigo,tipot,a) %>% filter(grepl('LPRO_02|LPRO_03',codigo)) %>% group_by(codigo) %>% summarise(n=n()) #CANTIDAD DE CLASTOS INDIFERENCIADOS | |
margasLPRO_02=27 | |
totalLPRO_02=100 | |
margasLPRO_03=33 | |
totalLPRO_03=100 | |
pma=(margasLPRO_02+margasLPRO_03)/(totalLPRO_02+totalLPRO_03) #PROPORCIÓN MUESTRAL DE ÉXITOS AGRUPADA | |
pma | |
qma=1-pma #PROPORCIÓN MUESTRAL DE FRACASOS AGRUPADA | |
qma | |
pr1=margasLPRO_02/totalLPRO_02 | |
pr2=margasLPRO_03/totalLPRO_03 | |
z=(pr1-pr2)/sqrt((pma*qma/totalLPRO_02)+(pma*qma/totalLPRO_03)) | |
z | |
pnorm(z,lower.tail = F) #NO SE RECHAZA LA HIPÓTESIS NULA, POR TRATARSE DE UN VALOR DE p MAYOR A 0.05 | |
#PROPORCIONES DE CALIZAS ENTRE LPRO_01 y LPRO_02 | |
d %>% dplyr::select(codigo,tipot,a) %>% filter(grepl('LPRO_01|LPRO_02',codigo),grepl('calizas',tipot)) %>% group_by(codigo) %>% summarise(n=n()) #CANTIDAD DE CALIZAS SEGÚN MUESTRAS | |
d %>% dplyr::select(codigo,tipot,a) %>% filter(grepl('LPRO_01|LPRO_02',codigo)) %>% group_by(codigo) %>% summarise(n=n()) #CANTIDAD DE CLASTOS INDIFERENCIADOS | |
calizasLPRO_01=80 | |
totalLPRO_01=100 | |
calizasLPRO_02=30 | |
totalLPRO_02=100 | |
pma=(calizasLPRO_01+calizasLPRO_02)/(totalLPRO_01+totalLPRO_01) #PROPORCIÓN MUESTRAL DE ÉXITOS AGRUPADA | |
pma | |
qma=1-pma #PROPORCIÓN MUESTRAL DE FRACASOS AGRUPADA | |
qma | |
pr1=calizasLPRO_01/totalLPRO_01 | |
pr2=calizasLPRO_02/totalLPRO_02 | |
z=(pr1-pr2)/sqrt((pma*qma/totalLPRO_01)+(pma*qma/totalLPRO_02)) | |
z | |
pnorm(z,lower.tail = F) #SE RECHAZA LA HIPÓTESIS NULA, POR TRATARSE DE UN VALOR DE p MENOR A 0.05 | |
#6) COMPARAR MEDIAS DEL EJE a DE MARGAS DE LAS MUESTRAS LPRO_02 Y LPRO_03. HACER LO MISMO CON LAS CALIZAS DE LAS MUESTRAS LPRO_01 Y LPRO_02. ¿QUÉ CONCLUSIONES SE OBTIENEN DE ESTE RESULTADO? | |
#AL IGUAL QUE EN CASOS ANTERIORES, SE PUEDEN SEPARAR LAS MUESTRAS, CREAR NUEVOS OBJETOS Y CORRER LA PRUEBA, PERO SE AÑADEN PASOS Y SE CREAN OBJETOS INNECESARIOS | |
#EJE a DE LAS MARGAS DE LPRO_02 Y LPRO_03 | |
#USANDO "[" | |
t.test(a~codigo,d[d$codigo %in% c('LPRO_02','LPRO_03')&d$tipot %in% c('margas/lutitas'),c('codigo','a')]) | |
#USANDO dplyr | |
t.test(a~codigo, d %>% filter(grepl('LPRO_02|LPRO_03',codigo),grepl('margas/lutitas',tipot))) | |
#EJE a DE LAS CALIZAS DE LPRO_01 Y LPRO_02 | |
#USANDO "[" | |
t.test(a~codigo,d[d$codigo %in% c('LPRO_01','LPRO_02')&d$tipot %in% c('calizas'),c('codigo','a')]) | |
#USANDO dplyr | |
t.test(a~codigo, d %>% filter(grepl('LPRO_01|LPRO_02',codigo),grepl('calizas',tipot))) | |
#7) COMPARAR MEDIAS DEL EJE a DE MARGAS ENTRE LAS POSICIONES GEOMORFOLÓGICAS "banco mediano (distal)" Y "banco mediano (proximal)". HACER LO MISMO CON LAS CALIZAS. ¿QUÉ CONCLUSIONES SE OBTIENEN DE ESTE RESULTADO? | |
#EJE a DE LAS MARGAS DE ENTRE UNIDADES GEOMORFOLÓGICAS | |
#USANDO dplyr | |
inner_join(d,d.env) %>% filter(grepl('margas/lutitas',tipot)) %>% group_by(unidad_geom) %>% summarise(n=n()) #CANTIDAD DE CALIZAS SEGÚN UNIDADES GEOMORFOLÓGICAS | |
t.test(a~unidad_geom, inner_join(d,d.env) %>% filter(grepl('margas/lutitas',tipot))) | |
#EJE a DE LAS CALIZAS DE ENTRE UNIDADES GEOMORFOLÓGICAS | |
#USANDO dplyr | |
inner_join(d,d.env) %>% filter(grepl('calizas',tipot)) %>% group_by(unidad_geom) %>% summarise(n=n()) #CANTIDAD DE CALIZAS SEGÚN UNIDADES GEOMORFOLÓGICAS. MUCHA ASIMETRÍA ENTRE LAS MUESTRAS, LO CUAL DESACONSEJA LA APLICACIÓN DE LA PRUEBA | |
t.test(a~unidad_geom, inner_join(d,d.env) %>% filter(grepl('calizas',tipot))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment