====== Códigos no R ======

===== IPMs, filogenias e sinal filogenético =====


Este é o código para geração dos IPMs, para as duas filogenias e para o teste de Mantel.
Aqui também encontram-se os arquivos necessários para gerar o output do código,  e a área de trabalho após rodar o código (um pouco demorado).

Os pacotes necessários são:
  * ade4
  * ape
  * IPMpack

{{peicdez2012_modif.txt|Planilha de dados PEIC}}

{{makesurvobj1.r|função para criação dos objetos de sobrevivencia modificada}}

{{harmon_glorpp.r|Modelo nulo para Mantel - Harmon & Glor}}

{{mantelr.txt|Teste de Mantel com o modelo nulo acima incorporado}}

{{corte5.r|Código abaixo}}

{{corte5.zip|Área de trabalho gerada a partir do código abaixo - Parte 1}}

{{corte5.z01.pdf|Área de trabalho gerada a partir do código abaixo - Parte 2 - retirar ".pdf"}}

{{corte5.z02.pdf|Área de trabalho gerada a partir do código abaixo - Parte 3 - retirar ".pdf"}}

{{corte5.z03.pdf|Área de trabalho gerada a partir do código abaixo - Parte 4 - retirar ".pdf"}}


<code>
setwd("C:/Users/Gabriel/Documents/Mestrado/Analises/Cortes")

peic=read.table("peicDez2012_modif.txt", header=T, sep="\t", as.is=TRUE)

#setwd("C:/Users/Gabriel/Documents/Mestrado/Analises/IPM")

#lista de espécies com recrutas, para poder ter estimativa de fecundidade
recrut.sp=sort(unique(peic[peic$status_05=="R", "specie"]))
adult.censo1=sort(unique(peic[peic$status_05=="A" & peic$dbh_04>50, "specie"]))
listasp=sort(recrut.sp[which(recrut.sp %in% adult.censo1)])#108 especies


#retirando as especies q nao estao nas filogenias
listasp=listasp[-which(listasp=="Annona neosericea")]
listasp=listasp[-which(listasp=="Malouetia cestroides")]
listasp=listasp[-which(listasp=="Bactris setosa")]
listasp=listasp[-which(listasp=="Croton sphaerogynus")]
listasp=listasp[-which(listasp=="Nectandra oppositifolia")]
listasp=listasp[-which(listasp=="Persea willdenovii")]
listasp=listasp[-which(listasp=="Tibouchina trichopoda")]
listasp=listasp[-which(listasp=="Cabralea canjerana")]
listasp=listasp[-which(listasp=="Ficus eximia")]
listasp=listasp[-which(listasp=="Myrcia glabra")]
listasp=listasp[-which(listasp=="Myrceugenia myrcioides")]
listasp=listasp[-which(listasp=="Ocotea dispersa")]
listasp=listasp[-which(listasp=="Psychotria mapourioides")]
listasp=listasp[-which(listasp=="Ecclinusa ramiflora")]
listasp=listasp[-which(listasp=="Symplocos laxiflora")]


#estimando a fecundidade

fecundity=data.frame(listasp)
for(i in 1:length(listasp))
{
fecundity[i,"N offspring50"]=length(peic[peic$status_05=="R" & peic$specie==listasp[i] & peic$dbh_09>50, "arv"])
fecundity[i,"A"]=length(peic[peic$status_05=="A" & peic$specie==listasp[i] & peic$dbh_04>50, "arv"])
fecundity[i,"AS"]=length(peic[peic$status_05=="AS" & peic$specie==listasp[i] & peic$dbh_09>50, "arv"])
}

fecundity[,"N adults"]=fecundity[,"A"]+fecundity[,"AS"]
fecundity[,"fec50"]=fecundity[,"N offspring50"]/fecundity[,"N adults"]


#selecionando colunas que interessam
peic.select=peic[peic$specie %in% listasp, c("specie","dbh_04", "dbh_09", "status_05", "status09")]

#criando a coluna survival
peic.select[,"surv"]=NA
peic.select[peic.select$status_05!="R" & peic.select$status09=="A","surv"]=1
peic.select[peic.select$status_05!="R" & peic.select$status09=="AS","surv"]=1
peic.select[peic.select$status09=="D","surv"]=0
peic.select=peic.select[peic.select$status09!="E" & peic.select$status09!="NE",]

###MODELO DE PROJECAO INTEGRAL
#criando as matrizes IPM
library(IPMpack)

IPMdata=peic.select[,c("specie","dbh_04","dbh_09","surv")]
names(IPMdata)=c("specie","size","sizeNext","surv")
IPMdata[is.na(IPMdata$size),"surv"]=NA
IPMdata[,"fec"]=NA
IPMdata[,"stage"]="continuous"
IPMdata[is.na(IPMdata$size),"stage"]=NA
IPMdata[,"stageNext"]="continuous"
IPMdata[is.na(IPMdata$sizeNext),"stageNext"]=NA
head(IPMdata)
IPMdata[,"size"]=as.numeric(IPMdata[,"size"])
IPMdata[,"stage"]=as.factor(IPMdata[,"stage"])
IPMdata[,"stageNext"]=as.factor(IPMdata[,"stageNext"])


#checando o numero de individuos para cada especie
numero=data.frame(specie=listasp)
for(i in 1:length(listasp))
{
numero[i,"n"]=nrow(IPMdata[IPMdata$specie==listasp[i] & !is.na(IPMdata$size) & !is.na(IPMdata$sizeNext),c("size","sizeNext")])
}

#espécies com n<5
excluir5=numero[numero$n<5,"specie"]#especies que vão sair de listasp
listasp5=listasp[-which(listasp %in% excluir5)]

#espécies com n<20
#excluir20=numero[numero$n<20,"specie"]#especies que vão sair de listasp
#listasp20=listasp[-which(listasp %in% excluir20)]

#espécies com n<50
#excluir50=numero[numero$n<50,"specie"]#especies que vão sair de listasp
#listasp50=listasp[-which(listasp %in% excluir50)]

#espécies com n<100
#excluir100=numero[numero$n<100,"specie"]#especies que vão sair de listasp
#listasp100=listasp[-which(listasp %in% excluir100)]


################################################################
################################################################

listasp=listasp5


#Crescimento
gr=list()
 for(i in 1:length(listasp))
{
gr[[i]] = makeGrowthObj(IPMdata[IPMdata$specie==listasp[i],])
}

#Sobrevivencia
#precisa aumentar as iteracoes do glm para os modelos convergirem
#portanto, modifiquei a funçao makeSurvObj:

source("makeSurvObj1.r")

sv=list() 
for(i in 1:length(listasp))
{
sv[[i]] = makeSurvObj1(IPMdata[IPMdata$specie==listasp[i],], explanatoryVariables = "size+size2")
}

##Create P Matrices
#definindo os tamanhos minimos e maximos das matrizes P
tamanho=data.frame(specie=listasp)

for(i in 1:length(listasp))
{
tamanho[i,"min"]=min(IPMdata[IPMdata$specie==listasp[i],"sizeNext"], na.rm=TRUE)
tamanho[i,"max"]=max(IPMdata[IPMdata$specie==listasp[i],"sizeNext"], na.rm=TRUE)
}

Pmatrix=list()#uma lista de matrizes com os mesmos parametros para todas as especies
for(i in 1:length(listasp))
{
Pmatrix[[i]] <- createIPMPmatrix(nBigMatrix = 100,minSize = tamanho[i,"min"]*0.9 , maxSize=tamanho[i,"max"]*1.1, growObj = gr[[i]], survObj = sv[[i]],correction = "constant")
}

#Fecundidade

IPMdata[, "fec"]=0
for(i in 1:length(listasp))
{
IPMdata[which(IPMdata$size>=50 & IPMdata$surv==1 & IPMdata$specie==listasp[i]), "fec"]=fecundity[i, "fec50"]
}
head(IPMdata)

fv50=list()
for(i in 1:length(listasp))
{
fv50[[i]] = makeFecObj(IPMdata[IPMdata$specie==listasp[i],])
}

Fmatrix50=list()
for(i in 1:length(listasp))
{
Fmatrix50[[i]] <- createIPMFmatrix(nBigMatrix = 100, tamanho[i,"min"]*0.9 ,maxSize = tamanho[i,"max"]*1.1 ,fecObj = fv50[[i]])
}

IPM50=list()
for(i in 1:length(listasp))
{
IPM50[[i]]=Pmatrix[[i]]+Fmatrix50[[i]]
}

lambdas=data.frame(specie=listasp)
for(i in 1:length(listasp))
{
lambdas[i,"lambda.fec50"]=eigen(IPM50[[i]])$value[1]
}

#sobrevivencia em cada classe s (sigma)
s=list()
for(p in 1:length(listasp))
{
s[[p]]=numeric(0)
for(j in 1:ncol(Pmatrix[[p]]))
{
s[[p]][j]=sum(Pmatrix[[p]][,j])
}
}

#crescimento positivo (da classe j para classe i) e crescimento negativo (classe i para classe j)
gr=list()
for(p in 1:length(listasp))
{
gr[[p]]=Pmatrix[[p]]/s[[p]]
}

#crescimento positivo apenas (gama ij)
gij=list()
for(p in 1:length(listasp))
{
gij[[p]]=gr[[p]]
gij[[p]][upper.tri(gij[[p]])]=0
diag(gij[[p]])=0
}

#crescimento negativo apenas (ro ij)
rij=list()
for(p in 1:length(listasp))
{
rij[[p]]=gr[[p]]
rij[[p]][lower.tri(rij[[p]])]=0
diag(rij[[p]])=0
}

#fecundidade individual de j para i (fi ij)

fij50=list()
for(p in 1:length(listasp))
{
fij50[[p]]=Fmatrix50[[p]]/s[[p]]
}

#Survival elasticity

Es50=list()
for(p in 1:length(listasp))
{
Es50[[p]]=(s[[p]]/lambdas[p,"lambda.fec50"])*(diag(sens(IPM50[[p]]))*(1-(apply((gij[[p]]), MARGIN=1, FUN=sum))-(apply((rij[[p]]), MARGIN=1, FUN=sum))))+apply(sens(IPM50[[p]])*gij[[p]], MARGIN=1, FUN=sum)+apply(sens(IPM50[[p]])*rij[[p]], MARGIN=1, FUN=sum)+apply(sens(IPM50[[p]])*fij50[[p]], MARGIN=1, FUN=sum)
}

#Positive growth elasticity

Eg50=list()
for(p in 1:length(listasp))
{
Eg50[[p]]=numeric(0)
Eg50[[p]]=(gij[[p]]/lambdas[p,"lambda.fec50"])*((diag(sens(IPM50[[p]]))*(-s[[p]]))+ ( sens(IPM50[[p]])* s[[p]] ))
}


#negative growth elasticities

Er50=list()
for(p in 1:length(listasp))
{
Er50[[p]]=numeric(0)
Er50[[p]]=(rij[[p]]/lambdas[p,"lambda.fec50"])*((diag(sens(IPM50[[p]]))*(-s[[p]]))+ ( sens(IPM50[[p]])* s[[p]] ))
}

#fecundity elasticities

Ef50=list()
for(p in 1:length(listasp))
{
Ef50[[p]]=numeric(0)
Ef50[[p]]=(fij50[[p]]/lambdas[p,"lambda.fec50"])*( sens(IPM50[[p]])* s[[p]] )
}

#Elasticidades das taxas vitais

E50=numeric(0)
S50=numeric(0)
G50=numeric(0)
F50=numeric(0)
for(p in 1:length(listasp))
{
E50[p]=sum(abs(Es50[[p]]))+sum(abs(Eg50[[p]]))+sum(abs(Er50[[p]]))+sum(abs(Ef50[[p]]))
S50[p]=sum(abs(Es50[[p]]))/E50[p]
G50[p]=(sum(abs(Eg50[[p]]))+sum(abs(Er50[[p]])))/E50[p]
F50[p]=sum(abs(Ef50[[p]]))/E50[p]
}

S50=as.numeric(S50)
G50=as.numeric(G50)
F50=as.numeric(F50)

#distancias ecologicas
dist.f50=as.matrix(dist(data.frame(F50,S50,G50)))
colnames(dist.f50)=listasp
rownames(dist.f50)=listasp



###############
## FILOGENIA ##
###############
library(ape)

##Arvore Wikstrom
WIK="((((((((Rubiaceae_Cordiera_myrciifolia:71,Rubiaceae_Amaioua_intermedia:71,Rubiaceae_Posoqueria_latifolia:71)Rubiaceae:36,((Bignoniaceae_Jacaranda_puberula:38,(Bignoniaceae_Tabebuia_alba:38,Bignoniaceae_Tabebuia_cassinoides:38)Tabebuia:0)Bignoniaceae:25,Oleaceae_Chionanthus_filiformis:63)Lamiales:44)EuasteridsI:5,(Araliaceae_Schefflera_angustissima:107,(Aquifoliaceae_Ilex_dumosa:49,Aquifoliaceae_Ilex_pseudobuxus:49,Aquifoliaceae_Ilex_theezans:49)Aquifoliaceae:58)EuasteridsII:5)25:5,((Primulaceae_Cybianthus_brasiliensis:15,(Primulaceae_Myrsine_guianensis:15,Primulaceae_Myrsine_umbellata:15,Primulaceae_Myrsine_venosa:15)Rapanea:0)Primulaceae:99,((Styracaceae_Styrax_glabratus:67,((Sapotaceae_Manilkara_subsericea:54,Sapotaceae_Pouteria_beaurepairei:54)Sapotaceae:0,(Theaceae_Laplacea_fructicosa:54,Pentaphylacaceae_Ternstroemia_brasiliensis:54)Theaceae:0)23:13)22:7,Clethraceae_Clethra_scabra:74)18:40)Ericales:3)Asterids:5,(Nyctaginaceae_Guapira_opposita:113,Olacaceae_Heisteria_silvianii:113)15:9)14:33,((Thymelaeaceae_Daphnopsis_racemosa:107,(Anacardiaceae_Tapirira_guianensis:61,(Sapindaceae_Matayba_guianensis:56,Meliaceae_Guarea_macrophylla:56)37:5)Sapindales:46)35:14,(((Melastomataceae_Miconia_chartacea:41,Melastomataceae_Miconia_cubatanensis:41,Melastomataceae_Miconia_saldanhaei:41)Melastomataceae:47,(Myrtaceae_Blepharocalyx_salicifolius:68,Myrtaceae_Calyptranthes_concinna:68,(Myrtaceae_Eugenia_neoglomerata:68,Myrtaceae_Eugenia_stigmatosa:68,Myrtaceae_Eugenia_sulcata:68,Myrtaceae_Eugenia_umbelliflora:68)Eugenia:0,(Myrtaceae_Myrcia_hebepetala:68,Myrtaceae_Myrcia_ilheosensis:68,Myrtaceae_Myrcia_brasiliensis:68)Gomidesia:0,Myrtaceae_Marlierea_racemosa:68,(Myrtaceae_Myrcia_pulchra:68,Myrtaceae_Myrcia_pubipetala:68,Myrtaceae_Myrcia_insularis:68,Myrtaceae_Myrcia_multiflora:68,Myrtaceae_Myrcia_racemosa:68,Myrtaceae_Myrcia_splendens:68,Myrtaceae_Myrcia_sp.1:68)Myrcia:0,Myrtaceae_Neomitranthes_glomerata:68,Myrtaceae_Pimenta_pseudocaryophyllus:68,Myrtaceae_Psidium_cattleyanum:68,Myrtaceae_Siphoneugena_guilfoyleiana:68)Myrtaceae:20)Myrtales:26,(((Urticaceae_Cecropia_pachystachya:76,Urticaceae_Coussapoa_microcarpa:76)Urticaceae:3,((Fabaceae_Abarema_brachystachya:56,Fabaceae_Abarema_langsdorffii:56)Abarema:0,Fabaceae_Andira_anthelmia:56,Fabaceae_Albizia_pedicellaris:56,Fabaceae_Hymenolobium_janeirense:56,Fabaceae_Ormosia_arborea:56)Fabaceae:23)Fabales:22,(Celastraceae_Maytenus_robusta:94,(((Euphorbiaceae_Alchornea_triplinervia:69,Euphorbiaceae_Aparisthmium_cordatum:69,Euphorbiaceae_Maprounea_guianensis:69,Peraceae_Pera_glabrata:69)Euphorbiaceae:12,(Chrysobalanaceae_Hirtella_hebeclada:78,(((Clusiaceae_Calophyllum_brasiliense:45,Clusiaceae_Clusia_criuva:45,Clusiaceae_Garcinia_gardneriana:45)Clusiaceae:23,Malpighiaceae_Byrsonima_ligustrifolia:68)52:9,(Humiriaceae_Humiriastrum_dentatum:76,(Ochnaceae_Ouratea_parviflora:74,Erythroxylaceae_Erythroxylum_amplifolium:74)55:2)53:1)51:1)50:3)Malpighiales:10,(Cunoniaceae_Weinmannia_paulliniifolia:77,Elaeocarpaceae_Sloanea_guianensis:77)Oxalidales:14)47:3)46:7)EurosidsI:13)17:7)Rosids:34)Eudicots:24,(((((Lauraceae_Aiouea_saligna:34,Lauraceae_Aniba_viridis:34,Lauraceae_Endlicheria_paniculata:34,Lauraceae_Nectandra_grandiflora:34,(Lauraceae_Ocotea_aciphylla:34,Lauraceae_Ocotea_dispersa:34,Lauraceae_Ocotea_pulchella:34,Lauraceae_Ocotea_pulchra:34,Lauraceae_Ocotea_venulosa:34)Ocotea:0)Lauraceae:51,(Monimiaceae_Mollinedia_boracensis:71,Monimiaceae_Mollinedia_schottiana:71)Monimiaceae:14)Laurales:57,(Annonaceae_Guatteria_australis:26,(Annonaceae_Xylopia_brasiliensis:26,Annonaceae_Xylopia_langsdorffiana:26)Xylopia:0)Annonaceae:116)5:12,(Arecaceae_Astrocaryum_aculeatissimum:73,Arecaceae_Euterpe_edulis:73,Arecaceae_Geonoma_schottiana:73,Arecaceae_Syagrus_romanzoffiana:73)Arecaceae:81)3:4,Chloranthaceae_Hedyosmum_brasiliense:158)Magnoliids:21)Angiosperms:189,Podocarpaceae_Podocarpus_sellowii:368)Seed Plants;"
wiktree= read.tree(text=WIK)


#retirando da arvore especies que nao possuem demografia
listasp1=character(0)
for(i in 1:length(listasp))
{
listasp1[i]=paste(sort(unique(peic[peic$specie==listasp[i], "Familia"])),"_",sort(unique(peic[peic$specie==listasp[i], "genero"]))[1],"_",sort(unique(peic[peic$specie==listasp[i], "epitsp"]))[1], sep="")
}

listaspA=listasp1
remvA=wiktree$tip.label[-which(wiktree$tip.label %in% listaspA)]
wiktree=drop.tip(wiktree,remvA)

##Arvore molecular
load("C:\\Users\\Gabriel\\Documents\\Mestrado\\Analises\\Filogenia\\Alexandre\\newtrees.Rdata")

#adequar os nomes
tree_chrono$tip.label[tree_chrono$tip.label=="Aquifoliaceae_Ilex_amara"]="Aquifoliaceae_Ilex_dumosa"
tree_chrono$tip.label[tree_chrono$tip.label=="Cecropiaceae_Cecropia_pachystachya"]="Urticaceae_Cecropia_pachystachya"
tree_chrono$tip.label[tree_chrono$tip.label=="Cecropiaceae_Coussapoa_microcarpa"]="Urticaceae_Coussapoa_microcarpa"
tree_chrono$tip.label[tree_chrono$tip.label=="Euphorbiaceae_Pera_glabrata_NB11"]="Peraceae_Pera_glabrata_NB11"
tree_chrono$tip.label[tree_chrono$tip.label=="Fabaceae_Balizia_pedicelaris"]="Fabaceae_Albizia_pedicellaris"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrsinaceae_Cybianthus_peruvianus"]="Primulaceae_Cybianthus_brasiliensis"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrsinaceae_Rapanea_ferruginea"]="Primulaceae_Myrsine_ferruginea"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrsinaceae_Rapanea_guianensis"]="Primulaceae_Myrsine_guianensis"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrsinaceae_Rapanea_umbellata"]="Primulaceae_Myrsine_umbellata"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrsinaceae_Rapanea_venosa"]="Primulaceae_Myrsine_venosa"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrtaceae_Gomidesia_affinis"]="Myrtaceae_Myrcia_hebepetala"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrtaceae_Gomidesia_fenzliana"]="Myrtaceae_Myrcia_ilheosensis"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrtaceae_Gomidesia_schaueriana"]="Myrtaceae_Myrcia_brasiliensis"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrtaceae_Myrcia_rostrata"]="Myrtaceae_Myrcia_splendens"
tree_chrono$tip.label[tree_chrono$tip.label=="Myrtaceae_Myrcia_sp"]="Myrtaceae_Myrcia_sp.1"
tree_chrono$tip.label[tree_chrono$tip.label=="Podocarpus_sellowii"]="Podocarpaceae_Podocarpus_sellowii"
tree_chrono$tip.label[tree_chrono$tip.label=="Rubiaceae_Alibertia_myrciifolia"]="Rubiaceae_Cordiera_myrciifolia"
tree_chrono$tip.label[tree_chrono$tip.label=="Tabebuia_cassinoides_rbcL"]="Bignoniaceae_Tabebuia_cassinoides"
tree_chrono$tip.label[tree_chrono$tip.label=="Theaceae_Gordonia_fruticosa"]="Theaceae_Laplacea_fructicosa"
tree_chrono$tip.label[tree_chrono$tip.label=="Theaceae_Ternstroemia_brasiliensis"]="Pentaphylacaceae_Ternstroemia_brasiliensis"
tree_chrono$tip.label[tree_chrono$tip.label=="Thymeleaceae_Daphnopsis_racemosa"]="Thymelaeaceae_Daphnopsis_racemosa"

listaspB=listasp1
listaspB[listaspB=="Arecaceae_Euterpe_edulis"]="Arecaceae_Euterpe_edulis_NB1"
listaspB[listaspB=="Chloranthaceae_Hedyosmum_brasiliense"]="Chloranthaceae_Hedyosmum_brasiliense_NB3"
listaspB[listaspB=="Clethraceae_Clethra_scabra"]="Clethraceae_Clethra_scabra_NB4"
listaspB[listaspB=="Clusiaceae_Calophyllum_brasiliense"]="Clusiaceae_Calophyllum_brasiliense_NB5"
listaspB[listaspB=="Clusiaceae_Clusia_criuva"]="Clusiaceae_Clusia_criuva_NB21"
listaspB[listaspB=="Clusiaceae_Garcinia_gardneriana"]="Clusiaceae_Garcinia_gardneriana_NB6"
listaspB[listaspB=="Erythroxylaceae_Erythroxylum_amplifolium"]="Erythroxylaceae_Erythroxylum_amplifolium_NB7"
listaspB[listaspB=="Euphorbiaceae_Alchornea_triplinervia"]="Euphorbiaceae_Alchornea_triplinervia_NB8"
listaspB[listaspB=="Euphorbiaceae_Maprounea_guianensis"]="Euphorbiaceae_Maprounea_guianensis_NB10"
listaspB[listaspB=="Peraceae_Pera_glabrata"]="Peraceae_Pera_glabrata_NB11"
listaspB[listaspB=="Fabaceae_Andira_anthelmia"]="Fabaceae_Andira_anthelmia_NB12"
listaspB[listaspB=="Meliaceae_Guarea_macrophylla"]="Meliaceae_Guarea_macrophylla_NB14"
listaspB[listaspB=="Ochnaceae_Ouratea_parviflora"]="Ochnaceae_Ouratea_parviflora_NB16"
listaspB[listaspB=="Sapindaceae_Matayba_guianensis"]="Sapindaceae_Matayba_guianensis_NB17"
listaspB[listaspB=="Sapotaceae_Pouteria_beaurepairei"]="Sapotaceae_Pouteria_beaurepairei_NB23"


#retirar da arvore especies que nao ha demografia

remvB=tree_chrono$tip.label[-which(tree_chrono$tip.label %in% listaspB)]
tree_chrono=drop.tip(tree_chrono, remvB)


##########################
## 3 SINAL FILOGENETICO ##
##########################

##Mantel

# árvore Wikstrom
dist.wik=cophenetic(wiktree)

#árvore molecular
dist.molec=cophenetic(chronoMPL(tree_chrono))

###distancia ecológica
#a)reordenando as especies na ordem da arvore wikstrom
elasticidades=data.frame(listasp1,F50,S50,G50)
elasWIK=elasticidades[match(wiktree$tip.label,elasticidades$listasp1),]

#b)reordenando as especies na ordem da arvore molecular
elas2=elasticidades
elas2$listasp1=listaspB
elasMOL=elas2[match(tree_chrono$tip.label,elas2$listasp1),]

#Distancia ecologica para arvore WIK
dist.ecolA=as.matrix(dist(elasWIK[,c(2,3,4)]))
colnames(dist.ecolA)=wiktree$tip.label
rownames(dist.ecolA)=wiktree$tip.label

#Distancia ecologica para arvore molecular
dist.ecolB=as.matrix(dist(elasMOL[,c(2,3,4)]))
colnames(dist.ecolB)=tree_chrono$tip.label
rownames(dist.ecolB)=tree_chrono$tip.label

#MANTEL Lapointe and Garland 2001
library(ade4)

source("harmon&glorPP.r")
source("mantelR.txt")

mantelWIK=phyloMantelr(phy=wiktree,m1=dist.ecolA,m2=dist.wik,nperm=10000,graph=TRUE)
mantelMOLEC=phyloMantelr(phy=tree_chrono,m1=dist.ecolB,m2=dist.molec,nperm=10000,graph=TRUE,rpt=F)



################
###resultados###
################

#elast
setwd("C:/Users/Gabriel/Documents/Mestrado/Dissertação/CORTE5")
elast=elasticidades
elast[,2]=round(elast[,2],4)
elast[,3]=round(elast[,3],4)
elast[,4]=round(elast[,4],4)
elast$listasp1=as.character(elast$listasp1)

#lambdas
lamb.round=lambdas
lamb.round[,2]=as.numeric(round(lamb.round[,2],3))
lamb.round$specie=as.character(lamb.round$specie)
for(i in 1:length(listasp))
{
lamb.round$specie[i]=paste(sort(unique(peic[peic$specie==listasp[i], "Familia"])),"_",sort(unique(peic[peic$specie==listasp[i], "genero"]))[1],"_",sort(unique(peic[peic$specie==listasp[i], "epitsp"]))[1], sep="")
}


</code>