Lichen

Setup

Load packages

##paquetes
library(readxl)
library(ggplot2)
library(data.table)
library(tidyr)
library(dplyr)
library(purrr)
library(lattice)
library(boot)
library(nlme)
library(MASS)
library(lme4)
library(lsmeans)
library(vegan)
library(rjags)

Load data

##Datos
morfo <- read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vTpy_uV7A-fV7cgZ_uJ0o3_m9nrC_hTrljJCPgrzZKYq6nYV-51XyI79IpXr8GhaBpvQl-zl1qjoYZX/pub?output=csv")
foro <- read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vRB8FvUGkqvWsdylid0DTXZcRbWtkWG7xZcEurreoXMnRPg3AZ_rggTiHZwAYKCK6rKGfqu9ADAIYgk/pub?output=csv")
zona <- read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vQvJDYjxqSXT3URQcsREYKBMreV2fxiKhTNksEziaMMIxSxCO3rQxrFwq32dlk33dPvK9VTM3k2o_yO/pub?output=csv")

## Warning in read.table(file = file, header = header, sep
## = sep, quote = quote, : incomplete final line found by
## readTableHeader on 'https://docs.google.com/spreadsheets/d/e/
## 2PACX-1vQvJDYjxqSXT3URQcsREYKBMreV2fxiKhTNksEziaMMIxSxCO3rQxrFwq32dlk33dPvK9VTM3k2o_yO/
## pub?output=csv'

dat <- merge(morfo, foro)
data <- merge(zona, dat) #datos
data <- data.table(data)
data[, pre := rep(1, nrow(data))]
morfo1 <- with(data, unique(Morfo))
label <- 1:length(morfo1)
morfo2 <- data.table("Morfo" = morfo1, "Label" = label)
data <- merge(data, morfo2)

---

Part 1: Data

los datos muestran la presencia de morfos en diferentes forofitos

dat <- data %>% group_by(Nombre, Forofito, ID_Forofito) %>% summarise(Morfos = sum(pre), CAP = mean(CAP), .groups = "drop")
dat1 <- data %>% group_by(Nombre, Forofito, ID_Forofito) %>% summarise(Riqueza = sum(pre), CAP = mean(CAP), .groups = "drop")
dat2 <- dat1 %>% group_by(Nombre, Forofito) %>% summarise(Promedio = mean(Riqueza), Desviacion = sd(Riqueza), CAP = mean(CAP), .groups = "drop")
dat3 <- dat1 %>% group_by(Nombre, Forofito) %>% summarise(Promedio = mean(Riqueza), Desviacion = sd(Riqueza), CAP = mean(CAP), .groups = "drop")
dat10 <- data[ , lapply(.SD, sum), .SDcols = "pre", by = c("Nombre","Forofito", "Label", "Morfo")]
dat11 <- with(dat10, split(dat10, Forofito))
datt <- data[ , lapply(.SD, sum), .SDcols = "pre", by = "Morfo"]
todo <- datt$pre

---

Part 2: Research question

factores como la especie de forofito o la zona influyen en el numero de morfos encontrados por forofito?

se muestrearon dos zonas un colina con mucho viento y el valle con poco viento en general en l colina los forofitos de la misma especie tenian mayor numero de morfos

d <- rep(dat1$Nombre, dat1$Riqueza)
d2 <- rep(dat1$Forofito, dat1$Riqueza)
full <- data.frame(d, d2)
colnames(full) <- c("Zona", "Forofito")

ggplot(data = full, mapping = aes(x = Forofito, fill = Zona)) +
  geom_bar(alpha = 0.5, position = "fill", width = 0.5) +
  scale_fill_manual(values = c("khaki2", "olivedrab")) +#honeydew3, khaki2
  labs(y = "Proporción") +
  theme_classic()

resumen de datos por Zona

dat2

## # A tibble: 10 x 5
##    Nombre Forofito     Promedio Desviacion   CAP
##    <chr>  <chr>           <dbl>      <dbl> <dbl>
##  1 Colina Opuntia           4.3       2.26  26.8
##  2 Colina Pilosocereus      5.1       1.85  42.7
##  3 Colina Stenocereus       5         1.76  47.8
##  4 Colina Vachellia         8.8       2.74  46.3
##  5 Colina Zanthoxylum       7.5       1.65  25.9
##  6 Valle  Opuntia           3.8       1.23  23.3
##  7 Valle  Pilosocereus      4.1       1.66  44.6
##  8 Valle  Stenocereus       4.3       1.25  25.3
##  9 Valle  Vachellia         6.2       2.39  44.4
## 10 Valle  Zanthoxylum       6.5       2.22  20.5

---

Part 3: Exploratory data analysis

boxplot de los datos

dat <- dat %>%
  mutate(Tipo = ifelse(Forofito %in% c("Opuntia", "Pilosocereus", "Stenocereus"), "Cactus", ifelse(Forofito %in% c("Vachellia", "Zanthoxylum"), "Árbol", "Trambo")))

ggplot(dat, aes(Morfos, Forofito, col = Tipo, fill = Tipo)) + 
  geom_violin(scale = "area", alpha = 0.5 ) +
  scale_fill_manual(values = c("gold3", "forestgreen")) +
  scale_colour_manual(values = c("gold3", "forestgreen")) +
  theme_minimal()

densidad de circunferencias

with(dat, densityplot(~CAP|Forofito, groups = Tipo, par.settings = list(superpose.symbol = list(alpha = 1.0, col = c("gold3", "forestgreen"), pch = 19),
                                                                        superpose.line = list(alpha = 0.5, col = c("gold3", "forestgreen")))))

con exepcion de Vachellia la densidad por especie es similar

---

Part 4: Modeling

se hicieron varios modelos

primero uno simple

m1 <- lm(Promedio ~ Forofito + Nombre, data = dat2)#el modelo con mas soporte tanton en BIC como en AIC es el que no tiene CAP
shapiro.test(resid(m1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.98417, p-value = 0.9836

qqnorm(resid(m1))
qqline(resid(m1))

pero el modelo simple presume normalidad y deben ser continuos pero los datos no son continuos por este motivo y como son un conteo tambien hicimos un modelo generaliado poisson

m2 <- glm(formula = Riqueza ~ Forofito + Nombre, family = poisson, 
    data = dat1)

m2$deviance/m2$df.residual#datos subdispersos

## [1] 0.670312

tabien hicimos un modelo mixto con el forofito como variable aleatoria no la especie del forofito si ni el forofito elegido para muestrear

m3 <- glmer(Riqueza ~ (1|ID_Forofito) + Forofito + Nombre, family = "poisson", data = dat1)

## boundary (singular) fit: see ?isSingular

para elegir el modelo se evaluo el AIC y el BIC

AIC(m2, m3)

##    df      AIC
## m2  6 423.9348
## m3  7 425.9348

BIC(m2, m3)

##    df      BIC
## m2  6 439.5658
## m3  7 444.1710

summary(m2)

## 
## Call:
## glm(formula = Riqueza ~ Forofito + Nombre, family = poisson, 
##     data = dat1)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.98723  -0.52586  -0.09002   0.41927   1.87937  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           1.49794    0.11749  12.749  < 2e-16 ***
## ForofitoPilosocereus  0.12734    0.15237   0.836   0.4033    
## ForofitoStenocereus   0.13815    0.15198   0.909   0.3634    
## ForofitoVachellia     0.61619    0.13788   4.469 7.86e-06 ***
## ForofitoZanthoxylum   0.54719    0.13960   3.920 8.87e-05 ***
## NombreValle          -0.20939    0.08528  -2.455   0.0141 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 103.895  on 99  degrees of freedom
## Residual deviance:  63.009  on 94  degrees of freedom
## AIC: 423.93
## 
## Number of Fisher Scoring iterations: 4

con <- emmeans(m2, ~Forofito, method = "tukey")
contrast(con)

##  contrast            estimate     SE  df z.ratio p.value
##  Opuntia effect        -0.286 0.0965 Inf -2.960  0.0051 
##  Pilosocereus effect   -0.158 0.0918 Inf -1.725  0.1056 
##  Stenocereus effect    -0.148 0.0914 Inf -1.614  0.1065 
##  Vachellia effect       0.330 0.0769 Inf  4.298  0.0001 
##  Zanthoxylum effect     0.261 0.0787 Inf  3.321  0.0022 
## 
## Results are averaged over the levels of: Nombre 
## Results are given on the log (not the response) scale. 
## P value adjustment: fdr method for 5 tests

coefficients(m2)

##          (Intercept) ForofitoPilosocereus  ForofitoStenocereus 
##            1.4979435            0.1273394            0.1381503 
##    ForofitoVachellia  ForofitoZanthoxylum          NombreValle 
##            0.6161861            0.5471933           -0.2093949

el modelo uno no lo pruebo debido a basarse en promedios y por esto tener otros datos solo lo hicimos por curiosidad el modelo con mas soporte es el fijo y este muestra diferencias significativas

###bootraping

tambien se hicieron boots del numero de morfos promedio por forofito

set.seed(950728)
pilo <- function(data,i){
        d<-data[i,]
        d1 <- d %>% filter(Forofito == "Pilosocereus")
        return(with(d1, c(mean(Riqueza), sd(Riqueza))))
}

bootpi <- boot(dat1, pilo, R = 10000)
pilo1 <- boot.ci(bootpi)

va <- function(data,i){
        d<-data[i,]
        d1 <- d %>% filter(Forofito == "Vachellia")
        return(with(d1, c(mean(Riqueza), sd(Riqueza))))
}

bootva <- boot(dat1, va, R = 10000)
vai1 <- boot.ci(bootva)

opu <- function(data,i){
        d<-data[i,]
        d1 <- d %>% filter(Forofito == "Opuntia")
        return(with(d1, c(mean(Riqueza), sd(Riqueza))))
}

bootopu <- boot(dat1, opu, R = 10000)
opu1 <- boot.ci(bootopu)

ste <- function(data,i){
        d<-data[i,]
        d1 <- d %>% filter(Forofito == "Stenocereus")
        return(with(d1, c(mean(Riqueza), sd(Riqueza))))
}

bootste <- boot(dat1, ste, R = 10000)
ste1 <- boot.ci(bootste)

zan <- function(data,i){
        d<-data[i,]
        d1 <- d %>% filter(Forofito == "Zanthoxylum")
        return(with(d1, c(mean(Riqueza), sd(Riqueza))))
}
bootzan <- boot(dat1, zan, R = 10000)
zan1 <- boot.ci(bootzan)

valor p del boot entre las comparaciones en los boots

po <- sum(bootpi$t[,1]>=bootopu$t[,1])/bootopu$R #.871
ps <- sum(bootpi$t[,1]>=bootste$t[,1])/bootopu$R #.2837
pv <- sum(bootpi$t[,1]>=bootva$t[,1])/bootopu$R #0
pz <- sum(bootpi$t[,1]>=bootzan$t[,1])/bootopu$R #2e-04
so <- sum(bootste$t[,1]>=bootopu$t[,1])/bootopu$R #.9653
sv <- sum(bootste$t[,1]>=bootva$t[,1])/bootopu$R #0
sz <- sum(bootste$t[,1]>=bootzan$t[,1])/bootopu$R #6e-04
ov <- sum(bootopu$t[,1]>=bootva$t[,1])/bootopu$R #0
oz <- sum(bootopu$t[,1]>=bootzan$t[,1])/bootopu$R #0
vz <- sum(bootva$t[,1]>=bootzan$t[,1])/bootopu$R #.533
po#pilo vs opu

## [1] 0.8463

ps#pilo vs ste

## [1] 0.4736

pv#pilo vs vac

## [1] 0

so#ste vs opu

## [1] 0.8799

sv#ste vs vac

## [1] 0

sz#ste vs zan

## [1] 1e-04

ov#opu vs vac

## [1] 0

oz#opu vs zan

## [1] 0

vz#vac vs zan

## [1] 0.7402

igualmente se encontro valores menores entre la comparacion entre cactus vs arboles

###Herachical modeling

hicimos un modelo jerarquico beyesiano similar al modelo generalizado

mod_string1 = "model {
for (i in 1:length(Riqueza)) {
        Riqueza[i] ~ dpois(lam[Forofito[i]])
        
    }
    
    for (j in 1:max(Forofito)) {
        lam[j] ~ dgamma(alfa, beta)
    }

    mu ~ dgamma(2.0, 1.0/2.6)
    sig ~ dexp(1.0)

    alfa = mu^2 / sig^2
    beta = mu / sig^2
}"



params1 <- c("lam", "mu", "sig")
Forofito <- as.numeric(as.factor(dat1$Forofito))
Riqueza <- dat1$Riqueza
data_jags1 <- data.frame(Riqueza, Forofito)
data_jags1 <- as.list(data_jags1)
mod1 <- jags.model(textConnection(mod_string1), data = data_jags1, n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 7
##    Total graph size: 215
## 
## Initializing model

update(mod1, 1e3)

mod_sim <- coda.samples(model = mod1,
                        variable.names = params1,
                        n.iter = 2e4)

mod_csim <- as.mcmc(do.call(rbind, mod_sim)) # combined chains
summary(mod_csim)

## 
## Iterations = 1:60000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 60000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean     SD Naive SE Time-series SE
## lam[1] 4.263 0.4517 0.001844       0.002161
## lam[2] 4.732 0.4647 0.001897       0.002085
## lam[3] 4.778 0.4671 0.001907       0.002057
## lam[4] 7.214 0.5910 0.002413       0.002947
## lam[5] 6.785 0.5634 0.002300       0.002653
## mu     5.542 0.7125 0.002909       0.004768
## sig    1.485 0.5775 0.002358       0.005184
## 
## 2. Quantiles for each variable:
## 
##          2.5%   25%   50%   75% 97.5%
## lam[1] 3.4180 3.951 4.250 4.562 5.185
## lam[2] 3.8523 4.412 4.722 5.040 5.669
## lam[3] 3.8970 4.457 4.766 5.086 5.733
## lam[4] 6.1120 6.805 7.196 7.601 8.426
## lam[5] 5.7314 6.397 6.766 7.157 7.938
## mu     4.1963 5.094 5.514 5.954 7.062
## sig    0.6894 1.086 1.377 1.759 2.918

los cuantiles bayesianos muestran los mismo en generl en los cactus(lam1, lam2, lam3) se espera contar menos morfos

---

##UNTB

por ultimo hicimos un modelo neutral con nuestros datos

##Simulacion UNTB

neutral3 = function(data, datam, x, theta, nsim = 10, D = 1, cycles = 1e4) {
  d <- sort(data[[x]]$pre)
  label <- sort(data[[x]]$Label)
  J <- sum(d)
  S <- length(d)
  
  dm <- sort(datam)
  Jm <- sum(dm)
  Sm <- length(dm)
  nu <- theta/(theta + sum(dm)*10 - 1)
  m <- 1.0
  comu <- matrix(nrow = J, ncol = nsim)
  wc <- matrix(nrow = J, ncol = nsim)
  Jm <- Jm*10
  meta <- matrix(nrow = Jm, ncol = nsim)
  wm <- matrix(nrow = Jm, ncol = nsim)
  comu[,1] <- rep(label, d)
  meta[,1] <- rep(1:Sm, dm*10)
  
  for ( i in 2:nsim){   
    k = i - 1
    
    cod.sp <- comu[,k]
    meta.sp <- meta[,k] 
    
    deathk <- sample(1:J, D)
    meta.deathk <- sample(1:Jm, D)
    ##Births                    
    outside <- sample(c(T, F), size = D, replace = T, prob = c(m, 1-m))
    
    if(runif(1) < nu){
      meta.outside <- T
    } else {meta.outside <- F}
    
    birthkin <- sample(1:J, D - sum(outside), replace = T)
    meta.birthkin <- sample(1:Jm, D - sum(meta.outside))
    
    birthkout <- sample(1:Jm, sum(outside), replace = T)
    meta.birthkout <- sample(1:Jm, sum(meta.outside), replace = T)              
    
    ##Subtitutions
    if(length(birthkin) > 0){
      cod.sp[deathk[!outside]] <- cod.sp[birthkin]
    }
    
    if(length(birthkout) > 0){
      cod.sp[deathk[outside]] <- meta.sp[birthkout]
    }
    if(length(meta.birthkin) > 0){
      meta.sp[meta.deathk[!meta.outside]] <- meta.sp[meta.birthkin]
    }
    
    if(length(meta.birthkout) > 0){
      meta.sp[meta.deathk[meta.outside]] <- max(meta.sp) + 1
      
    }
    cod.sp -> comu[,i]
    meta.sp -> meta[,i]
    
  }
  return(comu)}

j <- c(81, 92, 93, 150, 140, 556)
masi <- c(8.869934, 16.01548, 10.15438, 11.58415, 16.07991, 14.96903)
set.seed(950728)
dat10 <- data[ , lapply(.SD, sum), .SDcols = "pre", by = c("Nombre","Forofito", "Label", "Morfo")]
dat11 <- with(dat10, split(dat10, Forofito))
datt <- data[ , lapply(.SD, sum), .SDcols = "pre", by = "Morfo"]
todo <- datt$pre
nsim <- 1e3
pri <- neutral3(data = dat11, datam = todo, x = 1, theta = 14.96903, nsim = nsim)
tab <- apply(pri, 2, table)

opu <- tab[nsim:(nsim-99)]
largo <- c()
for(i in 1:100){
  largo[i] <- length(opu[[i]])
}
labs <- paste0(rep("Opuntia", 100), rep(1:100, largo))
opu1 <- unlist(opu)
spp <- names(opu1)
prueba <- cbind(labs, spp, opu1)
colnames(prueba) <- c("Simulacion", "Especie", "Abundancia")

pri <- neutral3(data = dat11, datam = todo, x = 2, theta = 14.96903, nsim = nsim)
tab <- apply(pri, 2, table)

pilo <- tab[nsim:(nsim-99)]
largo <- c()
for(i in 1:100){
  largo[i] <- length(pilo[[i]])
}
labs <- paste0(rep("Pilosocereus", 100), rep(1:100, largo))
pilo1 <- unlist(pilo)
spp <- names(pilo1)
prueba1 <- cbind(labs, spp, pilo1)
colnames(prueba1) <- c("Simulacion", "Especie", "Abundancia")

pri <- neutral3(data = dat11, datam = todo, x = 3, theta = 14.96903, nsim = nsim)
tab <- apply(pri, 2, table)

ste <- tab[nsim:(nsim-99)]
largo <- c()
for(i in 1:100){
  largo[i] <- length(ste[[i]])
}
labs <- paste0(rep("Stenocereus", 100), rep(1:100, largo))
ste1 <- unlist(ste)
spp <- names(ste1)
prueba2 <- cbind(labs, spp, ste1)
colnames(prueba2) <- c("Simulacion", "Especie", "Abundancia")

pri <- neutral3(data = dat11, datam = todo, x = 4, theta = 14.96903, nsim = nsim)
tab <- apply(pri, 2, table)

vac <- tab[nsim:(nsim-99)]
largo <- c()
for(i in 1:100){
  largo[i] <- length(vac[[i]])
}
labs <- paste0(rep("Vachellia", 100), rep(1:100, largo))
vac1 <- unlist(vac)
spp <- names(vac1)
prueba3 <- cbind(labs, spp, vac1)
colnames(prueba3) <- c("Simulacion", "Especie", "Abundancia")

pri <- neutral3(data = dat11, datam = todo, x = 5, theta = 14.96903, nsim = nsim)
tab <- apply(pri, 2, table)

zan <- tab[nsim:(nsim-99)]
largo <- c()
for(i in 1:100){
  largo[i] <- length(zan[[i]])
}
labs <- paste0(rep("Zanthoxyllum", 100), rep(1:100, largo))
zan1 <- unlist(zan)
spp <- names(zan1)
prueba4 <- cbind(labs, spp, zan1)
colnames(prueba4) <- c("Simulacion", "Especie", "Abundancia")

combinado <- rbind(prueba, prueba1, prueba2, prueba3, prueba4)

combinado <- data.frame(Simulacion = combinado[,1], Especie = combinado[,2], Abundancia = combinado[,3])
combinado[,3] <- as.numeric(as.character(combinado[,3]))

prese <- combinado %>% group_by(Simulacion) %>% spread(Especie, Abundancia)

prese[is.na(prese)] <- 0
rownames(prese) <- prese[,1][[1]]

## Warning: Setting row names on a tibble is deprecated.

diste <- vegdist(prese[,-1], "bray")
z1 <- as.matrix(diste)

op1 <- z1[1:100,101:200]#opu-pilo
os1 <- z1[1:100,201:300]#opu-ste
ov1 <- z1[1:100,301:400]#opu-vac
oz1 <- z1[1:100,401:500]#opu-zan
ps1 <- z1[101:200,201:300]#pilo-ste
pv1 <- z1[101:200,301:400]#pilo-vac
pz1 <- z1[101:200,401:500]#pilo-zan
sv1 <- z1[201:300,301:400]#ste-vac
sz1 <- z1[201:300,401:500]#ste-zan
vz1 <- z1[301:400,401:500]#vac-zan

foros <- with(data, unique(Forofito))

op2 <- cbind(rep(foros[4], 100), rep(foros[5], 100), colMeans(op1))
os2 <- cbind(rep(foros[4], 100), rep(foros[2], 100), colMeans(os1))
ov2 <- cbind(rep(foros[4], 100), rep(foros[1], 100), colMeans(ov1))
oz2 <- cbind(rep(foros[4], 100), rep(foros[3], 100), colMeans(oz1))
ps2 <- cbind(rep(foros[5], 100), rep(foros[2], 100), colMeans(ps1))
pv2 <- cbind(rep(foros[5], 100), rep(foros[1], 100), colMeans(pv1))
pz2 <- cbind(rep(foros[5], 100), rep(foros[3], 100), colMeans(pz1))
sv2 <- cbind(rep(foros[2], 100), rep(foros[1], 100), colMeans(sv1))
sz2 <- cbind(rep(foros[2], 100), rep(foros[3], 100), colMeans(sz1))
vz2 <- cbind(rep(foros[1], 100), rep(foros[3], 100), colMeans(vz1))


simu <- rbind(op2, os2, ov2, oz2, ps2, pv2, pz2, sv2, sz2, vz2)
simu <- data.frame(Forofito1 = simu[,1], Forofito2 = simu[,2], BrayCurtis = simu[,3])
simu[,3] <- as.numeric(as.character(simu[,3]))

lines <- c(0.4682081, 0.4597701, 0.4891775, 0.5022624, 0.3297297, 0.4628099, 0.4137931, 0.4320988, 0.5450644, 0.4068966)
obser <- rep(lines, each = 100)
inter2 <- cbind(simu, obser)

perdon por no explicarlo tanto pero la idea era ver si asumiendo neutralidad podria encontrar diferencias en compocion similares a las observadas hay otras formas de hacer la simulacion pero lo decidimos hacer con solo la info de los datos, en este caso las diferencias en la compocion observada son mayores, sin embargo falta hacerlo con datos totalmente simulados con lo hacen en los art generalmente

ggplot(inter2, aes(BrayCurtis)) + geom_histogram() + facet_wrap(Forofito1~Forofito2) + geom_vline(aes(xintercept = obser, col = "red"), linetype = "dashed") + labs(colour = "Observado", x = "Bray-Curtis")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.