Introdução

Na semana passada, eu publiquei esse post sobre o desmatamento da Amazônia. Um dos meus objetivos era fazer uma previsão do desmatamento para o ano de 2019 (período de Agosto de 2018 a Julho 2019). Apesar de conseguir fazer uma previsão, eu não consegui reproduzir as taxas anuais oficiais de desmatamento publicadas pelo Instituto Nacional de Pesquisas Espaciais (INPE), disponíveis para os anos de 1988 a 2018.

Não conseguir reproduzir as taxas me deixou muito incomodado. A reprodutibilidade é um dos pilares de pesquisas científicas, e serve como uma receita básica para garantir que resultados são confiáveis e transparentes, como também para garantir que compreendemos exatamente o que foi feito. Mas não consegui. Quando publiquei o post, recebi alguns comentários com orientações sobre como as estimativas são calculadas, e finalmente foi possível replicar as taxas anuais. Porém, deu muito trabalho, muito mais trabalho do que deveria ser necessário. Foi difícil encontrar os dados corretos, mas pior ainda foi reproduzir as análises como são feitas pelo INPE para obter as estimativas finais.

Essas dificuldades me fizeram pensar sobre o que quer dizer, de fato, reprodutibilidade. Pesquisando na internet sobre o tema, encontrei esse artigo (Fidler and Wilcox 2018), intitulado Reprodutibilidade de Resultados Científicos. Uma das questões que estava me incomodando foi esclarecida pelo artigo. Apesar de muitas vezes as palavras repetição, replicabilidade e reprodutibilidade serem utilizadas como sinônimos, elas não necessariamente significam a mesma coisa. Em geral, em disciplinas relacionadas a computação e a minha área de estatística, reprodutibilidade se refere à capacidade de reproduzir as computações realizadas, ou seja, o termo se refere exclusivamente a compartilhar e documentar suficientemente os dados e os códigos. Nessas disciplinas, replicabilidade, por outro lado, se refere ao refazer completamente o experimento.

Caso essa distinção não tenha ficado clara, tentarei esclarecer com um exemplo simplista. Suponha que você, leitor, jogue uma moeda três vezes e observe três coroas. Ao fazer as contas você obtêm como resultado que 100% dos arremessos foram coroas, o qual poderia até levar você a concluir que moedas sempre caem com a coroa virada para cima. Como um teste independente, seu amigo joga a moeda mais cinco vezes, e observa quatro caras e uma coroa. Ao fazer as contas e constatar como resultado que apenas 20% dos arremessos foram coroas, esse amigo pode concluir que os seus resultados estavam incorretos. Além de não conseguir replicar os resultados, para ele a moeda geralmente cai com a cara virada para cima.

Seu amigo replicou o seu experimento, mas os resultados foram diferentes. Na vida real, replicabilidade está associada a repetir o experimento todo mesmo esperando que os resultados não sejam necessariamente iguais. Pode tanto levar a correções de erros quanto levar ao aumento de nosso conhecimento. No nosso exemplo, a moeda, na verdade, tem a mesma chance de cair com a cara ou com a coroa virada para cima.

Já reprodutibilidade não se refere a refazer o experimento, mas sim a usar os dados obervados num experimento e refazer as contas para ver se é possível chegar ao mesmo resultado. No exemplo acima, seria o caso do amigo apenas contar novamente as suas três coroas e refazer seus cálculos para constatar que os resultados de fato mostram que 100% dos arremessos resultaram em coroa. Mesmo que a sua conclusão em si estivesse errada, seu amigo conseguiu reproduzir os seus resultados.

Reproduzindo as taxas anuais de desmatamento

Nesta seção discutirei como finalmente consegui reproduzir as taxas anuais de desmatamento da Amazônia produzidas pelo INPE. As informações sobre desmatamento se baseiam nos dados obtidos por um dos seus sistemas de apuração via satélite (PRODES). É a partir dessas imagens, que são gerados os dados utilizados para o cálculo das estatísticas do desmatamento. Mais informações sobre esse sistema podem ser encontradas nesse post.

Esse cálculo não é simples. Um dos problemas é que existem várias bancos de dados diferentes gerados pelos PRODES. É difícil para alguém sem experiência nesses bancos de dados entender qual deve ser utilizado. Ao fazer uma estimativa do desmatamento no post anterior, utilizei os dados georeferenciados. Com exceção das informações de geo-localização das áreas desmatadas (incremento com relação ao ano anterior), não existem muitas outras informações disponíveis nesse banco de dados. O dicionário das informações contidas nesse banco de dados não é muito detalhado. A descrição sobre o seu conteúdo é essa:

“Áreas desmatadas à partir de 2008 discretizadas por ano. O mapeamento utiliza imagens do satélite Landsat ou similares, para registrar e quantificar as áreas desmatadas maiores que 6,25 hectares. O PRODES considera como desmatamento a remoção completa da cobertura florestal primária por corte raso, independentemente da futura utilização destas áreas.”

Minha experiência na primeira tentativa de reproduzir as estimativas foi frustrante. Posteriormente descobri um arquivo em pdf com a metodologia de cálculo da taxa anual de desmatamento do PRODES. Esse relatório foi fundamental para entender como essas taxas são de fato calculadas. De maneira simplificada, depois de determinado o incremento de desmatamento observado em cada foto do satélite, são necessárias três etapas para se chegar às estimativas finais:

1. Ajuste para cobertura de nuvens: Em muitas imagens geradas pelo satélite, algumas áreas são cobertas por núvens. Nesse passo estima-se a quantidade de desmatamento que ocorreu embaixo das nuvens. Para isso presupõem-se que a taxa de desmatamento observada nas áreas descobertas se manterá nas áreas cobertas por núvens.
   
2. Estimativa proporcional para data de referência: A data de referência para as estimativas anuais é 1º de agosto. Para cada área da Amazônia, busca-se uma imagem com a menor cobertura de núvens possível. Para isso, imagens de diferentes períodos podem ser utilizadas. Para garantir que estimativas feitas a partir de imagens de diferentes datas mantenham a mesma data de referência, é necessário ajustá-las.
   
3. Detecção de Outliers: Com as estimativas em mãos, descartam-se ajustes muito grandes feitos por causa da cobertura das núvens. Nesses casos, utiliza-se apenas os incrementos visíveis de desmatamento.

Existem mais dois bancos de dados do PRODES. Um contém os dados de entrada, na qual existem informações sobre a cobertura de nuvens. Com esse banco é possível replicar as taxas anuais, porém é preciso calcular as estimativas das áreas sob as nuvens. Optei por trabalhar com o outro banco, que já contém as taxas calculadas. Nesse banco, todas as informações necessárias para o cálculo das estimativas estão disponíveis (para os anos de 2002 a 2018). Utilizando esse banco foi mais fácil verificar os cálculos feitos para reproduzir os resultados, pois os benchmarks estão no próprio banco de dados. Nesse caso, não é necessário reproduzir o cálculo do efeito de cobertura de núvens.

No relatório com a metodologia para o cálculo das taxas anuais, o segundo passo é descrito de forma bastante superficial. Mostra-se apenas um caso específico, onde todas as datas de referência e as imagens de satélites utilizadas estão dentro dos períodos de seca. Não entrarei em detalhes aqui, porém existem vários outros cenários possíveis. Especificamente, o relatório diz na página 30 que:

“O procedimento detalhado do cálculo da taxa dependerá das datas de aquisição da imagem no ano em análise e nos dois anos anteriores. Para ilustrar o procedimento, tomou-se o exemplo no qual as duas datas de aquisição das imagens estão dentro da estação seca (vide Figura 19). Os demais casos são calculados de maneira similar.”

Na imagem abaixo (figura 20 do relatório metodológico), ilustra-se como calcular a taxa anual para uma área específica no cenário descrito acima. Como pode-se ver na imagem, algumas medidas auxiliares são utilizadas. Destacarei aqui apenas Td2: taxa de desmatamento diária entre a imagem do ano analisado e a imagem do ano anterior e nd1: número de dias de estação seca entre a imagem do ano anterior e o final da estação seca.

Para calcular todos os outros casos ao qual a metodologia se refere levou muito tempo. E foi muito mais difícil do que o relatório dá a entender. Foram necessárias muitas tentativas e erros até conseguir reproduzir todas as medidas auxiliares. Ao final desse longo processo, descobri o que parece ser um erro, não um erro de cálculo do INPE, mas do manual sobre a metodologia. Em alguns casos, a quantidade nd1 utilizada para o cálculo da taxa anual tem que ser diferente do que a utilizada para calcular a taxa diária Td2. No manual, a mesma quantidade nd1 é sempre utilizada nos dois casos.

Na tabela 1 abaixo calculei essa quantidade conforme o manual sugere. Das 4837 áreas entre 2002 e 2018, em 3939 delas não há problema em utilizar a quantidade nd1 como descrita no manual para calcular tanto a taxa anual quanto a taxa diária. Porém, como é possível ver na linha destacada em laranja, em 898 áreas, se utilizarmos o mesmo nd1, sem alterações, para calcular a taxa diária Td2, essa estimativa estará errada em média -0,084 \(Km^2\) , e consequentemente a taxa anual em média por -8,3 \(Km^2\). A quantidade nd1 para o cálculo da taxa diária Td2 não pode depender da data de referência, definida como sendo o dia 211 do calendário juliano. Ao modificar isso, todas as estimativas são reproduzídas corretamente.

#carregando bibliotecas
require(tidyverse)
require(rvest)

#definindo função

get_data_saida <- function(url){

  print(url)
  pag <- read_html(url)
  df <- pag %>% html_table(header=TRUE)
  df <- df[[1]]
  if ("" %in% names(df) & !("Outliers" %in% names(df))){
    pos <- grep("^$",names(df))
      names(df)[pos] <- "Outliers"
  }
  if ("" %in% names(df)){
    df <- df[,!(names(df) == "")]
  }
  if ("Outliers" %in% names(df)){
    df <- df %>% mutate_at(vars(-Sta,-Outliers),list(~as.numeric(.)))
    df$Outliers <- str_replace(df$Outliers,"rate$","rates")
    df$Outliers <- str_replace(df$Outliers,"(cloud$|clds|clould)","clouds")
  } else {
    df <- df %>% mutate_at(vars(-Sta),list(~as.numeric(.)))
  }
  df$ano <- as.numeric(str_extract(url,"[0-9]{4}"))
  df <- df %>% filter(str_detect(PathRow,"^[0-9]{3,5}$"))
  df <- df %>% rename(pathrow=PathRow,uf=Sta)


  return(df)

}

#executando código

urls <- "http://www.obt.inpe.br/OBT/assuntos/programas/amazonia/prodes/taxas-calculadas-por-imagem-de-satelite-pos-2002"
pag <- read_html(urls)

links = pag %>% html_nodes(xpath="//a[@class='external-link']") %>% html_attr("data-mce-href")
links = links[!is.na(links)]

df.saida <- map_df(links,~get_data_saida(.))

#ajuste manual para um caso com informações faltantes no site
df.saida$Outliers[df.saida$pathrow == "22263" & df.saida$uf == "MA" & df.saida$ano==2005] <- "rates"

#alterando nomes das variáveis que serão reproduzidas
df.fim <- df.saida
df.fim <- df.fim %>% rename(rate.inpe=Rate,PercRate.inpe=PercRate,DRate2.inpe=DRate2,nd2r.inpe=nd2r,nd1.inpe=nd1,DRate1.inpe=DRate1,nd1r.inpe=nd1r,ndnxtY.inpe=ndnxtY,Outliers.inpe=Outliers)

#ordenando a base de dados.
df.fim <- df.fim %>% arrange(pathrow,uf,ano)

#reproduzindo resultados
df.fim <- df.fim %>% mutate(
  nd1a = case_when(
    Jul0 > EndClim ~ 0,
    Jul0 < StClim ~ EndClim - StClim + 1,
    TRUE ~ EndClim -  Jul0 + 1
  ),
  nd2a = case_when(
    Jul1 > EndClim ~ EndClim - StClim + 1,
    Jul1 < StClim ~ 0,
    TRUE ~ Jul1 - StClim + 1
  ),
  nd2 = case_when(
    Jul2 > EndClim ~ EndClim - StClim + 1,
    Jul2 < StClim ~ 0,
    TRUE ~ Jul2 - StClim + 1
  ),
  nd1_ = case_when(
    Jul1 > EndClim ~ 0,
    Jul1 < StClim ~ EndClim - StClim + 1,
    TRUE ~ pmax(EndClim - Jul1 + 1,0)
  ),
  nd1 = case_when(
    Jul1 > EndClim ~ 0,
    Jul1 <= StClim & EndClim < 211 ~ 211 - StClim + 1,
    Jul1 <= StClim & StClim < 211 ~ EndClim - 211 + 1,
    Jul1 <= StClim & EndClim >= 211 ~ EndClim - StClim + 1,
    Jul1 > StClim & Jul1 <= 211 ~ pmax(EndClim - 211 + 1,0),
    TRUE ~ pmax(EndClim - Jul1 + 1,0)
  ),
  nd2r = case_when(
    EndClim < 211 ~ pmax(EndClim - StClim + 1,0),
    EndClim >= 211 ~ pmax(211 - StClim + 1,0)
  ),
  nd1r = case_when(
    Jul1 <= 211 ~ 0,
    211 < StClim & Jul1 < StClim ~ 0,
    211 < StClim & Jul1 == StClim ~ 1,
    Jul1 > StClim & Jul1 <= EndClim & StClim > 211 ~ Jul1 - StClim + 1,
    Jul1 > StClim & Jul1 <= EndClim & StClim <= 211 ~ Jul1 - 211 + 1,
    Jul1 > EndClim & EndClim >= 211 & StClim > 211 ~ EndClim - StClim + 1,
    Jul1 > EndClim & EndClim >= 211 & StClim <= 211 ~ EndClim - 211 + 1,
    Jul1 > EndClim & EndClim < 211 ~ 0,
    Jul1 <= StClim ~ EndClim - StClim + 1
  ),
  DRate2 = ifelse(is.infinite(CorrInc / (nd1_ + nd2)),0, CorrInc / (nd1_ + nd2)),
  DRate2 = ifelse(nd1_ + nd2 == 0,0,DRate2),
  DRate2.w = ifelse(is.infinite(CorrInc / (nd1 + nd2)),0, CorrInc / (nd1 + nd2)),
  DRate2.w = ifelse(nd1 + nd2 == 0,0,DRate2.w),
  DRate1 = ifelse(is.infinite(CorrLstY / (nd1a + nd2a)),0,CorrLstY / (nd1a + nd2a)),
  DRate1 = ifelse(nd1a + nd2a == 0,0,DRate1),
  RATE = DRate2 * (nd1 + nd2r) + DRate1 * (nd1r),
  RATE.w = DRate2.w * (nd1 + nd2r) + DRate1 * (nd1r),
  PercRate =    round(100*((RATE / CorrInc)-1)),
  outliers = case_when(
    (PercClds >100 & Increm > 50) | (lag(PercClds) >100 & lag(Increm) > 50) ~ "cloulds",
    100 * (((RATE - DRate1*nd1r)  - CorrInc) / CorrInc) > 50 ~ "rates"
  ),
  final = ifelse(Outliers.inpe %in% c("rates","clouds"),Increm,rate.inpe),
  status.check = ifelse(abs(rate.inpe - RATE) <= 0.02,"OK","ERRADO"),
  status.w.check = ifelse(abs(rate.inpe - RATE.w) <= 0.02,"OK","ERRADO"),
  rate.check = round(rate.inpe - RATE,3),
  rate1.check = round(DRate1.inpe - DRate1,3),
  rate2.check = round(DRate2.inpe - DRate2,3),
  rate2w.check = round(DRate2.inpe - DRate2.w,3),
  nd1.check = nd1.inpe - nd1,
  nd1r.check = nd1r.inpe - nd1r,
  nd2r.check = nd2r.inpe - nd2r,
  out.check = Outliers.inpe ==outliers
)

xtabs(~df.fim$outliers + df.fim$Outliers.inpe ,addNA = TRUE)

##                df.fim$Outliers.inpe
## df.fim$outliers      clouds rates <NA>
##         cloulds    5     16     1    0
##         rates      9      0    41    0
##         <NA>    4480      3     0  282

Na tabela 2 abaixo reproduzimos as taxas anuais do PRODES, com a alteração acima mencionada da quantidade nd1 para o cálculo da taxa diária Td2. Todas as estatísticas intermediárias foram recalculadas para chegar a esse resultado, com exceção dos ajustes relacionados à cobertura das nuvens e da detecção de outliers. Utilizando a explicação sobre os outliers contida na metodologia, consegui reproduzir corretamente 95% dos outliers. Por causa disso, na reprodução mencionada acima, utilizamos os outliers previamente identificados na base de dados. Algumas pequenas diferenças podem ser observadas, porém acredito que elas se referem a erros de arredondamento.

df.ano <- df.fim %>% group_by(ano) %>% summarise(
  replicado = round(sum(final,na.rm = TRUE),0)
)

df.prodes <- tibble(
  ano = c(1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997,
          1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007,
          2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017,
          2018),
  prodes = c(21050, 17770, 13730, 11030, 13786, 14896, 14896, 29059,
             18161, 13227, 17383, 17259, 18226, 18165, 21651, 25396,
             27772, 19014, 14286, 11651, 12911, 7464, 7000, 6418, 4571, 5891,
             5012, 6207, 7893, 6947, 7536)
)
df.ano <- df.ano %>% left_join(df.prodes)
df.ano$dif <- df.ano$replicado - df.ano$prodes
df.ano$perc <- round(100 * (df.ano$dif / df.ano$prodes),1)

No contexto dos dados de desmatamento do PRODES, do meu ponto de vista reprodutibilidade significa facilitar ao máximo que qualquer pessoa leiga consiga reproduzir as estimativas anuais divulgadas pelo INPE, ou seja, que todos os arquivos importantes, como bases de dados, relatórios técnicos, dicionários de dados e código, estejam no mesmo local. Além disso, os dicionários de dados deveriam ser bem mais detalhados, e discutir com mais detalhe qual a finalidade de cada banco de dados como também a diferença entre eles.

Mas mais importante que tudo, seria disponibilizar um arquivo com o código necessário para reproduzir as estimativas das taxas anuais. Quando aprendemos a programar, é comum ouvir pessoas dizendo que o código é tão importante quanto aos bancos de dados, porque o código te mostra, exatamente, como saímos do ponto A e chegamos ao ponto B. Nesse caso, o ponto A é um dos bancos de dados do PRODES, e o ponto B seria as estimativas das taxas anuais do desmatamento. Um exemplo de um arquivo de código para reproduzir as taxas anuais está nesse link. Nesse arquivo, incluí todo o código necessário em R¹ para baixar os dados do banco com as taxas calculadas e reproduzir todas as taxas anuais de 2002 a 2018.

Conclusão

Reprodutibilidade não é algo desafiador. Inclusive pode ser bem chato de fazer, com foco nos detalhes. Dá uma sensação de ser algo redundante, quase que perda de tempo. Talvez esse post sobre reprodutibilidade das taxas anuais do PRODES tenha sido o mais chato que já fiz. Levei dias para conseguir reproduzir as taxas, o que me deixou bastante frustrado.

Entretando, eu também já estive do outro lado. Entre os anos de 2011 e 2013 trabalhei nos EUA, e num dos projetos no qual fui alocado o cliente final era o IRS americano, que é o equivalente a Receita Federal. Reproduzir todas as análises estatísticas é obrigatório não só para projetos do IRS, como também para qualquer projeto do Governo Americano. Nesse projeto, existia uma pessoa cujo trabalho era reproduzir, exatamente, minhas análises estatísticas. Não só isso, eu tinha que realizar os cálculos no mesmo escritório onde essa pessoa trabalhava. Para piorar, as análises eram em grande parte automatizadas. Não é preciso dizer o quanto isso me frustrava. Algumas vezes me dava a sensação que estavam me perseguindo, parecia que ele tinha algo pessoal contra mim.

Porém, em mais de uma oportunidade, o meu “reprodutor” encontrou erros nas minhas contas. Com o tempo, percebi que era muito melhor para mim, e para o projeto, que quando os erros ocorressem fossem descobertos rapidamente, antes de que se perpetuassem para o cliente final. E acredite, erros sempre vão ocorrer. Quando estamos muito envolvidos com algo, às vezes é dificil enxergar erros que são evidentes para outras pessoas menos envolvidas no projeto. Nos melhores casos, são erros novos, que não haviam ocorrido antes. Mas muitas vezes são erros bobos, facilmente evitáveis para quem está com a mente limpa e descansada. Eu defendi minha tese de doutorado em 2010. Para minha frustação, quase toda vez que leio algum trecho da mesma, encontro um erro novo!

Meu ponto é: reprodutibilidade pode parecer algo chato, algo pessoal de alguém que está pegando no seu pé. Porém é, de fato, algo importante, não só para qualquer projeto, mas também para a Ciência. É uma forma de entender melhor o que foi feito, além de corrigir erros. Por mais que não tenha tanto glamour quanto criar algo novo, é um passo importante para que o conhecimento seja passado para futuras gerações. Como disse Júlio Verne no livro “Viagem ao Centro da Terra”:

“Ciência…é feita de erros, mas erros que são importantes de se fazer, porque de pouco em pouco eles nos levam à verdade!”