Curso R-Biodiversidad, UASD-US, de enero de 2016

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)

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *