Hipercubo Latino

Chalom, 2012

pse

Resumo do que eu entendi

O que é o Hipercubo Latino?

• É uma técnica utilizada para amostrar, do espaço paramétrico de um modelo, valores que serão utilizados como input das simulações
• O Hipercubo Latino parece ser uma técnica que permite que todas as sub-áreas das distribuições dos parâmetros sejam bem representadas e, em contrapartida, não é tão custosa quanto outras técnicas
• O Hipercubo Latino não exige que utilizemos técnicas de análise quantitativa específicas após sua aplicação

Antes de usar o Hipercubo Latino: conhecendo seus parâmetros

1. Quais são os parâmetros do modelo?
2. Definindo o espaço paramétrico: quais são os limites de valores que cada parâmetro pode assumir?
3. Identificar a distribuição de probabilidade de cada parâmetro, a partir do nosso conhecimento sobre a biologia do sistema
4. Construir a “joint probability distribution”. Esta função leva em conta não só a distribuição individual dos parâmetros, mas também a correlação entre eles

Usando o Hipercubo Latino!

Ver dúvidas abaixo…

Depois de usar o Hipercubo Latino: análise quantitativa do output

1. Análise de incerteza: quanto da variação nos valores utilizados para os parâmetros (input) é traduzida em variação nos resultados obtidos (output)
2. Análise de sensibilidade: quanto da variação nos resultados (output) pode ser atribuída à variação dos valores utilizados para cada parâmetro (input). Para isso, utilizar coeficientes de correlação parciais, que descontam o efeito de todos os outros parâmetros no resultado para analisar o efeito de um de cada vez. A relação entre os parâmetros e os resultados podem ser classificadas em: 1. Relação Linear (coeficiente de Pearson); 2. Relação Monotônica (coeficiente de Spearman); 3. Medidas de Tendência Central (?) (teste de Kruskal-Wallis); 4. Tendências em Variabilidade (?) (método FAST ou índices Sobol). Para uma boa análise, devemos aplicar todos os métodos.

Última coisa a fazer

• Depois de rodar o modelo com os parâmetros iniciais selecionados, talvez seja necessário aumentar o número de amostras. Mas como saber o número ideal de amostras? Ao compararmos duas ou mais “baterias” de simulações, em que cada uma tenha contido diversas combinações de parâmetros (amostras), podemos concluir que chegamos a um bom N amostral quando estas baterias apresentarem resultados similares nos testes de sensibilidade (i.e., apresentarem rankings de importância dos parâmetros semelhantes). E como definir se dois conjuntos de amostras são similares? SBMA (Symmetrized Blest Measure of Association).

E no final, o que nos resta?

• Dado que temos a “joint probability distribution” e que realizamos as análises quantitativas, devemos obter a distribuição de probabilidade de nossa variável resposta. Agora, sabemos qual a importância relativa de cada parâmetro para os resultados obtidos com o nosso modelo.

Observações

• Ao rodar as primeiras simulações, lembrar que, como se trata de um modelo estocástico, há variação nos resultados obtidos com diferentes amostras não só devido ao peso relativo de cada parâmetro que o compõe (incerteza epistêmica), mas a efeitos estocásticos inatos ao modelo (incerteza estocástica)

Dúvidas

• Como construir a “joint probability distribution”? E por que?
• Como identificar as possíveis correlações entre os parâmetros?
• Como realizar a análise de incerteza? Ela é a própria “joint probability distribution” ou a distribuição de probabilidade da variável resposta?
• Usando o Hipercubo Latino em si: para fazer a amostragem, devemos dividir as distribuições de probabilidade de todos os parâmetros em partes equiprováveis ou dividir diretamente a “joint probability distribution” (e retirar um valor - aleatório ou o ponto médio - de cada divisória)? E depois, como usar os valores selecionados para formar as combinações de valores que vamos utilizar de fato?

Bom, vamos lá, começar pelo começo.

— Estrutura do output — É importante decidir qual é a estrutura do output antes de rodar o LHS - porque essa estrutura vai decidir quais são as perguntas que você pode fazer.

O LHS foi pensado, lá pelos seus criadores, pra modelos nos quais a cada vez que você roda, produzem um único número. Então a gente agrega todas essas saídas do modelo (ou seja, uma lista de números), e compara com os parâmetros pra saber como a variação de cada parâmetro afeta a variação da saída. Por exemplo, será que aumentar o número de espécies inicial aumenta ou diminui o número de espécies final?

“Mas Chalom, você disse que dá pra fazer análise de vários outputs ao mesmo tempo!”

É verdade. Mas “vários outputs ao mesmo tempo” quer dizer: o programa que vc escreve é o mesmo, ele vai rodar as análises ao mesmo tempo, mas a análise de cada output é conceitualmente diferente. Então a gente analisa se “aumentar o número inicial de espécies aumenta o número final de espécies” usando o mesmo comando em que a gente analisa se “aumentar o número inicial de espécies aumenta o número total de propágulos”, mas isso é diferente de juntar as coisas e perguntar “aumentar o número inicial de espécies muda a saída do modelo como um todo”. Fez sentido?

(Tem um jeito de juntar tudo e perguntar “qual parâmetro influencia mais o resultado do modelo COMO UM TODO”, mas ele é meio experimental e eu não tenho muita certeza de como funcionaria. Se vc topar ser cobaia de laboratório junto comigo a gente pode tentar!)

E a análise de incerteza vai ser feita só com a distribuição da variável resposta. Ela não tem mais nenhuma relação com a joint probability distribution dos parâmetros

— Socorrooo, precisa reescrever meu código? —

Não! Dá pra gente escrever uma função que pega o resultado do modelo e extrai só o que a gente quer! Vamos supor que a hubbels.game está devolvendo uma lista na qual a primeira posição é um data frame, a segunda é uma lista de vetores e …. e a posição 27 é o número total de propágulos. A gente escreve uma função em volta da hubbels.game assim:

hubbels.game.externa ← function (s, a, nc, dc, nsim) { x ← hubbels.game(s, a, nc, dc, nsim) return (x[27]) }

Ou melhor, essa função externa pode ser usada pra calcular qualquer coisa em cima do resultado da hubbels.game, como usar um “apply”, “aggregate”, etc etc etc!

— Que diabos é a joint probability distribution? —

OK. A gente sabe a distribuição de cada parâmetro (x e y, variando uniformemente entre 0 e 1, por exemplo). Então a gente sabe distribuir em uma reta cada um dos parâmetros: por exemplo, metade dos x vai estar acima do 0.5, metade abaixo. Será que a gente sabe como é a distribuição 2D desses pontos? Quantos pontos vão ter AO MESMO TEMPO x < 0.5 e y < 0.5? A resposta mais intuitiva é que devem ser 25% dos pontos (mais 25% com x <0.5 e y > 0.5, etc etc). Mas a resposta certa pra isso éeeeee… não sei. Depende.

O ponto é resposta depende de como x e y se relacionam entre si. Vamos fazer o exemplo intuitivo:

library(pse)
LH1 ← LHS(model=NULL, factors=2, N=100, opts=list(COR=matrix(c(1,0,0,1), nc=2)))

Isso cria um objeto LH1 com dois parâmetros, variando uniformemente entre 0 e 1 (default), e FORÇA esses parâmetros a serem não correlacionados (também é default, mas eu decidi explicitar isso aqui ;). Vamos ver se isso é verdade mesmo:

cor(LH1$data)

Isso vai mostrar uma matriz de correlação com 1 na diagonal (x está correlacionado com x, ufa!) e 0 fora da diagonal. E vamos ver onde esses pontos gerados foram parar:

plot(LH1$data)

Bonitinho, né? Mas eu posso gerar outro hipercubo com AS MESMAS distribuições marginais e uma figura bem diferente!

LH2 ← LHS(model=NULL, factors=2, N=100, opts=list(COR=matrix(c(1,0.7,0.7,1), nc=2)))
cor(LH2$data)
plot(LH2$data)

Viu como a figura muda? Agora todos os pontos estão concentrados perto da linha x = y. Essa figura 2D é a “joint probability distribution” - e ela é afetada tanto pelas distribuições marginais que a gente escolheu quanto pela matriz de correlação dada. “Tá, mas COMO eu decido qual é a melhor matriz de correlação?” Pra uma análise exploratória, deixa no 0 (default). Isso só precisa ser levado em conta se a gente tem dados empíricos muito específicos (por exemplo, amostras sobre uma determinada espécie na qual a gente sabe que o tamanho da ninhada correlaciona fortemente com o tamanho corporal).

— Última dúvida! —

Usando o Hipercubo Latino em si: para fazer a amostragem, devemos dividir as distribuições de probabilidade
de todos os parâmetros em partes equiprováveis ou dividir diretamente a “joint probability distribution”
(e retirar um valor - aleatório ou o ponto médio - de cada divisória)?

Pra usar o LHS, a gente “picota” e embaralha cada distribuição marginal. Pegando uma distribuição normal de média 12 e sd 5, picotando em 50 partes e depois embaralhando tudo em um único comando:

N = 50
sample(qnorm(p = 1:N/N - 1/N/2, mean = 12, sd = 5))

(Isso é feito “por você” dentro da função LHS)

E depois, como usar os valores selecionados para formar as combinações de valores que vamos utilizar de fato?

Jeito “comprido”: crie um vetor representando cada um dos parâmetros:

N = 50
x ← sample(qnorm(p = 1:N/N - 1/N/2, mean = 12, sd = 5))
y ← sample(qnorm(p = 1:N/N - 1/N/2, mean = 54, sd = 3))

Defina sua função (vou fazer uma função bobinha) e execute um mapply dela nos seus parâmetros:

model ← function(x, y) x+y
result ← mapply(model, x, y)
plot(ecdf(result))

Análise de incerteza completa!

Jeito rápido: use a função LHS que “faz tudo por você”!

q ← “qnorm”
q.arg ← list(list(12, 5), list(54, 3))
model ← function(x) x[1] + x[2]
result ← LHS(model=model, N=50, factors=2, q=q, q.arg=q.arg)
plotecdf(result)

1. Eu entendi direitinho agora o que é a joint probability distribution, mas não entendi direito quando ela é usada no método do hipercubo (ou nas análises que se seguem). Ou ela não é usada e é “só” para nós compreendermos melhor nossos próprios parâmetros?

2. Devo interpretar o gráfico resultante da análise de incerteza (plotecdf) como: no eixo x, os possíveis resultados que obtenho com meu modelo; no eixo y, a proporção em que eles aparecem nas simulações que rodei? Entendi assim, mas não tenho total certeza.

1. A J.P.D. não é usada explicitamente em nenhum momento. Preciso melhorar isso no texto :P O que a gente usa é cada marginal distribution e, *se for o caso*, a matriz de correlação. As M.D. + a matriz de correlação são matematicamente equivalentes a uma J.P.D. ;)

2. Isso, só tomando cuidado que a linha é cumulativa. Ou seja, a cada ponto x, o y correspondente é a proporção das simulações que deram esse x OU MENOR.

← Voltar para Artigos