################################ GUANAJUATO ################################ # Análisis EIC 2015 Guanajuato. # Análisis Exploratorio de Datos remove(list=ls()) # Fijamos directorio de trabajo setwd("~/Documents/Difusion/Curso Prepa 2018") #Cargamos la base de datos sobre población, es un csv datpob = read.csv("TR_PERSONA11.csv",header=TRUE) # 723,696 86 names(datpob) # [1] "ID_VIV" "ID_PERSONA" "ENT" # [4] "NOM_ENT" "MUN" "NOM_MUN" # [7] "LOC50K" "NOM_LOC" "COBERTURA" #[10] "ESTRATO" "UPM" "FACTOR" #[13] "NUMPER" "SEXO" "EDAD" #[16] "PARENT" "PARENT_OTRO_C" "IDENT_MADRE" #[19] "IDENT_PADRE" "SERSALUD" "AFRODES" #[22] "ACTA_NAC" "DHSERSAL1" "DHSERSAL2" #[25] "PERTE_INDIGENA" "ENT_PAIS_NAC" "NACIONALIDAD" #[28] "HLENGUA" "QDIALECT_C" "QDIALECT_INALI" #[31] "HESPANOL" "ELENGUA" "ASISTEN" #[34] "MUN_ASI" "NOM_MUN_ASI" "ENT_PAIS_ASI" #[37] "TIE_TRASLADO_ESCU" "MED_TRASLADO_ESC1" "MED_TRASLADO_ESC2" #[40] "MED_TRASLADO_ESC3" "ESCOLARI" "NIVACAD" #[43] "ALFABET" "ESCOACUM" "MUN_RES10" #[46] "NOM_MUN_RES10" "ENT_PAIS_RES10" "SITUA_CONYUGAL" #[49] "IDENT_PAREJA" "CONACT" "OCUPACION_C" #[52] "SITUACION_TRAB" "AGUINALDO" "VACACIONES" #[55] "SERVICIO_MEDICO" "UTILIDADES" "INCAP_SUELDO" #[58] "SAR_AFORE" "CREDITO_VIVIENDA" "INGTRMEN" #[61] "ACTIVIDADES_C" "MUN_TRAB" "NOM_MUN_TRAB" #[64] "ENT_PAIS_TRAB" "TIE_TRASLADO_TRAB" "MED_TRASLADO_TRAB1" #[67] "MED_TRASLADO_TRAB2" "MED_TRASLADO_TRAB3" "ACTI_SIN_PAGO1" #[70] "ACTI_SIN_PAGO2" "ACTI_SIN_PAGO3" "ACTI_SIN_PAGO4" #[73] "ACTI_SIN_PAGO5" "ACTI_SIN_PAGO6" "ACTI_SIN_PAGO7" #[76] "ACTI_SIN_PAGO8" "HIJOS_NAC_VIVOS" "HIJOS_FALLECIDOS" #[79] "HIJOS_SOBREVIV" "FECHA_NAC_M" "FECHA_NAC_A" #[82] "SOBREVIVENCIA" "EDAD_MORIR_D" "EDAD_MORIR_M" #[85] "EDAD_MORIR_A" "TAMLOC" nomMunAbrev = c( "ABA","ACA","SMA","APA","APG","ATA","CEL","DOB","COM","COR","CRT", "CUE","DRM","DHG","GTO","HUA","IRA","JAR","JER","LEO","MOR","OCA", "PEN","PBN","PUR","ROM","SAL","SLV","SDU","SFE","SFR","SJI","SLP", "SCT","JUR","MAR","SIL","TCU","TAR","TBL","URI","VDS","VIC","VGR", "XCH","YUR") # Cuántos habitantes hay en el estado? # La estimación del total de la población de Guanajuato es total = sum(datpob$FACTOR) # 5853677 # Cuántos hay por municipio? sum(datpob[ datpob$NOM_MUN=="León",]$FACTOR) # 1578626 sum(datpob[ datpob$NOM_MUN=="Celaya",]$FACTOR) # 494304 sum(datpob[ datpob$MUN==14,]$FACTOR) # 152113 # FACTOR = Factor de expansión aa = unique(datpob$MUN) nm = length(aa) nhab = rep(0,nm) for(i in 1:nm){ nhab[i] = sum(datpob[ datpob$MUN==aa[i],]$FACTOR) } uu = data.frame(nomMunAbrev,nhab,aa) zz = unique(datpob$NOM_MUN) # En el presente análisis, no haremos inferencias globales. Sólo # exploraremos la muestra, como si fuera toda la población competa # Esto es, no haremos más uso de la variable FACTOR # Ingreso por trabajo (mensual) aa = datpob$INGTRMEN aa = aa[!is.na(aa)] length(aa) aa = aa[aa<999998] summary(aa) #summary(datpob$FACTOR) hist(aa[aa<6000]) hist(aa[aa<100000]) hist(aa[aa<50000]) hist(aa[aa<20000]) hist(aa) # Qué proporción de personas ganan más de 8000 al mes? naa = length(aa) 100 * length(aa[aa>8000]) / naa a #Cargamos la base de datos sobre viviendas, es un csv datviv = read.csv("TR_VIVIENDA11.csv",header=TRUE) # 178,100 88 # [1] "ID_VIV" "ENT" "NOM_ENT" "MUN" # [5] "NOM_MUN" "LOC50K" "NOM_LOC" "COBERTURA" # [9] "ESTRATO" "UPM" "FACTOR" "CLAVIVP" #[13] "PAREDES" "TECHOS" "PISOS" "CUADORM" #[17] "TOTCUART" "COCINA" "LUGAR_COCINA" "COMBUSTIBLE" #[21] "ESTUFA" "ELECTRICIDAD" "FOCOS" "FOCOS_AHORRA" #[25] "AGUA_ENTUBADA" "ABA_AGUA_ENTU" "ABA_AGUA_NO_ENTU" "TINACO" #[29] "CISTERNA" "BOMBA_AGUA" "REGADERA" "BOILER" #[33] "CALENTADOR_SOLAR" "AIRE_ACON" "PANEL_SOLAR" "SERSAN" #[37] "CONAGUA" "USOEXC" "DRENAJE" "DESTINO_BASURA" #[41] "SEPARA_BASURA" "SEPARA_ALI" "SEPARA_ABO" "SEPARA_VEN" #[45] "REFRIGERADOR" "LAVADORA" "HORNO" "AUTOPROP" #[49] "RADIO" "TELEVISOR" "TELEVISOR_PP" "COMPUTADORA" #[53] "TELEFONO" "CELULAR" "INTERNET" "SERV_TV_PAGA" #[57] "TENENCIA" "NUM_DUE_VIV1" "NUM_DUE_VIV2" "ESCRITURAS_VIV" #[61] "FORMA_ADQUI" "FINANCIAMIENTO" "DEUDA" "NUMPERS" #[65] "INGR_PEROTROPAIS" "INGR_PERDENTPAIS" "INGR_AYUGOB" "INGR_JUBPEN" #[69] "TERRENO_AGROPE" "TERRENO_VIVERO" "NUM_DUE_TERR" "ALIM_ADU1" #[73] "ALIM_ADU2" "ALIM_ADU3" "ING_ALIM_ADU1" "ING_ALIM_ADU2" #[77] "ING_ALIM_ADU3" "ALIM_MEN1" "ALIM_MEN2" "ALIM_MEN3" #[81] "ING_ALIM_MEN1" "ING_ALIM_MEN2" "ING_ALIM_MEN3" "TAMLOC" #[85] "TIPOHOG" "JEFE_SEXO" "JEFE_EDAD" "INGTRHOG" aa = datviv$SERSAN aa = aa[!is.na(aa)] sum(aa==1) / (sum(aa==1)+sum(aa==2)+sum(aa==3)) table(aa) sum(is.na(aa)) aa = unique(datpob$MUN) nm = length(aa) propsan = rep(0,nm) for(i in 1:nm){ aa = datviv$SERSAN[datviv$MUN==i] aa = aa[!is.na(aa)] propsan[i] = sum(aa==1) / (sum(aa==1)+sum(aa==2)+sum(aa==3)) } data.frame(nomMunAbrev,round(100*propsan)) datviv$ID_VIV a aa = datpob$EDAD aa = aa[aa<200] hist(aa) boxplot(aa) summary(aa) datpob[datpob$EDAD==110,] # Cuál es el municipio más joven? edadM = rep(0,nm) for(i in 1:nm){ aa = datpob$EDAD[datpob$MUN==i] aa = aa[!is.na(aa)] edadM[i] = median(aa) } data.frame(nomMunAbrev,edadM) i = 1 sum(is.na(aa)) a # Qué proporción de viviendas cuenta con internet? aa = datviv$INTERNET length(aa) table(aa) sel = which(!(is.na(aa))) aa = aa[sel] naa = length(aa) sum(aa == 5) / naa # Qué proporción de viviendas cuenta con internet? (por municipio) internetmun = rep(0,46) for(i in 1:46){ aa = datviv$INTERNET[datviv$MUN==i] sel = which(!(is.na(aa))) aa = aa[sel] naa = length(aa) internetmun[i] = round(100*sum(aa == 5) / naa,1) } data.frame(nomMunAbrev,internetmun) a # un mapa de esto # Indicador de mujeres (no NA's, ok) selMuj = (datpob$SEXO==3) # Indicador de rango de edades (en toda la población) (no NA's, ok) rango = seq(15,50,by=5) zz = findInterval( datpob$EDAD, rango ) # 0<15, 1=15-19, 2=20-24, 3=25-29, 4=30-34, 5=35-39, 6=40-44, 7=45-49, 8>49 uu = (datpob$HIJOS_NAC_VIVOS < 26) # para evitar códigos 98 y 99 # 98 = Omisión del tema # 99 = No especificado # Nombre completo de municipios nommuni = as.character(unique(datpob$NOM_MUN)) # iris setosa = iris[1:50,] versic = iris[51:100,] virgin = iris[101:150,] mean(setosa[,1]) mean(versic[,1]) mean(virgin[,1]) mean(setosa[,2]) mean(versic[,2]) mean(virgin[,2]) tapply(iris$Sepal.Length,iris$Species,sd) summary(iris) cor(setosa[,1],setosa[,2]) plot(setosa[,1],setosa[,2]) x = seq(-1,1,len=100) y1 = sqrt(1-x^2) y2 = -y1 plot(x,y1,ylim=c(-1,1)) points(x,y2) xc = c(x,x) yc = c(y1,y2) cbind(xc,yc) cor(xc,yc) cor(setosa[,1:4]) cor(versic[,1:4]) cor(virgin[,1:4]) x = rnorm(100,mean=10,sd=1) mean(x) x = runif(100,min=-1,max=1) y = runif(100,min=-1,max=1) plot(x,y) cor(x,y) y = 3*x cor(x,-y) x = c(1,1,-1,-1) y = c(1,-1,1,-1) cor(x,y) library("aplpack") iris_sample = iris[sample(1:dim(iris)[1],size=49,replace=F),] faces(iris_sample[1:4],face.type=1,labels=iris_sample$Species) faces(iris_sample[1:4],face.type=1) iris # Prube de hipótesis 1- pbinom(7,size=10,prob=.5) # error tipo I = 0.0546875 # error tipo II pp = seq(.5,1,len=100) pbinom(7,size=10,prob=pp) # 10.2 Wackerly pbinom(12,size=20,prob=.8) 1-pbinom(12,20,.6) 1-pbinom(12,20,.4) qbinom(.01,20,.8) # 10.3 pbinom(11,size=20,prob=.8) 1-pbinom(11,20,.6) 1-pbinom(11,20,.4) # Comparación de dos algoritmos para elegir # u punto al azar en una circunferencia par(mfrow=c(2,1),mar=c(1,1,1,1)) # Algoritmo 1 N = 200 x = runif(N,min=-1,max=1) y = sample(c(-1,1),size=N,replace = TRUE)*sqrt(1-x^2) plot(x,y,pch=20,col="blue") abline(h=0,v=0) # Algoritmo 2 N = 200 tet = runif(N,min=0,max=2*pi) x = cos(tet) y = sin(tet) plot(x,y,pch=20,col="blue") abline(h=0,v=0) # Elijo al azar una secante # Qué longitud tiene? N = 200000 tet1 = runif(N,min=0,max=2*pi) x1 = cos(tet1) y1 = sin(tet1) tet2 = runif(N,min=0,max=2*pi) x2 = cos(tet2) y2 = sin(tet2) sec = sqrt((y2-y1)^2+(x2-x1)^2) hist(sec) sample(1:6,size=1) # Elegir al azar un punto en la superficie de una esfera N = 200 lat = runif(N,min=-90,max=90) lon = runif(N,min=0,max=2*pi) plot(lon,lat) pbinom(12,20,.8) 1 - pbinom(12,20,.6) 1 - pbinom(12,20,.4) # Función de potencia # Quiero una grafica de beta vs p ps = seq(0,.8,len=100) betas = 1 - pbinom(12,20,ps) plot(ps,betas,type="l",col="blue",lwd=2) 1 - ( pbinom(21,36,.5)-pbinom(14,36,.5) )