Calculando as transições
Este script e ainda experimental para o calculo das transicoes.
### Preparando os dados o cálculo das transições ###
#Estes procedimentos foram preparados para uma planilha típica de dados de demografia de árvores, onde as linhas são os indivíduos e as colunas as medidas tomadas, sendo que repetições temporais da mesma medida representam colunas diferentes.
#O indivíduo pode ter uma das seguintes medidas: DAP(diâmetro à altura do peito), para plantas com DAP até 70mm (com paquímetro), PAP (perímetro à altura do peito - fita métrica) ou DAS (diâmetro à altura do solo), quando à altura do peito ele apresenta DAP<10mm ou tem altura inferior a 1,3m.
#As medidas em cada tempo estão identificadas por um número: 0, 1 e 2
#Preparando a planilha a ser salva como txt:
#1. Espaços em branco (principalmente nos outros tipos de medida): o R coloca NA automaticamente, mas lembrar de trocar o código "NA" (de não se aplica) no status, pois pode haver confusão. Trocar para "NE" (de não existe).
#2. Excluir observações
#3. Simplificar títulos das colunas!
#4. Não esquecer de arrumar os códigos das colunas de status
teste <- read.table("demo.gua.txt", header=T, sep="\t")
#Inserção das colunas com as áreas basais (para a somatória dos que tem caules múltiplos e dos que tem mais de um tipo de medida)
teste$ab.dap0 <- round((pi*(teste$dap0^2))/4,2) ## valores arredondados para dois dígitos
teste$ab.dap1 <- round((pi*(teste$dap1^2))/4,2)
teste$ab.das0 <- round((pi*(teste$das0^2))/4,2)
teste$ab.das1 <- round((pi*(teste$das1^2))/4,2)
teste$ab.pap0 <- round((teste$pap0^2)/(4*pi),2)
teste$ab.pap1 <- round((teste$pap1^2)/(4*pi),2)
#Somatória das áreas basais por placa
ab.final <- aggregate(teste[,c("ab.dap0", "ab.dap1", "ab.das0", "ab.das1", "ab.pap0", "ab.pap1")], by=list(teste$placa), FUN=sum, na.rm=T) ##indexação de todas as linhas do objeto teste e as colunas especificadas
colnames(ab.final)=c("placa", "dap0", "dap1", "das0", "das1", "pap0", "pap1") ##já que o arquivo já se chama ab, simplificamos os nomes das colunas
#Transformando as medidas dap e pap numa só: somatória dos valores dessas duas colunas por linha:
ab.final$dap1=rowSums(ab.final[,c("dap1","pap1")],na.rm=T) ##essa função rowSums aparece no help da sum(); essa indexação é muito importante de saber!
ab.final$dap0=rowSums(ab.final[,c("dap0","pap0")],na.rm=T)
head(ab.final) ##conferida a soma pelas plantas que tinham só pap
#Gerando um valor de dap final:
diam.final=sqrt(ab.final[,-1]*(4/pi)) ##transformando a área basal em um valor de dap final, exceto pela coluna 1 de ab.final, que é a coluna de placa!
diam.final[diam.final==0]=NA ## prestar atenção porque essa indexação é diferente(sem a vírgula!) é para matriz ou data.frame; aqui não queremos manter os zeros inseridos quando fizemos as somatórias, por isso, passamos tudo que é zero para NA
#Agora retornamos a coluna placa para esse novo objeto diam.final
diam.final$placa=ab.final$placa
#Só que a coluna placa foi adicionada no final e não precisamos mais das colunas de pap0 e pap1
diam.final=diam.final[,c(7,1:4)] ##reordenamos as colunas por indexação e ainda removemos as colunas 5 (pap0) e 6 (pap1)
diam.final ##ok, agora esse arquivo já tem os diâmetros finais para os dois tempos (dap e das)
### Preparando os dados para as transições ###
#Criando Status para o diam.final: precisamos adicionar colunas com os status0 e status1:
#Ao invés de buscar essas colunas no objeto teste, onde cada linha é um fuste, vamos criar as colunas status, pois são três apenas os estados possíveis.
#Status possíveis:
#V=viva
#M=morta
#NE=não existe status para as plantas recrutadas num tempo anterior ao recrutamento
#Primeiro criamos as duas novas colunas desses status com NA's:
diam.final$status0=NA
diam.final$status1=NA
diam.final ##conferindo...
#Definindo as condições para cada um dos status:
#Para o tempo 0:
diam.final$status0[is.na(diam.final$dap0) & is.na(diam.final$das0) ]="NE" ##na coluna status0 todos os indivíduos que tiverem dap0 e das0 iguais a "NA" atribuímos "NE" (ou seja, todos que não existiam no 0 mas existem no 1)
diam.final$status0[diam.final$dap0>0 | diam.final$das0>0]="V" ##na coluna status0 todos os indivíduos com dap0 ou das0 maior do que 0
#Para o tempo 1:
diam.final$status1[diam.final$status0=="V" & is.na(diam.final$dap1) & is.na(diam.final$das1)]="M" ##na coluna status1 todos os indivíduos que tem status0=V e que são NA no dap1 e no das1 atribuímos "M"
diam.final$status1[diam.final$dap1>0 | diam.final$das1>0]="V" ##por outro lado, todos em que o dap1 ou o das1 são maiores do que zero atribuímos "V"
#Uma questão importante: nunca fazer comandos redundantes, pois dessa forma não podemos enxergar possíveis erros.
#Por isso para a atribuição de cada status0 e 1 na primeira expressão fazemos uma indexação e na segunda outra e não no segundo comando simplesmente mandamos chamar de outro status tudo que não estava contemplado no primeiro
#table(is.na(diam.final$dap1)) ##apenas conferindo se havia mortas (quantos NA's no dap1 que não são das...)
#diam.final[is.na(diam.final$dap1),] ##quem eles são
diam.final$status0
diam.final$status1
diam.final
#placas.prob=diam.final[is.na(diam.final$status1),"placa"] ##apenas uns testes para os NA's que tinham aparecido estranhamente... mas já foi solucionado, só mantive aqui para guardar o comando match
#match(placas.prob,teste$placa)
#teste[match(placas.prob,teste$placa),]
diam.final[is.na(diam.final$status1),] ##checando os NA's após a correção mencionada acima, resta apenas a placa 4185, que está sem nenhum valor de medida
#Uma vez que os diâmetros finais estão prontos e os status também, podemos estabelecer as classes de tamanho:
### Estabelecendo as classes ###
#Checagem dos mínimos e máximos:
min(diam.final$das0,na.rm=T)
max(diam.final$das0,na.rm=T)
hist(diam.final$das0)
#Criação de uma coluna para a classe de tamanho do tempo zero (primeiro com NA's)
diam.final$classe0=NA
#Atribuição de classes a essa coluna:
diam.final$classe0[diam.final$das0>0]="p" ##classe "plântula" para todos que tem das
diam.final$classe0[diam.final$dap0>=10 & diam.final$dap0<40]="j1" ##classe "juvenil 1"; notar que a indexação é feita em dois momentos!
diam.final$classe0[diam.final$dap0>=40 & diam.final$dap0<70]="j2" ##classe "juvenil 2"
diam.final$classe0[diam.final$dap0>=70 & diam.final$dap0<100]="a1" ##classe "adulto 1"
diam.final$classe0[diam.final$dap0>=100]="a2" ##classe "adulto 2"
diam.final
#Criação de uma coluna para a classe de tamanho do tempo 1 (primeiro com NA's)
diam.final$classe1=NA
diam.final$classe1[diam.final$das1>0]="p"
diam.final$classe1[diam.final$dap1>=10 & diam.final$dap1<40]="j1"
diam.final$classe1[diam.final$dap1>=40 & diam.final$dap1<70]="j2"
diam.final$classe1[diam.final$dap1>=70 & diam.final$dap1<100]="a1"
diam.final$classe1[diam.final$dap1>=100]="a2"
diam.final
### As transições, finalmente! ###
#Preciso contar quantos permaneceram na mesma classe e quantos mudaram:
#Começando pelos que permaneceram na mesma classe:
table(diam.final$classe0==diam.final$classe1)
#FALSE TRUE
# 9 181 ##só nove transições!
#Mas são 192 observações, cadê as outras 2?
#Ah! são NA's nas classes, problemas a conferir na planilha original...
#Quais são então as transições que acontecem NESSE CONJUNTO DE DADOS?
diam.final[diam.final$classe0!=diam.final$classe1,]
# p -> j1
# j1 -> j2
# j2 -> a1
# a1 -> a2
table(diam.final$classe0)
table(diam.final$classe1)
#Achei que um comando desse tipo faria a contagem das combinações entre classes, mas não deu certo!
#xtabs(cbind(classe0, classe1) ~ ., data = diam.final)) ##não funciona!
#Calculando as transições uma a uma:
# p -> p: permanência
diam.final[diam.final$classe0=="p" & diam.final$classe1=="p",] ##indexando linha e coluna ele mostra quais são as linhas de dados que satisfazem as condições
table(diam.final$classe0=="p" & diam.final$classe1=="p") ##fazendo um teste lógico combinado com o table ele tabula o número de verdadeiros e falsos
#FALSE TRUE
# 145 46
table(diam.final$classe0=="p") ##quantos indivíduos pertenciam à classe p no tempo 0
#FALSE TRUE
# 143 48
tr.pp <- round(46/48,4) ##não consegui automatizar isso!
tr.pp
[1] 0.9583
####################TENTEI CRIAR UM COMANDO GENÉRICO, MAS NÃO DESCOBRI ######################################################################
#tr.pp <- table(table(diam.final$classe0=="p" & diam.final$classe1=="p")=="F") ##como faço para ele dizer o número de falsos?
#tr.pp <- round(table(table(diam.final$classe0=="p" & diam.final$classe1=="p")=="F")/192, 4) ##como faço para ele dizer o número de falsos?
#tr.pp
#dim(diam.final[1]) ##como faço para dizer para ele pegar o número de linhas do objeto e utilizar no cálculo?
##############################################################################################################################################
# p -> j1: crescimento
diam.final[diam.final$classe0=="p" & diam.final$classe1=="j1",]
table(diam.final$classe0=="p" & diam.final$classe1=="j1")
#FALSE TRUE
# 189 2
table(diam.final$classe0=="p")
#FALSE TRUE
# 143 48
tr.pj1 <- round(2/48,4)
tr.pj1
#[1] 0.0417
# p -> j2: crescimento
diam.final[diam.final$classe0=="p" & diam.final$classe1=="j2",]
table(diam.final$classe0=="p" & diam.final$classe1=="j2")
# FALSE
# 191
tr.pj2 <- 0/48
tr.pj2
#[1] 0
# p -> a1: crescimento - é zero, mas preciso calcular aqui?
diam.final[diam.final$classe0=="p" & diam.final$classe1=="a1",]
table(diam.final$classe0=="p" & diam.final$classe1=="a1")
# FALSE
# 191
tr.pa1 <- 0/48
tr.pa1
#[1] 0
# p -> a2: crescimento - é zero, mas preciso calcular aqui?
diam.final[diam.final$classe0=="p" & diam.final$classe1=="a2",]
table(diam.final$classe0=="p" & diam.final$classe1=="a2")
# FALSE
# 191
tr.pa2 <- 0/48
tr.pa2
#[1] 0
# j1 -> p: retração
diam.final[diam.final$classe0=="j1" & diam.final$classe1=="p",]
table(diam.final$classe0=="j1" & diam.final$classe1=="p")
# FALSE
# 191
table(diam.final$classe0=="j1")
#FALSE TRUE
# 125 66
tr.j1p <- 0/66
tr.j1p
#[1] 0
# j1 -> j1: permanência
diam.final[diam.final$classe0=="j1" & diam.final$classe1=="j1",]
table(diam.final$classe0=="j1" & diam.final$classe1=="j1")
#FALSE TRUE
# 126 65
tr.j1j1 <- round(65/66,4)
tr.j1j1
#[1] 0.9848
# j1 -> j2: crescimento
diam.final[diam.final$classe0=="j1" & diam.final$classe1=="j2",]
table(diam.final$classe0=="j1" & diam.final$classe1=="j2")
#FALSE TRUE
# 190 1
tr.j1j2 <- round(1/66,4)
tr.j1j2
#[1] 0.0152
# j1 -> a1: crescimento
diam.final[diam.final$classe0=="j1" & diam.final$classe1=="a1",]
table(diam.final$classe0=="j1" & diam.final$classe1=="a1")
#FALSE
# 191
tr.j1a1 <- round(0/66,4)
tr.j1a1
#[1] 0
# j1 -> a2: crescimento
diam.final[diam.final$classe0=="j1" & diam.final$classe1=="a2",]
table(diam.final$classe0=="j1" & diam.final$classe1=="a2")
#FALSE
# 191
tr.j1a2 <- round(0/66,4)
tr.j1a2
#[1] 0
# j2 -> p: retração
diam.final[diam.final$classe0=="j2" & diam.final$classe1=="p",]
table(diam.final$classe0=="j2" & diam.final$classe1=="p")
#FALSE
# 191
table(diam.final$classe0=="j2")
#FALSE TRUE
# 161 30
tr.j2p <- round(0/30,4)
tr.j2p
#[1] 0
# j2 -> j1: retração
diam.final[diam.final$classe0=="j2" & diam.final$classe1=="j1",]
table(diam.final$classe0=="j2" & diam.final$classe1=="j1")
#FALSE
# 191
tr.j2j1 <- round(0/30,4)
tr.j2j1
#[1] 0
# j2 -> j2: permanência
diam.final[diam.final$classe0=="j2" & diam.final$classe1=="j2",]
table(diam.final$classe0=="j2" & diam.final$classe1=="j2")
#FALSE TRUE
# 163 28
tr.j2j2 <- round(28/30,4)
tr.j2j2
#[1] 0.9333
# j2 -> a1: crescimento
diam.final[diam.final$classe0=="j2" & diam.final$classe1=="a1",]
table(diam.final$classe0=="j2" & diam.final$classe1=="a1")
#FALSE TRUE
# 189 2
tr.j2a1 <- round(2/30,4)
tr.j2a1
#[1] 0.0667
# j2 -> a2: crescimento
diam.final[diam.final$classe0=="j2" & diam.final$classe1=="a2",]
table(diam.final$classe0=="j2" & diam.final$classe1=="a2")
#FALSE
# 191
tr.j2a2 <- round(0/30,4)
tr.j2a2
#[1] 0
# a1 -> p: retração
diam.final[diam.final$classe0=="a1" & diam.final$classe1=="p",]
table(diam.final$classe0=="a1" & diam.final$classe1=="p")
#FALSE
# 190
table(diam.final$classe0=="a1")
#FALSE TRUE
# 170 21
tr.a1p <- round(0/21,4)
tr.a1p
#[1] 0
# a1 -> j1: retração
diam.final[diam.final$classe0=="a1" & diam.final$classe1=="j1",]
table(diam.final$classe0=="a1" & diam.final$classe1=="j1")
#FALSE
# 190
tr.a1j1 <- round(0/21,4)
tr.a1j1
#[1] 0
# a1 -> j2: retração
diam.final[diam.final$classe0=="a1" & diam.final$classe1=="j2",]
table(diam.final$classe0=="a1" & diam.final$classe1=="j2")
#FALSE
# 190
tr.a1j2 <- round(0/21,4)
tr.a1j2
#[1] 0
# a1 -> a1: permanência
diam.final[diam.final$classe0=="a1" & diam.final$classe1=="a1",]
table(diam.final$classe0=="a1" & diam.final$classe1=="a1")
#FALSE TRUE
# 174 16
tr.a1a1 <- round(16/21,4)
tr.a1a1
#[1] 0.7619
# a1 -> a2: crescimento
diam.final[diam.final$classe0=="a1" & diam.final$classe1=="a2",]
table(diam.final$classe0=="a1" & diam.final$classe1=="a2")
#FALSE TRUE
# 186 4
tr.a1a2 <- round(4/21,4)
tr.a1a2
#[1] 0.1905
##AQUI NÃO ESTÁ SOMANDO 1!ACHO QUE SÃO OS NA'S
# a2 -> p: retração
diam.final[diam.final$classe0=="a2" & diam.final$classe1=="p",]
table(diam.final$classe0=="a2" & diam.final$classe1=="p")
#FALSE
# 191
table(diam.final$classe0=="a2")
#FALSE TRUE
# 165 26
tr.a2p <- round(0/26,4)
tr.a2p
#[1] 0
# a2 -> j1: retração
diam.final[diam.final$classe0=="a2" & diam.final$classe1=="j1",]
table(diam.final$classe0=="a2" & diam.final$classe1=="j1")
#FALSE
# 191
tr.a2j1 <- round(0/26,4)
tr.a2j1
#[1] 0
# a2 -> j2: retração
diam.final[diam.final$classe0=="a2" & diam.final$classe1=="j2",]
table(diam.final$classe0=="a2" & diam.final$classe1=="j2")
#FALSE
# 191
tr.a2j2 <- round(0/26,4)
tr.a2j2
#[1] 0
# a2 -> a1: retração
diam.final[diam.final$classe0=="a2" & diam.final$classe1=="a1",]
table(diam.final$classe0=="a2" & diam.final$classe1=="a1")
#FALSE
# 191
tr.a2a1 <- round(0/26,4)
tr.a2a1
#[1] 0
# a2 -> a2: permanência
diam.final[diam.final$classe0=="a2" & diam.final$classe1=="a2",]
table(diam.final$classe0=="a2" & diam.final$classe1=="a2")
#FALSE TRUE
# 165 26
tr.a2a2 <- round(26/26,4)
tr.a2a2
#[1] 1
### AGORA JUNTANDO TUDO NUMA MATRIZ ###
dados <- c(tr.pp, tr.pj1, tr.pj2, tr.pa1, tr.pa2, tr.j1p, tr.j1j1, tr.j1j2, tr.j1a1, tr.j1a2, tr.j2p, tr.j2j1, tr.j2j2, tr.j2a1, tr.j2a2, tr.a1p, tr.a1j1, tr.a1j2, tr.a1a1, tr.a1a2, tr.a2p, tr.a2j1, tr.a2j2, tr.a2a1, tr.a2a2)
trans.matrix <- matrix(data = dados, nrow = 5, ncol = 5, byrow = FALSE)
trans.matrix
colnames (trans.matrix) <- c("p", "j1", "j2", "a1", "a2")
rownames (trans.matrix) <- c("p", "j1", "j2", "a1", "a2")
trans.matrix