Crescimento Populacional e a Seleção Natural

A teoria da evolução é talvez a teoria mais importante para a biologia. Ela unificou muitos ramos científicos atribuídos a biologia que antes dela eram praticamente áreas distintas, como fisiologia, taxonomia, comportamento, etc.

Temos a frase celebre do nosso colega Theodosius Dobzhansky.

Nada em biologia faz sentido se não for à luz da evolução

Mas será que a gente realmente entende como a evolução funciona?

Quero dizer, provavelmente, todo mundo tem uma noção básica. O que a gente aprende na escola é tranquilo, mas o problema é que muitas vezes nos lembramos do exemplo básico e fazemos generalizações muito amplas, ou pior ainda, generalizações curtas e não conseguimos ver o processo como ele realmente ocorre.

Na verdade a teoria da evolução é muito mais do que vamos ver aqui, mas vamos olhar como funciona o negócio de variação no número de descendentes produzidos. Quem não lembra do exemplo da mariposa e da revolução industrial. A Biston betularia. Essa ai é aquela mariposa que existia de dois tipo, ou dois fenótipos, como preferir.

Tinhamos a versão, ou fenótipo mais claro.
claro

E a versão, ou fenótipo, mais escuro.
escuro

Ai antes da poluição, a mais clarinha era comum, e a escura era muito rara, com o acumulo de fuligem nos troncos das arvores, a clarinha começou a não ser muito boa para se camuflar, a escura, da cor da fuligem, começou a se dar bem e o cenário se inverteu, com a clarinhas passando a ser a mais rara e a escura sendo a mais comum.

Daqui pra frente vamos chamar de fenótipos sempre. Se a gente fosse la no campo e conta-se todas as mariposas existentes, poderíamos dizer que existe uma quantidade N de mariposas, essa quantidade sendo independente do fenótipo, que é a cor para as mariposas.

E elas vão vivendo, se reproduzindo a uma certa taxa, e bem isso nos faz lembrar da equação de crescimento populacional.

N_t = N_0 \cdot r^t

Lembre-se que estamos falando do crescimento discreto ai em cima, não confunda com a formula N_t = N_0 \cdot e^{r \cdot t} , a gente já discutiu sobre crescimento discreto versus crescimento contínuo aqui, de uma olhada para relembrar e ver a diferença.

Mas continuando, então temos nosso crescimento populacional, vamos ficar com ele por enquanto para facilitar o raciocinio, se começar a meter aquele número e de Euler vai ficar tudo mais difícil.

N_t = N_0 \cdot r^t

So que essa equação é valida para toda a população dessa determinada espécie, todos os indivíduos de uma região que compõem uma população de uma certa espécie, nossa mariposinha por exemplo.

Vamos pensar nos fenótipos claro e escuro como fenótipo A e fenótipo a, cada um deles então vai ter uma frequência na população, que chamaremos de f_A e f_a e cada um deles vai ter sua taxa de crescimento R, r_A e r_a, que também é conhecido como o fitness absoluto (Absolute fitness), que a princípio vamos considerar como igual a taxa de crescimento da população como um todo, escrevendo isso de forma mais clara, vamos assumir que r_A = r_a = R.

Primeiro, vamos pensar que o número inicial de indivíduos A seria igual a sua frequência no total inicial da população.

N_{A,0} = N_0 \cdot f_A

Certo, se a população inicial é N_0, sabemos que uma fração dela é de mariposas claras, a fração é f_A, por isso se multiplicar o total pela fração N_0 \cdot f_A da N_{A,0}

O mesmo vai ser valido para a população de a,

N_{a,0} = N_0 \cdot f_a

Lembrando que f_A + f_a = 1, tem que ser o total da população.

Vamos em frente, A e a são a população total dessa espécie.

N_0 = N_0 \cdot f_A + N_0 \cdot f_a

Sem muita mágica até aqui.
Certo, agora como o crescimento vai funcionar aqui?
Muito simples, para a próxima geração (vamos pensar apenas na diferença de uma geração para a próxima) temos a quantidade inicial vezes a taxa de crescimento, so que por cada um dos fenótipos.

Para A temos

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a

Mas como r_A = r_a = r, o fitness dos dois são iguais, então poderíamos simplesmente dizer.

N_{A,1} = N_0 \cdot f_A \cdot r

e

N_{a,1} = N_0 \cdot f_a \cdot r

Então podemos ver que a população total no tempo 1 N_1 vai ser

N_1 = N_0 \cdot f_A \cdot r + N_0 \cdot f_a \cdot r

Veja que o termo N_0, que é a população total inicial e o termo rocorrem nas duas multiplicações, então podemos simplificar isso.

N_1 = N_0 \cdot r \cdot (f_A + f_a)

Lembre-se que f_A + f_a = 1 já que são frequência, e esses dois A e a tem que ser o todo.

Então apos uma geração a frequência de A será a quantidade de indivíduos A pelo total de indivíduos da população, a frequência f_A será o seguinte.

{f_{A,1}} = {{N_0 \cdot f_{A,0} \cdot r}\over{N_0 \cdot r}}

Como é tudo uma multiplicação podemos cortar os N_0 e r e assim ficamos somente com.

f_{A,1}= f_{A,0}

A frequência na próxima geração não mudou nada.

Mas é isso que esperamos não? Porque nesse exemplo r_A = r_a = r, ou seja os dois fenótipos eram igualmente bons, mas e se r_A \ne r_a, um deles pode ser melhor que o outro, sei la r_A > r_a ou r_A < r_a, mas ambos os casos r_A \ne r_a, como fica esse mesmo esquema?

A principio, a próxima geração vai ter a mesma formula do crescimento.

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a,

Ou ainda, generalizando para o tempo t, temos

N_{A,t} = N_0 \cdot f_A \cdot {r_A}^t

E para a temos

N_{a,t} = N_0 \cdot f_a \cdot {r_a}^t

Ou seja a população N_t será

N_{t} = N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t

Agora vamos fazer aquela conta denovo, após t gerações, a frequência f_A vai ser a quantidade de individuos A sobre o total N

{f_{A,t}} = {N_0 \cdot f_A \cdot {r_A}^t \over{N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t}}

Agora nos podemos avaliar a frequência de a depois de t gerações assim como fizemos antes, mas podemos simplificar um pouco essa equação. Podemos começar tirando esse N_0 em cima e em baixo na divisão.

A frequência f_A vai ser o total dela sobre o tamanho total da população.
{f_{A,t}} = {f_A \cdot {r_A}^t \over{f_A \cdot {r_A}^t + f_a \cdot {r_a}^t}}

Agora podemos dividir em cima e em baixo por r_A e ficamos com o seguinte

{f_{A,t}} = {{f_A}\over{{f_A} + {f_a \cdot {{r_a}^t\over{{r_A}^t}}}}}

E finalmente temos a nossa equação para acompanhar como a frequência gênica muda ao longo das gerações.
Se a gente pensar um pouco, quando a razão {r_a}^t\over{{r_A}^t} é 1, ambos os fenótipos são igualmente bons, essa razão da um, e teremos  {f_A}\over{f_A + f_a \cdot 1}, como  f_A + f_a = 1, teremos  f_A dividido por 1 e nada nunca muda. Agora qualquer desvio na razão {r_a}^t\over{{r_A}^t} para mais ou para menos que 1, teremos mudanças.

Mas de qualquer forma, vamos jogar alguns números e fazer uma simulação no R.

Precisamos estabelecer algumas condições iniciais, e pensar por quantas gerações queremos simular.

# ************************************************** # Condições iniciais # ************************************************** N0 <- 1000 # População inicial rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a" fA <- 0.3 # Frequência do alelo "A" max_gen <- 20 # Número de gerações a simular

Veja que baseado naquelas quantidades, o resto dos ingredientes são derivados

# ************************************************** # Calculando variáveis derivadas # ************************************************** fa <- 1.0 - fA # Frequência do alelo "a" NA_0 <- N0 * fA # População inicial do alelo "A" Na_0 <- N0 * fa # População inicial do alelo "a"

Mas vamos la, criamos uma função para com aquelas condições iniciais e usando as equações que estabelecemos acima para ver o que acontece.

evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}

E testamos nossa funçãozinha.

> evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10) Nt NA_t Na_t fA_t fa_t [1,] 1000.000 300.0000 700.000 0.3 0.7 [2,] 1440.000 432.0000 1008.000 0.3 0.7 [3,] 1728.000 518.4000 1209.600 0.3 0.7 [4,] 2073.600 622.0800 1451.520 0.3 0.7 [5,] 2488.320 746.4960 1741.824 0.3 0.7 [6,] 2985.984 895.7952 2090.189 0.3 0.7 [7,] 3583.181 1074.9542 2508.227 0.3 0.7 [8,] 4299.817 1289.9451 3009.872 0.3 0.7 [9,] 5159.780 1547.9341 3611.846 0.3 0.7 [10,] 6191.736 1857.5209 4334.215 0.3 0.7 [11,] 7430.084 2229.0251 5201.059 0.3 0.7

Temos como entradas as quantidades iniciais, que repetimos na primeira linha, a geração zero, e temos na saída o tamanho da população total como primeira coluna, como ela vai mudando sendo cada geração uma linha, como cada fenótipo muda, bem como suas frequências ao longo das gerações, mas como olhar para uma matriz de números é meio chato, vamos representar esses dados em dois gráficos, um de como está a população e outro de como estão as frequências.

01

Certo, bem simples não, a população está crescendo exponencialmente, e as frequências genicas estão constantes, como previmos, porque as taxas de crescimento iniciais que colocamos para essa simulação são iguais.

rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a"

Agora lembre-se que nem tudo na vida são flores, se tivermos taxas de crescimento abaixo de 1, a população estará diminuindo.

Vamos pegar esse caso aqui.

#N0 = 1000 rate_A = 0.7 rate_a = 0.7 fA = 0.3

02

Mas ainda temos as frequências dos fenótipos iguais, aqui ta todo mundo se ferrando igual. Agora quando a gente diz que está havendo evolução é quando as frequência começam a mudar, quando aparece um novo mutante que derrepente se da melhor, se alguém tem a taxa de crescimento maior que sua contraparte, vemos o seguinte.

#N0 = 1000 rate_A = 2.0 rate_a = 1.5 fA = 0.02

03

Veja que aqui, nenhum dos fenótipos está necessáriamente se ferrando, só temos um fenótipo que é melhor que o outro, mas ninguém está diminuindo em quantidade ao longo do tempo, a população está crescendo como um todo, todos os fenótipos estão crescendo, apenas um mais rápido que o outro, mas isso faz com que a frequência dos fenótipos mude. Quando eu digo que as vezes a evolução pode não ser tão simples, porque como no exemplo das mariposas, a gente não pensa nesse caso assim, a gente sempre pensa que se alguém é pior, ele tem que estar se ferrando, como no próximo exemplo

#N0 = 1000 rate_A = 1.2 rate_a = 0.9 fA = 0.02

04

Essa é a situação que a gente tem na cabeça para as mariposas, pelo menos a primeira que vinha na minha, se a gente pensar na linha azul como as mariposas brancas e vermelho as mariposas escuras, a gente quer que alguém seja muito comum e alguém raro, ai esse cara mais comum, as mariposas brancas, começam a se dar mal por causa da revolução industrial, as mariposas escuras, que antes eram raras e ruim, começam a se dar bem, e temos uma inversão de quem é mais raro e quem é mais comum, mas vemos que a população decaiu um pouco e depois voltou a crescer. Mas como vimos, a evolução pode acontecer mesmo quando todo mundo está se dando bem, com a população crescendo, e além disso ainda, podemos pensar em outra situação.

#N0 = 10000 rate_A = 0.8 rate_a = 0.6 fA = 0.02

05

E é isso ai, ambos os fenótipos estão se ferrando, ambos estão diminuindo sua frequência. Cada vez que você olha a mata, você ve menos e menos daquela espécia que você observa, mas isso não quer dizer que a evolução não possa estar agindo ali. Veja que aqui não é quem se ferra e quem se da bem, aqui a situação é quem se ferra menos, quem pode perdurar aos tempos difíceis, talvez essa espécie seja extinta, mas talvez não, e quando os tempos mudarem, ela pode voltar a crescer, mas as frequências dos fenótipos estão alteradas já na próxima largada para o crescimento populacional.

Mas, não é a toa que nos humanos somos péssimos geradores de números aleatórios, pelo menos eu, sempre tendia a pensar em evolução, no caso de uma espécie que tem um fenótipo se tornando mais raro, porque a população dele está diminuindo, e outro se tornando mais comum porque a população dele está aumentando, mas esteja a população estacionada, aumentando ou diminuindo, as frequências de alelos, fenótipos, snp, como queria pensar, podem estar mudando, e todos esses cenários tem que estar na nossas cabeças.

So que a primeira coisa que vem na minha cabeça é um fenótipo se dando bem e outro se dando mal. Mas existem muitos mais possibilidades. Então é preciso se policiar para não tomar conclusões precipitadas, e conduzir uma análise de dados que de possibilidade de todos os casos serem avaliados corretamente.

Alias, para teste em sistemas mais complexos, sempre esses teste envolvem simulações com parâmetros aleatórios, porque muitas vezes o que podemos pensar não contempla todas as possibilidades, pensa em quantas árvores podemos desenhar com n espécies, eu nunca pensei que existissem tantas possibilidades assim, e é fácil se perder nas contas. Mas ao tentar algumas simulações, passar tudo para um modelo e ir tentando várias possibilidades de parâmetros, podemos talvez ver aquilo que nunca tinha passado pela nossa cabeça.

Bem é isso ai, as equações la do começo do post eu peguei de um curso do coursera, o Computational Molecular Evolution, ótimo curso alias, uma excelente seleção de tópicos para a ementa e conceitos bem complicados são explicados de maneira bem simples, pena que a maior parte pratica é usando o paup, mas agora é ir adaptando tudo pro R com as ferramentas que tivermos a mãos. Scripts no repositório recologia e aqui em baixo, e até o próximo post.

# **************************************************
# Condições iniciais
# **************************************************
 
N0      <- 1000 # População inicial
rate_A  <- 1.2  # Taxa de crescimento do alelo "A"
rate_a  <- 1.2  # Taxa de crescimento do alelo "a"
fA      <- 0.3  # Frequência do alelo "A"
max_gen <- 20   # Número de gerações a simular
 
# **************************************************
# Calculando variáveis derivadas
# **************************************************
 
fa   <- 1.0 - fA  # Frequência do alelo "a"
NA_0 <- N0 * fA   # População inicial do alelo "A"
Na_0 <- N0 * fa   # Polulação inicial do alelo "a"
 
# **************************************************
# Simulação
# **************************************************
 
evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}
 
evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10)
 
saida<-evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 50     rate_A = 1.2  rate_a = 1.2  fA = 0.3
saida<-evosimu(N0=50,rate_A=1.2,rate_a=1.2,fA=0.3,max_gen=20)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 0.7  rate_a = 0.7  fA = 0.3
saida<-evosimu(N0=100,rate_A=0.7,rate_a=0.7,fA=0.3,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topright",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 2.0  rate_a = 1.5  fA = 0.02
saida<-evosimu(N0=1000,rate_A=2.0,rate_a=1.5,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 1.2  rate_a = 0.9  fA = 0.02
saida<-evosimu(N0=1000,rate_A=1.2,rate_a=0.9,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 10000  rate_A = 0.8  rate_a = 0.6  fA = 0.02
saida<-evosimu(N0=10000,rate_A=0.8,rate_a=0.6,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")

Simulando sequências de DNA com Cadeias de Markov (Markov Chain)

03

Opa, uma forma excelente de entender melhor analise estatística, é simular dados. Eu sempre digo isso por aqui, mas aprendi nos livros do Marc Kéry.
Sempre que estamos simulando dados, estamos fazendo o processo inverso da análise de dados. Na análise a gente olha amostras e tenta entender como elas são geradas, quais as regras, para biologia as regras impostas pela natureza, que produzem essas amostras.

Agora para aprender melhor a construir árvores filogenéticas, usar elas para corrigir filogenia em análises como aqui, a gente precisa da capacidade de simular sequâncias de DNA, para testar nossas análises, se elas funcionam.Mas como sempre, nem tudo é tão simples como parece.

O DNA é composto de um alfabeto de quatro letras, A de adenina, C de citosina, T de timina e G de guanina.

O jeito mais simples de fazer isso, é simular o sorteio repetido dentro desse alfabeto.
Vamos supor um alfabeto mais simples, uma moeda, a moeda só cara ou coroa, então uma moeda honesta a gente simularia assim no R

> sample(c("Cara","Coroa"),30,rep=TRUE,prob=c(0.5,0.5)) [1] "Coroa" "Coroa" "Coroa" "Cara" "Coroa" "Cara" "Cara" "Cara" "Coroa" "Coroa" "Cara" "Cara" [13] "Cara" "Coroa" "Cara" "Coroa" "Cara" "Coroa" "Cara" "Cara" "Coroa" "Cara" "Coroa" "Coroa" [25] "Cara" "Coroa" "Coroa" "Coroa" "Cara" "Cara"

Usamos a função sample, que o primeiro argumento é um vetor com o alfabeto que queremos retirar as amostras, depois informamos que o sorteio é com repetição, o default é sem repetição. Temos que informar também quantas amostras queremos, aqui 30, e por ultimo um vetor de probabilidades, que para uma moeda honesta é 50% para cara e 50% para coroa.

Agora a gente pode seguir a mesma linha e simular dessa forma uma sequencia de DNA. Basta trocar o alfabeto usado no sorteio e o vetor de probabilidades, que tem que ser do tamanho do alfabeto, claro.

> sample(c("A","G","C","T"),30,rep=TRUE,prob=c(0.25,0.25,0.25,0.25)) [1] "G" "C" "G" "C" "G" "A" "T" "T" "A" "G" "G" "A" "T" "G" "C" "T" "C" "C" "A" "A" "C" "T" "A" "A" [25] "T" "G" "T" "A" "G" "C"

Mas talvez sequências geradas assim sejam simplesmente um monte de letras, e sequências de DNA não são tão aleatórias assim. Como assim?
Bem por exemplo, podemos citar o conteúdo de GC na sequência. A se liga com T e C com G, so que quando A se liga com T, essa ligação possui duas pontes de hidrogênio, enquanto C com G formam três pontes de hidrogênio, assim quanto mais o conteúdo de CG na cadeia, mais estável ela fica, do jeito que estamos gerando sequências, a tendência é elas sempre terem quantidades iguais de AT e CG, já que estamos sorteando cada base de forma independente.
Poderíamos gerar sequências com maior conteúdo de CG simplesmente aumentando a probabilidade dessas letras do nosso alfabeto de DNA, mas podemos também condicionar a simulação a quando tiver acumulando muito AT, termos uma chance maior de tender a um estado com mais CG usando cadeias de markov para simular o DNA.

Vamos tentar fazer um exemplo pra ficar mais fácil.

A gente tem uma sequência com variáveis \{X_1, X_2, X_3, X_4 , \dots, X_n \} e cada variável tem um estado E, estado pra gente é qual letra do do alfabeto um X_n qualquer está, sendo o alfabeto ATCG, bem E = \{A,\ C,\ ,T,\ G\}X_n.

Mas para pegar o jeito primeiro, vamos assumir um alfabeto mais simples, pense que E = \{1,\ 2\} e temos a sequência \{X_1, X_2, X_3, X_4 , \dots, X_n \}.
Quando olhamos para X_n = i qualquer, o nosso processo está no estado i no tempo n.

A expressão P(X_1=i) denota a probabilidade que o processo esteja no estado i no tempo 1.

O evento no qual o processo muda seu estado de i para j (transição) entre os tempos um para o dois corresponde ao evento (X_2 = j \vert X_1 = i), sendo que essa barrinha significa “dado que”, lembre-se de probabilidade condicional que a gente comentou aqui.
Essa probabilidade é denotada como P(X_2 = j \vert X_1 = i), mas de modo geral, a probabilidade de transição de i para j do ponto n ao n+1 pode ser descrita de forma mais geral como P(X_{n+1} = j \vert X_n = i). Dessa forma, podemos pensar ja nas probabilidades relacionadas a toda a sequencia que vamos gerar e probabilidades ligadas a todo o alfabeto, todos os estados.
Todas essas probabilidades pode ser visualizadas de uma forma mais simples numa matriz, a matriz de probabilidades de transição, onde cada elemento da matriz é P_{ij}(X_{n+1} = j \vert X_n = i).

A matrix de probabilidade de transição contem a probabilidade discreta (condicional) em cada uma de suas linhas (que somam 1 sempre). Para um processo de Markov, o estado no ponto n+1 depende do estado do ponto n, mas não nos estados dos tempos anteriores.

Então no nosso processo com o alfabeto E = \{1,\ 2\}, temos o seguinte.

  \begin{array}{ccc}   & 1 & 2 \\  1 & P_{11} & P_{12} \\  2 & P_{21} & P_{22} \end{array}

E como a gente viu la em cima, isso implica no seguinte.

  \begin{array}{ccc}   & 1 & 2 \\  1 & (X_{n+1} = 1 \vert X_n = 1) & (X_{n+1} = 2 \vert X_n = 1) \\  2 & (X_{n+1} = 1 \vert X_n = 2) & (X_{n+1} = 2 \vert X_n = 2) \end{array}

Certo, agora a gente pode colocar isso na pratica, a gente pode fazer uma função para simular sequencias de dna usando cadeias de Markov

markov1 <- function(x,P,n) {
    seq <- x
    for(k in 1:(n-1)){
        seq[k+1] <- sample(x, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}

Lembrando que nosso exemplo ainda não é o dna, os estados são E = \{1,\ 2\}.
Na função a gente recebe x, que é o alfabeto, dai podemos usar a função sample para sortear o próximo estado, e usar as probabilidades da linha da matriz que do estado atual.
Bem temos que iniciar também a cadeia, aqui vamos simplesmente começar com x.

Temos que estabelecer uma matriz de transição, vamos inventar alguns valores aqui.

  \begin{array}{ccc}   & 1 & 2 \\  1 & {5\over{6}} & {1\over{6}} \\  2 & {1\over{2}} & {1\over{2}} \end{array}

No R é so a gente fazer uma matriz.

> P <- matrix(c(1/6,5/6,0.5,0.5), 2, 2, byrow=TRUE) > P [,1] [,2] [1,] 0.1666667 0.8333333 [2,] 0.5000000 0.5000000

Ai temos que arrumar os nomes de linhas e colunas da matriz, para ela funcionar corretamente no samples, estabelecemos os estado atual e é só simular nossa sequencia.

> rownames(P)<-c(1,2) > colnames(P)<-c(1,2) > x<-c(1,2) > markov1(x,P,30) [1] 1 2 2 2 1 2 1 2 2 1 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 1 2 2 1 2

Certo, agora como isso funciona para sequências de DNA?
Vamos tentar adicionar a propriedade que falamos la em cima, vamos colocar uma maior quantidade de C e G na sequência.

P <- matrix(c(
1/6,5/6,0,0,
1/8,2/4,1/4,1/8,
0,2/6,3/6,1/6,
0,1/6,3/6,2/6),4,4,byrow=TRUE)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")

O que a gente precisa é pensar na matriz de transição.
Olha la, depois de um A, a gente poe uma grande chance de ir para C
Se a gente tiver em C, a gente tem a maior chance de ficar em C, segundo maior em G.
E assim vai.

> P a c g t a 0.1666667 0.8333333 0.00 0.0000000 c 0.1250000 0.5000000 0.25 0.1250000 g 0.0000000 0.3333333 0.50 0.1666667 t 0.0000000 0.1666667 0.50 0.3333333

A gente ainda tem que fazer uma pequena modificação, vamos mandar a parte uma chance de cada uma das bases para iniciar a cadeia, chamando ela de pi0. No mais temos o mesmo sistema anterior, e vamos guardar os resultado num vetor de character.

markov2 <- function(StateSpace,P,pi0,n) {
    seq <- character(n)
    seq[1] <- sample(StateSpace, 1, replace=TRUE, pi0)
    for(k in 1:(n-1)) {
        seq[k+1] <- sample(StateSpace, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}

Então geramos um vetor para iniciar a sequência.

pi0 <- c(1/4,1/4,1/4,1/4)

E bang, usamos nossa função, podemos ver o começo da cadeia com a função head, e podemos usar a função table e ver quanto temos de cada base.

> x <- markov2(StateSpace,P,pi0,1000) > head(x) [1] "c" "t" "t" "t" "g" "g" > table(x) x a c g t 54 399 364 183

E olha ae, temos bastante C e G, pouquinho A a um pouco mais de T.

Mudando um pouco a matriz de transição, podemos ainda aumentar a chance de alguns aminoácidos, lembre-se que três bases vão ser codificados em um aminoácido, As sequencia que produzem phenylalanine são TTT e TTC. No começo da tabelinha de codon a gente ve esses dois grupinhos.

Geramos uma nova matriz de transição

P <- matrix(c(.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.28,.01,0.70),4,4,byrow=T)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
x <- markov2(StateSpace,P,pi0,30000)

E geramos nossa sequência de DNA.

> library(seqinr) > table(getTrans(x)) * A C D F H I L M N P Q R S T V W Y 2 2 62 1 5157 20 105 2255 2 1 15 1 23 2192 1 93 1 67

Podemos usar a função getTrans para traduzir a sequencia em aminoácidos, mas a gente ja fez essa função aqui.
E pimba, temos 5157 phenylalanina, enquanto o aminoácido que tem a segunda maior quantidade é a Leucina que tem 2255.

Agora você pode pensar, poxa, agora como diabos eu vou ficar inventando essas matrizes de transição? Eu não sou criativo, como eu vou tirar esse monte de número da minha cabeça?
Mas você não precisa ficar inventado, a na verdade essa matriz de transição já existe, nas sequencias que a gente já sequenciou. Você pode pegar uma sequência qualquer no genbank e estimar a matriz de transição da qual ela veio.

A gente pode contar o número de transições e dividir pelo total da linha. Podemos assim estimar as probabilidades de transição de maior probabilidade, mais verossimilhantes. Alias qualquer probabilidade zero será exato nesse caso, porque serão casos que nunca existem.

Primeiro a gente conta as transições.

> nr <- count(x,2) > nr aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt 4 1 2 295 66 56 61 6419 1 4 5 267 231 6542 208 15837

Temos que fazer uma matriz para receber os resultados e colocar as contagens nele

A <- matrix(NA,4,4)
A[1,1]<-nr["aa"]; A[1,2]<-nr["ag"]; A[1,3]<-nr["ac"]; A[1,4]<-nr["at"]
A[2,1]<-nr["ga"]; A[2,2]<-nr["gg"]; A[2,3]<-nr["gc"]; A[2,4]<-nr["gt"]
A[3,1]<-nr["ca"]; A[3,2]<-nr["cg"]; A[3,3]<-nr["cc"]; A[3,4]<-nr["ct"]
A[4,1]<-nr["ta"]; A[4,2]<-nr["tg"]; A[4,3]<-nr["tc"]; A[4,4]<-nr["tt"]

E vemos o seguinte

> A a c g t a 4 2 1 295 c 1 5 4 267 g 66 61 56 6419 t 231 208 6542 15837

Igual a matriz de probabilidade de transições, mas temos contagens aqui.
ai a gente transforma essas contagens em probabilidade, dividindo pela soma da linha.

rowsumA <- apply(A, 1, sum)
Phat <- sweep(A, 1, rowsumA, FUN="/")
rownames(Phat) <- colnames(Phat) <- c("a","g","c","t")

E podemos agora comparar os valores originais

> P a c g t a 0.01 0.01 0.01 0.97 c 0.01 0.01 0.01 0.97 g 0.01 0.01 0.01 0.97 t 0.01 0.28 0.01 0.70

Contra os valores estimados

> round(Phat,3) a g c t a 0.013 0.007 0.003 0.977 g 0.004 0.018 0.014 0.964 c 0.010 0.009 0.008 0.972 t 0.010 0.009 0.287 0.694

Fica bem próximo não, claro que esperamos melhores resultado com sequências maiores, mas é um resultado expressivo. Pensando desse jeito, a gente poderia pegar qualquer sequência conhecida, estimar a matriz de transição dela e simular outras sequências para testar analises. Como a gente disse, a gente faz analise de dados para estimar um modelo, mas com o modelo a gente simula dados.

Ficamos por aqui, por agora, porque ainda vamos melhorar isso ae. Os scripts sempre no repositório recologia e é isso, até o próximo post.

Referências:
Wim P. Krijnen – Applied Statistics for Bioinformatics using R

sample(c("Cara","Coroa"),30,rep=TRUE,prob=c(0.5,0.5))
 
sample(c("A","G","C","T"),30,rep=TRUE,prob=c(0.25,0.25,0.25,0.25))
 
markov1 <- function(x,P,n) {
    seq <- x
    for(k in 1:(n-1)){
        seq[k+1] <- sample(x, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}
 
P <- matrix(c(1/6,5/6,0.5,0.5), 2, 2, byrow=TRUE)
 
rownames(P)<-c(1,2)
colnames(P)<-c(1,2)
x<-c(1,2)
 
markov1(x,P,30)
 
markov2 <- function(StateSpace,P,pi0,n) {
    seq <- character(n)
    seq[1] <- sample(StateSpace, 1, replace=TRUE, pi0)
    for(k in 1:(n-1)) {
        seq[k+1] <- sample(StateSpace, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}
 
P <- matrix(c(
1/6,5/6,0,0,
1/8,2/4,1/4,1/8,
0,2/6,3/6,1/6,
0,1/6,3/6,2/6),4,4,byrow=TRUE)
 
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
 
pi0 <- c(1/4,1/4,1/4,1/4)
 
x <- markov2(StateSpace,P,pi0,1000)
 
head(x)
table(x)
 
P <- matrix(c(.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.28,.01,0.70),4,4,byrow=T)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
x <- markov2(StateSpace,P,pi0,30000)
 
#install.packages("seqinr")
library(seqinr)
table(getTrans(x))
 
nr <- count(x,2)
nr
 
A <- matrix(NA,4,4)
rownames(A) <- colnames(A) <- c("a","c","g","t")
A[1,1]<-nr["aa"]; A[1,2]<-nr["ag"]; A[1,3]<-nr["ac"]; A[1,4]<-nr["at"]
A[2,1]<-nr["ga"]; A[2,2]<-nr["gg"]; A[2,3]<-nr["gc"]; A[2,4]<-nr["gt"]
A[3,1]<-nr["ca"]; A[3,2]<-nr["cg"]; A[3,3]<-nr["cc"]; A[3,4]<-nr["ct"]
A[4,1]<-nr["ta"]; A[4,2]<-nr["tg"]; A[4,3]<-nr["tc"]; A[4,4]<-nr["tt"]
 
A
 
rowsumA <- apply(A, 1, sum)
Phat <- sweep(A, 1, rowsumA, FUN="/")
rownames(Phat) <- colnames(Phat) <- c("a","g","c","t")
 
P
round(Phat,3)

A distribuição Chi-quadrado

Ola, vamos dar uma olhada aqui na distribuição chi-quadrado e tentar entender melhor porque da para usar ela para selecionar modelos.

A gente já falou da distribuição f aqui e da distribuição t aqui, e agora adicionando a distribuição chi-quadrado ao repertório.

A distribuição chi-quadrado sempre tem um papel importante no teste de hipóteses sobre frequências.

Basicamente, se a gente tem um conjunto de n elementos \{Z_1,Z_2,\cdots,Z_n\} independentes e normalizados, a soma dos quadrados deles vai ter uma distribuição aleatória chi-quadrado com n graus de liberdades.

Então:

\chi^2=\sum\limits_{i=1}^n\{Z_i^2\}

Mas o que quer dizer normalizado? No post anterior, aqui, a gente realizou essa operação, é só diminuir o valor pela média e dividir pelo desvio.

Na distribuição normal a gente faz:

{x-\alpha} \over \sigma

Agora no teste de qui-quadrado a gente faz o famoso:

{Observado\ -\ Esperado} \over {Esperado}

Certo, agora vamos contextualizar isso.
Na graduação, muita gente deve ter participado em algum momento de um trabalho, pelo menos eu participei, de reproduzir moscas Drosophila melanogaster. Famosas nos trabalhinhos de genética.

Chega uma hora que a gente deixa elas cruzarem e depois ia ver se havia realmente uma segregação independente dos cromossomos. No caso do cromossomos ligados ao sexo, temos fêmeas XX e machos XY.

Ai lembrando as leis do nosso colega Gregor Mendel , podemos fazer a quadrado de punnet.

Fêmea X X M X XX XX a c Y XY XY h O

Ou seja esperamos 50% de machos e 50% de fêmeas, uma razão de 1:1.

Então se nasceu 100 moscas, dessas 50 machos e 50 fêmeas, isso é perfeitamente o que esperamos.
Mas se nascer 49 macho e 51 fêmeas, não é exatamente o esperado, mas esta bem próximo, agora se a gente ver 5 machos e 95 fêmeas, isso é bem longe do que esperamos. Mas como transformar esse próximo ou longe em um número?
Bang, chegamos ao chi-quadrado. Vamos fazer uma função de chi-quadrado no R para o nosso caso e fazer alguns teste.

chi2_moscas<-function(observado) {
    if(is.vector(observado) && length(observado)==2) {
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
}

então essa função recebe como entrada um vetor de dois números, a nossa contagem de machos e fêmeas, e vai calcular a estatística chi-quadrado, dado que nossa hipótese para as mosquinhas realmente seja uma razão sexual de 1:1, como o senhor Mendel propôs, não estamos tirando da cartola esse valor de 50% para cada sexo.

Veja que na função, a gente adicionou um teste para ver que a gente está entrando com os dados corretamente, caso alguém mande os dados errados, é comum a gente ver verificações em funções, fizemos uma aqui.
Usando nossa função a gente vê.

> chi2_moscas(c(200)) [1] "Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!"

Opa, entramos com um dado que só tem uma contagem, isso não faz o menor sentido, então não calculamos nada, a gente só faz conta quando faz sentido.

Agora vamos supor que pegamos nossa garrafinha la, e contamos 221 moscas macho e 475 moscas fêmeas, nossa isso ta bem fora da casinha, ta muito diferente de 1:1, se a gente calcula o chi-quadrado a gente vê…

> chi2_moscas(c(221,475)) [1] 92.6954

Um número grandão. Mas é grande mesmo? Vamos pensar que em outra garrafa, a gente pegou 332 machos e 375 fêmeas, isso é bem mais razoável segundo a hipótese não? Ta mais com cara de 1:1, e quando a gente calcula o valor de chi-quadrado a gente vê…

> chi2_moscas(c(332,375)) [1] 2.615276

Um valor baixinho. Legal, então um valor grandão quer dizer uma discrepância muito grande da nossa hipótese, e um valor baixo quer dizer que estamos diante de algo bem próximo da hipótese.
Isso não é difícil de enxergar com a formula que estamos fazendo.

Estamos calculando na função chi-quadrado o seguinte

{\chi^2}={\sum\limits_{i=1}^2} { {(Observado\ -\ Esperado)^2} \over {Esperado} }

Veja que vai até 2 a somatória, porque temos machos e fêmeas aqui nesse problema. O quadrado é para não ter problemas com números negativos, mas o que acontece é que se o observado é próximo do número esperado, vai dar um número pequenininho ali no numerador, se for exatamente o mesmo número, da zero, e zero dividido por qualquer coisa da zero, se for um número pequeno, vai ser algo pequeno dividido por um número que vai dar um número pequeno, agora conforme essa diferença aumenta, o numerador aumenta, e cada vez o número fica maior.

Continuando com esse número total de 100 moscas, para ficar fácil fazer contas, podemos calcular o chi-quadrado para os casos de 49-51, 48-52 … até 1-99 e ver o que acontece.

01

Olha ai, conforme a discrepância em relação a nossa hipótese aumenta, o valor de chi-quadrado aumenta também, exponencialmente.

Mas num é qualquer discrepância que a gente vai para e assumir que o senhor Mendel tava errado. A gente tem que lembrar como as coisas estão acontecendo.

Cada mosca tem 50% de chance de ser macho o fêmea, independente da contagem atual, pode já ter nascido 99 moscas fêmeas, a próxima, se Mendel está certo, tem 50% de ser macho. Nos, seres humanos, somo péssimos geradores de números aleatórias por isso, se você contar 99 fêmeas, você vai primeiro pensar que algo acontece, as vezes sem levar em conta que esse resultado é perfeitamente possível dentro da teoria de Mendel, ele é raro, mas perfeitamente possível.

Agora a gente pode fazer um simulação no R e ver o quão raro é um caso desse, alias, a gente pode ao invés de guardar quantas moscas de qual sexo nasceram, simular os nascimentos e guardar somente o valor de chi-quadrado para o resultado.

Então é so a gente ir simulando os nascimento usando a função sample

> sample(c("Macho","Femea"),100,replace=T) [1] "Femea" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Femea" [12] "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Macho" "Macho" "Femea" "Femea" [23] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" [34] "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [45] "Femea" "Macho" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" [56] "Macho" "Femea" "Femea" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" "Macho" "Macho" [67] "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" [78] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [89] "Macho" "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" [100] "Macho"

E ir guardando o resultado num vetor, já que é somente um número.

for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}

Agora para facilitar a, a gente olha o resultado num histograma.

02

Humm, então existem poucas chances de valores extremos, mas valores próximos de zero são bem comuns.
Bem quando alguém descreve uma distribuição, o cara faz uma probability density function ou PDF, que nada mais é que uma função que descreve esse resultado. O R ja tem pronta um monte desses PDF, normalmente é d grudado no nome da distribuição, aqui é o dchisq.

Vamos dar uma olhada, se é parecido com o que encontramos.

03

Olha ai que legal, fica bem parecido com o que encontramos, ou seja antes da gente fazer a simulação, ja tinha como saber como seria o resultado com essa distribuição. Eu marquei o quantile de 95% também. O que esse risco azul ta mostrando, que 95% dos casos estão a esquerda desse número, e 5% a direita.

É fantástico pensar que isso funciona, podemos comparar esse numero que calculamos com o quantile com os valores de chi-quadrado que calculamos e ver se ele funciona de verdadinha.

> sum(estatistica<qchisq(0.95,1))/length(estatistica)
[1] 0.9456

E num é que funciona, 0.9456 são menor que ele, isso é praticamente 95%, muito boa predição.

Agora se o nosso valor que encontramos na garrifinha for muito extremo, podemos começar a suspeitar da teoria do Mendel, não que isso implique que ela é errada, mas podemos suspeitar, isso é o tal do p ser significativo, é a gente começar a suspeitar que algo estranho ta acontecendo, esse algo é permitido acontecer, mas incomum, muito incomum.

Vamos pensar em dois casos aqui, vamos pensar que contamos 43 machos e 57 femeas, se a gente calcular o valor de qui-quadrado

> experimento1<-c(43,57) > chi2_moscas(experimento1) [1] 1.96

Da 1.96, isso é bem baixinho, se a gente compara com os resultados da nossa simulação

> sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
[1] 0.1277

tem pelo menos 13% dos valores tão baixos quanto 1.96, não tem muito porque levantar suspeitas sobre a teoria de Mendel, e podemos comparar o procedimento que estamos fazendo com o teste de chi-quadrado já pronto do R, usando a função chisq.test

> chisq.test(experimento1,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento1 X-squared = 1.96, df = 1, p-value = 0.1615

E a gente tem o mesmo valor da estatística calculada, 1.96 (estamos fazendo a mesma continha), e um valor de p que traz a mesma conclusão, 0.16.

Agora vamos pensar que em outra experimento vemos um resultado bem diferente, 21 machos e 79 fêmeas.
Aqui começa a ter cara de um desvio na razão sexual nervoso. A gente pode calcular o valor de chi-quadrado e vemos o seguinte.

> experimento2<-c(21,79) > chi2_moscas(experimento2) [1] 33.64

Nossa, 33.64, se a gente comprar com os valores que obtemos ao acaso, fazendo a simulação onde machos e fêmeas tem 50% de chance cada de nascer.

> sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
[1] 0

Nenhum, absolutamente nenhuma das 10 mil simulações foi tão extrema assim. Ou seja, algo deve estar acontecendo, ou a teoria do Mendel esta errada aqui, ou alguém esta pregando uma peça no laboratório, não sabemos o que, varias coisas podem influenciar esse resultado, mas sabemos que se for para confiar no azar para ver um resultado assim, esse resultado é muitoooo raro.

Podemos fazer o teste de chi-quadrado com a função do R aqui também.

> chisq.test(experimento2,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento2 X-squared = 33.64, df = 1, p-value = 6.631e-09

E olha ai, um valor de p de 6.631e-09, lembrando que isso está em notação cientifica, esse e-09 significa que esse número é…

> format(6.631e-09,scientific=F) [1] "0.000000006631"

Muito legal né, então podemos usar essa distribuição para compara se algo está próximo da hipótese ou longe. Agora outro coisa, além de razão sexual, que tem uma distribuição que segue uma distribuição chi-quadrado são valores de verosimilhança, ou likelihood, e AIC também (AIC é loglikelihood corrigido meu, claro que tem que ser igual o likelihood quanto a distribuição).

E onde entra o chi-quadrado aqui?

Imagine que temos duas variáveis, um preditor e uma resposta, mas que nao tem relação nenhuma, os dois são amostras ao acaso.

> preditor<-rnorm(30) > resposta<-rnorm(30)

A gente pode tentar modelar eles tanto como variáveis que tem uma relação entre si, como um modelo onde a resposta simplesmente vem de uma distribuição normal, ou seja, assumir que ela simplesmente vem de uma distribuição normal é mais simples que assumir que ela tem relação com algo.

> modelo1<-lm(resposta~preditor) > modelo2<-lm(resposta~1)

A gente pode avaliar esses dois modelos e seus parâmetros.

> summary(modelo1) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.7239 -0.6622 -0.1236 0.5255 2.1173 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06622 0.17136 0.386 0.702 preditor 0.00845 0.19107 0.044 0.965 Residual standard error: 0.9379 on 28 degrees of freedom Multiple R-squared: 6.984e-05, Adjusted R-squared: -0.03564 F-statistic: 0.001956 on 1 and 28 DF, p-value: 0.965 > summary(modelo2) Call: lm(formula = resposta ~ 1) Residuals: Min 1Q Median 3Q Max -1.7339 -0.6600 -0.1162 0.5272 2.1214 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06592 0.16826 0.392 0.698 Residual standard error: 0.9216 on 29 degrees of freedom

Mas e ai, qual é melhor? Ue o que tiver o melhor likelihood, o menor custo, aquele que acumula menos resíduo. Mas essa diferença pode ser ao acaso, a diferença de likelihood, e podemos testar isso usando a distribuição chi-quadrado.

> anova(modelo1,modelo2,test="Chisq") Analysis of Variance Table Model 1: resposta ~ preditor Model 2: resposta ~ 1 Res.Df RSS Df Sum of Sq Pr(>Chi) 1 28 24.628 2 29 24.630 -1 -0.0017202 0.9647

Olha ai, comparamos esses modelos e vemos que eles são iguais, tem a mesma capacidade de predição. Então por parcimônia, é melhor eu escolher o mais simples uai, porque ficar com algo mais complexo, se a complexidade nao ajuda em nada, essa é a idea da navalha de occan, e vamos seguir ela, mas tendo um teste para sustentar essa escolha, pode diminuir a discussão entre nos e nossos pares, afinal se todo mundo concorda com esses métodos, vamos todos chegar as mesmas conclusões.

Ficamos por aqui, agora com a distribuição chi-quadrado melhor compreendida, espero eu, os scripts sempre no repositório recologia e é isso, até o próximo post.

set.seed(123)
chi2_moscas<-function(observado) {
 
    if(is.vector(observado) && length(observado)==2) {
 
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
 
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
 
}
 
 
chi2_moscas(c(200))
chi2_moscas(c(221,475))
chi2_moscas(c(332,375))
 
 
dados<-vector()
 
for(i in 1:50) {
    dados[i]<-chi2_moscas(c(50-i,50+i))
}
 
#figura 1
plot(dados,frame=F,type="p",pch=19,cex=0.5,xaxt="n",ylab="Estatisthica Chi-quadrado",xlab="Razão de Machos para Femeas")
axis(1,at=c(1,13,25,38,49),label=c("1:1","1:13","1:25","1:38","1:49"))
 
 
sample(c("Macho","Femea"),100,replace=T)
 
estatistica<-vector()
 
for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}
 
#figura 2
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
 
 
 
#figura 3
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
curve(dchisq(x,1),0,15,add=T,lwd=3,lty=2,col="red")
abline(v=qchisq(0.95,1),col="blue",lty=3)
legend("topright",lwd=c(3,1),col=c("red","blue"),lty=c(2,3),legend=c("Distribuição Chi-Quadrado","Quantile de 95%"),
       bty="n")
 
 
sum(estatistica<qchisq(0.95,1))/length(estatistica)
 
 
experimento1<-c(43,57)
chi2_moscas(experimento1)
sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
chisq.test(experimento1,p=c(0.5,0.5))
 
 
experimento2<-c(21,79)
chi2_moscas(experimento2)
sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
chisq.test(experimento2,p=c(0.5,0.5))
 
format(6.631e-09,scientific=F)
 
preditor<-rnorm(30)
resposta<-rnorm(30)
 
modelo1<-lm(resposta~preditor)
modelo2<-lm(resposta~1)
 
summary(modelo1)
summary(modelo2)
 
anova(modelo1,modelo2,test="Chisq")

Machine Learning – O algoritimo Gradient descent

Vamos denovo meter o nariz na area de Machine Learning, a gente tinha começado a falar disso aqui, foi complicado, mas agora vamos ver algo mais leve e legal.

Alias, na faculdade, em geral, acho que em quase todo curso que tem calculo, todo mundo odeia calculo.

math-rage

Mas um ponto ruim, é que a gente tem que aprender um monte de coisa primeiro, e só depois usar.
Calculo sempre tem aqueles exemplo de velocidade, aceleração, etc. Isso eu acho bem chato.

Mas calculo tem sua utilidade para optimização, que é mais legal que velocidade, na minha opinião. Mas ai você deve pergunta, como isso funciona?

A gente ja comentou, por exemplo aqui, que sempre o que a gente ta fazendo, é ajustando um modelo, que tem uma parte determinística, e uma parte estocástica.

Acontece que o gradiant descent é um algoritmo de optimização, que é baseado em derivada e serve para otimizar.

Mas o que é otimizar?

Quando a gente fala em forrageamento ótimo, como aqui e aqui, a gente tem em mente coisas como ganhar o máximo de energia no menor tempo possível, ou ficar o menor tempo possível exposto, procurando comida.

Aqui é a mesma coisa, queremos normalmente achar o máximo, ou o minimo de alguma função. Mas máximo ou minimo do que? Aqui, para a regressão linear, da verossimilhança, ou likelihood.

Existe um custo, que é o tamanho do erro que vamos ter de acordo com os valores de parametros que assumimos. Custo aqui é o tamanho do desvio que temos para um dado set de parâmetros.

Está confuso certo? Vamos tentar ver numa figura. Suponha que temos uma amostra de 30 valores, e para cada amostra temos o valor de um preditor, que queremos usar ele da melhor forma possível.

Primeiro, vamos determinar aqui, qual o nosso modelo?

Bem, para a regressão linear é:

y=\alpha+\beta \cdot X

Veja que essas letras são whatever, podemos usar a mesma letra com subscritos também:

y=\theta_1+\theta_2 \cdot X

Então a gente tem essa formula, e quer achar os melhores valores de \theta_n para um conjunto de dados de resposta, dado nossos preditores.

Vamos simular, como sempre alguns valores, e olhar como podemos escolher dois valores de \theta_n e dar um parecer, se eles são bons ou ruins, qual o custo deles.

Se a gente assumir \theta_1 = 1 e \theta_2 = 1, temos a formula y=1+1 \cdot X e vemos o seguinte em relação os pontos.

01

Agora essa escolha dos valores \theta_1 = 1 e \theta_2 = 1 é ruim? Não sei, mas podemos calcular a distancia de y até o que esperaríamos de acordo com nossa formula.

> (1+1*X[1])-y[1] [1] -7.841167

Esse -7 é o quão ruim nossa previsão é, e podemos repetir isso para todas as amostras.
Mas dependendo do modelo, esses valores vão ser positivos ou negativos, para evitar isso elevamos essa operação de diminuição ao quadrado, e é só somar.

Em termos matemáticos, o que estamos falando que esse custo é:

Temos nossa hipótese h que é:

h=\theta_1+\theta_2 \cdot X^i

e o custo será:

J(\theta_1,\theta_2)={1\over{2\cdot m}}\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2

Porque esse custo está assim?
Bem o 1\over{2\cdot m} a gente tem que usar, para ponderar o número de amostras, senão, como temos uma somatoria, esse número so vai crescer, crescer, crescer. Como a gente ta dividindo por algo relacionado ao número de amostras, ele ta ficando como uma média.
Aqui, a gente ta dizendo, h_\theta(x^{(i)}), se usarmos essa hipótese h=\theta_1+\theta_2 \cdot X^i, teremos para um valor de X, tal previsão para y.
E no final de tudo, estamos elevando tudo ao quadrado, porque se não for assim, quando a linha estiver passando acima do ponto, da um valor positivo, quando ta abaixo um valor negativo, como tudo elevado ao quadrado é positivo, a gente consegue ver um número relativo a distancia como no gráfico.

Bem se está difícil de entender, podemos escrever isso como uma função no R.

computeCost<- function(X, y, theta) {
    m <- length(y)
    J <- 0
    J = (1/(2*m)) * sum((X%*%theta-y)^2)
    return(J)
}

A gente entra com os valores de X, y e os valores de theta que queremos e calculamos o custo. Que intuitivamente, podemos pensar como algo proporcional a somar todos os pauzinhos vermelhos no gráfico que fizemos e dividir essa soma pelo número de amostras.

Aqui a função deve estar meio confusa, mas porque estamos usando álgebra de matrizes. No R %*% multiplica matrizes. De uma olhada aqui para entender melhor porque isso da certo, mas da, e se eu explicar aqui vou enrolar demais.

Mas o ponto que quero chegar, podemos calcular esse custo, se os parâmetros estão bons ou não, mas isso so faze sentido quando comparamos uma escolha de parâmetros com outra. Aqui se a gente escolher \theta_1 = 5 e \theta_2 = 2, vemos o seguinte:

02

Essa é uma escolha muito melhor. A gente ja tentou uma busca em grid aqui, que consiste em testar sistematicamente todos os valores possiveis. Mas isso é ruim, porque a gente acaba por computar a verossimilhança para um monte de conjunto de parâmetros, valores de \theta que a gente sabe mais ou menos que vão ser ruim.

Ai que entre o gradient descent, e se a toda escolha de parâmetro, você souber pra onde ir, souber como melhorar, se vai melhorar ou não mudar?

Até agora a gente tem uma hipótese.

h=\theta_1+\theta_2 \cdot X

E a gente tem como calcular se estamos escolhendo bons valores ou não para \theta

J(\theta_1,\theta_2)={1\over{2\cdot m}}\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2

Quanto melhor os valores de parâmetros, menor o valor de J. E nosso objetivo e conseguir deixar J o menor valor possível.

Lembra-se que com derivada a gente consegue traçar aquelas restas tangentes, de uma olhada aqui.
Aqui esta o pulo do gato, se ao passar de um ponto para o próximo, a derivada está ficando menos negativa, menos inclinada, a gente ta indo pra baixo, se ela esta ficando positiva, a gente ta subindo, aumentado o J, indo na direção errada.

Então usando derivadas, a gente sabe se deve aumentar ou diminuir os valores de \theta para conseguir um valor menor de J.

Como colocar isso em pratica? Precisamos calcular a derivada do custo para cada parâmetro.

Mas o algoritimo vai funcionar assim:

Passo 1:
Começamos com algum valor arbitrário para \theta_1 e \theta_2.

Passo 2:
Calculamos a derivada para cada parâmetro e atual e fazemos um update do parâmetro atual para um novo valor até não conseguir mais minimizar J(\theta_1,\theta_2).

{\theta_j }:=\theta_j - \alpha \cdot \frac{\partial}{\partial \theta_j} J(\theta_1,\theta_2)

Veja que a gente multiplica por esse \alpha que é o nosso “learning rate”, quanto maior, maior o passo, ou mais brusca a mudança no \theta, e dele ta multiplicando a nossa derivada naquele ponto.

Mas antes de disso a gente vai normalizar os dados para ficar mais bonitos os gráficos depois.
Normalizar os dados significa que vamos simplesmente centralizar eles no zero. A forma dos dados não muda. Veja só.

03

Mudamos os valores no eixos, ele tem agora o centro de X e Y no zero, mas continua a mesma forma geral. Isso implica que o intercepto é zero agora.

Mas tudo bem, o nosso algoritmo vai ficar assim no R.

gradientDescent<- function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=3)
    for(iter in 1:num_iters) {
        h <- ((X %*% theta)-y)
        temp1 <- theta[1,] - alpha * (1/m) * sum( h * X[,1] )
        temp2 <- theta[2,] - alpha * (1/m) * sum( h * X[,2] )
        theta[1,] = temp1
        theta[2,] = temp2
        J_history[iter,1] = temp1
        J_history[iter,2] = temp2
        J_history[iter,3] = computeCost(X, y, theta);
    }
    return(J_history)
}

Então recebemos os dados de X e y, os parâmetros \theta (isso se chama theta) iniciais, e um número de iterações que queremos realizar, quantos updates queremos fazer.

Dai a gente vai ver quantos dados temos e salvar em m, preparar uma matriz J_history para guardar os valores de \theta conforme vamos mudando, para olhar de pois, bem como os valores de J. Calculamos a derivada e fazemos simultaneamente o update para os dois valores de \theta e guardamos eles e fazemos isso até rodar o número num_iters de iterações que pedimos.

Vamos começar com uns valores de \theta bem ruins.

> theta<-matrix(c(-5,-5),ncol=1,nrow=2) > theta [,1] [1,] -5 [2,] -5 > computeCost(X, y, theta) [1] 29.82811

A gente aplica o algoritimo gradient descent

> iterations = 1500 > alpha = 0.01 > resultado<-gradientDescent(X, y, theta, alpha, iterations) > resultado [,1] [,2] [,3] [1,] -4.950000e+00 -4.942143785 29.24624141 [2,] -4.900500e+00 -4.884846846 28.67573437 [3,] -4.851495e+00 -4.828103778 28.11636476 . . . [1498,] -1.447077e-06 0.985122862 0.01427153 [1499,] -1.432606e-06 0.985122890 0.01427153 [1500,] -1.418280e-06 0.985122918 0.01427153

Aqui, a coluna 1 é o intercepto, a coluna 2 é a inclinação e a coluna 3 é o custo.
A gente começou em 29, e veja que no final, estava em 0.0143 o custo e não dava mais para sair disso.

Se a gente fizer uma regressão linear, que é outra forma de caçar os melhores parâmetros, como visto aqui.

lm(y~X[,2]) Call: lm(formula = y ~ X[, 2]) Coefficients: (Intercept) X[, 2] -2.238e-17 9.851e-01

Vemos praticamente a mesma inclinação e intercepto. Não confunda os número por causa da notação cientifica. O intercepto é um 0,000…2238, um numero muito próximo de zero, e a inclinação é
0.985, bem próximo do 0.985122918 que achamos com o gradient descent.

Se a gente for olhar como os valores do custo foram mudando, veremos o seguinte:

04

Ele foi diminuindo, mais rápido no começo, quando a inclinação da derivada era maior, depois mais devagar quando a inclinação da derivada era menor, até não mudar, quando a inclinação chegou próximo de zero.

Se comparar o modelo de regressão e os resultados do gradient descent, vemos que chegamos a mesma definição de parâmetros.

05

Sendo que a gente foi mudando os parâmetros de forma certeira, sempre minimizando o custo.

06

A gente pode fazer uma analogia dessas mudança como a escalada de uma montanha, o topo da montanha (menor custo) é onde a gente quer chegar, e dado que saímos de um ponto, vamos subindo sempre pra cima, nunca para baixo, até chegar ao cume, mas a gente começa indo rápido e diminuindo o passo ao longo da subida para não passar o cume.

07

A gente pode pegar e olhar, nessa montanha, como foram os nossos 200 primeiros passos, para ter uma ideia.

08

Muito legal certo? E um uso bem legal de derivada, conseguir andar certeiro até o maximum likelihood.
Tudo bem que para regressão linear, ninguém vai usar isso, mas a gente viu que isso funciona, e no caso de analises como o MDS (Multidimensional scaling), a gente ta fazendo algo nesse sentido, a gente tem uma matrix de distancia e esta tentando posicionar os pontos num plano, ou reta, ou sei la, mas estamos tentado posicionar e medindo se está ficando bom de acordo com uma função de custo, uma loss function.

Veja que é fácil modificar para outras funções(desde que sejam derivaveis) o gradient descent.
Peguemos o caso de uma regressão multiplica. Basta analisar para n valores de \theta.

gradientDescentMulti<-function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=1+nrow(theta))
    for(iter in 1:num_iters) {
        for(j in 1:length(theta)) {
            theta[j,]<-theta[j,] - alpha * (1/m) * sum( ((X %*% theta)-y) * X[,j] )
            J_history[iter,j] = theta[j,]
        }
        J_history[iter,nrow(theta)+1] = computeCost(X, y, theta);
    }
    return(J_history)
}

Veja que eu adicionei um loop ali dentro para cada valor de \theta, ao invés das duas linhas no caso anterior, mas a ideia continua a mesma. So com mais \theta.

Aqui a nossa hipótese se torna:

h=\theta_1+\theta_2 \cdot X1^i+\theta_3 \cdot X2^i

Um modelo aditivo, com uma soma a mais, iniciamos os \theta e rodamos nosso algoritmo gradient descent

> theta<-matrix(c(0,0,0),ncol=1,nrow=3) > theta [,1] [1,] 0 [2,] 0 [3,] 0 > computeCost(X, y, theta) [1] 0.4833333 > gradientDescentMulti(X, y, theta, alpha, 1500) [,1] [,2] [,3] [,4] [1,] -1.846902e-18 0.005242762 0.005947532 0.47707775 [2,] -3.745557e-18 0.010452588 0.011853113 0.47090590 [3,] -5.621950e-18 0.015629671 0.017717052 0.46481668 . . . [1498,] -1.664308e-16 0.808733080 0.863209808 0.01529028 [1499,] -1.664447e-16 0.808733320 0.863210047 0.01529028 [1500,] -1.665603e-16 0.808733558 0.863210284 0.01529028

E temos nas três primeiras colunas os parâmetros e a última sendo o valor do custo J.

> lm(y~X[,2]+X[,3]) Call: lm(formula = y ~ X[, 2] + X[, 3]) Coefficients: (Intercept) X[, 2] X[, 3] -1.730e-16 8.088e-01 8.632e-01

Comparamos com o resultado da regressão linear e vemos que ele está funcionando perfeito novamente.
Podemos ainda avaliar a interação, sem mexer em nada na função.

A hipótese fica assim:

h=\theta_1+\theta_2 \cdot X1^i+\theta_3 \cdot X2^i+ \theta_4 \cdot X1 \cdot X2

Rodando a função vemos o seguinte.

> theta [,1] [1,] 0 [2,] 0 [3,] 0 [4,] 0 > computeCost(X, y, theta) [1] 0.4833333 > gradientDescentMulti(X, y, theta, alpha, 1500) [,1] [,2] [,3] [,4] [,5] [1,] 1.170360e-18 0.005229905 0.007172164 0.004086374 0.4738325917 [2,] 3.297842e-06 0.010411918 0.014276559 0.008114284 0.4645249067 [3,] 9.813365e-06 0.015546491 0.021313839 0.012084525 0.4554062184 . . . [1498,] 2.092109e-02 0.585063765 0.773159318 0.259232115 0.0001656415 [1499,] 2.092109e-02 0.585063777 0.773159330 0.259232111 0.0001656415 [1500,] 2.092109e-02 0.585063789 0.773159342 0.259232106 0.0001656415

E denovo conseguimos os mesmos valores de parâmetros da regressão.

> lm(y~X[,2]*X[,3]) Call: lm(formula = y ~ X[, 2] * X[, 3]) Coefficients: (Intercept) X[, 2] X[, 3] X[, 2]:X[, 3] 0.02092 0.58507 0.77316 0.25923

Bem, é isso. As vezes ver uma utilidade diferente para derivadas pode tornar menos massante estudar calculo. Eu achei bem legal quando vi esse algoritimo na aula de aprendizado de máquina nesta matéria do coursera aqui. Ela é muito legal, mas infelizmente tudo nela está na linguagem do matlab e octave, mas não é difícil passar as coisas para o R, se eu consegui, qualquer um consegue.

Bem temos o script aqui em baixo, esse foi o tópico das duas primeiras semanas, se eu conseguir, vou tentar passar a parte de modelos de classificação para o R também, e quem sabe, com sorte, eu dou conta de passar a parte de rede neural que é super legal, mas façam o curso la, que é mil vezes melhor explicado que aqui e está no começo de uma turma.

Ficamos por aqui, os scripts sempre no repositório recologia e é isso, até o próximo post.

computeCost<- function(X, y, theta) {
    m <- length(y)
    J <- 0
    J = (1/(2*m)) * sum((X%*%theta-y)^2)
    return(J)
}
 
gradientDescent<- function(X, y, theta, alpha, num_iters) {
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=3)
 
    for(iter in 1:num_iters) {
 
        h <- ((X %*% theta)-y)
 
        temp1 <- theta[1,] - alpha * (1/m) * sum( h * X[,1] )
        temp2 <- theta[2,] - alpha * (1/m) * sum( h * X[,2] )
 
        theta[1,] = temp1
        theta[2,] = temp2
 
        J_history[iter,1] = temp1
        J_history[iter,2] = temp2
 
        J_history[iter,3] = computeCost(X, y, theta);
 
    }
 
    return(J_history)
}
 
gradientDescentMulti<-function(X, y, theta, alpha, num_iters) {
 
    m = length(y)
    J_history = matrix(0,nrow=num_iters,ncol=1+nrow(theta))
 
    for(iter in 1:num_iters) {
 
        for(j in 1:length(theta)) {
            theta[j,]<-theta[j,] - alpha * (1/m) * sum( ((X %*% theta)-y) * X[,j] )
            J_history[iter,j] = theta[j,]
        }
 
        J_history[iter,nrow(theta)+1] = computeCost(X, y, theta);
    }
 
    return(J_history)
}
 
#gerando  dados
X<-runif(30,1,10)
y<-rnorm(30,5+2*X)
 
#figura 1
plot(y~X,main="Modelo 1+1*X",frame=F,xlim=c(0,10),ylim=c(0,30))
abline(a=1,b=1)
for(i in 1:length(X)) {
    points(c(X[i],X[i]),c(y[i],1+1*X[i]),type="l",lty=1,col="red",lwd=2)
}
 
#figura 2
plot(y~X,main="Modelo 5+2*X",frame=F,xlim=c(0,10),ylim=c(0,30))
abline(a=5,b=2)
for(i in 1:length(X)) {
    points(c(X[i],X[i]),c(y[i],5+2*X[i]),type="l",lty=1,col="blue",lwd=2)
}
 
#figura 3
par(mfrow=c(2,1))
plot(y~X,main="Dados sem normalização",frame=F)
abline(v=mean(X),h=mean(y),lty=2)
plot(scale(y)~scale(X),main="Dados normalizados",frame=F)
abline(v=mean(scale(X)),h=mean(scale(y)),lty=2)
 
#transformando em matrix
X<-matrix(c(rep(1,length(X)),scale(X)),ncol=2)
y<-matrix(scale(y),ncol=1)
 
#gerando theta inicial
theta<-matrix(c(-5,-5),ncol=1,nrow=2)
theta
 
#custo inicial
computeCost(X, y, theta)
 
#outros parametros
iterations = 1500
alpha = 0.01
 
 
#gradient descent
 
resultado<-gradientDescent(X, y, theta, alpha, iterations)
 
resultado
lm(y~X[,2])
#figura 4
plot(resultado[1:200,3],type="b")
 
#figura 5
plot(y~X[,2],frame=F)
abline(lm(y~X[,2]),col="red",lty=2)
abline(a=resultado[1500,1],b=resultado[1500,2],col="blue",lty=3)
legend("topleft",lty=c(2,3),col=c("red","blue"),legend=c("Regressão Linear","Gradient Descent"))
 
#figura 6
plot(resultado[,1],resultado[,2],type="b")
 
a<-seq(-2,2,by=0.1)
b<-seq(-2,3,by=0.1)
J<-matrix(0,nrow=length(a),ncol=length(b))
 
for(i in 1:length(a)) {
    for(j in 1:length(b)) {
        J[i,j]<--1*computeCost(X, y, matrix(c(a[i],b[j]),ncol=1,nrow=2))
 
    }
}
 
#figura 7
persp(a,b,J,theta=45,phi = 30,xlab="Intercepto",ylab="Inclinação",zlab="Menor custo")
 
#figura 8
contour(a,b,J)
points(resultado[1:300,1],resultado[1:300,2],type="b",cex=0.3)
 
##Regressão Multipla
X1<-runif(30,1,10)
X2<-runif(30,1,10)
y<-rnorm(30,5+2*X1+2*X2)
 
X<-matrix(c(rep(1,length(X1)),scale(X1),scale(X2)),ncol=3)
y<-matrix(scale(y),ncol=1)
 
theta<-matrix(c(0,0,0),ncol=1,nrow=3)
theta
 
computeCost(X, y, theta)
 
gradientDescentMulti(X, y, theta, alpha, 1500)
lm(y~X[,2]+X[,3])
 
X1<-runif(30,1,10)
X2<-runif(30,1,10)
y<-rnorm(30,5+2*X1+2*X2+3*X1*X2)
 
X<-matrix(c(rep(1,length(X1)),scale(X1),scale(X2),scale(X1)*scale(X2)),ncol=4)
y<-matrix(scale(y),ncol=1)
 
theta<-matrix(c(0,0,0,0),ncol=1,nrow=4)
theta
 
computeCost(X, y, theta)
 
gradientDescentMulti(X, y, theta, alpha, 1500)
lm(y~X[,2]*X[,3])

Usando Offset em modelos lineares

A gente falou aqui sobre modelos com Overdispersion e aqui sobre modelos com Zeros-Inflados.

Usando a linguagem Bugs, por um lado a gente pode sofre um pouco por ser obrigado a definir o modelo, mas por outro lado, a gente ganha no melhor entendimento do que está acontecendo. Overdispersion é um exemplo, que somente depois que eu vi o modelo no Bugs que eu fui entender o que diabos era aquele número que saio no glm quando a gente usava a family=quasipoisson.

Agora ta na hora de entender o que significa offset, aquela opção que você pode usar nos modelos lineares da função lm, nos modelos gerais lineares na função glm, afinal, em quase todo lugar.

Se a gente olhar a documentação da função com ?lm, a gente vai ver a descrição o argumento offset:

“this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one are specified their sum is used. See model.offset.”

No meu inglês nórdico, podemos ler:

“Este pode ser usado para especificar um componente conhecido a priori para ser incluído como preditor linear no ajuste. Este deve ser NULL ou um vetor numérico de tamanho igual ao número de casos. Um ou mais termos de offset podem ser incluídos na formula ao invés ou tal como, e se mais de um for especificado, sua soma será usada, olhe model.offset.”

Se você olhar a função model.offset, vai cair nas funções que transformam formula em matrizes, álgebra linear e não da mais para entender nada.

Mas vamos pensar de maneira mais simples. No Glm de Poisson, como estamos olhando, nos assumimos que as contagens esperadas são adequadamente descritas pelos efeitos das covariáveis. Porém, frequentemente, nos temos algo como uma “janela de contagem” que não é constante. Por exemplo, áreas de estudos que não tem o mesmo tamanho ou amostras que tem duração diferentes, a gente pode não determinar antes o tempo que vai olhar uma lagoa, simplesmente vamos olhando e contando quantos sapinhos vemos la. Para adicionar essa variação no modelo de Poisson, por exemplo se temos contagens de coelhos em vários campos de tamanho diferentes, nos podemos definir a área como offset, mais especificamente, como o modelo de Poisson a gente usa função de ligação logarítmica, a gente define o log da área de cada amostras. Simplificando, a gente modela a densidade como a resposta, e não a contagem per se.

Vamos fazer nosso exemplo aqui, vamos pensar em coelhinhos. Estamos contando coelhinhos em pastos e plantações, vamos fazer varias amostras, mais ou menos assim:

01

Cada pontinho é uma amostra. Mas cada uma dessas amostras são de areas cercados, com curva de nível, varias coisas no campo que dividem elas, então não são contínuos iguais, temos tamanhos diferentes para cada amostras. Podemos representar esses tamanhos aqui com o tamanho da bolinhas sendo equivalente ao tamanho da amostra.

02

No modelo de Poisson a gente tem:

C_i \sim Poisson(\lambda_i)

Cada amostra vem de uma distribuição de Poisson com média \lambda.
Agora como cada amostra foi feita num pequeno pasto de área A, nosso modelo então é na verdade:

C_i \sim Poisson(A_i \cdot \lambda_i)

Dessa forma, o preditor linear fica assim:

log(A_i \cdot \lambda_i)

E pelas regras de logaritmo, a gente pode separar essa multiplicação dentro do log em dois logs.

log(A_i) + log(\lambda_i)

Agora esse \lambda pode vir de um preditor linear, por exemplo com a variável x sendo o tipo de pasto que coletamos. Então

log(A_i \cdot \lambda_i) = log(A_i) + \alpha + \beta \cdot x_i

Veja que ficou assim porque podemos fazer aquela separação da área e do \lambda, por causa da multiplicação no log.

De forma mais simples, estamos pensando nesse log da área com um coeficiente 1 multiplicando ele, ou seja, a área fazendo o efeito dela proporcionalmente, mas isso pode ser mudado, usando um coeficiente para a área.

log(A_i \cdot \lambda_i) = \beta_0 \cdot log(A_i) + \alpha + \beta \cdot x_i

Agora é so passar isso para a linguagem Bugs da seguinte forma.

 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }

São duas linhas, a linha que tem a ~ é a parte estocástica, e a linha com <- é a parte determinística, que está sendo somado, além do preditor linear, o offset da área, ja logaritimizado, a gente não ta mandando pro bugs a área real, a gente já calcula o log no R e só mando o valor do log. Então é juntar os dados, configurar o MCMC, fazer a função de inicialização das cadeias e rodar o modelo. E para frisar, veja que a gente ta mandando o termo logA que é o log(A).

# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))

Rodamos o modelo, mesmo com cadeias pequenas ele converge rápido, é um modelo relativamente simples e vemos o seguinte.

> print(out, dig = 3) Inference for Bugs model at "Offset.txt", Current: 3 chains, each with 1100 iterations (first 100 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.579 0.126 0.323 0.493 0.583 0.669 0.809 1.003 3000 beta 0.966 0.148 0.685 0.861 0.966 1.064 1.255 1.001 3000 deviance 97.335 1.974 95.430 95.960 96.745 98.080 102.902 1.001 3000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 2.0 and DIC = 99.3 DIC is an estimate of expected predictive error (lower deviance is better).

Todos os parâmetros com um Rhat perto de 1, e vemos a diferença entre os dois tipos de locais que observamos.
Lembrando que beta aqui é a diferença entre os dois locais, isso pela forma como escrevemos o preditor linear, de uma olhada aqui se não entendeu essa parte. Mas vemos que essa diferença tem média 0.966, mas mais importante que isso o intervalo de confiança de 95% dessa distribuição vai de 0.685 a 1.255, ou seja , não engloba 0 (nenhuma) diferença, alias, com 0.148 de desvio padrão para esse parâmetro, o valor de zero diferença está muitos desvios padrão de distancia da média, temos muita certeza da diferença.

Agora a titulo de curiosidade vamos explorar aqui o que acontece ao ignorarmos o offset, ou seja construir o modelo de forma errada.

> glm.fit.sem.offset <- glm(C ~ x, family = poisson) > summary(glm.fit.sem.offset) Call: glm(formula = C ~ x, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -1.4298 -0.7816 -0.1419 0.8385 1.7360 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.8245 0.1270 14.367 < 2e-16 *** xPlantação 0.9355 0.1499 6.242 4.31e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 60.948 on 19 degrees of freedom Residual deviance: 17.615 on 18 degrees of freedom AIC: 103.45 Number of Fisher Scoring iterations: 4

Veja que no modelo bayesiano com offset, o alpha ou intercepto, como quiser chamar, ficou 0.579, sem o offset aqui ele ficou como 1.8245, ou seja muito maior, super-estimamos muito a contagem de coelhos, esse valor vai dar quantos coelhos esperamos achar no pasto e nas plantações, lembrando que nas plantações o lambda vai ser alpha mais beta, ou aqui intercept + xPlantação, então por mais que a diferença até que fica parecida, temos os valores de lambda bem super-estimado.

Mas não confunda com o problema ser modelo bayesiano ou não, se ajustarmos o modelo frequentista temos praticamente os mesmo parametros do modelo bayesiana, com o mesmo erro.

> glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A)) > summary(glm.fit.com.offset) Call: glm(formula = C ~ x, family = poisson, offset = log(A)) Deviance Residuals: Min 1Q Median 3Q Max -1.68068 -0.44061 -0.06695 0.31836 1.69544 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.5811 0.1270 4.575 4.75e-06 *** xPlantação 0.9658 0.1499 6.445 1.16e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 59.836 on 19 degrees of freedom Residual deviance: 13.535 on 18 degrees of freedom AIC: 99.376 Number of Fisher Scoring iterations: 4

Apesar da grande diferença nas estimativas dos parâmetros, o uso ou não de offset, não é facil detectar que ele esta faltando, temos que ter pensado isso a priori. Se a gente olhar os resíduos de ambos os modelos não vemos nenhuma tendencia, pelo menos eu não vejo.

03

Temos um desvio residual menor no modelo empregando o offset, 13.535 com offset contra 17.615 sem offset, mas a gente nunca saberia disso se não fizer o ajuste com o offset, a novamente, não pegaria a necessidade disso olhando somente um plot dos resíduos sem offset.

Bem, é isso, mais um conceito legal para se entender, e facilitar a vida ao ler artigos e quem saber resolver os problemas de alguém. O repositório recologia esta la. E é isso, até o próximo post.

Referência:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

### Gerando dados
set.seed(1234)
n.locais <- 10
A <- runif(n = 2*n.locais, 2,5)		# Area de 2 a 5 km2 
x <- gl(n = 2, k = n.locais, labels = c("Pasto", "Plantação"))
preditor.linear <- log(A) + 0.69 +(0.92*(as.numeric(x)-1))
lambda <- exp(preditor.linear)
C <- rpois(n = 2*n.locais, lambda = lambda)
 
#Figura 1
stripchart(C~x,vertical=T,pch=19,method="stack",offset=0.8,at=c(1.2,2.2),frame=F,xlab="Contagem",ylab="Local")
 
#Figura 2
stripchart(C~x,vertical=T,pch=19,method="jitter",offset=0.8,at=c(1,2),frame=F,xlab="Contagem",ylab="Local",cex=A/2)
 
### Análise com OpenBugs
# Definindo o modelo
sink("Offset.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- 1 * logA[i] + alpha + beta *x[i]
 }
}
",fill=TRUE)
sink()
 
# Juntando os dados numa lista para o OpenBugs
bugs.data <- list(C = C, x = as.numeric(x)-1, logA = log(A), n = length(x))
 
# Função para iniciar as cadeias com valores ao acaso
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros a estimar
params <- c("alpha","beta")
 
# MCMC settings
nc <- 3    # Número de cadeias
ni <- 1100 # Tamanho de cada cadeia
nb <- 100  # Número de amostras a descartar como burn-in
nt <- 2    # Thinning rate
 
# Inicie o Gibbs-sampler
library(R2OpenBUGS)
 
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params, model.file="Offset.txt",
            n.thin=nt,n.chains=nc,n.burnin=nb, n.iter=ni)
 
#Resultado
print(out, dig = 3)
 
### Analise usando R
 
glm.fit.sem.offset <- glm(C ~ x, family = poisson)
glm.fit.com.offset <- glm(C ~ x, family = poisson, offset = log(A))
summary(glm.fit.sem.offset)
summary(glm.fit.com.offset)
 
#figura 3
par(mfrow=c(2,1))
plot(resid(glm.fit.sem.offset),ylim=c(-2.5,2),main="Sem Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)
plot(resid(glm.fit.com.offset),ylim=c(-2.5,2),main="Com Offset",pch=19,ylab="Resíduos",xlab="Amostras",cex=1.3)