Estrutura da variação morfológica no espaço

Se simpatria entre grupos morfológicos pode ser usada como evidência indireta de isolamento reprodutivo, ainda que parcial, se numa análise da distribuição espacial da variação morfológica houver indicação de locais onde a variacao observada é maior do que o esperado ao acaso, isso iria sugerir que no grupo há mais de uma espécie. Inversamente, se a variação morfológica observa for menor que o esperado ao acaso, não haveria indicação de que diferentes morfologias coexistem.
O seguinte script testa se há locais na distribuição geográfica de um grupo de interesse onde a variação morfológica seja maior ou menor que o esperado ao acaso, aleatorizando os próprios dados para gerar uma distribuição de valores esperados (aleatórios) com a qual o valor observado pode ser comparado.

#limpa os objetos
rm(list=ls())
#habilita pacotes necessarios
library(maps)
library(cluster)
 
#le os dados 
#dados com latitude  e longitude das amostras
specs = read.table("formularioMestradoFINAL.csv",sep="\t",as.is=T,header=T,na.strings = c("NA",""))
cls = c("WikiEspecimenID", "NOME",  "LATITUDE", "LONGITUDE")
vl = !is.na(specs$WikiEspecimenID)
specs = specs[vl,cls]
specs = unique(specs)
rownames(specs) = specs$WikiEspecimenID
 
#dados morfologicos das amostras
morfo = read.table("dados_medias_varsel.csv",sep='\t',as.is=T,row.names = 1,header=T, na.strings = c("NA",""))
#calcula a distancia morfologica entre as amostras
colnames(morfo)  #deve ser apenas as variaveis
morfodist = as.matrix(daisy(morfo, metric = "gower", stand = F,warnBin = F))
#nome das colunas e linhas dessa matrix deve ser o mesmo identificaor de specs
 
#filtra os dados da especie a ser testada e que tem dados geograficos
spps = unique(specs$NOME)
#cria o filtro
vl = specs$NOME=="Ocotea longifolia" & !is.na(specs$NOME) & !is.na(specs$LONGITUDE) & !is.na(specs$LATITUDE)
#filtra os dados
data.df = specs[vl,]
 
#calcula o numero de amostras em cada celula de 1x1 graus decimais
ji <- function(xy, origin=c(0,0), cellsize=c(1,1)) {
  t(apply(xy, 1, function(z) cellsize/2+origin+cellsize*(floor((z - origin)/cellsize))))
}
JI <- ji(cbind(data.df$LONGITUDE, data.df$LATITUDE))
data.df$X <- JI[, 1]
data.df$Y <- JI[, 2]
data.df$Cell <- paste(data.df$X, data.df$Y)
#conta o numero de amostras em cada celula
counts <- by(data.df, data.df$Cell, function(d) c(d$X[1], d$Y[1], nrow(d)))
counts.m <- matrix(unlist(counts), nrow=3)
rownames(counts.m) <- c("X", "Y", "Count")
 
#plota um mapa mostrando essa contagem
count.max <- max(counts.m["Count",])
xl = range(data.df$LONGITUDE)+c(-0.5,0.5)
yl = range(data.df$LATITUDE)+c(-0.5,0.5)
pdf(file='mapaDensidadeColetaOLongifolia.pdf',height = 5,width=6)
par(mgp=c(1.2,0.5,0),cex.axis=0.8,tck=-0.01,mar=c(1,2,1,2))
#plota o mapa de todos os pontos
map(xlim=xl,yl=yl)
#tamanho dos pontos proporcional ao numero de amostras
ccex = sqrt(counts.m["Count",]*10/100)
points(counts.m["X",] + 1/2, counts.m["Y",] + 1/2, cex=ccex,pch = 19, col="black")
box()
axis(side=1)
axis(side=2,las=2)
axis(side=4,las=2)
axis(side=3)
 
#quais pontos tem mais do que 3 amostras?
td = t(counts.m)
vl = td[,3]>=3
#colore esses pontos no mapa
points(counts.m["X",vl] + 1/2, counts.m["Y",vl] + 1/2, cex=ccex[vl],pch = 19, col="red")
dev.off()
 
#para cada um desse locais, testa quem tem maior ou menor variacao que o esperado
cells = paste(td[vl,1],td[vl,2])
ospvs = NULL
for(l in 1:length(cells)) {
  print(paste("testando a celula", cells[l]))
  #pega as plantas que estao na na celula
  vl2 = data.df$Cell%in%cells[l]
  ids = as.character(data.df$WikiEspecimenID[vl2])
  #distancias morfologicas entre essas plantas
  mf = morfodist[ids,ids]
  #distancia media morfologica OBSERVADA na celula
  obsdist = mean(mf[lower.tri(mf)])
 
  #gera uma distribuicao aleatoria de distancias morfologicas (ESPERADO)
  asdist = NULL
  for(i in 1:1000) {
    #aleatoriza a matriz morfologica
    nn = sample(colnames(morfodist))
    omd = morfodist
    colnames(omd) = nn
    rownames(omd) = nn
    mf = omd[ids,ids]
    esdist = mean(mf[lower.tri(mf)],na.rm=T)
    asdist= c(asdist,esdist)
  }
  #o valor observado e significativamente maior que o esperado?
  #ou seja quantos valores esperados sao maiores que o observado?
  pv = sum(asdist>obsdist,na.rm=T)/length(asdist)
  if (pv<=0.05) {
   print(paste("celula", cells[l], "tem mais  variacao morfologica do que o esperado"))
  }
  ospvs = c(ospvs,pv)
}
res = data.frame(CELULA=cells,PV.MAIOR.ESPERADO=ospvs,PV.MENOR.ESPERADO=(1-ospvs))
#quantos pontos foram testados
nrow(res)
 
#quantos tem mais variacao que o esperado?
sum(res$PV.MAIOR.ESPERADO<0.05)
 
#quantos tem menos variacao que o esperado?
sum(res$PV.MENOR.ESPERADO<0.05)
 
#salva o resultado
write.table(res,file='resultadoTesteVariacaoEspacia.csv',sep='\t',na='',row.names = F)