#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(ggplot2) #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)))