Rosalind – Introduction to Pattern Matching

Bem, é chegada a hora de voltar a ver os problemas do Rosalind. Faz um bom tempo que não resolvia nenhum exercício de la.
Então vamos começar com um fácil, ou não.

No exercício Introduction to Pattern Matching, a gente tem que implementar uma arvore de prefixos ou trie para as sequências.
Basicamente essa é uma arvore que permite guardar dados sem ter que repetir prefixos, e tem vários usos, por exemplo, se a gente falar um prefixo, a partir dela é fácil achar quais palavras tem aquele prefixo.

Então a diferença aqui é que ao invés de um alfabeto completo, a gente vai ter apenas o alfabeto das sequencias de dna, que é o ATCG.

A entrada então é uma serie de sequencias

ATAGA ATC GAT

E devemos retornar a árvore trie formado por ela, falando as arestas e qual a letra que ela liga na árvore, sengo que vamos anotar as arestas na ordem que são escritas.

1 2 A 2 3 T 3 4 A 4 5 G 5 6 A 3 7 C 1 8 G 8 9 A 9 10 T

Então a gente vai ser por exemplo, em python, apos um parse, uma lista de strings que são as cadeias

['ATAGA', 'ATC', 'GAT']

E a gente vai ter que construir a arvore trie. Quando eu estudei isso com minha querida professora de estrutura de dados, a gente implementava tudo na linguagem C, usando ponteiros, que de certa forma facilitava as coisas. Isso é possível de ser feito em Python também, usando classes, já que tudo é um ponteiro para objetos, assim como Java, eu acho, pelo menos aqui tem um exemplo legal pelo Robert Sedgewick, que aparece no curso dele de algorítimos do coursera.

Mas lendo um pouco eu vi essa implementação bem legal usando lista e dicionários. Então ao invés de usar ponteiros, a gente usa dicionários.

Bem primeiro a gente cria uma classe em python, para segurar essa estrutura

1
2
3
4
class Trie(object):
    def __init__(self):
        self.contagem = 1
        self.raiz = [self.contagem, {}]

Veja que atributo contagem é para manter a contagem dos nós que estão sendo criados, a gente começa em 1 que é a raiz, dai vamos incrementar ao inserir um nó, e como não temos nada, começamos com um dicionario vazio.
Em python, um dicionário vai ser uma estrutura que a gente cria com essa chaves, que vai ter uma chave e um valor para essa chave, separado por :.
Por exemplo:

{'A':adenina, 'T':timina, 'C':citosina, 'G': guanina}

A chave tem que ser um identificador, como uma letra ou número, pra gente vai ser as letras, e o valor associado com a chave vai ser…

Outra chave, ou o próximo nível da árvore. Legal né.
Poderíamos só usar o dicionario ser não tivéssemos que manter o registro da entrada de um novo nó na árvore, mas como temos que fazer isso, vamos usar uma lista, com o primeiro elemento da lista sendo o valor da ordem, e o segundo o dicionário.

1
2
3
4
5
6
7
    def inserir(self, sequencia):
        node = self.raiz
        for base in sequencia:
            if base not in node[1]:
                self.contagem = self.contagem +1
                node[1][base] = [self.contagem, {}]
            node = node[1][base]

Então agora vamos ver como vamos inserir algo.
Para inserir, criamos a função inserir, que recebe a sequência que desejamos inserir.
Primeiro pegamos a lista que é a raiz, veja que essa lista tem dois elementos, no elemento 0 temos a contagem que não vamos mexer, e no elemento 1 o dicionário.
Dai para cada base na sequência, usando o loop do for, vamos ver se a base é uma chave na lista, veja que usamos o not in node[1], que é o dicionário, então ele ta comparando essa base com o conjunto de chaves, se ele não encontrar, ele aumenta 1 a contagem, e cria uma nova lista, que é um nível abaixo da na árvore, e la adiciona uma lista com a contagem e um dicionário vazio, mas ele adiciona isso no dicionário da raiz, usando a base como chave, e pega esse nova lista feita como raiz para a próxima base.
Se base já estiver no dicionário, a gente não entra nesse if e simplesmente vamos para o dicionário que tem como chave a base em questão

Então para a primeira letra, no inicio a raiz vai ser

[1,{}]

E vai entrar um A, então isso vai ser adicionado ali no dicionario

[1,{A:[2,{}]}]

E assim sucessivamente.

Adicionando todas as sequências de exemplos a gente fica com o seguinte:

[1, {'A': [2, {'T': [3, {'A': [4, {'G': [5, {'A': [6, {}]}]}], 'C': [7, {}]}]}], 'G': [8, {'A': [9, {'T': [10, {}]}]}]}]

E dai para frente, é só formatar a saída, que basicamente é como imprimir a árvore em pré ordem, assim como é feito em árvores primarias. Imprimindo a aresta, e qual a base que gerou ela.

1 2 A 2 3 T 3 4 A 4 5 G 5 6 A 3 7 C 1 8 G 8 9 A 9 10 T

As vezes é meio difícil visualizar essa árvore, mas se a gente olhar o dicionário com uma indentação, da para começar a se ligar no que esta acontecendo.

[1, {'A': [2, {'T': [3, {'A': [4, {'G': [5, {'A': [6, {}] } ] } ], 'C': [7, {} ] } ] } ], 'G': [8, {'A': [9, {'T': [10, {} ]} ]} ]} ]

E para que nunca ouviu falar de arvore TRIE, basicamente a gente está escrevendo nos dados a seguinte árvore.

01

Fiz essa representação usando o pacote igraph, so para ajudar a visualizar, segue o código em R.

1
2
3
4
5
6
library(igraph)
exemplo<-graph.formula( 1--2,2--3,3--4,4--5,5--6,3--7,1--8,8--9,9--10)
V(exemplo)$name<-c("Raiz","A","T","A","G","A","C","G","A","T")
local<-structure(c(215, 184, 187, 190, 191, 195, 148, 242, 243, 246, 
236, 191, 145, 91, 43, 0, 97, 194, 146, 94), .Dim = c(10L, 2L))
plot(exemplo,layout=local)

Lembrando que para criar o local, primeiro eu so plotei o grafo usando o tkplot, depois que eu arrumei os vértices, eu peguei as coordenadas com a função get.coordinates(), dei um dput na saída dela para ter algo que eu pudesse escrever num objeto para colocar aqui. Assim eu pude omitir a parte do tkplot e ainda deixar meu grafo reproduzível.

Bem é isso ai, o script vai estar la no repositório recologia, ele so insere, porque não precisamos remover nem modificar nada para esse exercício, acho que ele vai ser mais de introdução para comparação de sequencias, que não não cheguei em exercícios de alinhamento no rosalind. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
class Trie(object):
    def __init__(self):
        self.contagem = 1
        self.raiz = [self.contagem, {}]
 
    def inserir(self, sequencia):
        node = self.raiz
        for base in sequencia:
            if base not in node[1]:
                self.contagem = self.contagem +1
                node[1][base] = [self.contagem, {}]
            node = node[1][base]
 
    def imprime(self):
        print self.raiz
 
def formata_recursivo(node):
    for base, node2 in node[1].iteritems():
        print node[0], node2[0], base
        formata_recursivo(node2)
 
if __name__ == '__main__':
    entrada = open("trie.in").readlines()
    dados=[]
 
    for linha in entrada:
        dados.append(linha.strip())
 
    print dados
 
    trie_sequencias = Trie()
 
    for sequencia in dados:
        trie_sequencias.inserir(sequencia)
 
    trie_sequencias.imprime()
 
    formata_recursivo(trie_sequencias.raiz)

Evolução – Fisherian Optimality Models – Parte 3

E vamos continuar olhando os modelos ótimos baseados na ideia de Fisher. A gente já viu dois exemplos aqui, na parte 1 e parte 2, agora vamos ver um caso um pouco mais complicado, do ponto de vista do cálculo pelo menos, ao mudar somente um pressuposto matemática em relação a parte 2, mais especificamente o pressuposto 3.

Pressupostos gerais:

1. O organismo é iteroparo
2. Fecundidade, F, aumenta com o tamanho do corpo x, que não muda depois da maturidade atingida.
3. Sobrevivência, S, diminui com o tamanho do corpo x.
4. O fitness, W, é uma função da fecundidade e da sobrevivência.

Pressupostos matemáticos:

1. Maturidade ocorre na idade 1 depois da qual não há mais crescimento.

2. A Fecundidade aumenta linearmente com o tamanho na maturidade, resultando na fecundidade sendo uma função uniforme da idade, t:

F_t=a_F+b_{F}x

3. A taxa instantânea de mortalidade aumenta linearmente com o tamanho do corpo obtido na idade 1 e é constante por unidade de tempo. A partir desse pressuposto, a sobrevivência na idade t é dada por

S_t = e^{-(a_S+b_Sx)t}

Note primeiro que a sobrevivência não é mais independente do tamanho do corpo e que para fazer a sobrevivência uma função que diminui em relação ao tamanho do corpo, dado que a função é exponencial, nos alteramos o a_S-b_Sx com a_S+b_Sx.

4. Fitness, W, é o sucesso reprodutivo esperado ao longo da vida, R_0, dado como o produto acumulado da sobrevivência e da fecundidade

W = R_0 = \sum\limits_{t=1}^\infty F_{t} S_t = (a_F + b_Fx)e^{-(a_S+b_Sx)t}

Agora a gente não pode retirar do somatório os efeitos dependentes da idade e assim o tamanho ótimo para o corpo não será como o visto na parte 2.

Gráfico da função de fitness

Antes de tentar encontrar um valor ótimo, nos primeiro plotamos a função do fitness para verificar que há um ponto máximo, que a função não tem uma forma esquisita de forma que faça as rotinas para encontrar o ponto máximo falharem, independente do ponto de inicio. Nos vamos assumir os mesmo parâmetros do exemplo da parte 2 (a_F=0, b_F=4,a_S=1, b_S=0.5), assim a função

W = \sum\limits_{t=1}^\infty 4xe^{-(1+0.5x)t}

O somatório acima vai até o infinito, o que geralmente não é uma opção para a análise numérica. Assim, primeiro, nos temos que decidir para quantas gerações, n, nos precisamos considerar esse somatório. Para isso nos setamos o tamanho do corpo, x, em algum valor arbitrário mas razoável, vamos dizer x=1, e fazemos o seguinte:

Geramos uma sequencia de inteiros de 1 a n e atribuímos isso a um vetor chamado idade.

A partir dai nos criamos outro vetor chamado
Wt, que é o componente especifico da idade da equação, 4xe^{-(1+0.5x)t}

Por último nos calculamos a soma do vetor Wt usando a função sum.

1
2
3
4
5
6
Somatorio <- function(n,x=1) {
    #print(x)
    Idade <- seq(from=1, to=n)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}

Então para o peso = 1, temos a seguinte figura usando este código:

1
2
3
4
nmax <- 20
n <- matrix(seq(from=1, to=nmax))
W <- apply(n,1,Somatorio)
plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3)

01

O somatório rapidamente se aproxima de sua assintota, seu valor assintótico, assim se fizer uma somatório até a idade 20, esse valor deve ser adequado para todos os valores de x.

Apenas para ilustrar um pouco mais, o que aconteceria se escolhêssemos outros valores de x, vamos fazer mais algumas tentativas.

Usamos o seguinte código:

1
2
3
4
5
6
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
for(i in 1:4) {
    W <- apply(n,1,function(n) {Somatorio(n,x=i)})
    plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3,
         main=paste("Tamanho do corpo",i),frame=F)
}

Um pequeno detalhe, se alguém ficou curioso o porque do print na função somatório, quando eu fiz esse for, eu fiquei preocupado com o escopo do código, se ele ia pegar o i certo, dai imprimi o i para garantir, veja que dentro do for o R pegou o valor de i que estava neste escopo e o n do apply.

02

Apesar de diferente o valor da assintota, a gente sempre tende pra ela. Por curiosidade, a gente pode se perguntar, porque diabos isso acontece? Bem basta olhar que a função que está no pressuposto 3, veja que a mortalidade aumenta com o tempo, ou seja a sobrevivência diminui ao longo do tempo, então a cada tempo, a cada tempo que se passa, a contribuição da idade para a sobrevivência é menor, um ano a mais para alguém “velho” não faz mas muita diferença, porque ele já vai morrer e não vai mais contribuir para mais descendentes. Podemos olhar como a função que estamos somando, como ela é para cada tempo.

03

Então conseguimos nos convencer que não é preciso somar a função de fitness até o infinito, até o 20 ta bom, ja que a cada novo valor do somatório, a gente soma um numero menorzinho, menorzinho, cada vez mais insignificante, então procedemos para tanto.

04

E olhando assim, da para ver que temos um ponto máximo para o tamanho do corpo, agora é só encontrá-lo.

Encontrando o máximo usando cálculo.

Agora que nos temos certeza que há um fitness máximo, nos sabemos seu valor aproximado pelo olhômetro da figura anterior e nos podemos procurar esse valor. O somatório nesse caso é resolvível, e para tanto nos vamos tentar chegar a solução exata usando cálculo.
Uma série que ocorre frequentemente em modelos de historia de vida, como o que temos aqui, são as séries geométricas.

\sum\limits_{t=1}^\infty = 1+a+a^2+a^3+ \dots = \frac{1}{1-a}

onde \left| a \right| < 1[/latex] (valor absoluto é menor que 1). No nosso caso, a função de mortalidade é uma série geométrica. Para simplificar, vamos estabelecer que [latex]a_F + b_F x = A[/latex] e [latex]a_S+b_S x = B[/latex], o que nos da [latex]W= \sum\limits_{t=1}^\infty F_tS_t = \sum\limits_{t=1}^\infty (a_F+b_Fx)e^{(a_S+b_Sx)t} = \sum\limits_{t=1}^\infty (A)e^{(B)t}[/latex] E vendo dessa forma, como é uma multiplicação ali dentro do somatório, podemos arrancar ele fora, como fizemos na parte 2 [latex]W = A \sum\limits_{t=1}^\infty e^{(B)t}[/latex] Agora para mudar essa equação, para a formula fechada do série geométrica, note que [latex]e^{-Bt}[/latex] pode ser escrito como [latex]e^Be^{B(t-1)}[/latex] e assim podemos arrancar esse cara também para fora do somatório [latex]W = A \sum\limits_{t=1}^\infty e^{(B)t} = Ae^{-B} \sum\limits_{t=1}^\infty e^{-B(t-1)}[/latex] E agora mudamos o somatório para a formula fechada [latex]W = Ae^{-B} (\sum\limits_{t=1}^\infty e^{-B(t-1)}) = Ae^{-B} (\frac{Ae^{-B}}{(1-e^{-B})})[/latex] Agora nos temos uma função "simples", simples no sentido que ela é diferenciável, porque ta complicado esse monte de equação. Mas é possível seguir em frente e diferenciar para achar o 0 da derivada. Podemos usar a regra da cadeia para encontrar a derivada, que fica. [latex]\frac{dW}{dx} = \frac{(b_F)e^{-(a_S+b_Sx)}}{1-e^{-(a_S+b_S x)}} + \frac{(a_F+b_F x)(-b_S e^{-(a_S+b_S x)})}{1-e^{-(a_S+b_S x)}} + \frac{(a_F+b_Fx)e^{-(a_s+b_Sx)}(-1)(b_S e^{-(a_S+b_Sx)})}{(1-e^{-(a_S+b_S x)})^2} [/latex] Da para destrinchar passo a passo, mas vamos deixar isso para outro post. Agora nos precisamos achar o valor de x no qual [latex]\frac{dW}{dx}=0[/latex], de forma que podemos simplificar essa equação final da derivada um pouco já que [latex]\frac{e^{-(a_S+b_Sx)} }{ 1-e^{-(a_s+b_Sx)}}[/latex] é um elemento comum aos três termos que temos, assim a equação que temos que resolver é [latex](b_F) + (a_F+b_Fx)(-b_S) + \frac{(a_F+b_Fx)(-1)(b_Se^{a_S+b_Sx})}{1-e^{a_S+b_Sx}} = 0[/latex] E para resolver isso no R podemos usar a função uniroot. Agora é simples, é só escrevermos essa função no R, usando os valores que definimos para os parâmetros e chamar a função uniroot. Então definimos a função

1
2
3
4
5
funcao <- function(x){
    (4)+
        (0+4*x)*(-0.5)+
            (0+4*x)*(-1)*(0.5*exp(-(1+0.5*x)))/ (1-exp(-(1+0.5*x)))
}
> solucao <- uniroot(funcao, interval=c(0,4)) > solucao$root [1] 1.682812

Como nos estamos apenas interessado nos valores positivos, pesos negativos não fariam nenhum sentido biológico, o limite para a busca da raiz é colocado como 0, excluindo possíveis raízes negativas, Se houver duas raízes positivas, a raiz mais próxima do menor limite sera encontrada e devolvida pela função. Nos sabemos que para esse caso, a partir da figura do fitness, que existe apenas uma raiz positiva, e assim não precisamos procurar por mais valores.

Aproximação numérica

Agora, vamos dar uma olhada como seria para encontrar o valor máximo a partir de uma aproximação numérica.

Nos vamos usar o mesmo somatório anterior, exceto pelo fato que vamos tomar o resultado como um valor negativo, porque a função nlm do R acha o mínimo de uma função, e ao colocar o resultado negativo, nos deixamos a função de ponta cabeça, mas ainda acharemos o valor que nos interessa. Precisamos de uma valor inicial também para iniciar a busca, e setamos esse valor como 1, já que sabemos, no olhômetro do gráfico, que esse é um valor razoável.

Então alterando a função para retornar menos o somatório

1
2
3
4
5
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(-sum(Wt))
}

E encontrando a estimativa de valor ótimo para o peso

> nlm(Somatorio, p=1) $minimum [1] -1.268755 $estimate [1] 1.68281 $gradient [1] -2.111179e-09 $code [1] 1 $iterations [1] 6

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, principalmente nas equações, que é fácil passar error, ou mande um e-mail e até mais. Olhando a solução por métodos numéricos, realmente da uma tentação em ignorar o uso do cálculo, mas não custa nada olhar como faz. O livro do Roff de evolução é realmente muito legal e não terminamos nem o inicio, mas é devagar e sempre.

Referência:
Derek A. Roff 2010 – Modeling Evolution – an introduction to numerical methods. Oxford University Press 352pp.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
##################################
## Figuras
##################################
Somatorio <- function(n,x=1) {
    #print(x)
    Idade <- seq(from=1, to=n)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}
 
nmax <- 20
n <- matrix(seq(from=1, to=nmax))
W <- apply(n,1,Somatorio)
 
#Figura 1
plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3)
 
 
#Figura 2
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
for(i in 1:4) {
    W <- apply(n,1,function(n) {Somatorio(n,x=i)})
    plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3,
         main=paste("Tamanho do corpo",i),frame=F)
}
 
#Figura 3
curve(4*1*exp(-(1+0.5*1)*x),0,2,frame=F)
 
#
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}
 
 
x<- matrix(seq(from=0,to=5,length=100))
W <- apply(x,1,Somatorio)
 
#Figura 4
plot(x,W,type="l", xlab="Tamanho do corpo, x", ylab="Fitness, W",las=1,
     lwd=2,lty=2,col="red",frame=F)
 
 
##################################
## Usando Calculo
##################################
funcao <- function(x){
    (4)+
        (0+4*x)*(-0.5)+
            (0+4*x)*(-1)*(0.5*exp(-(1+0.5*x)))/ (1-exp(-(1+0.5*x)))
}
 
#Raiz
solucao <- uniroot(funcao, interval=c(0,4))
solucao$root
 
##################################
## Usando aproximação numérica
##################################
 
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(-sum(Wt))
}
 
nlm(Somatorio, p=1)

Gráficos usando pacote Lattice

Opa, após um tempinho sem postagens, estou de volta. Essas últimas semanas foram meio corridas, muita coisa para fazer e não sobrou muito tempo. Mas nunca podemos parar de estudar, então aqui estamos novamente.

Recentemente eu tive que fazer alguns gráficos, de variáveis com mais de um fator, e me vi apanhando para conseguir implementar tudo que queria na figura. Eu tenho tentado aprender um pouco mais sobre ggplot2, mas ainda não sai do sistema padrão de gráficos do R.

Bem, vamos ver um exemplo aqui, e o que da para fazer.

Vamos supor que tenhamos medidas do tamanho e massa de indivíduos de uma espécie de passarinho e queremos trabalhar com a alometria dessas medidas. Mas temos essas medidas para quatro populações distintas.

Os dados são mais ou menos assim:

> head(dados) pop tamanho massa 1 1 -0.9132259 248.1997 2 1 -0.5202886 323.5949 3 1 0.2195127 300.5737 4 1 1.4554830 311.6083 5 1 -1.1484635 253.0550 6 1 1.4192977 303.5195

Então queremos ver a relação entre massa e tamanho, logo a primeira ideia que vem a cabeça é fazer um gráfico de pontos, um scatterplot entre tamanho e medida.

Como os dados estão organizadinhos como tidy-data, isso fica fácil, é usar a função plot

1
plot(massa~tamanho,data=dados)

E temos a figura.
Figura 01
01

Mas olhando os dados ali em cima, sabemos que são quatro populações diferentes, logo olhar tudo junto é no mínimo perigoso, então precisamos de algum jeito de separar elas, bem a primeira coisa que vem na minha cabeça é que podemos usar cores diferentes, ou símbolos. Como não somos limitados ao preto e branco para a simples visualização, colocamos quatro cores diferentes

1
2
plot(massa~tamanho,pch=19,frame=F,xlab="Tamanho",ylab="Massa",col=dados$pop,data=dados)
legend("topleft",pch=19,col=1:4,legend=paste("População",1:4),bty="n")

Mas como ninguém é adivinho, vamos usamos a função legend para fazer uma legenda, para ter uma associação com qual cor representa qual população, isso é importante, porque para quem sabe da onde vem os dados, isso pode trazer muito sentido, além disso podemos caprichar colocando nomes melhores nos eixos com xlab e ylab e colocar um bolinha fechada mudando o símbolo com o pch, como as populações são representadas com números aqui, usamos o col para mudar ela, apenas apontando ele para a linha do data.frame, mas se fossem palavras, o mais rápido seria usar as.numeric no vetor com os nomes.

Figura 02
02

Bem a figura ja fica mais bonitinha, mas nesse caso onde estamos colocando a massa no eixo y, temos como preditores o tamanho e a população, 2 preditores, isso pode começar a ficar confuso, então podemos na verdade quebrar a figura em quadros, ai que entra o tal das figuras do tipo lattice, o mais simples de fazer isso é alterando o argumento do par mfrow, que vai construir os gráficos por linha, ou mfcol para construir por coluna, acho que o mfrow é o mais intuitivo de usar, e mais comum nos exemplos, mas uma outra possibilidade, que eu acho interessante é usar a função layout.

No layout, a gente faz uma matriz, por exemplo:

> matrix(1:4,ncol=2,nrow=2,byrow=T) [,1] [,2] [1,] 1 2 [2,] 3 4

E isso quer dizer que vamos dividir a tela em 4 partes (número de elementos da matriz), a primeira parte vai ser do primeiro gráfico que fizermos, a segunda do segundo, e assim sucessivamente.

Então usamos informamos essa matriz para o layout:

1
layout(matrix(1:4,ncol=2,nrow=2,byrow=T))

E podemos fazer nossa figura

1
2
3
4
5
6
7
8
9
10
for(i in 1:4) {
    plot(massa~tamanho,pch=19,frame=F,xlab="Tamanho",ylab="Massa",
         col=i,data=dados[dados$pop==i,],main=paste("População",i),
         xlim=c(min(dados$tamanho),max(dados$tamanho)),
         ylim=c(min(dados$massa),max(dados$massa)) )
    ordem <- order(dados[dados$pop==i,"tamanho"])
    ajuste<-loess(massa~tamanho,dados[dados$pop==i,])
    points(ajuste$fitted[ordem]~dados[dados$pop==i,"tamanho"][ordem],
           type="l",col="red",lty=3,lwd=2)
}

Agora usamos um loop para fazer a figura, porque temos que fazer 4 gráficos, e não um so. Vamos adicionar uma linha de tendencia também usando a função loess.

Figura 03
03

E o resultado ja fica bem legal, agora so para ilustrar como o layout pode fazer coisas legais, veja que a gente pode definir a matrix de layout assim

> matrix(c(1:3,rep(4,3)),ncol=3,nrow=2,byrow=T) [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 4 4

Isso quer dizer que usaremos a parte de cima da tela para os três primeiros gráficos, e a parte de baixo para o quarto, talvez não o melhor resultado, mas a ideia é que podemos repetir números para dar mais espaço a algum gráfico em detrimento do outro, basta repetir ele, usando uma matriz de quantas linhas e colunas quisermos, basta colocar os números que queremos repetindo nos espaços que devem ser dedicados ao gráfico.

Figura 04
04

Mas de modo geral, estamos reinventando a roda, como é comum a necessidade de fazer figuras desse tipo, e pode se tornar bem complicado e longo o código, levando ao erro. Mas o pacote lattice está ai para resolver, a ideia dele é ele fazer o layout e fazer os gráficos por separação de um segundo preditor.

1
2
library(lattice)
xyplot(massa~tamanho|pop,data=dados)

Então é so abrir o pacote lattice, que vem junto com o R por padrão, e separar o segundo fator por uma barrinha reta “|”, se tivermos um terceiro fator preditor, tipo espécie, será colocar “|pop+ espécie”. Mas simples assim.

Figura 05
05

Agora o que achei mais complicado no começo, foi como mexer em características dentro do plot, veja que la em cima mudávamos cor, símbolo, etc, adicionávamos linha de tendência.

Veja que coisas como

1
xyplot(massa~tamanho|pop,data=dados,pch=19,col=dados$pop)

e

1
xyplot(massa~tamanho|pop,data=dados,pch=19,col=1:4)

Vão falhar.

Figura 06
06

Bem uma diferença do lattice para o sistema de gráfico padrão do R é que no sistema padrão os gráficos a gente vai fazendo as coisas iterativamente, primeiro a gente faz o gráfico com bolinhas, depois adiciona linha de tendencia, depois faz uma legenda, depois adiciona mais pontos.

O lattice, assim como o ggplot2, é pa pum, você tem que planejar tudo, fazer o comando e mandar. Dai ele tem padrões muito legais, mas por exemplo para mudar a cor, uma possibilidade é a seguinte.

1
2
3
meu.panel = function(x,y,subscripts) {
    panel.xyplot(x,y,col=dados$pop[subscripts],pch=19)
}

Você tem que fazer sua própria função panel, a função panel é a que vai fazer a figura dentro de cada painel. Essa é a forma mais direta de você mudar o que bem entender dentro de cada painel. Existe vários panels prontos, e até argumentos que mandamos no plot que são passados para os panel padrão, mas fazer uma função panel nem é o fim do mundo. Basicamente fazemos o mesmo que fizemos dentro do loop que fizemos la em cima, mas o lattice tem copias de quase todas as funções do plot, então ao invés de usar plot, usamos panel.xyplot, e veja que o panel está recebendo três argumentos, os valores de x, de y e o subscript referentes aos dados originais, assim podemos pegar qual cor ele deve ser do vetor dados$pop.
Existem panel.points, panel.abline, acho que panel.qualquer coisa que funciona no plot, só pesquisar no help do lattice. Mas então definimos nossa função de panel e mandamos ela para o xyplot.

1
xyplot(massa~tamanho|pop,data=dados,panel=meu.panel)

Figura 07
07

Outra coisa, ao invés de usar loess, usamos panel.loess, ta tudo prontinho, na mão. Sendo que o loess já faz a linha direto, não precisamos nem ajustar o modelo e predizer os pontos com ele.

1
2
3
4
meu.panel = function(x,y,subscripts) {
    panel.xyplot(x,y,col=dados$pop[subscripts],pch=19)
    panel.loess(x, y,col="red",lty=2,lwd=2)
}

Fora isso, outras alterações comuns são querer mudar a cor, estilo do quadradinho que fica escrito o nome dos paneis, bem eu coloquei preto e branco uma vez para uma publicação, e sempre reutilizo o mesmo set, já que quase nunca quero mudar muito mais coisas mesmo.

1
trellis.par.set(strip.background=list(col="white"),box.rectangle=list(col="black"),add.text=list(font=3),box.umbrella=list(col="black"),box.dot=list(col="black"),plot.symbol=list(col="black",pch=1,cex=1))

Veja que tem parâmetros ai para mudar os boxplot feito com o irmão do xyplot que é o bwplot, que faz boxplot tipo lattice, mas isso não tem problema.

Depois de rodar ele, é so usar meu novo panel com o loess adicionado, e podemos dar configurações do strip pelo xyplot também, no caso aqui eu so vou alterar os nomes via argumento da função xyplot, isso é muito legal, para pegar configurações de gráficos bonitas que outras pessoas disponibilizam e so mudar coisas como os nomes, que são específicos.

1
2
xyplot(massa~tamanho|pop,data=dados,panel=meu.panel,xlab="Tamanho",ylab="Massa",
       strip = strip.custom(factor.levels =paste("População",1:4)))

Figura 08
08

Bem é isso ai, o script vai estar la no repositório recologia, espero do fundo do coração retormar os posts de evolução, para chegar logo a parte de invasão de espécies/fenótipos. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
library(lattice)
set.seed(1)
 
##################################
## Gerando os dados
##################################
n.grupos <- 4 # Número de populações
n.amostras <- 30 # Número de Cobras por amostras
n <- n.grupos * n.amostras # Total de amostras (Total de passarinhos)
pop <- gl(n = n.grupos, k = n.amostras) # Indicador da população
# Tamanho do corpo (cm)
tamanho.original <- runif(n, 45, 70)
mn <- mean(tamanho.original)
sd <- sd(tamanho.original)
#Mesmo que usar o comando scale(), isso faz a média ser zero
tamanho <- (tamanho.original - mn) / sd
Xmat <- model.matrix(~pop*tamanho-1-tamanho) #Aqui estamos fazendo o
#contrario que o comando formula faz, para simular uai
intercepito.mean <- 230	 # intercepto médio
intercepito.sd <- 20 # desvio dos interceptos
inclinação.mean <- 60 # inclinação média
inclinação.sd <- 30 # desvio das inclinações
intercepitos<-rnorm(n = n.grupos, mean = intercepito.mean, sd = intercepito.sd)
inclinações <- rnorm(n = n.grupos, mean = inclinação.mean, sd = inclinação.sd)
todos.os.efeitos <- c(intercepitos,inclinações) # Juntando tudo
lin.pred <- Xmat[,] %*% todos.os.efeitos # Preditores lineares
eps <- rnorm(n = n, mean = 0, sd = 30) # residuos
# resposta = preditor linear + residuo
massa <- lin.pred + eps
dados<-data.frame(pop,tamanho,massa)
#deixando somente os dados em dataframe, para acessar certo
rm(list=ls()[-1])
 
##################################
##Figuras tipo lattice
##################################
head(dados)
 
#figura 1
plot(massa~tamanho,data=dados)
 
#figura 2
plot(massa~tamanho,pch=19,frame=F,xlab="Tamanho",ylab="Massa",col=dados$pop,data=dados)
legend("topleft",pch=19,col=1:4,legend=paste("População",1:4),bty="n")
 
#figura 3
layout(matrix(1:4,ncol=2,nrow=2,byrow=T))
for(i in 1:4) {
    plot(massa~tamanho,pch=19,frame=F,xlab="Tamanho",ylab="Massa",
         col=i,data=dados[dados$pop==i,],main=paste("População",i),
         xlim=c(min(dados$tamanho),max(dados$tamanho)),
         ylim=c(min(dados$massa),max(dados$massa)) )
    ordem <- order(dados[dados$pop==i,"tamanho"])
    ajuste<-loess(massa~tamanho,dados[dados$pop==i,])
    points(ajuste$fitted[ordem]~dados[dados$pop==i,"tamanho"][ordem],
           type="l",col="red",lty=3,lwd=2)
}
 
#figura 4
layout(matrix(c(1:3,rep(4,3)),ncol=3,nrow=2,byrow=T))
for(i in 1:4) {
    plot(massa~tamanho,pch=19,frame=F,xlab="Tamanho",ylab="Massa",
         col=i,data=dados[dados$pop==i,],main=paste("População",i),
         xlim=c(min(dados$tamanho),max(dados$tamanho)),
         ylim=c(min(dados$massa),max(dados$massa)) )
    ordem <- order(dados[dados$pop==i,"tamanho"])
    ajuste<-loess(massa~tamanho,dados[dados$pop==i,])
    points(ajuste$fitted[ordem]~dados[dados$pop==i,"tamanho"][ordem],
           type="l",col="red",lty=3,lwd=2)
}
 
#figura 5
xyplot(massa~tamanho|pop,data=dados)
 
#figura 6
xyplot(massa~tamanho|pop,data=dados,pch=19,col=dados$pop)
 
#figura 7
meu.panel = function(x,y,subscripts) {
    panel.xyplot(x,y,col=dados$pop[subscripts],pch=19)
}
 
xyplot(massa~tamanho|pop,data=dados,panel=meu.panel)
 
                                        #figura 8
meu.panel = function(x,y,subscripts) {
    panel.xyplot(x,y,col=dados$pop[subscripts],pch=19)
    panel.loess(x, y,col="red",lty=2,lwd=2)
}
 
trellis.par.set(strip.background=list(col="white"),box.rectangle=list(col="black"),add.text=list(font=3),
box.umbrella=list(col="black"),box.dot=list(col="black"),plot.symbol=list(col="black",pch=1,cex=1))
 
xyplot(massa~tamanho|pop,data=dados,panel=meu.panel,xlab="Tamanho",ylab="Massa",
       strip = strip.custom(factor.levels =paste("População",1:4)))