setwd("~/Dropbox/Proyecto_Invasión_Neusa/data")
# Librerias necesarias
library("aod", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
library("arm", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
##
## arm (Version 1.10-1, built: 2018-4-12)
## Working directory is /Users/mazinger/Dropbox/Proyecto_Invasión_Neusa/data
library("ggplot2", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
library("tidyverse", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## ── Attaching packages ──────────────────────────────── tidyverse 1.2.1 ──
## ✔ tibble 2.0.1 ✔ purrr 0.3.1
## ✔ tidyr 0.8.3 ✔ dplyr 0.8.0.1
## ✔ readr 1.3.1 ✔ stringr 1.4.0
## ✔ tibble 2.0.1 ✔ forcats 0.3.0
## ── Conflicts ─────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
library("caret", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
library("rms", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: Hmisc
## Loading required package: survival
##
## Attaching package: 'survival'
## The following object is masked from 'package:caret':
##
## cluster
## The following object is masked from 'package:aod':
##
## rats
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
##
## backsolve
library("car", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: carData
##
## Attaching package: 'car'
## The following objects are masked from 'package:rms':
##
## Predict, vif
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:arm':
##
## logit
library("rcompanion", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
library("lmtest", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'lmtest'
## The following object is masked from 'package:rms':
##
## lrtest
# Lectura base de datos para Retamo liso
invasion3<-read.table("retlis.csv", sep = ",", header = T)
attach(invasion3)
# Conversión a factor para variables categóricas
via<-factor(via)
lado.cam<-factor(lado.cam)
cobertura<-factor(cobertura)
pendiente<-factor(pendiente)
nombre<-factor(nombre)
rep<-factor(rep)
invas.tipo<-factor(invas.tipo)
tipol<-factor(tipol)
# Modelo final
invaslogit.low3<-glm(invas.low~alt.prom+borde+msnm.prom+msnm.rango+via+rep+invas.tipo+tipol, family = "binomial")
summary(invaslogit.low3)
##
## Call:
## glm(formula = invas.low ~ alt.prom + borde + msnm.prom + msnm.rango +
## via + rep + invas.tipo + tipol, family = "binomial")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5045 -1.0556 0.5699 0.9000 2.2810
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.287e+01 8.827e+02 -0.015 0.988
## alt.prom 1.023e+00 2.296e-01 4.453 8.46e-06 ***
## borde 3.403e-02 3.008e-02 1.131 0.258
## msnm.prom -9.496e-04 7.148e-04 -1.328 0.184
## msnm.rango 4.608e-03 4.611e-03 1.000 0.318
## via2 -3.930e-01 3.523e-01 -1.115 0.265
## via3 5.040e-01 4.024e-01 1.253 0.210
## rep2 -3.613e-01 1.072e+03 0.000 1.000
## rep3 1.167e+01 8.827e+02 0.013 0.989
## rep5 1.453e+01 8.827e+02 0.016 0.987
## invas.tipo2 -3.073e-01 3.671e-01 -0.837 0.403
## tipol2 -6.120e-02 3.881e-01 -0.158 0.875
## tipol3 8.972e-01 5.868e-01 1.529 0.126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 609.21 on 457 degrees of freedom
## Residual deviance: 510.97 on 445 degrees of freedom
## AIC: 536.97
##
## Number of Fisher Scoring iterations: 13
Anova(invaslogit.low3, type="II", test="Wald")
## Analysis of Deviance Table (Type II tests)
##
## Response: invas.low
## Df Chisq Pr(>Chisq)
## alt.prom 1 19.8312 8.459e-06 ***
## borde 1 1.2799 0.257913
## msnm.prom 1 1.7646 0.184054
## msnm.rango 1 0.9990 0.317549
## via 2 11.9736 0.002512 **
## rep 3 7.3630 0.061185 .
## invas.tipo 1 0.7006 0.402578
## tipol 2 2.4287 0.296908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Odds ratios
exp(cbind(OR = coef(invaslogit.low3), confint(invaslogit.low3)))
## Waiting for profiling to be done...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## OR 2.5 % 97.5 %
## (Intercept) 2.584459e-06 NA 3.156731e+71
## alt.prom 2.780564e+00 1.793980e+00 4.419473e+00
## borde 1.034616e+00 9.872263e-01 1.086581e+00
## msnm.prom 9.990509e-01 9.979293e-01 1.000445e+00
## msnm.rango 1.004619e+00 9.967279e-01 1.014735e+00
## via2 6.750267e-01 3.333582e-01 1.335448e+00
## via3 1.655401e+00 7.459356e-01 3.634837e+00
## rep2 6.967725e-01 5.181001e-05 7.897878e+03
## rep3 1.165814e+05 1.342143e-61 NA
## rep5 2.045951e+06 1.646176e-72 NA
## invas.tipo2 7.354568e-01 3.478806e-01 1.481536e+00
## tipol2 9.406391e-01 4.426889e-01 2.045196e+00
## tipol3 2.452782e+00 8.470507e-01 8.940355e+00
# Pseudo R2
nagelkerke(invaslogit.low3)
## $Models
##
## Model: "glm, invas.low ~ alt.prom + borde + msnm.prom + msnm.rango + via + rep + invas.tipo + tipol, binomial"
## Null: "glm, invas.low ~ 1, binomial"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.161268
## Cox and Snell (ML) 0.193065
## Nagelkerke (Cragg and Uhler) 0.262472
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -12 -49.123 98.247 1.2269e-15
##
## $Number.of.observations
##
## Model: 458
## Null: 458
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"
# Significancia del modelo (siempre y cuando significativa, la adición de las variables en comparación a un modelo nulo puede explicar la magnitud de la invasión)
lrtest(invaslogit.low3)
## Likelihood ratio test
##
## Model 1: invas.low ~ alt.prom + borde + msnm.prom + msnm.rango + via +
## rep + invas.tipo + tipol
## Model 2: invas.low ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 13 -255.48
## 2 1 -304.61 -12 98.247 1.227e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Sobredispersión (no debe ser superior a 1.5)
summary(invaslogit.low3)$deviance / summary(invaslogit.low3)$df.residual
## [1] 1.148241
#examen gráfico de los residuos
plot(predict(invaslogit.low3),residuals(invaslogit.low3),col=c("blue","red")[1+invas.low])
abline(h=0,lty=2,col="grey")
lines(lowess(predict(invaslogit.low3),residuals(invaslogit.low3)),col="black",lwd=2)
# El modelo muestra que las variables seleccionadas en el conjunto de datos completo (con las dos especies incluidas)
#no son adecuadas para predecir la invasion de retamo liso, el pseudo R2 del Retamo espinoso es menor al R2 del modelo con las dos especies
#además, el modelo de Retamo liso muestra predicción perfecta, por tanto se analiza el n disponible en las variables categoricas para hacer la predicción
xtabs(~invas.low+via)
## via
## invas.low 1 2 3
## 0 20 117 38
## 1 46 127 110
xtabs(~invas.low+rep)
## rep
## invas.low 1 2 3 5
## 0 1 2 19 153
## 1 0 0 1 282
# No existen o existen muy pocas observaciones para los niveles 1,2 y 3 de rep (estados reproductivos) para hacer una predicción adecuada esta es la razón por la cual el modelo
# muestra una predicción perfecta
xtabs(~invas.low+invas.tipo)
## invas.tipo
## invas.low 1 2
## 0 14 161
## 1 38 245
xtabs(~invas.low+tipol)
## tipol
## invas.low 1 2 3
## 0 156 15 4
## 1 230 29 24
#Conclusión: El modelo de Retamo liso muestra que la presencia de esta especie no conlleva a magnitudes de invasión superiores a 10 m2. El modelo usando las dos especies puede predecir mejor la magnitud de la invasión ya que no presenta prediciones perfectas.