#viene copiado de "ER_quemadol_pcor_smoutri.r", en carpeta ER
bichoss=read.table(file="insect.csv",sep=" ",header=T)
sitio.plantas=read.table(file="plants.csv",sep=" ",header=T)
## para cada sp de abejas
source("pcor.r")
mod1=matrix(rep(0,7*18),nrow=18,byrow=T)
dimnames(mod1)=list(c("a.pla","a.abund","a.mean","a.celd","a.nidos","pla.s","abund.mean","abund.celd","abund.nidos","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos","anio.mean","anio.celd","anio.nidos"),c("meg.A","meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig"))
AIC=rep(0,3)
coeff=rep(0,18)


#estiamtion of path coeficients, AIC and C of the model 1 for the 7 bee species
for(n in 1:7){
  sp=c("meg.A", "meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig")
  
  b=subset(bichoss,codigo.i==sp[n])# meg.A meg.cte,mou.tri,ant.dec,ant.rub,xyl.ord,tri.lat,ant.vig
  
  
  b=merge(sitio.plantas,b,by.x="sitio",by.y="sitio",all=T)
  for (i in 1:dim(b)[1]){
    #if (is.na(b$codigo.i[i])==TRUE){b$codigo.i[i]="b"}
    if (is.na(b$mean.no.celd[i])==TRUE){b$mean.no.celd[i]=NA}
    if (is.na(b$celd.tot[i])==TRUE){b$celd.tot[i]=0}
    if (is.na(b$nidos.tot[i])==TRUE){b$nidos.tot[i]=0}
  }
  
  
  # variables
  s=b$S
  a=b$altitud
  mean=b$mean.no.celd
  pla=b$raref
  abund=b$abund.flores
  celd=b$celd.tot
  nidos=b$nidos.tot
  anio=b$anio.q
  #modelo 1, 
  
  var=data.frame(mean,celd,nidos)
  #var=("mean","celd","nidos")
  
  #AIC estimation
  for(m in 1:3){
    v=var[,m]
    data=data.frame(anio,abund,a,s,pla,v)
    
       
    C=-2*(log(cor.test(anio,a, method = c("pearson"))[3]$p.value)
          +log(pcor.test(s,a, data[,c("pla")], method = c("pearson"))[2]$p.value)
          +log(pcor.test(pla,abund, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(s,abund, data[,c("pla","a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(anio,abund, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(pla,anio, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(anio,s, data[,c("pla")], method = c("pearson"), na.rm = T)[2]$p.value))
    #C=12.11525
    #gl=14
    #p=
    #corrección para calcular AICc = C + 2K (n/(n-K-1))
    AICc=C +2*11 *(14/(14-11-1))
    #AICc=166.1152
    AIC[m]=AICc#AIC for diferent species are the same as the response variable is not included in the C formula
  }
  
  
  #path coefficients
  #armar matriz para guardar los coeff. coeff.labels=c("a.pla","a.abund","a.v","pla.s","abund.v","pla.v","s.v","anio.v")
  coeff[1]=cor.test(a,pla, method = c("pearson"), na.rm = T)[4]$estimate[1]
  coeff[2]=cor.test(a,abund, method = c("pearson"), na.rm = T)[4]$estimate[1]
  coeff[3]=pcor.test(a,mean, data[,c("anio","pla","s","abund")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[4]=pcor.test(a,celd, data[,c("anio","pla","s","abund")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[5]=pcor.test(a,nidos, data[,c("anio","pla","s","abund")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[6]=cor.test(pla,s, method = c("pearson"), na.rm = T)[4]$estimate[1]
  coeff[7]=pcor.test(abund,mean, data[,c("anio","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[8]=pcor.test(abund,celd, data[,c("anio","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[9]=pcor.test(abund,nidos, data[,c("anio","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[10]=pcor.test(pla,mean, data[,c("anio","abund","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[11]=pcor.test(pla,celd, data[,c("anio","abund","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[12]=pcor.test(pla,nidos, data[,c("anio","abund","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[13]=pcor.test(s,mean, data[,c("anio","pla","abund","a")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[14]=pcor.test(s,celd, data[,c("anio","pla","abund","a")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[15]=pcor.test(s,nidos, data[,c("anio","pla","abund","a")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[16]=pcor.test(anio,mean, data[,c("abund","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[17]=pcor.test(anio,celd, data[,c("abund","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff[18]=pcor.test(anio,nidos, data[,c("abund","pla","a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  
  
  mod1[,n]=coeff
  
}

# total effects
alt.tot.mean=mod1[3,]+mod1[1,]*mod1[10,]+mod1[1,]*mod1[6,]*mod1[13,]+mod1[2,]*mod1[7,]
alt.tot.celd=mod1[4,]+mod1[1,]*mod1[11,]+mod1[1,]*mod1[6,]*mod1[14,]+mod1[2,]*mod1[8,]
alt.tot.nidos=mod1[5,]+mod1[1,]*mod1[12,]+mod1[1,]*mod1[6,]*mod1[15,]+mod1[2,]*mod1[9,]
div.tot.mean=mod1[10,]+mod1[6,]*mod1[13,]
div.tot.celd=mod1[11,]+mod1[6,]*mod1[14,]
div.tot.nidos=mod1[12,]+mod1[6,]*mod1[15,]

## transformación z ##
mod1.z=(0.5*log((1+mod1)/(1-mod1)))*(14-3)## acá vá la inversa de la varianza = N-3, es para todos igual podría no usarla en este caso
alt.tot.mean.z=(0.5*log((1+alt.tot.mean)/(1-alt.tot.mean)))*(14-3)
alt.tot.celd.z=(0.5*log((1+alt.tot.celd)/(1-alt.tot.celd)))*(14-3)
alt.tot.nidos.z=(0.5*log((1+alt.tot.nidos)/(1-alt.tot.nidos)))*(14-3)
div.tot.mean.z=(0.5*log((1+div.tot.mean)/(1-div.tot.mean)))*(14-3)
div.tot.celd.z=(0.5*log((1+div.tot.celd)/(1-div.tot.celd)))*(14-3)
div.tot.nidos.z=(0.5*log((1+div.tot.nidos)/(1-div.tot.nidos)))*(14-3)


# grafico s/z
names=c("A","B","C")

split.screen(c(2,4))
title(main='a', outer=F)
screen(1)
boxplot(mod1[10,],mod1[11,],mod1[12,],main="Flower diversity",names=names)
abline(0,0,col="red")
arrows(1,-0.1,1,0.47,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.26,2,0.23,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.23,3,0.26,code=3,angle=90,length=0.1,col="blue") 
screen(2)
boxplot(div.tot.mean,div.tot.celd,div.tot.nidos,main="Flower diversity total effect",names=names)
abline(0,0,col="red")
arrows(1,-0.2,1,0.24,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.37,2,0.3,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.35,3,0.36,code=3,angle=90,length=0.1,col="blue") 
screen(3)
boxplot(mod1[13,],mod1[14,],mod1[15,],main="Temporal stability",names=names)
abline(0,0,col="red")
arrows(1,-0.46,1,-0.1,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.27,2,0.18,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.23,3,0.21,code=3,angle=90,length=0.1,col="blue") 
screen(4)
boxplot(mod1[16,],mod1[17,],mod1[18,],main="Fire",names=names)
abline(0,0,col="red")
arrows(1,-0.25,1,0.42,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.26,2,0.31,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.23,3,0.34,code=3,angle=90,length=0.1,col="blue") 
screen(5)
boxplot(mod1[3,],mod1[4,],mod1[5,],main="Elevation",names=names)
abline(0,0,col="red")
arrows(1,-0.54,1,0.16,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.46,2,0.06,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.45,3,0.04,code=3,angle=90,length=0.1,col="blue") 
screen(6)
boxplot(alt.tot.mean,alt.tot.celd,alt.tot.nidos,main="Elevation total effect",names=names)
abline(0,0,col="red")
arrows(1,-0.66,1,0.09,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.45,2,0.02,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.48,3,0.02,code=3,angle=90,length=0.1,col="blue") 
screen(7)
boxplot(mod1[7,],mod1[8,],mod1[9,],main="Flower abundance",names=names)
abline(0,0,col="red")
arrows(1,-0.19,1,0.41,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.26,2,0.11,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.14,3,0.14,code=3,angle=90,length=0.1,col="blue") 

#dev.copy2pdf(file="coeff_quemado.pcor_smoutri.pdf")

# bootstrap (confidence intervals)
eff.tot=rbind(alt.tot.mean.z,alt.tot.celd.z,alt.tot.nidos.z,div.tot.mean.z,div.tot.celd.z,div.tot.nidos.z)
todo=rbind(mod1.z,eff.tot)

xbar.obs <- mean(mean)
xbar.res <- rep(0, 100000)
ci=matrix(rep(0,2*24),nrow=24,byrow=T)
dimnames(ci)=list(c("a.pla","a.abund","a.mean","a.celd","a.nidos","pla.s","abund.mean","abund.celd","abund.nidos","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos","anio.mean","anio.celd","anio.nidos","alt.tot.mean.z","alt.tot.celd.z","alt.tot.nidos.z","div.tot.mean.z","div.tot.celd.z","div.tot.nidos.z"),c("0.025","0.975"))
for (a in 1:24){
  for ( i in 1:100000 )
  {
    x.res <- sample(todo[a,], replace =TRUE)
    xbar.res[i] <- mean(x.res)
  }
  
  range(xbar.res)
  summary(xbar.res)
  xbar.res <- sort(xbar.res)
  ci[a,] <- c(xbar.res[2500], xbar.res[97500])
  
}

#destransformo los datos ci
ci=as.data.frame(ci)
ci$des.0.025=((exp(0.1818182*ci[,1]))-1)/((exp(0.1818182*ci[,1]))+1)
ci$des.0.095=((exp(0.1818182*ci[2]))-1)/((exp(0.1818182*ci[,2]))+1)


##### model 2 ########
mod2=matrix(rep(0,7*11),nrow=11,byrow=T)
dimnames(mod2)=list(c("a.pla","a.mean","a.celd","a.nidos","pla.s","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos"),c("meg.A","meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig"))
AIC2=rep(0,3)
coeff2=rep(0,11)



for(n in 1:7){
  sp=c("meg.A", "meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig")
  
  b=subset(bichoss,codigo.i==sp[n])# meg.A meg.cte,mou.tri,ant.dec,ant.rub,xyl.ord,tri.lat,ant.vig
  
  
  b=merge(sitio.plantas,b,by.x="sitio",by.y="sitio",all=T)
  for (i in 1:dim(b)[1]){
    #if (is.na(b$codigo.i[i])==TRUE){b$codigo.i[i]="b"}
    if (is.na(b$mean.no.celd[i])==TRUE){b$mean.no.celd[i]=NA}
    if (is.na(b$celd.tot[i])==TRUE){b$celd.tot[i]=0}
    if (is.na(b$nidos.tot[i])==TRUE){b$nidos.tot[i]=0}
  }
  
  
  #standardize variables
  s=(b$S-mean(b$S))/sd(b$S)
  a=(b$altitud-mean(b$altitud))/sd(b$altitud)
  mean=(b$mean.no.celd-mean(b$mean.no.celd,na.rm=T))/sd(b$mean.no.celd,na.rm=T)
  pla=(b$raref-mean(b$raref))/sd(b$raref)
  celd=(b$celd.tot-mean(b$celd.tot))/sd(b$celd.tot)
  nidos=(b$nidos.tot-mean(b$nidos.tot))/sd(b$nidos.tot)
  
  
  
  var=data.frame(mean,celd,nidos)
  #var=("mean","celd","nidos")
  
  for(m in 1:3){
    v=var[,m]
    data=data.frame(a,s,pla,v)
    
    #################################################################################################
    
    
    pcor.test(s,a, data[,c("pla")], method = c("pearson"), na.rm = T)
    
    
    
    C=-2*(log(pcor.test(s,a, data[,c("pla")], method = c("pearson"))[2]$p.value))
    #C=3.687032
    #gl=2
    #p=0.84
    #corrección para calcular AICc = C + 2K (n/(n-K-1))
    AICc=C +2*11 *(14/(14-11-1))
    #AICc=
    AIC2[m]=AICc
  }
  
  
  #calculo indices
  #armar matriz para guardar los coeff. coeff.labels=c("a.pla","a.v","pla.s","pla.v","s.v")
  coeff2[1]=cor.test(a,pla, method = c("pearson"), na.rm = T)[4]$estimate[1]
  coeff2[2]=pcor.test(a,mean, data[,c("pla","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[3]=pcor.test(a,celd, data[,c("pla","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[4]=pcor.test(a,nidos, data[,c("pla","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[5]=cor.test(pla,s, method = c("pearson"), na.rm = T)[4]$estimate[1]
  coeff2[6]=pcor.test(pla,mean, data[,c("a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[7]=pcor.test(pla,celd, data[,c("a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[8]=pcor.test(pla,nidos, data[,c("a","s")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[9]=pcor.test(s,mean, data[,c("pla","a")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[10]=pcor.test(s,celd, data[,c("pla","a")], method = c("pearson"), na.rm = T)[1]$estimate
  coeff2[11]=pcor.test(s,nidos, data[,c("pla","a")], method = c("pearson"), na.rm = T)[1]$estimate
  
  
  mod2[,n]=coeff2
  
}


# calculo los efectos totales
alt.tot.mean.2=mod2[2,]+mod2[1,]*mod2[6,]+mod2[1,]*mod2[5,]*mod2[9,]
alt.tot.celd.2=mod2[3,]+mod2[1,]*mod2[7,]+mod2[1,]*mod2[5,]*mod2[10,]
alt.tot.nidos.2=mod2[4,]+mod2[1,]*mod2[8,]+mod2[1,]*mod2[5,]*mod2[11,]
div.tot.mean.2=mod2[6,]+mod2[5,]*mod2[9,]
div.tot.celd.2=mod2[7,]+mod2[5,]*mod2[10,]
div.tot.nidos.2=mod2[8,]+mod2[5,]*mod2[11,]

## transformación z ##
mod2.z=(0.5*log((1+mod2)/(1-mod2)))*(14-3)## acá vá la inversa de la varianza = N-3, es para todos igual podría no usarla en este caso
alt.tot.mean.z2=(0.5*log((1+alt.tot.mean.2)/(1-alt.tot.mean.2)))*(14-3)
alt.tot.celd.z2=(0.5*log((1+alt.tot.celd.2)/(1-alt.tot.celd.2)))*(14-3)
alt.tot.nidos.z2=(0.5*log((1+alt.tot.nidos.2)/(1-alt.tot.nidos.2)))*(14-3)
div.tot.mean.z2=(0.5*log((1+div.tot.mean.2)/(1-div.tot.mean.2)))*(14-3)
div.tot.celd.z2=(0.5*log((1+div.tot.celd.2)/(1-div.tot.celd.2)))*(14-3)
div.tot.nidos.z2=(0.5*log((1+div.tot.nidos.2)/(1-div.tot.nidos.2)))*(14-3)


# grafico s/z
names=c("A","B","C")
split.screen(c(2,3))
title(main='b', outer=F)
screen(1)
boxplot(mod2[6,],mod2[7,],mod2[8,],main="Flower diversity",names=names)
abline(0,0,col="red")
arrows(1,-0.12,1,0.32,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.21,2,0.15,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.21,3,0.17,code=3,angle=90,length=0.1,col="blue") 
screen(2)
boxplot(div.tot.mean.2,div.tot.celd.2,div.tot.nidos.2,main="Flower diversity total effect",names=names)
abline(0,0,col="red")
arrows(1,-0.20,1,0.27,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.30,2,0.19,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.34,3,0.23,code=3,angle=90,length=0.1,col="blue") 
screen(3)
boxplot(mod2[9,],mod2[10,],mod2[11,],main="Temporal stability",names=names)
abline(0,0,col="red")
arrows(1,-0.26,1,-0.004,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.22,2,0.15,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.24,3,0.17,code=3,angle=90,length=0.1,col="blue") 
screen(4)
boxplot(mod2[2,],mod2[3,],mod2[4,],main="Elevation",names=names)
abline(0,0,col="red")
arrows(1,-0.44,1,0.2,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.41,2,0.03,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.41,3,0.01,code=3,angle=90,length=0.1,col="blue") 
screen(5)
boxplot(alt.tot.mean.2,alt.tot.celd.2,alt.tot.nidos.2,main="Elevation total effect",names=names)
abline(0,0,col="red")
arrows(1,-0.41,1,0.13,code=3,angle=90,length=0.1,col="blue")
arrows(2,-0.46,2,0.03,code=3,angle=90,length=0.1,col="blue")
arrows(3,-0.47,3,0.02,code=3,angle=90,length=0.1,col="blue") 

#dev.copy2pdf(file="coeff_quemado.z2.pcor_modsimpl_smoutri.pdf")

# bootstrap
eff.tot2=rbind(alt.tot.mean.z2,alt.tot.celd.z2,alt.tot.nidos.z2,div.tot.mean.z2,div.tot.celd.z2,div.tot.nidos.z2)
todo2=rbind(mod2.z,eff.tot2)
s.mean.ci2.s=matrix(rep(0,2*100),nrow=100,byrow=T)
s.mean.ci2.a=matrix(rep(0,2*100),nrow=100,byrow=T)



xbar.obs <- mean(mean)
xbar.res <- rep(0, 100000)
ci2=matrix(rep(0,2*17),nrow=17,byrow=T)
dimnames(ci2)=list(c("a.pla","a.mean","a.celd","a.nidos","pla.s","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos","alt.tot.mean.z","alt.tot.celd.z","alt.tot.nidos.z","div.tot.mean.z","div.tot.celd.z","div.tot.nidos.z"),c("0.025","0.975"))
for (a in 1:17){
  for ( i in 1:100000 )
  {
    x.res <- sample(todo2[a,], replace =TRUE)
    xbar.res[i] <- mean(x.res)
  }
  
  range(xbar.res)
  summary(xbar.res)
  xbar.res <- sort(xbar.res)
  ci2[a,] <- c(xbar.res[2500], xbar.res[97500])
  
}

xbar.res2=sort(s.mean.ci2.s[,2])
ci2.s.mean.s=c(xbar.res2[2], xbar.res2[95])

xbar.res2=sort(s.mean.ci2.a[,2])
ci2.s.mean.a=c(xbar.res2[2], xbar.res2[95])

#destransformo los datos
ci2=as.data.frame(ci2)
ci2$des.0.025=((exp(0.1818182*ci2[,1]))-1)/((exp(0.1818182*ci2[,1]))+1)
ci2$des.0.095=((exp(0.1818182*ci2[2]))-1)/((exp(0.1818182*ci2[,2]))+1)

#calculo las medianas de los dos modelos
med1=rep(0,dim(mod1)[1])
for (i in 1:dim(mod1)[1]) { 
  med1[i]=median(mod1[i,])
}

med2=rep(0,dim(mod2)[1])
for (i in 1:dim(mod2)[1]) { 
  med2[i]=median(mod2[i,])
}



########## grado #########

dat08=read.table(file="datos.grado.csv",sep=" ",header=TRUE)
## hago bootstrap
riq=data.frame(matrix(NA,nrow=7,ncol=2))
names(riq)=c("codigo.i","r")
riq[1,1]="meg.A"
riq[2,1]="meg.cte"
riq[3,1]="ant.dec"
riq[4,1]="ant.rub"
riq[5,1]="xyl.ord"
riq[6,1]="tri.lat"
riq[7,1]="ant.vig"
riq.m=data.frame(matrix(NA,nrow=7,ncol=100))
s.index=data.frame(matrix(NA,nrow=7,ncol=2))
names(s.index)=c("codigo.i","s.ind")
s.index[1,1]="meg.A"
s.index[2,1]="meg.cte"
s.index[3,1]="ant.vig"
s.index[4,1]="ant.rub"
s.index[5,1]="ant.dec"
s.index[6,1]="xyl.ord"
s.index[7,1]="tri.lat"
s.index.m=data.frame(matrix(NA,nrow=7,ncol=100))

# richness and simpson index for 30 individual of each bee species
for(z in 1:100){
  
  sp=c("meg.A", "meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig")
  a=array(0,2)
  
  # sampling of 30 individuals
  for (i in 1:7){ 
    b=subset(unique(dat08[,c(1,2)]),codigo.i==sp[i])
    ci=sample(b$id.trampa,size=30,replace=T)
    a=c(a,ci)
  }
  a=a[-c(1,2)]## a tiene los id que necesito para calcular los indices, con esto debo hacer un montecarlo
  
  d=subset(dat08,id.trampa==a[1])
  names(d)=c("id.trampa","codigo.i","sitio","codigo.p","sum.sp","sum.tot","prop")    
  
  
  for (m in 1:length(a)){
    e=subset(dat08,id.trampa==a[m])
    d=rbind(d,e)
  }
  d=d[-(dim(subset(dat08,id.trampa==a[1]))[1]),]#acá tengo el subset
  
  # simpson index and richness 
  riqueza=unique(d[,c(2,4)])
  riqueza$i=rep(1,dim(riqueza)[1])                    
  r=aggregate(riqueza$i,by=list(codigo.i=riqueza$codigo.i), sum)# r tiene la riqueza por especie calculada con 30 individuos para cada especie
  riq=merge(riq,r,by.x="codigo.i",by.y="codigo.i") 
  riq=riq[,-2]
  riq.m[z]=riq[,2]
  
  #calculo indice simpson
  simpson.prop=aggregate(d$prop,by=list(codigo.i=d$codigo.i,codigo.p=d$codigo.p),na.rm=T,sum)
  names(simpson.prop)[3]="prop.sum"
  simpson.sum=aggregate(simpson.prop$prop.sum,by= list(codigo.i=simpson.prop$codigo.i),sum)
  names(simpson.sum)[2]="sum"
  simpson.ind=merge(simpson.prop,simpson.sum)
  simpson.ind$prop=simpson.ind$prop.sum/simpson.ind$sum
  simpson.ind$prop2=simpson.ind$prop^2
  ind.simp=aggregate(simpson.ind$prop2,by=list(codigo.i=simpson.ind$codigo.i),sum)
  names(ind.simp)[2]="sum.prop"
  ind.simp$index=1/ind.simp$sum.prop
  simpson.index=ind.simp[,c(1,3)]
  s.index.m[z]=simpson.index[,2]
}
#estimate means of the simpson index and ricness
riq.meg.A=mean(t(riq.m)[,1])
riq.meg.cte=mean(t(riq.m)[,2])
riq.ant.dec=mean(t(riq.m)[,3])
riq.ant.rub=mean(t(riq.m)[,4])
riq.xyl.ord=mean(t(riq.m)[,5])
riq.tri.lat=mean(t(riq.m)[,6])
riq.ant.vig=mean(t(riq.m)[,7])

s.index.meg.A=mean(t(s.index.m)[,1])
s.index.meg.cte=mean(t(s.index.m)[,2])
s.index.ant.dec=mean(t(s.index.m)[,3])
s.index.ant.rub=mean(t(s.index.m)[,4])
s.index.xyl.ord=mean(t(s.index.m)[,5])
s.index.tri.lat=mean(t(s.index.m)[,6])
s.index.ant.vig=mean(t(s.index.m)[,7])
#


#order species as in the rest of the analysis: "ant.dec", "ant.rub" ,"ant.vig" ,"meg.A" ,"meg.cte" ,"tri.lat", "xyl.ord"
s.ind.ord=c(s.index.meg.A,s.index.meg.cte,s.index.ant.dec,s.index.ant.rub,s.index.xyl.ord,s.index.tri.lat,s.index.ant.vig)
ri.ord=c(riq.meg.A,riq.meg.cte,riq.ant.dec,riq.ant.rub,riq.xyl.ord,riq.tri.lat,riq.ant.vig)


###### correlations for table 2

cor.test(ri.ord, mod2[6,],method="spearman")
cor.test(ri.ord, mod2[7,],method="spearman")
cor.test(ri.ord, mod2[8,],method="spearman")
cor.test(s.ind.ord, mod2[6,],method="spearman")
cor.test(s.ind.ord, mod2[7,],method="spearman")
cor.test(s.ind.ord, mod2[8,],method="spearman")


### explore effects of elevation
cor.test(ri.ord,mod2[2,],method="spearman")
cor.test(ri.ord,mod2[3,],method="spearman")
cor.test(ri.ord,mod2[4,],method="spearman")
cor.test(s.ind.ord,mod2[2,],method="spearman")
cor.test(s.ind.ord,mod2[3,],method="spearman")
cor.test(s.ind.ord,mod2[4,],method="spearman")

#####################
#### todos juntos####
#####################
celdillas=subset(bichoss,codigo.i!="lepid8")
celdillas=subset(celdillas,codigo.i!="the.sp2")
cled.sitio=aggregate(celdillas$celd.tot,by=list(sitio=celdillas$sitio),sum)
names(cled.sitio)[2]="no.celd"
nest.sitio=aggregate(celdillas$nidos.tot,by=list(sitio=celdillas$sitio),sum)
names(nest.sitio)[2]="no.nests"

#estiamtion of path coeficients, AIC and C of the model 1 for the pool of species
model1=matrix(rep(0,7*18),nrow=18,byrow=T)
dimnames(mod1)=list(c("a.pla","a.abund","a.mean","a.celd","a.nidos","pla.s","abund.mean","abund.celd","abund.nidos","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos","anio.mean","anio.celd","anio.nidos"),c("meg.A","meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig"))
AIC=rep(0,3)
coeff=rep(0,18)

  
  
  bbb=merge(sitio.plantas,cled.sitio,by.x="sitio",by.y="sitio",all=T)
  bbb=merge(bbb,nest.sitio,by.x="sitio",by.y="sitio",all=T)
  
 
  
  # variables
  s=bbb$S
  a=bbb$altitud
  pla=bbb$raref
  abund=bbb$abund.flores
  celd=bbb$no.celd
  nidos=bbb$no.nests
  anio=bbb$anio.q
  #modelo 1, 
  
  var=data.frame(celd,nidos)

  
  #AIC estimation
  for(m in 1:2){
    v=var[,m]
    data=data.frame(anio,abund,a,s,pla,v)
    
    
    C=-2*(log(cor.test(anio,a, method = c("pearson"))[3]$p.value)
          +log(pcor.test(s,a, data[,c("pla")], method = c("pearson"))[2]$p.value)
          +log(pcor.test(pla,abund, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(s,abund, data[,c("pla","a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(anio,abund, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(pla,anio, data[,c("a")], method = c("pearson"), na.rm = T)[2]$p.value)
          +log(pcor.test(anio,s, data[,c("pla")], method = c("pearson"), na.rm = T)[2]$p.value))
    #C=12.11525
    #gl=14
    #p=
    #corrección para calcular AICc = C + 2K (n/(n-K-1))
    AICc=C +2*11 *(14/(14-11-1))
    #AICc=166.1152
    AIC[m]=AICc#AIC 
  }
  
  
  #path coefficients
  
 cor.test(a,pla, method = c("pearson"), na.rm = T)
 cor.test(a,abund, method = c("pearson"), na.rm = T)
pcor.test(a,celd, data[,c("anio","pla","s","abund")], method = c("pearson"), na.rm = T)
  pcor.test(a,nidos, data[,c("anio","pla","s","abund")], method = c("pearson"), na.rm = T)
  cor.test(pla,s, method = c("pearson"), na.rm = T)
  pcor.test(abund,celd, data[,c("anio","pla","a","s")], method = c("pearson"), na.rm = T)
pcor.test(abund,nidos, data[,c("anio","pla","a","s")], method = c("pearson"), na.rm = T)
 pcor.test(pla,celd, data[,c("anio","abund","a","s")], method = c("pearson"), na.rm = T)
 pcor.test(pla,nidos, data[,c("anio","abund","a","s")], method = c("pearson"), na.rm = T)
  pcor.test(s,celd, data[,c("anio","pla","abund","a")], method = c("pearson"), na.rm = T)
  pcor.test(s,nidos, data[,c("anio","pla","abund","a")], method = c("pearson"), na.rm = T)
  pcor.test(anio,celd, data[,c("abund","pla","a","s")], method = c("pearson"), na.rm = T)
 pcor.test(anio,nidos, data[,c("abund","pla","a","s")], method = c("pearson"), na.rm = T)
  
  


# total effects #####faltan modificarlos!!!!!!!
alt.tot.mean=mod1[3,]+mod1[1,]*mod1[10,]+mod1[1,]*mod1[6,]*mod1[13,]+mod1[2,]*mod1[7,]
alt.tot.celd=mod1[4,]+mod1[1,]*mod1[11,]+mod1[1,]*mod1[6,]*mod1[14,]+mod1[2,]*mod1[8,]
alt.tot.nidos=mod1[5,]+mod1[1,]*mod1[12,]+mod1[1,]*mod1[6,]*mod1[15,]+mod1[2,]*mod1[9,]
div.tot.mean=mod1[10,]+mod1[6,]*mod1[13,]
div.tot.celd=mod1[11,]+mod1[6,]*mod1[14,]
div.tot.nidos=mod1[12,]+mod1[6,]*mod1[15,]
#hasta acá hay que modificar



##### model 2 ########
mod2=matrix(rep(0,7*11),nrow=11,byrow=T)
dimnames(mod2)=list(c("a.pla","a.mean","a.celd","a.nidos","pla.s","pla.mean","pla.celd","pla.nidos","s.mean","s.celd","s.nidos"),c("meg.A","meg.cte","ant.dec","ant.rub","xyl.ord","tri.lat","ant.vig"))
AIC2=rep(0,3)
coeff2=rep(0,11)



  #standardize variables
  s=bbb$S
  a=bbb$altitud
  pla=bbb$raref
  celd=bbb$no.celd
  nidos=bbb$no.nests
  
  
  
  var=data.frame(celd,nidos)
  #var=("mean","celd","nidos")
  
  for(m in 1:2){
    v=var[,m]
    data=data.frame(a,s,pla,v)
    
    #################################################################################################
    
    
    pcor.test(s,a, data[,c("pla")], method = c("pearson"), na.rm = T)
    
    
    
    C=-2*(log(pcor.test(s,a, data[,c("pla")], method = c("pearson"))[2]$p.value))
    #C=3.687032
    #gl=2
    #p=0.84
    #corrección para calcular AICc = C + 2K (n/(n-K-1))
    AICc=C +2*11 *(14/(14-11-1))
    #AICc=
    AIC2[m]=AICc
  }
  
  
  #calculo indices
  #armar matriz para guardar los coeff. coeff.labels=c("a.pla","a.v","pla.s","pla.v","s.v")
  cor.test(a,pla, method = c("pearson"), na.rm = T)
  pcor.test(a,celd, data[,c("pla","s")], method = c("pearson"), na.rm = T)
  pcor.test(a,nidos, data[,c("pla","s")], method = c("pearson"), na.rm = T)
  cor.test(pla,s, method = c("pearson"), na.rm = T)
  pcor.test(pla,celd, data[,c("a","s")], method = c("pearson"), na.rm = T)
  pcor.test(pla,nidos, data[,c("a","s")], method = c("pearson"), na.rm = T)
  pcor.test(s,celd, data[,c("pla","a")], method = c("pearson"), na.rm = T)
  pcor.test(s,nidos, data[,c("pla","a")], method = c("pearson"), na.rm = T)
  
  
  
  
  


######falta modificar#####
# calculo los efectos totales
alt.tot.mean.2=mod2[2,]+mod2[1,]*mod2[6,]+mod2[1,]*mod2[5,]*mod2[9,]
alt.tot.celd.2=mod2[3,]+mod2[1,]*mod2[7,]+mod2[1,]*mod2[5,]*mod2[10,]
alt.tot.nidos.2=mod2[4,]+mod2[1,]*mod2[8,]+mod2[1,]*mod2[5,]*mod2[11,]
div.tot.mean.2=mod2[6,]+mod2[5,]*mod2[9,]
div.tot.celd.2=mod2[7,]+mod2[5,]*mod2[10,]
div.tot.nidos.2=mod2[8,]+mod2[5,]*mod2[11,]
#####hasta acá######