Esta guía se hizo en el marco de un curso/taller impartido en Málaga en julio de 2015, y se fue actualizando con los avances de los paquetes. En letra monospace están los comandos que hay que escribir (o copiar y pegar) en la consola de R (y después presionar enter para ejecutarlos). Nota que todos los comandos son sensibles a mayúsculas y minúsculas; que hay que utilizar siempre apóstrofes y comillas rectos (ej. “texto” y no “texto”); y que R sólo está listo para recibir un nuevo comando cuando hay un signo de prompt (>) al final de la consola (atención a cierres de paréntesis y demás).
Descargar R (http://cran.r-project.org) e instalarlo.
Descargar RStudio (interfaz amigable para R; http://www.rstudio.com/products/rstudio/download), instalarlo y abrirlo.
Instalar los paquetes fuzzySim (http://fuzzysim.r-forge.r-project.org) y modEvA (http://modeva.r-forge.r-project.org/), que de momento no están disponibles en el CRAN (el repositorio oficial de R), sino en R-Forge (una plataforma más experimental). Si tienes instalada la última versión de R, deberían funcionar los siguientes comandos:
install.packages("fuzzySim", repos = "http://R-Forge.R-project.org")
install.packages("modEvA", repos = "http://R-Forge.R-project.org")Si sale un mensaje de error, o para cerciorarte de que instalas la última versión de los paquetes, puedes descargarlos (fuzzySim de https://r-forge.r-project.org/R/?group_id=1853 y modEvA de https://r-forge.r-project.org/R/?group_id=1876; versión .tar.gz para Linux y Mac, versión .zip para Windows), guardarlos en el disco e instalarlos desde allí. En RStudio, esto se hace con el menú Tools - Install packages - Install from: Package Archive File, buscando después en el disco dónde se ha guardado cada paquete.
Para importar y exportar tablas en formato .dbf, también vas a necesitar el paquete foreign. Éste está disponible en el CRAN, por tanto sólo tienes que darle al menú Tools - Install packages, Install from: Repository (CRAN), escribir foreign, mantener Install dependencies seleccionado y darle a Install. Alternativamente, puedes instalar el paquete con el comando:
install.packages("foreign")Se puede crear un proyecto en RStudio que permite guardar, y más tarde recuperar, el trabajo realizado en un determinado tema. Darle a menú File - New project, elegir New Directory, Empty project. En Directory name, escribir el nombre de la carpeta que se va a crear (‘Proyecto_curso’); en Create project as subdirectory of, elegir la carpeta donde ésta se va a guardar - en este caso, el escritorio. Darle a Create project. En el escritorio habrá ahora una carpeta llamada ‘Proyecto_curso’ y, dentro de ésta, un archivo llamado ‘Proyecto_curso.Rproj’. En esta carpeta vas a guardar los documentos del curso.
Para cargar (activar) los paquetes externos a utilizar en esta sesión de R, puedes marcarlos en la ventana inferior derecha de RStudio (pestaña Packages) o utilizar los siguientes comandos:
library(fuzzySim)
library(modEvA)
library(foreign)Para establecer el directorio de trabajo – es decir, la carpeta donde, por defecto, R va a buscar y guardar archivos – puedes utilizar el menú Session - Set Working Directory - Choose Directory y localizar la carpeta donde vas a trabajar, o utilizar el comando:
setwd("C:/Documents and settings/User/Escritorio/Proyecto_curso")(Reemplazar por la ruta de acceso a la carpeta en tu ordenador, que puedes obtener pinchando con el ratón derecho el nombre de un archivo de esa carpeta y dándole a Properties; utilizar barras normales (/) y no barras hacia atrás como las que usa Windows.) Por defecto, si has creado un proyecto en RStudio, el directorio de trabajo ya se ha definido automáticamente como el mismo donde se guardó el proyecto.
Importar la tabla con los datos a R. La tabla puede estar, por ejemplo, en formato .csv o .txt, con los valores (columnas) separados por comas, espacios, punto y coma, o tabuladores. Los valores decimales se recomienda que vayan separados por punto y no coma (ej. 0.999 y no 0,999). Nuestra tabla se llama ‘especiesmont10X10.csv’. En RStudio se puede utilizar el menú Tools - Import Dataset - From Text File y buscar la tabla en el disco. Otra opción es utilizar el comando:
datos <- read.table(file = "especiesmont10X10.csv", header = TRUE, sep = ",", dec = ".")El archivo (file) debe estar en la carpeta de trabajo (o incluir la ruta de acceso, separada por barras ‘/’) y escribirse entre comillas rectas y respetando las mayúsculas y minúsculas. El argumento header = TRUE indica que la primera fila de la tabla contiene los nombres de las columnas. El argumento sep indica el caracter que separa las columnas en la tabla a importar – en este caso son comas, pero podemos especificar otro separador, como espacio (sep = " "), punto y coma (sep = ";") o tabulador (sep = "\t"), cuando sea apropiado. El argumento dec indica el caracter que separa los decimales en la tabla a importar – en este caso son puntos, pero se pueden especificar comas cuando sea apropiado. Como para cualquer función de R, puedes obtener información con help(read.table).
La función read.table, las funciones relacionadas read.csv y read.delim (que difieren de ésta en los valores de los argumentos por defecto), leen tablas en formato de texto (.txt, .csv). Si la tabla a importar estuviera en formato .dbf, tendríamos que utilizar la función read.dbf del paquete foreign:
datos <- read.dbf(file = "especiesmont10X10.dbf")Una vez importada la tabla, puedes visualizarla, o algunas de sus características:
head(datos) # primeras filas de la tabla
str(datos) # estructura de la tabla
names(datos) # nombres de las columnasSi has creado un proyecto de RStudio y se guardan los cambios, no tienes que volver a importar estos datos en futuras sesiones. Cuando abres un proyecto de nuevo, puedes visualizar lo que has guardado pinchando en .Rhistory en la ventana inferior derecha de RStudio (pestaña Files).
Para modelar una sola especie con un modelo lineal generalizado (GLM), se puede utilizar la función glm de R. Ejemplo para un modelo univariante, i.e., con una sola variable (ej. la altitud) como predictora:
modelo.perdicera <- glm(HIE_FAS ~ ALTI, family = binomial, data = datos) # hace el modelo
modelo.perdicera # muestra el modeloPodemos crear también un modelo multivariante, con varias variables predictoras. Atención que al crear un objeto con el mismo nombre que otro que ya existía (en este caso, ‘modelo.perdicera’), éste se sobreescribe sin aviso:
names(datos) # ver en qué columnas están las variables predictoras
paste(names(datos)[21:45], collapse = " + ")
# copiar variables de la consola para utilizar en el comando siguiente:
modelo.perdicera <- glm(HIE_FAS ~ ORS + ORW + ALTI + ALTI2 + PEND + PEND2 +
DAUT + DIHE + DIPR + ESCO + ETP + ETR + HUEN + HUJU +
INSO + PERM + PM24 + PMR + PREC + RADS + TENE + TJUL +
TMED + U100 + U500, family = binomial, data = datos)
modelo.perdicera # muestra el modeloAntes de meter todas las variables en el modelo, se puede analizar la multicolinealidad con la función multicol de fuzzySim:
names(datos) # ver en qué columnas están las variables
mult <- multicol(subset(datos, select = 21:45)) # especificar sólo las columnas con variables!
mult # muestra la multicolinealidad
subset(mult, VIF > 10) # muestra cuáles tienen VIF elevadoTambién se puede hacer una pre-selección de variables que tengan una relación significativa con la distribución de la especie, utilizando el False Discovery Rate (FDR). Hay que especificar correctamente las columnas que contienen la especie y las variables en nuestra tabla:
names(datos) # ver en qué columnas están las especies y variables
FDR(datos, sp.cols = 8, var.cols = 21:45)Con la función corSelect se puede también seleccionar, de cada par de variables correlacionadas por en cima de determinado umbral (por ejemplo, 0.8), la que tenga una relación más significativa con la distribución de la especie:
names(datos) # ver en qué columnas están las especies y variables
corSelect(datos, sp.cols = 8, var.cols = 21:45, cor.thresh = 0.8) Luego se hace el modelo sin las variables excluidas (por ejemplo, en el caso del FDR, las variables ‘ETR’, ‘U100’ y ‘PERM’). Sobre éstas, se puede aplicar también un procedimiento de selección por pasos, basado en el criterio de información de Akaike (AIC), con la función step de R:
modelo.perdicera <- step(glm(HIE_FAS ~ ORS + ORW + ALTI + ALTI2 + PEND +
PEND2 + DAUT + DIHE + DIPR + ESCO + ETP + HUEN +
HUJU + INSO + PM24 + PMR + PREC + RADS + TENE +
TJUL + TMED + U500, family = binomial,
data = datos))
modelo.perdicera # muestra el modelo
summary(modelo.perdicera) # muestra los estadísticos de los coeficientesA menudo (como es el caso), la selección por AIC deja en el modelo variables no significativas. Éstas se pueden eliminar, también en un procedimiento por pasos hacia atrás, con la función modelTrim del paquete fuzzySim:
modelo.perdicera <- modelTrim(modelo.perdicera) # elimina las no significativas
modelo.perdicera # muestra el modelo
summary(modelo.perdicera) # muestra los estadísticos de los coeficientesPara obtener las predicciones del modelo y la favorabilidad:
Hiefas_P <- modelo.perdicera$fitted.values
Hiefas_F <- Fav(model = modelo.perdicera)
Hiefas_P
Hiefas_FLa función multGLM del paquete fuzzySim permite hacer modelos para una o más especies a la vez, con la posibilidad de selección de variables por correlación (corSelect), FDR, AIC (step) y/o significación (modelTrim), y guarda los valores predichos automáticamente, incluida la favorabilidad y (si lo especificamos) la Y. En esta función las variables se especifican de forma distinta que en la función glm:
names(datos) # ver en qué columnas están las especies y las variables
mis.modelos <- multGLM(data = na.omit(datos), sp.cols = 8:11, var.cols = 21:45,
id.col = 1, step = TRUE, FDR = TRUE, corSelect = FALSE,
trim = TRUE, Y = TRUE, Fav = TRUE)Si las variables y/o especies a modelar no estuvieran contiguas en la tabla, se pueden especificar por separado:
mis.modelos.2 <- multGLM(data = na.omit(datos), sp.cols = c(8, 10, 11),
var.cols = c(21, 24:30, 41, 45), id.col = 1,
step = TRUE, FDR = TRUE, corSelect = FALSE, trim = TRUE,
Y = TRUE, Fav = TRUE) El objeto resultante de multGLM (‘mis.modelos’) es una lista que contiene dos objetos: otra lista llamada ‘models’ con los modelos, y una tabla llamada ‘predictions’ con sus predicciones:
mis.modelos$models # muestra los modelos
names(mis.modelos$models) # muestra los nombres de los modelos
head(mis.modelos$predictions) # primeras filas de la tabla de predichosEl sumario del modelo contiene mucha información interesante:
summary(modelo.perdicera) # modelo creado por la función 'glm'
summary(mis.modelos$models$HIE_FAS) # modelo de la lista creada por 'multGLM'
summary(mis.modelos$models$CAP_PYR)
summary(mis.modelos$models$ALY_DIC)
summary(mis.modelos$models$VIP_LAT)Podemos atribuir nombres más cortos a los modelos de multGLM, para que sea más fácil utilizarlos más adelante:
mod.hiefas <- mis.modelos$models$HIE_FAS
mod.cappyr <- mis.modelos$models$CAP_PYR
mod.alydic <- mis.modelos$models$ALY_DIC
mod.viplat <- mis.modelos$models$VIP_LATEl AIC del modelo viene incluido en el summary, pero puedes obtener su valor individual (y más exacto) con summary(mod.hiefas)$aic.
Para obtener la ecuación final del modelo, puedes utilizar la función getModEqn del paquete modEvA. Por defecto, esta función nos da la ecuación de la Y (predictor lineal, en la escala de las variables), pero podemos pedirle la ecuación de la P (probabilidad) o de la F (favorabilidad, que conlleva un cambio en el intercepto):
getModEqn(modelo.perdicera)
getModEqn(mod.hiefas)
getModEqn(mod.hiefas, digits = 3) # para redondear a menos dígitos
getModEqn(mod.hiefas, digits = 3, type = "P")
getModEqn(mod.hiefas, digits = 3, type = "F")Se puede también visualizar la gráfica del modelo. El argumento ‘main’ define el título principal de la gráfica:
plotGLM(mod.hiefas)plotGLM(mod.hiefas, main = "Modelo perdicera")plotGLM(mod.cappyr, main = "Modelo cabra")plotGLM(mod.alydic, main = "Modelo partero")plotGLM(mod.viplat, main = "Modelo víbora")Calcular la cantidad de devianza explicada (D cuadrado):
Dsquared(mod.hiefas)
Dsquared(mod.hiefas, adjust = TRUE) # ajustado al nº de parámetrosCalcular distintos pseudo R cuadrados del modelo:
RsqGLM(mod.hiefas)Las medidas de discriminación evalúan la capacidad del modelo de distinguir las zonas de presencia de las zonas de ausencia de la especie. Dado que los valores predichos son continuos, estas medidas requieren definir un umbral a partir del cual se considera que el modelo predice que la especie estará presente. La función threshMeasures calcula varias medidas de discriminación para un valor umbral dado:
threshMeasures(model = mod.hiefas) # el umbral es 0.5 por defectothreshMeasures(model = mod.hiefas, ylim = c(0, 1)) # escala el eje ythreshMeasures(model = mod.hiefas, ylim = c(0, 1), thresh = 0.75)threshMeasures(model = mod.hiefas, ylim = c(0, 1), thresh = "preval")Para encontrar el valor umbral que optimiza cada medida de discriminación, se puede usar la función optiThresh:
optiThresh(model = mod.hiefas)También se puede optimizar el umbral para la mejor combinación de una pareja de medidas de evaluación relacionadas:
optiPair(model = mod.hiefas)# por defecto las dos medidas son sensibilidad y especificidad
# pero se puede especificar otras:
optiPair(model = mod.hiefas, measures = c("Omission", "Commission"))optiPair(model = mod.hiefas, measures = c("UPR", "OPR"))El área bajo la curva (AUC) resume la capacidad de discriminación del modelo sobre todos los umbrales posibles:
AUC(model = mod.hiefas)Las medidas de calibración miden el ajuste de las probabilidades predichas por el modelo a las frecuencias observadas en los datos. La función HLfit calcula la bondad de ajuste con el estadístico de Hosmer-Lemeshow. La significación se refiere a la diferencia entre observado y predicho; por tanto, cuanto más alta la p, mejor. Esta función aporta también la raiz cuadrada de la media de los cuadrados de los errores (RMSE). Atención que el método elegido para agrupar los datos (argumento bin.method) puede afectar bastante a los resultados:
HLfit(model = mod.hiefas, bin.method = "prob.bins")HLfit(model = mod.hiefas, bin.method = "n.bins", n.bins = 20)HLfit(model = mod.hiefas, bin.method = "quantiles")También se puede calcular los estadísticos de calibración de Miller (intercepto y pendiente, que idealmente deben ser respectivamente 0 y 1), pero estos no sirven para evaluar un modelo sobre los mismos datos utilizados en su construcción (el resultado siempre sale bueno), sino para ver cómo se aplica a otro conjunto de datos distinto:
MillerCalib(model = mod.hiefas)Consulta help(HLfit) y help(MillerCalib) para más información sobre estas funciones, las opciones disponibles y la interpretación de sus resultados.
La función multModEv permite evaluar a la vez todos los modelos de una lista dada (por ejemplo, la lista generada por la función multGLM, aunque contenga un sólo modelo), también con distintas opciones:
model_eval <- multModEv(models = mis.modelos$models)
model_eval <- multModEv(models = mis.modelos$models, thresh = "preval",
bin.method = "quantiles")
model_evalSe puede distinguir la cantidad de varianza del modelo atribuible a distintos factores y a sus combinaciones. La función varPart del paquete modEvA actualmente hace esto para hasta 3 factores - por ejemplo ambiental, humano y topográfico, pero hay que hacer alguna preparación antes. Primero, se define las variables relacionadas con cada factor:
names(datos) # ver en qué columna está cada variable, para elegirlas:
amb <- names(datos)[c(28:43)]
hum <- names(datos)[c(27, 44:45)]
top <- names(datos)[c(21:26)]
amb # variables del factor ambiental
hum # variables del factor humano
top # variables del factor topográficoDespués, ver qué variables están incluidas en el modelo y definir, de entre éstas, las que pertenecen a cada factor:
mod.vars <- names(mod.hiefas$coefficients)[-1]
amb.vars <- mod.vars[mod.vars %in% amb]
hum.vars <- mod.vars[mod.vars %in% hum]
top.vars <- mod.vars[mod.vars %in% top]
amb.vars # variables ambientales en el modelo
hum.vars # variables humanas en el modelo
top.vars # variables topográficas en el modeloAhora hay que hacer modelos parciales con las variables de cada factor (sin selección, incluirlas todas) y de cada par de factores:
paste(amb.vars, collapse = " + ")
paste(hum.vars, collapse = " + ")
paste(top.vars, collapse = " + ")
# copiar las variables de la consola para los comandos siguientes:
mod.amb <- glm(HIE_FAS ~ INSO + ESCO + DIHE + TMED + TJUL + ETP + PREC,
family = binomial, data = datos)
mod.hum <- glm(HIE_FAS ~ U500, family = binomial, data = datos)
mod.top <- glm(HIE_FAS ~ PEND + ALTI2 + PEND2 + ORW, family = binomial,
data = datos)
mod.amb.hum <- glm(HIE_FAS ~ INSO + ESCO + DIHE + TMED + TJUL + ETP + PREC +
U500, family = binomial, data = datos)
mod.amb.top <- glm(HIE_FAS ~ INSO + ESCO + DIHE + TMED + TJUL + ETP + PREC +
PEND + ALTI2 + PEND2 + ORW, family = binomial, data = datos)
mod.hum.top <- glm(HIE_FAS ~ U500 + PEND + ALTI2 + PEND2 + ORW,
family = binomial, data = datos)Para GLMs, hay varias formas de calcular la partición de varianza (ver Details en help(varPart)). Una de ellas implica obtener la favorabilidad dada por cada modelo parcial y calcular su correlación con la del modelo total:
Hiefas_F # calculada en la sección 3.1
mod.amb.F <- Fav(model = mod.amb)
mod.hum.F <- Fav(model = mod.hum)
mod.top.F <- Fav(model = mod.top)
mod.amb.hum.F <- Fav(model = mod.amb.hum)
mod.amb.top.F <- Fav(model = mod.amb.top)
mod.hum.top.F <- Fav(model = mod.hum.top)
favs <- cbind(Hiefas_F, mod.amb.F, mod.hum.F, mod.top.F, mod.amb.hum.F,
mod.amb.top.F, mod.hum.top.F)
corrs <- cor(favs, method = "spearman")[,1]Por fin, podemos calcular y representar la partición de varianza del modelo:
varPart(A = corrs["mod.amb.F"], B = corrs["mod.hum.F"], C = corrs["mod.top.F"],
AB = corrs["mod.amb.hum.F"], AC = corrs["mod.amb.top.F"],
BC = corrs["mod.hum.top.F"], ABC = 1, model.type = "GLM",
A.name = "ambiental", B.name = "humano", C.name = "topográfico")Para exportar las predicciones de los modelos y/o los resultados de su evaluación, se puede utilizar varios formatos:
write.table(mis.modelos$predictions, file = "predichos.csv", row.names = FALSE)
library(foreign) # necesario para formato dbf
write.dbf(mis.modelos$predictions, file = "predichos.dbf")El archivo que se va a guardar se puede nombrar como se quiera. Se recomienda usar el formato .dbf para representarlo en QGIS; los formatos .txt y .csv son de texto y no se pueden representar numéricamente de forma automática.
Los resultados de la evaluación múltiple se pueden exportar también:
write.csv(model_eval, file = "evaluacion.csv", row.names = FALSE)