Script de R del curso sobre Biodiversidad con R impartido el 20 de enero de 2016, seguir este vínculo.
#CARGA DE PAQUETES library(vegan) library(BiodiversityR) library(reshape) library(dplyr) library(tidyr) library(psych) library(ggplot2) library(fossil) library(plotKML) library(sp) library(plotGoogleMaps) library(maptools) #OPERACIONES BÁSICAS. LA GESTIÓN DE OBJETOS SE INTRODUCE A LO LARGO DEL SCRIPT. SE RECOMIENDA SEGUIR UN MANUAL INTRODUCTORIO DE R: https://cran.r-project.org/doc/contrib/R-intro-1.1.0-espanol.1.pdf 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 #LECTURA DE LOS DATOS ##MATRIZ DE ABUNDANCIA d.ma <- read.csv("http://geografiafisica.org/r/br/br.ma.csv",fileEncoding = 'latin1',comment.char = '#') #CREA LA MATRIZ DE ABUNDANCIA (COMO UN DATA FRAME DE R) A PARTIR DE UN CSV ALOJADO EN geografiafisica.org, DONDE CADA FILA ES UN INDIVIDUO. A PARTIR DE ESTA MATRIZ SE GENERARÁ LA DE COMUNIDAD str(d.ma) #MUESTRA LA ESTRUCTURA DE LA MATRIZ DE ABUNDANCIA head(d.ma) #MUESTRA LAS 6 PRIMERAS FILAS DE LA MATRIZ DE ABUNDANCIA. PARA MOSTRAR EL OBJETO COMPLETO, ESCRIBIR "d" (SIN LAS COMILLAS), PERO ESTO PROBABLEMENTE LLENE LA CONSOLA ##MATRIZ AMBIENTAL d.env <- read.csv("http://geografiafisica.org/r/br/br.env.csv",fileEncoding = 'latin1',comment.char = '#') #CREA MATRIZ "AMBIENTAL" (COMO UN DATA FRAME DE R) A PARTIR DE UN CSV ALOJADO EN geografiafisica.org, QUE CONTIENE ATRIBUTOS POR TRANSECTO. CADA FILA ES UN TRANSECTO head(d.env) #MUESTRA LAS 6 PRIMERAS FILAS DE LA MATRIZ AMBIENTAL. PARA MOSTRAR EL OBJETO COMPLETO, ESCRIBIR "d.env" (SIN LAS COMILLAS), PERO ESTO PROBABLEMENTE LLENE LA CONSOLA str(d.env) #MUESTRA LA ESTRUCTURA DE LA MATRIZ AMBIENTAL #VERIFICACIÓN DE INTEGRIDAD DE LOS DATOS table(d.ma$parcela) #FRECUENCIAS DE PARCELAS DESDE LA MATRIZ DE ABUNDANCIA. EL RESULTADO CORRESPONDE A LA ABUNDANCIA POR TRANSECTO, DADO QUE CADA FILA ES UN INDIVIDUO table(d.env$parcela) #CONTEO DE FRECUENCIAS DE PARCELAS DESDE LA MATRIZ AMBIENTAL. EL RESULTADO DEBERÍA SER 1 EN CADA CASO, DADO QUE CADA FILA ES UN TRANSECTO length(table(d.ma$parcela)) #NÚMERO DE PARCELAS ÚNICAS EN LA MATRIZ DE ABUNDANCIA. EL TOTAL DEBERÍA SER EQUIVALENTE AL NÚMERO DE TRANSECTOS length(table(d.env$parcela)) #NÚMERO DE PARCELAS ÚNICAS EN LA MATRIZ AMBIENTAL. EL TOTAL DEBERÍA SER EQUIVALENTE AL NÚMERO DE TRANSECTOS merge(as.data.frame(table(d.env$parcela)),as.data.frame(table(d.ma$parcela)),by='Var1',all=T) #ESTA UNIÓN DEBE MOSTRAR UN NÚMERO SIMILAR DE ENTRADAS ENTRE FUENTES. NO DEBEN APARECER "NA" EN NINGÚN CASO. LA COLUMNA "Freq.x" INDICA LA CANTIDAD DE VECES QUE APARECE EL NOMBRE DE UN MUESTREO EN LA MATRIZ AMBIENTAL. LA COLUMNA "Freq.y" INDICA LA CANTIDAD DE VECES QUE APARECE EL NOMBRE DE UN MUESTREO EN LA MATRIZ DE ABUNDANCIA #CONVERSIÓN DE LA MATRIZ DE ABUNDANCIA EN UNA MATRIZ DE COMUNIDAD d <- cast(d.ma,parcela~especie,fun.aggregate=length,value='id') #GENERACIÓN DE LA MATRIZ DE COMUNIDAD MEDIANTE LA TRASPOSICIÓN DE LA MATRIZ DE ABUNDANCIA AGREGANDO EL NÚMERO DE INDIVIDUOS DE CADA ESPECIE POR TRANSECTO. EL CAMPO DE ORDENACIÓN DE LA TABLA RESULTANTE ES "parcela" str(d) #MUESTRA LA ESTRUCTURA DE LA MATRIZ DE COMUNIDAD rownames(d) <- d$parcela #FIJACIÓN DE LOS NOMBRES DE TRANSECTOS COMO NOMBRES DE FILAS d <- d[,2:ncol(d)] #RECUPERACIÓN DE LAS COLUMNAS QUE CONTIENEN LOS DATOS DE ABUNDANCIA POR ESPECIES, DE MANERA QUE LA MATRIZ DE COMUNIDAD SOLO CONTENGA ESPECIES d <- as.data.frame(d) #CONVERSIÓN DEL OBJETO A UN DATA FRAME head(d[,1:4]) #MUESTRA LAS PRIMERAS 6 FILAS Y LAS PRIMERAS 4 COLUMNAS str(d) #MUESTRA LA ESTRUCTURA DE LA MATRIZ DE COMUNIDAD PREPARADA PARA ANÁLISIS DE DIVERSIDAD CON vegan #IMPORTANTE: COMPROBACIÓN DE QUE LAS MATRICES ESTÁN EN EL MISMO ORDEN pruebadeordenacion <- data.frame(enmatrizdeambiente=d.env$parcela,enmatrizdedecomunidad=rownames(d)) pruebadeordenacion$logica=pruebadeordenacion$enmatrizdeambiente==pruebadeordenacion$enmatrizdedecomunidad pruebadeordenacion #TODOS LOS RESULTADOS DE LA COLUMNA 'pruebadeordenacion' DEBEN APARECER COMO "TRUE" #NÚMERO DE TRANSECTOS POR... ##FORMACIÓN VEGETAL sort(table(d.env$formacion_vegetal)) ##RÉGIMEN BIOCLIMÁTICO sort(table(d.env$regimen_bioclimatico)) ##TOMADORES sort(table(d.env$tomadores)) ##MUNICIPIO sort(table(d.env$municipio)) #GENERACIÓN DE MATRICES DE COMUNIDAD POR... ##FORMACIÓN VEGETAL dfv <- t(sapply(lapply(unique(d.env$formacion_vegetal),function(x) d.env[d.env$formacion_vegetal==x,'parcela']),function(x) colSums(d[x,])));rownames(dfv)<-unique(d.env$formacion_vegetal) dfv <- as.data.frame(dfv) str(dfv) head(dfv[,1:4]) dfv.env <- data.frame(formacion_vegetal=rownames(dfv)) ##FORMACIÓN VEGETAL, PERO BASADO EN CÓDIGOS DEL CAMPO 'tipo' dfvc <- t(sapply(lapply(unique(d.env$tipo),function(x) d.env[d.env$tipo==x,'parcela']),function(x) colSums(d[x,])));rownames(dfvc)<-unique(d.env$tipo) dfvc <- as.data.frame(dfvc) str(dfvc) head(dfvc[,1:4]) ##RÉGIMEN BIOCLIMÁTICO drb <- t(sapply(lapply(unique(d.env$regimen_bioclimatico),function(x) d.env[d.env$regimen_bioclimatico==x,'parcela']),function(x) colSums(d[x,])));rownames(drb)<-unique(d.env$regimen_bioclimatico) drb <- as.data.frame(drb) str(drb) head(drb[,1:4]) drb.env <- data.frame(regimen_bioclimatico=rownames(drb)) #ANÁLISIS DE BIODIVERSIDAD: ÍNDICES, GRÁFICOS PROPIOS DEL ANÁLISIS DE BIODIVERSIDAD ##PARA EL CONJUNTO, RIQUEZA Y ABUNDANCIA specnumber(colSums(d)) #RIQUEZA TOTAL sum(rowSums(d)) #ABUNDANCIA TOTAL ##POR TRANSECTO, RIQUEZA Y ABUNDANCIA rowSums(d) #ABUNDANCIA POR TRANSECTO d.env$abundancia <- rowSums(d) #GUARDANDO EL VALOR DE ABUNDANCIA POR TRANSECTO EN LA MATRIZ AMBIENTAL. IMPORTANTE: SE REQUIERE LA COMPROBACIÓN DE QUE LAS MATRICES ESTÁN EN EL MISMO ORDEN. SE PUEDE COMPROBAR LA INTEGRIDAD DE ESTA UNIÓN MEDIANTE: rownames(d)==attr(rowSums(d),'names'); ESTE CÓDIGO VERIFICA SI LOS NOMBRES DE LAS FILAS ESTÁN EN EL MISMO ORDEN QUE LOS NOMBRES DE LOS ATRIBUTOS DE LA LISTA DE ABUNDANCIA DE rowSums(d). TODOS LOS RESULTADOS DEBEN APARECER COMO "TRUE" specnumber(d) #RIQUEZA POR TRANSECTO d.env$riqueza <- specnumber(d) #GUARDANDO EL VALOR DE ABUNDANCIA POR TRANSECTO EN LA MATRIZ AMBIENTAL. IMPORTANTE: SE REQUIERE LA COMPROBACIÓN DE QUE LAS MATRICES ESTÁN EN EL MISMO ORDEN. SE PUEDE COMPROBAR LA INTEGRIDAD DE ESTA UNIÓN MEDIANTE: rownames(d)==attr(specnumber(d),'names'); ESTE CÓDIGO VERIFICA SI LOS NOMBRES DE LAS FILAS ESTÁN EN EL MISMO ORDEN QUE LOS NOMBRES DE LOS ATRIBUTOS DE LA LISTA DE ABUNDANCIA DE specnumber(d). TODOS LOS RESULTADOS DEBEN APARECER COMO "TRUE" ##GRÁFICOS dev.new() #GENERA UNA VENTANA GRÁFICA INDEPENDIENTE. ESTA SENTENCIA SE REPITE MUCHO A LO LARGO DEL SCRIPT, Y EN LO ADELANTE APARECERÁ EN LA MISMA LÍNEA QUE EL CÓDIGO DE CREACIÓN DEL GRÁFICO CORRESPONDIENTE, SEPARADO POR UN ";" barplot(sort(rowSums(d),decreasing=T),las=3,ylim=c(0,100),ylab='abundancia (número de individuos)',main='Abundancia por transecto') #GRÁFICO DE BARRAS DE ABUNDANCIA POR TRANSECTO ##ÍNDICES DE BIODIVERSIDAD USANDO vegan PASO A PASO diversity(d) #CALCULA EL ÍNDICE DE SHANNON diversity(d,'simpson') #CALCULA LA EQUIDAD DE SIMPSON (1-D), DENOMINADO TAMBIÉN "GINI-SIMPSON" diversity(dfv) #CALCULA EL ÍNDICE DE SHANNON PARA LA MATRIZ DE COMUNIDAD POR FORMACIÓN VEGETAL diversity(dfv,'simpson') #CALCULA LA EQUIDAD DE SIMPSON (1-D), DENOMINADO TAMBIÉN "GINI-SIMPSON" PARA LA MATRIZ DE COMUNIDAD POR FORMACIÓN VEGETAL diversity(drb) #CALCULA EL ÍNDICE DE SHANNON PARA LA MATRIZ DE COMUNIDAD POR RÉGIMEN BIOCLIMÁTICO diversity(drb,'simpson') #CALCULA LA EQUIDAD DE SIMPSON (1-D), DENOMINADO TAMBIÉN "GINI-SIMPSON" PARA LA MATRIZ DE COMUNIDAD POR RÉGIMEN BIOCLIMÁTICO ##ÍNDICES CALCULADOS MASIVAMENTE Y ANEXADOS A MATRIZ AMBIENTAL indicesmasivo <- cbind(parcela=rownames(d),as.data.frame(sapply(c("Shannon","Simpson","inverseSimpson","Logalpha","Berger","Jevenness","Eevenness" ),function(x) diversityresult(d,index=x,method='s')))) #GENERA UNA TABLA CON VARIOS ÍNDICES DE FORMA MASIVA colnames(indicesmasivo)[2:ncol(indicesmasivo)]<-c('Shannon','Gini_Simpson','inverso_de_Simpson','Fisheralpha_Logalpha','Berger_Parker','equidad_de_Pielou_J_evenness','equidad_de_Buzas_Gibson_E_evenness') #COLOCANDO NOMBRES ADECUADOS A COLUMNAS indicesmasivo d.env <- merge(d.env,indicesmasivo) #ANEXA LOS ÍNDICES GENERADOS A LA MATRIZ AMBIENTAL ##LOGARITMOS, TAMBIÉN ANEXADOS logindicesmasivo <- cbind(parcela=rownames(d),as.data.frame(sapply(c("Shannon","Simpson","inverseSimpson","Logalpha","Berger","Jevenness","Eevenness" ),function(x) log(diversityresult(d,index=x,method='s'))))) #GENERA UNA TABLA CON LOS LOGARITMOS DE VARIOS ÍNDICES DE FORMA MASIVA colnames(logindicesmasivo)[2:ncol(logindicesmasivo)]<-paste('log',c('Shannon','Gini_Simpson','inverso_de_Simpson','Fisheralpha_Logalpha','Berger_Parker','equidad_de_Pielou_J_evenness','equidad_de_Buzas_Gibson_E_evenness'),sep='_') #COLOCANDO NOMBRES ADECUADOS A COLUMNAS logindicesmasivo d.env <- merge(d.env,logindicesmasivo) #ANEXA LOS ÍNDICES GENERADOS A LA MATRIZ AMBIENTAL d.env ##ESTIMADOR DE RIQUEZA CHAO1 estimadorchao1 <- data.frame(parcela=rownames(d),chao1=sapply(rownames(d),function(x) chao1(d[x,]))) #GENERA UNA TABLA CON EL ESTIMADOR DE CHAO1 d.env <- merge(d.env,estimadorchao1) #ANEXA LOS VALORES GENERADOS A LA MATRIZ AMBIENTAL ##DIVERSIDAD DE RENYI renyi(d) #VALORES DE DIVERSIDAD DE LA MATRIZ DE COMUNIDAD POR TRANSECTO PARA CADA NÚMERO DE HILL dev.new();plot(renyi(d)) #¡96 GRÁFICOS! dev.new();plot(renyi(d[1:2,])) #SOLO 2 dev.new();plot(renyi(d[1:24,])) #SELECCIÓN DE LOS PRIMEROS 24 MUESTREOS dev.new();plot(renyi(dfv)) #VALORES DE DIVERSIDAD DE LA MATRIZ DE COMUNIDAD POR FORMACIÓN VEGETAL PARA CADA NÚMERO DE HILL dev.new();plot(renyi(drb)) #VALORES DE DIVERSIDAD DE LA MATRIZ DE COMUNIDAD POR RÉGIMEN BIOCLIMÁTICO PARA CADA NÚMERO DE HILL ##CURVA DE ACUMULACIÓN DE ESPECIES. IMPORTANTE: SE REQUIERE HABER GENERADO LA COLUMNA DE 'abundancia' PREVIAMENTE EN LA MATRIZ AMBIENTAL d.env ###PARA EL CONJUNTO dev.new();accumplot(accumresult(d,y=d.env,scale='abundancia',method='rarefaction'),xlab='individuos',ylab='especies',ci.type='polygon',ci.col='yellow') title(main='Curva de acumulación de especies por individuos usando método de rarefacción') ###PARA MATA ATLÁNTICA dev.new();accumplot(accumresult(d,y=d.env,factor='formacion_vegetal',level='Mata atlantica',scale='abundancia',method='rarefaction'),xlab='individuos',ylab='especies',ci.type='polygon',ci.col='yellow') title(main='Curva de acumulación de especies por individuos usando método de rarefacción') ##RAREFACCIÓN Y CURVAS DE RAREFACCIÓN dev.new();rarecurve(d) #¡96 CURVAS, ILEGIBLE! dev.new();rarecurve(dfv) #POR FORMACIÓN VEGETAL dev.new();rarecurve(dfvc) #POR FORMACIÓN VEGETAL BASADO EN EL CAMPO 'tipo', PARA FACILITAR LA LEGIBILIDAD DEL GRÁFICO dev.new();rarecurve(drb) #POR RÉGIMEN BIOCLIMÁTICO ##TABLA Y CURVAS DE RANGO-ABUNDANCIA ###TABLAS, PARA EL CONJUNTO, POOLED rankabundance(d) #TABLA DE RANGO-ABUNDANCIA DE TODOS LOS MUESTREOS COMO SI FUESEN UNO SOLO (COMBINADOS, POOLED) rankabundance(d)[,2] #TABLA DE RANGO-ABUNDANCIA DE TODOS LOS MUESTREOS COMO SI FUESEN UNO SOLO (COMBINADOS, POOLED), MOSTRANDO SOLO LAS ESPECIES Y SU FRECUENCIA ###TABLAS, POR TRANSECTO tra <- sapply(d.env$parcela,function(x) rankabundance(d,y=d.env,factor='parcela',level=x)); names(tra) <- d.env$parcela #GENERACION DE LAS 96 TABLAS DE RANGO-ABUNDANCIA tra ###TABLAS, POR FORMACIÓN VEGETAL tra.fv <- sapply(dfv.env$formacion_vegetal,function(x) rankabundance(d,y=dfv.env,factor='formacion_vegetal',level=x)); names(tra.fv)<-dfv.env$formacion_vegetal #GENERACION DE LAS 8 TABLAS DE RANGO-ABUNDANCIA tra.fv ###TABLAS, POR RÉGIMEN BIOCLIMÁTICO tra.rb <- sapply(drb.env$regimen_bioclimatico,function(x) rankabundance(d,y=drb.env,factor='regimen_bioclimatico',level=x)); names(tra.rb)<-drb.env$regimen_bioclimatico #GENERACION DE LAS 8 TABLAS DE RANGO-ABUNDANCIA tra.rb ###GRÁFICOS DE RANGO-ABUNDANCIA dev.new();rankabuncomp(d,y=d.env,factor='parcela',scale='logabun') #ILEGIBLE. A CONTINUACIÓN, GRÁFICOS POR GRUPOS DE 24 TRANSECTOS dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[1:24,]),function(x) rankabunplot(rankabundance(d,y=d.env,factor='parcela',level=x),scale='logabun',main=x,xlim=c(0,20),ylim=c(0,20),specnames=1));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[25:48,]),function(x) rankabunplot(rankabundance(d,y=d.env,factor='parcela',level=x),scale='logabun',main=x,xlim=c(0,20),ylim=c(0,20),specnames=1));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[49:72,]),function(x) rankabunplot(rankabundance(d,y=d.env,factor='parcela',level=x),scale='logabun',main=x,xlim=c(0,20),ylim=c(0,20),specnames=1));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[73:96,]),function(x) rankabunplot(rankabundance(d,y=d.env,factor='parcela',level=x),scale='logabun',main=x,xlim=c(0,20),ylim=c(0,20),specnames=1));par(mfrow=c(1,1)) ##CURVA DE DOMINANCIA K ###POR TRANSECTO dev.new();rankabuncomp(d,y=d.env,factor='parcela',scale='accumfreq', log='x',main='Dominancia k de cada transecto') #¡96 LÍNEAS! ###POR FORMACIÓN VEGETAL dev.new();rankabuncomp(dfv,y=dfv.env,factor='formacion_vegetal',scale='accumfreq', log='x',main='Dominancia k de cada formación vegetal') ###POR RÉGIMEN BIOCLIMÁTICO dev.new();rankabuncomp(drb,y=drb.env,factor='regimen_bioclimatico',scale='accumfreq', log='x',main='Dominancia k de cada régimen bioclimático') ##AJUSTE POR AIC/BIC A VARIOS MODELOS DE ABUNDANCIA ###POR TRANSECTO radfit(d) dev.new() plot(radfit(d)) ###AGRUPADOS POR FORMACIÓN VEGETAL dev.new() plot(radfit(dfv)) ###AGRUPADOS POR RÉGIMEN BIOCLIMÁTICO dev.new() plot(radfit(drb)) ###AJUSTE A DISTRIBUCIÓN DE PRESTON ####AGRUPADOS POR TRANSECTO prestondistr(colSums(d)) veiledspec(colSums(d)) dev.new();plot(prestondistr(colSums(d))) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[1:24,]),function(x) plot(prestondistr(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[25:48,]),function(x) plot(prestondistr(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[49:72,]),function(x) plot(prestondistr(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[73:96,]),function(x) plot(prestondistr(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) ####AGRUPADOS POR FORMACIÓN VEGETAL dev.new();par(mfrow=c(2,4));par(mar=c(4,4,2,2));lapply(rownames(dfv),function(x) plot(prestondistr(dfv[rownames(dfv)==x,]),main=x));par(mfrow=c(1,1)) ####AGRUPADOS POR RÉGIMEN BIOCLIMÁTICO dev.new();par(mfrow=c(1,2));par(mar=c(4,4,2,2));lapply(rownames(drb),function(x) plot(prestondistr(drb[rownames(drb)==x,]),main=x));par(mfrow=c(1,1)) ###AJUSTE A DISTRIBUCIÓN FISHER fisherfit(colSums(d)) dev.new();plot(fisherfit(colSums(d))) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[1:24,]),function(x) plot(fisherfit(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[25:48,]),function(x) plot(fisherfit(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[49:72,]),function(x) plot(fisherfit(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) dev.new();par(mfrow=c(4,6));par(mar=c(4,4,2,2));lapply(rownames(d[73:96,]),function(x) plot(fisherfit(d[rownames(d)==x,]),main=x));par(mfrow=c(1,1)) ####AGRUPADOS POR FORMACIÓN VEGETAL dev.new();par(mfrow=c(2,4));par(mar=c(4,4,2,2));lapply(rownames(dfv),function(x) plot(fisherfit(dfv[rownames(dfv)==x,]),main=x));par(mfrow=c(1,1)) ####AGRUPADOS POR RÉGIMEN BIOCLIMÁTICO dev.new();par(mfrow=c(1,2));par(mar=c(4,4,2,2));lapply(rownames(drb),function(x) plot(fisherfit(drb[rownames(drb)==x,]),main=x));par(mfrow=c(1,1)) #ESTADÍSTICA DESCRIPTIVA DE ÍNDICES CON APOYO GRÁFICO ##GRÁFICOS DESCRIPTIVOS PARA LOS ÍNDICES dev.new() #GENERA UNA VENTANA GRÁFICA INDEPENDIENTE hist(rowSums(d)/100,col=100,xlab='densidad de individuos (individuos por metro cuadrado)',ylab='número de transectos',main='Número de transectos según densidad de individuos') #HISTOGRAMA DE LA DENSIDAD DE INDIVIDUOS dev.new() #GENERA UNA VENTANA GRÁFICA INDEPENDIENTE hist(log(rowSums(d))/100,col=100,xlab='logaritmo de la densidad de individuos (individuos por metro cuadrado)',ylab='número de transectos',main='Número de transectos según el logaritmo de la densidad de individuos') #HISTOGRAMA DEL LOGARITMO DE LA DENSIDAD DE INDIVIDUOS mean(rowSums(d)) #MEDIA DE LA ABUNDANCIA median(rowSums(d)) #MEDIANA DE LA ABUNDANCIA sd(rowSums(d)) #DESVIACIÓN ESTANDAR DE LA ABUNDANCIA var(rowSums(d)) #VARIANZA DE LA ABUNDANCIA ##ESTADÍSTICOS DESCRIPTIVOS PARA LOS ÍNDICES ###UNO A UNO, POCO PRÁCTICO mean(d.env$Shannon) median(d.env$Shannon) sd(d.env$Shannon) var(d.env$Shannon) min(d.env$Shannon) max(d.env$Shannon) ###MEDIANTE sapply sapply(d.env[sapply(d.env,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(d.env[sapply(d.env,is.numeric)],function(x) {c(media=mean(x,na.rm=T),mediana=median(x,na.rm=T))})) #RESULTADO TRASPUESTO ###MEDIANTE dplyr Y psych. USO DEL OPERADOR PIPE (%>%) CLAVE EN EL PROCESO. EL PAQUETE dplyr ES BASTANTE POTENTE, PERO REQUIERE EL DOMINICO DE SU SINTAXIS d.env %>% select(parcela,Shannon:log_equidad_de_Buzas_Gibson_E_evenness) %>% gather(indice,valor,-parcela) %>% group_by(indice) %>% do(describe(.$valor)) %>% select(-vars) options(dplyr.width = Inf) #PARA MOSTRAR TODAS LAS COLUMNAS d.env %>% select(parcela,Shannon:log_equidad_de_Buzas_Gibson_E_evenness) %>% gather(indice,valor,-parcela) %>% group_by(indice) %>% do(describe(.$valor)) %>% select(-vars) #AHORA SE VEN TODAS ##ALGUNOS GRÁFICOS ###UNO A UNO dev.new();par(mar=c(12,5,1,1));boxplot(Shannon~formacion_vegetal,d.env,las=2, main='Índice de Shannon según formación vegetal', ylab='nats') #LÍMITES DE Y AUTOMÁTICOS dev.new();par(mar=c(12,5,1,1));boxplot(Shannon~formacion_vegetal,d.env,las=2, main='Índice de Shannon según formación vegetal', ylab='nats',ylim=c(0,3)) #LÍMITES DE Y FIJOS dev.new();par(mar=c(12,5,1,1));boxplot(Gini_Simpson~formacion_vegetal,d.env,las=2, main='Índice de Gini-Simpson según formación vegetal', ylab='valor del índice') #LÍMITES DE Y AUTOMÁTICOS dev.new();par(mar=c(12,5,1,1));boxplot(Gini_Simpson~formacion_vegetal,d.env,las=2, main='Índice de Gini-Simpson según formación vegetal', ylab='valor del índice',ylim=c(0,1)) #LÍMITES DE Y FIJOS ###MASIVAMENTE pbppfv <- d.env %>% select(`formación vegetal`=formacion_vegetal,Shannon:log_equidad_de_Buzas_Gibson_E_evenness) %>% gather(indice,valor,-`formación vegetal`) dev.new();ggplot(pbppfv, aes(x=`formación vegetal`,y=valor,fill=`formación vegetal`)) + geom_boxplot() + theme(text = element_text(size = 18)) + facet_wrap(~indice,scales='free_y') + ggtitle('Biodiversidad medida según distintos índices por formación vegetal') + ylab('valores de los índices') + theme(axis.ticks = element_blank(), axis.text.x = element_blank()) #ESTADÍSTICA INFERENCIAL POR CRITERIOS DE AGRUPAMIENTO DE LA MATRIZ AMBIENTAL shapiro.test(log(rowSums(d))/100) #PRUEBA DE NORMALIDAD DE LA DENSIDAD DE INDIVIDUOS options(scipen = 9999) data.frame(valor_de_p_ST=t(d.env %>% select(Shannon:log_equidad_de_Pielou_J_evenness) %>% summarise_each(funs(shapiro.test(.)$p.value)))) d.env %>% select(formacion_vegetal,Shannon:log_equidad_de_Pielou_J_evenness) %>% gather(indice,valor,-formacion_vegetal) %>% group_by(indice,formacion_vegetal) %>% summarise(n=n()) %>% print(n=300) d.env %>% select(formacion_vegetal,Shannon:log_equidad_de_Pielou_J_evenness) %>% gather(indice,valor,-formacion_vegetal) %>% group_by(indice,formacion_vegetal) %>% summarise(p=shapiro.test(valor)$p.value) %>% arrange(desc(p)) %>% print(n=300) paov <- as.data.frame(d.env %>% select(formacion_vegetal,Shannon:log_equidad_de_Pielou_J_evenness) %>% gather(indice,valor,-formacion_vegetal)) levels(paov$indice) oneway.test(valor~formacion_vegetal,paov[paov$indice=='Shannon',]) oneway.test(valor~formacion_vegetal,paov[paov$indice=='Gini_Simpson',]) oneway.test(valor~formacion_vegetal,paov[paov$indice=='Fisheralpha_Logalpha',]) oneway.test(valor~formacion_vegetal,paov[paov$indice=='Berger_Parker',]) kruskal.test(valor~formacion_vegetal,paov[paov$indice=='Berger_Parker',]) ##CORRELACIONES cor(d.env[sapply(d.env,is.numeric)]) cor.test(d.env$x1,d.env$abundancia) #DIVERSIDAD BETA. BASADO EN TRABAJO DE PATRICIA KOLEFF. VER: http://geografiafisica.org/consigna/visitante/de_gd/!midiendo_la_diversidad_beta_para_datos_de_presencia_ausencia_vegan_usa_denominaciones_de_aqui_para_betadiver_Koleff_MUY_BUENO.pdf ##GENERACIÓN DE NOMBRES EN UN data.frame nbetad<- as.data.frame( cbind( 'abreviatura'= names( list(w = "(b+c)/(2*a+b+c)", '-1' = "(b+c)/(2*a+b+c)", c = "(b+c)/2", wb = "b+c", r = "2*b*c/((a+b+c)^2-2*b*c)", I = "log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) + (a+c)*log(a+c)) / (2*a+b+c)", e = "exp(log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) + (a+c)*log(a+c)) / (2*a+b+c))-1", t = "(b+c)/(2*a+b+c)", me = "(b+c)/(2*a+b+c)", j = "a/(a+b+c)", sor = "2*a/(2*a+b+c)", m = "(2*a+b+c)*(b+c)/(a+b+c)", '-2' = "pmin(b,c)/(pmax(b,c)+a)", co = "(a*c+a*b+2*b*c)/(2*(a+b)*(a+c))", cc = "(b+c)/(a+b+c)", g = "(b+c)/(a+b+c)", '-3' = "pmin(b,c)/(a+b+c)", l = "(b+c)/2", '19' = "2*(b*c+1)/((a+b+c)^2+(a+b+c))", hk = "(b+c)/(2*a+b+c)", rlb = "a/(a+c)", sim = "pmin(b,c)/(pmin(b,c)+a)", gl = "2*abs(b-c)/(2*a+b+c)", z = "(log(2)-log(2*a+b+c)+log(a+b+c))/log(2)") ) , 'fórmula'= list(w = "(b+c)/(2*a+b+c)", '-1' = "(b+c)/(2*a+b+c)", c = "(b+c)/2", wb = "b+c", r = "2*b*c/((a+b+c)^2-2*b*c)", I = "log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) + (a+c)*log(a+c)) / (2*a+b+c)", e = "exp(log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) + (a+c)*log(a+c)) / (2*a+b+c))-1", t = "(b+c)/(2*a+b+c)", me = "(b+c)/(2*a+b+c)", j = "a/(a+b+c)", sor = "2*a/(2*a+b+c)", m = "(2*a+b+c)*(b+c)/(a+b+c)", '-2' = "pmin(b,c)/(pmax(b,c)+a)", co = "(a*c+a*b+2*b*c)/(2*(a+b)*(a+c))", cc = "(b+c)/(a+b+c)", g = "(b+c)/(a+b+c)", '-3' = "pmin(b,c)/(a+b+c)", l = "(b+c)/2", '19' = "2*(b*c+1)/((a+b+c)^2+(a+b+c))", hk = "(b+c)/(2*a+b+c)", rlb = "a/(a+c)", sim = "pmin(b,c)/(pmin(b,c)+a)", gl = "2*abs(b-c)/(2*a+b+c)", z = "(log(2)-log(2*a+b+c)+log(a+b+c))/log(2)") ) ) nbetad$completo<-paste(nbetad$abreviatura,nbetad$fórmula) ##CÁLCULO DE LOS 24 ÍNDICES DE DIVERSIDAD BETA PARA LOS MUESTREOS DEL OBJETO d, Y ASIGNACIÓN DE NOMBRES A CADA ÍNDICE. betadiver(d,'j') #MATRIZ INMANEJABLE betadiver(d[16:25,])$a + betadiver(d[1:10,])$b + betadiver(d[1:10,])$c #ATRIBUTOS COMPARABLES betadiver(d[16:25,],'j') #COMPARAR CON LA TABLA DE VALORES SIGNIFICATIVOS DEL ÍNDICE DE JACCARD, SEGÚN R. REAL, IDENTIFICANDO LOS UMBRALES EN BASE A LOS ATRIBUTOS COMPARABLES. VER http://geografiafisica.org/consigna/visitante/de_gd/!tablas_de_valores_de_significancia_del_indice_de_similaridad_de_jaccard_beta.pdf betad<-lapply(1:24, function(x) betadiver(d,x));names(betad)<-nbetad$completo #REPRESENTACIÓN EN GOOGLEMAPS Y KML d.env.z24.sp<-d.env[d.env$x1>400000,] #CREA UN OBJETO ESPACIAL CON LOS TRANSECTOS DE LA ZONA UTM 24 d.env.z25.sp<-d.env[d.env$x1<400000,] #CREA UN OBJETO ESPACIAL CON LOS TRANSECTOS DE LA ZONA UTM 25 str(d.env.z24.sp) str(d.env.z25.sp) coordinates(d.env.z24.sp) <- ~x1 + y1 #INDICA AL PROGRAMA QUE LAS COORDENADAS DEL OBJETO d.env.sp ESTÁN EN LOS CAMPOS x1 E y1. ESTO CONVERTIRÁ A d.env.sp EN UN OBJETO ESPACIAL coordinates(d.env.z25.sp) <- ~x1 + y1 #INDICA AL PROGRAMA QUE LAS COORDENADAS DEL OBJETO d.env.sp ESTÁN EN LOS CAMPOS x1 E y1. ESTO CONVERTIRÁ A d.env.sp EN UN OBJETO ESPACIAL proj4string(d.env.z24.sp) <- CRS("+init=epsg:32724") #LE ASIGNA EL SISTEMA DE REFERENCIA DE COORDENADAS CORRESPONDIENTE AL OBJETO proj4string(d.env.z25.sp) <- CRS("+init=epsg:32725") #LE ASIGNA EL SISTEMA DE REFERENCIA DE COORDENADAS CORRESPONDIENTE AL OBJETO d.env.sp.z24.g <- spTransform(d.env.z24.sp, CRS("+proj=longlat +ellps=sphere +no_defs")) #TRANSFORMA A LAT/LON d.env.sp.z25.g <- spTransform(d.env.z25.sp, CRS("+proj=longlat +ellps=sphere +no_defs")) #TRANSFORMA A LAT/LON d.env.sp.g <- do.call("rbind", c(d.env.sp.z24.g,d.env.sp.z25.g)) #UNE AMBOS OBJETOS ESPACIALES plotGoogleMaps(d.env.sp.g,iconMarker='',zcol="parcela","br.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") kmlPoints(d.env.sp.g,kmlfile="br.kml",name=d.env.sp.g$parcela) #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 "parcela"
Dr. José Ramón Martínez Batlle (Ph.D)