Visualizando a distribuição normal multivariada com elipses

Uma distribuição interessante é a distribuição normal multivariada, que é um grupo de distribuições normais.

A parte interessante vem do fato que muito das coisas que a gente olha, são multivariadas, vem de várias distribuições, então entender como ela funciona pode ser útil, para captar melhor os modelos que fazemos para descrever o mundo.

Mas vamos ser mais específico, alguém já ficou pensando o que diabos é essa matriz que sai junto com alguns modelos?

Relembrando esse post aqui

Correlation of Fixed Effects: levsl1 levsl2 levsl3 levsl4 levsl2 0.990 levsl3 0.990 0.990 levsl4 0.990 0.990 0.990 levsl5 0.990 0.990 0.990 0.990

Veja que no final do modelo misto, sai essa matriz ai, que a gente, eu pelo menos, ignorava por um bom tempo, por não saber como diabos funciona uma matriz, para um monte de distribuição normal que são os parâmetros calculados ali em cima no modelo, opa, olha algo repetido aqui, monte de distribuição normal, e a distribuição normal multivariada é um monte de distribuição normal junta. Bem acho que agora já não parece que a gente ta so olhando mais uma distribuição maluca a esmo, e cedo ou tarde, isso vai aparecer.

Uma distribuição normal multivariada, é então como a distribuição normal, possui média e variância, com a diferença que ela é multi, então são várias médias e várias variâncias, so que não somente isso, ela introduz o conceito de correlação entre as amostras. Sem isso, a gente poderia trata-las de forma independente, mas ela é uma distribuição multivariada porque não se trata somente das médias e das variâncias, mas da correlação entre as amostras.

Na página da wikipedia sobre a distribuição, que tem aquele quadrinho legal do lado que sumariza várias coisas legais, como quais são os parâmetros da distribuição, ele fala que entre os parâmetros temos as médias, que são a localização no gráfico, e a matriz de correlação.

Vamos pensar num exemplo que é simples de representar graficamente, que é o de duas distribuições normais, um exemplo 2d.

Para gerar dados de uma distribuição, a gente pode usar a função rnorm, que vai perguntar quantas amostras queremos, e alguns parâmetros, que vão ser a média e o desvio-padrão (que ao quadrado é a variância)

Description Density, distribution function, quantile function and random generation for the normal distribution with mean equal to mean and standard deviation equal to sd. Usage dnorm(x, mean = 0, sd = 1, log = FALSE) pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) rnorm(n, mean = 0, sd = 1)

Mas existe, no pacote MASS, uma função equivalente ao rnorm para a distribuição normal multivariada que é o mvrnorm.

Ele vai precisar do número de amostras que você quer e os parâmetros, que aqui vai ser um vetor de médias, que é a localização das amostras, e o vetor de covariância, que é como se fosse a variância, mas tem a correlação entre as duas distribuições.

Vamos fazer aqui o seguinte exemplo, teremos primeiro como os parâmetros

   Parametros = \left\{       \begin{array}{lr}         \mu = {0,0}\\         \Sigma = \left(\begin{array}{cc} 1&0 \\ 0&1 \end{array} \right)       \end{array}     \right.

Veja então que é um vetor de médias e uma matriz de variância

mu<-c(0,0)
mu
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma

E podemos gerar uma amostra

amostras<-mvrnorm(n=1000, mu, Sigma)

E ver como fica

03

Veja que basicamente é um círculo, ou seja a gente tem o centro do círculo mais ou menos no ponto [0,0] e o raio vai ser mais ou menos a variância, mais ou menos porque é uma distribuição normal, então temos variação.

Agora so como exercício, olha que legal. Olhando o circulo no wikipedia, mais especificamente a equação paramétrica, a gente consegue reproduzir ela no R e desenha um círculo.

plot(amostras)
t<-seq(0, 2*pi, length = 100)
r<-1
a<-0
b<-0
points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col="red")

04

Muito legal né, inútil mas legal. So que com esse círculo ali, a gente pode ver quantos porcentos desses monte de ponto estão dentro do círculo. Como?
Por isso tem esse monte de formula de círculo, várias formas de equação, que cada uma serve pra uma coisa, a forma paramétrica a gente usa para desenhar, como ali em cima. A equação do círculo, que tem ali em cima da equação paramétrica no wikipedia, a gente faz uma função para verificar quem está dentro, quem está fora.

> ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE) [1] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE . . . [991] TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE

Agora a gente pode contar, quantos pontos estão dentro desse círculo.

> sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras) [1] 0.387

So contar com sum, lembrando que TRUE vale 1 e FALSE 0, então a gente soma esse vetor e divide pelo total de amostras e temos a nossa contagem.

Mas a gente pode tentar achar o valor de raio que englobe 95% dos pontos, que desenhe um circulo que englobe 95% dos pontos. Como? Repita isso para vários valores de raio.

> for(r in seq(0,4,0.5)) { + cat("Raio ",r," cobre ",(sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras))*100,"% dos pontos \n",sep="") + } Raio 0 cobre 0% dos pontos Raio 0.5 cobre 10.8% dos pontos Raio 1 cobre 38.7% dos pontos Raio 1.5 cobre 68.7% dos pontos Raio 2 cobre 86.6% dos pontos Raio 2.5 cobre 94.8% dos pontos Raio 3 cobre 98.7% dos pontos Raio 3.5 cobre 99.8% dos pontos Raio 4 cobre 100% dos pontos

Ha, olha ai, um raio de 2.5 cobriria praticamente 95% dos pontos. E como usar essa informação de um jeito legal, a gente pode fazer tipo um kernel pros pontos, ir colocando uns círculos pra gente ter uma noção melhor da distribuição dos pontos nesse plano.

05

Massa né, mas agora esse esquema usando o círculo vai pifar se tivermos covarianças. Vamos alterar a matriz de covariância.

Vamos usar os seguintes parâmetros agora para a distribuição

   Parametros = \left\{       \begin{array}{lr}         \mu = {0,0}\\         \Sigma = \left(\begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array} \right)       \end{array}     \right.

E olha o que vira.

06

Agora ja da pra compreender melhor essa matriz de covariância.

Veja que as variâncias tanto n eixo x como no eixo y são próximos de 1. Que é o que ta na diagonal central da matriz.
E o 0.8 ali é o quanto um eixo esta correlacionado um com o outro.

> var(amostras[,1]) [1] 0.9343743 > var(amostras[,2]) [1] 0.9513073 > cor(amostras) [,1] [,2] [1,] 1.0000000 0.7826158 [2,] 0.7826158 1.0000000

E a gente pode ir mexendo nesses parametros para ver o como vai ficando a figura.

07

E o circulo acabou de ir pro bebeleu, que ele não serve mais para tentar representar isso ai, mas a gente tem a elipse, que pode servir.

Olha no wikipedia la, a formula da elipse, é diferente, mas podemos simplesmente transcrever ela pro R a plotar para ver.

t<-seq(0, 2*pi, length = 100)
Xc<-0
Yc<-0
alpha<-pi/3
a<-1
b<-1
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")

E agora ao plotar no R vemos a elipse.

08

So que a elipse tem bem mais parâmetros que o circulo, que só tem o raio, e o ponto de origem. Aqui a gente tem o ponto de origem, um a e b e angulo alpha que vão controlar a forma da elipse. E o negócio é sempre sair jogando valor e ver no que da para pegar o feeling, mudando somente o alpha, veja o que da para fazer

08

Agora não vai dar mais para sair usando força bruta, como fizemos com o círculo, num vai dar muito certo aqui, porque são muitas combinações possíveis de parâmetros. Mas no pacote MASS do R, existe uma função chamada cov.mve que é um estimador pronto dos parâmetros da distribuição multivariada, que é o que precisamos para desenhar a elipse. Então podemos usar ela, e assim como no caso do circulo, podemos traçar uma elipse que cobre uma % desejada dos pontos.

09

Agora é uma questão de cobertura, assim como fizemos com o círculo, podemos desejar outros tamanhos de coberturas, ou %.

10

Mas a elipse so vai servir para a distribuição multivariada normal, as vezes quando a gente ta olhando esse tipo de dados, em 2d lembra muito dados espaciais, como quando olhando mapas, como no caso do John Snow, aqui a gente precisa de outras abordagens, como a distribuição posterior, que tem uma função pronta para o calculo no pacote de nosso querido Ben Bolker (O que ele faria?).

12

Outras soluções, talvez para uma primeira visualização, é usar imagens, onde cada pixel tem a densidade local expressa na intensidade de cor, por exemplo usando o hex do ggplot2.

13

E é isso ai, a gente so deu uma olhada numa distribuição, pra matar a curiosidade, e viu algumas formas de representação dela, bem como outras formas de representação de dados em 2d, para sei la, facilitar a visualização. Qualquer mapinha você vai precisar dessas coisas. A, outra coisa interessante.

01

02

Veja como a média, que é a localização, pode ser mantida, mesmo mudando os parâmetros de covariância. Veja que em cima e do lado do gráfico, a gente vê a distribuição de uma dimensão, mas no meio a gente ve a distribuição em duas dimensões.

Fim de ano ta difícil ter tempo para mexer com esses posts, mas já já acaba. Por enquanto é isso ai, apenas uma olhada em no mundo 2d. Os scripts sempre no repositório recologia e até o próximo post.

#Muito do codigo veio desse post
#http://www.sumsar.net/blog/2014/11/how-to-summarize-a-2d-posterior-using-a-highest-density-ellipse/
 
library(MASS)
library(cluster)
library(plotrix)
library(emdbook)
library(ggplot2)
library(hexbin)
#install.packages("devtools")
#library(devtools)
#install_github("rasmusab/bayesian_first_aid")
library(BayesianFirstAid)
 
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma
 
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
 
#
plot(amostras)
t<-seq(0, 2*pi, length = 100)
r<-1
a<-0
b<-0
points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col="red")
 
#Testando quem está dentro do circulo de raio r
sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras)
 
for(r in seq(0,4,0.2)) {
    cat("Raio ",r," cobre ",(sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras))*100,"% dos pontos \n",sep="")
}
 
#
plot(amostras,frame=F)
i=2
for(r in c(3.1,2.5,1.2)) {
    points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col=i)
    i=i+1
}
legend("topright",title="cobertura",legend=c("99%","95%","50%"),bty="n",lwd=4,lty=3,col=c(2,3,4))
 
#Distribuição não redondinha
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
Sigma
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
#
plot(amostras)
 
var(amostras[,1])
var(amostras[,2])
cor(amostras)
 
 
#
par(mfrow=c(2,2))
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=-2,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(1,-0.8,-0.8,1),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=2,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(1,0.8,0.8,2),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=-3,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(2,0.8,0.8,2),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=3,y=-3,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
 
 
#Desenhando Elipse
t<-seq(0, 2*pi, length = 100)
Xc<-0
Yc<-0
alpha<-pi/3
a<-1
b<-1
 
#
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
 
 
#
par(mfrow=c(2,2))
alpha<-pi/3.5
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topleft",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/3
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topleft",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/2
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("center",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/1.5
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topright",legend=paste("alpha=",round(alpha,2)),bty="n")
 
 
#Ajuste de elipse
ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * 0.50)
pontos_na_ellipse <- amostras[ajuste$best, ]
limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
plot(amostras)
points(limite_ellipse[,1],limite_ellipse[,2],type="l",lwd=5,lty=3,col="red")
legend("topleft", "50%",lwd=5,lty=3,col="red")
 
 
 
#Varios ajustes
plot(amostras)
i<-2
for(cobertura in c(0.99, 0.95, 0.5)) {
    ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * cobertura)
    pontos_na_ellipse <- amostras[ajuste$best, ]
    limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
    points(limite_ellipse[,1],limite_ellipse[,2],type="l",lwd=5,lty=3,col=i)
    i<-i+1
}
legend("topleft",legend=c("99%","95%","50%"),lwd=5,lty=3,col=2:4)
 
 
#Plot com poligonos
plot(amostras,type="n")
i<-2
for(cobertura in c(0.99, 0.95, 0.5)) {
    ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * cobertura)
    pontos_na_ellipse <- amostras[ajuste$best, ]
    limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
    polygon(limite_ellipse, col = i, border = NA)
    i<-i+1
}
points(amostras[,1],amostras[,2],pch=19)
legend("topleft",legend=c("99%","95%","50%"),pch=1,col=2:4)
 
#Distribuição posterior
plot(amostras)
HPDregionplot(amostras, prob = c(0.95, 0.75, 0.5), col=2:4, lwd=3,lty=3, add=TRUE)
legend("topleft", legend = c("95%", "75%", "50%"), col = 2:4, lty=3, lwd=3)
 
 
 
#ggplot2 hexagonos
qplot(amostras[,1], amostras[,2], geom=c("hex"))
 
#Bayes plot
amostras<-mvrnorm(n=50, c(0,0), Sigma)
ajuste <- bayes.cor.test(amostras[,1], amostras[,2])
plot(ajuste)

R Orientado a Objetos e um exemplo S3

Orientação a Objetos é um estilo de programação que pode simplificar muitos problemas.

A orientação a objetos se baseia em definir classes, criar e manipular objetos dessas classes. Um objeto nada mais é que uma variável com uma estrutura particular, como um Dataframe tem que a cara de uma planilha ou uma matrix de inteiros, que é uma lista de inteiros com uma estrutura de duas dimensões.

No R, é possível definir nossas próprias classes, se a gente quiser. Existem dois estilos de orientação a objeto no R, um mais antigo chamado de S3 e um mais novo chamado de S4.

Seguindo um desses estilos a gente pode criar classes e definir seus métodos, que são funções que sabem lidar com objetos daquela classe. Orientação a objetos não é necessariamente essencial, mas pode deixar as coisas bem organizadas.

Por exemplo a gente pode olhar a função mean, quando a gente digita mean, a gente ve o seguinte:

> mean function (x, ...) UseMethod("mean")

No R, a mean é uma função genérica que examina qual a classe do objeto do primeiro argumento e escolhe outra função apropriada para a classe.
Quando a gente digita mean, ela meramente para uma string chamada “mean” para uma função chamada UseMethod. A função UseMethod que chama a versão de mean que é designada para a classe do objeto do primeiro argumento.

Podemos ver quais métodos mean tem usando a função methods

> methods(mean) [1] mean.Date mean.default mean.difftime mean.POSIXct mean.POSIXlt

Se a classe é numeric, a gente cai em uma função mean, se a classe for Date, a gente usa outra função mean, e assim por diante.

Qualquer função genérica pode ser estendida para outras classes, se a gente criar nossa classe, para representar algo especial, poderíamos ter um mean.NossaClasse.

Mas ilustrando o caso acima, quando temos números fazemos a média aritmética

> objeto1<-c(30,1,2,3) > class(objeto1) [1] "numeric" > mean(objeto1) [1] 9

Mas se tivermos datas, cada mês acaba num dia diferente, então se acima ali fossem dias entre janeiro e fevereiro, teríamos:

> datas<- c("30/01/14", "01/02/14", "02/02/14", "03/02/14") > objeto2<-as.Date(datas, "%d/%m/%y") > class(objeto2) [1] "Date" > mean(objeto2) [1] "2014-02-01"

A média ta dia 1, que é o meio entre dia 30 -> 1 -> 2 -> 3. E não tem como ter fração de dia, fração de dias é medida em horas. Ou seja, cada classe é de um jeito, e apesar de tudo ser representado com números, o comportamento precisa ser diferente.

Basicamente, para escrever uma função que vai ser chamada automaticamente para substituir uma função genérica fu qualquer (fu pode ser a mean por exemplo, em computação, sempre o povo da exemplo de função como fu, sei la) para um tipo de classe, a gente so precisa definir ela como fu.NomeDaClasse. Claro que funções não precisam sempre de uma função genérica, a gente pode escrever uma função que so aceita objetos de uma classe a cabou, mas muitos nomes comuns são usados em múltiplos objetos, como mean, a gente calcula media de um monte de coisa. Ai é útil redefinir seu próprio método.

Outra função que tem varios métodos é a plot

> methods(plot) [1] plot.acf* plot.data.frame* plot.decomposed.ts* plot.default [5] plot.dendrogram* plot.density* plot.ecdf plot.factor* [9] plot.formula* plot.function plot.hclust* plot.histogram* [13] plot.HoltWinters* plot.isoreg* plot.lm* plot.medpolish* [17] plot.mlm* plot.ppr* plot.prcomp* plot.princomp* [21] plot.profile.nls* plot.spec* plot.stepfun plot.stl* [25] plot.table* plot.ts plot.tskernel* plot.TukeyHSD* Non-visible functions are asterisked

É avaliando a classe do objeto, que a plot decide o que fazer, por isso quando a gente da plot, as vezes sai boxplot, as vezes scatterplot, ou várias gráficos quando damos plot num modelo, o R olha a classe e alguém pensou qual deveria ser a escolha obvia de plot para aquela classe.

Mas ja olhamos bastante o que ta pronto, agora vamos tentar criar nossa própria classe so para ter o gostinho de fazer uma classe.

Vamos pegar um exemplo com plantinhas. Para entender como as plantas ocupam uma área, é importante entender sua dispersão de sementes. A dispersão de semente pode ser medida usando armadilhas para sementes, que são caixas ou panos de um tamanho conhecido que ficam colocadas a uma distância conhecida de planta que vai dispersar sementes por um período de tempo conhecido. O número de sementes que chega a cada armadilha pode ser então contadas, e assim podemos pensar como funciona a dispersão naquela planta.
Aqui, vamos assumir que as armadilhas de sementes são deixadas em uma linha reta em relação a planta mãe, e depois de um tempo as sementes de cada armadilha são contadas.

01

Aqui o pontinho verde é o centro da nossa plantinha, o circulo é a copa dela, aqui com raio 1, e os quadradinhos são nossas armadilhas, coletando sementes.

Agora a gente pode imaginar que os dados disponíveis para analisar de cada amostra aqui é um conjunto de distancias, de cada armadilha até a planta, a área de cada armadilha e sua contagem de sementes em cada armadilha.

Assim a gente pode montar uma classe que guarde exatamente esses dados, que nos interessam aqui, e também dar um nome pra ela, por exemplo, ArmadilhaDeSementes.

Veja que não existem métodos para essa classe.

> methods(class = "ArmadilhaDeSementes") no methods were found

Agora pelo descrição ali, a essa classe tem que ter as distâncias das armadilhas, as contagens de sementes de cada uma, e uma área de armadilha.

Para objetos S3, a gente cria primeiro uma função que cria um objeto daquele tipo.

ArmadilhaDeSementes <- function(distancias, contagem.sementes, area.armadilha = 0.0001) {
    if (length(distancias) != length(contagem.sementes))
        stop("Quantidade de distancias e contagens é diferente.")
    if (length(area.armadilha) != 1) stop("Area da armadilha é ambigua.")
    ArmadilhaDeSementes <- list(distancias = distancias,
                         contagem.semetes = contagem.sementes,
                         area.armadilha = area.armadilha)
    class(ArmadilhaDeSementes) <- "ArmadilhaDeSementes"
    return(ArmadilhaDeSementes)
}

Veja que essa função, é uma função, mas não é uma do tipo que faz alguma operação, como normalmente montamos, elas é um construtor de classe. Ela vai receber as informações e vai guardar elas de forma organizada, a organização da nossa classe para ser mais específico. Então a gente vai receber as distâncias e as contagens de sementes, podemos assumir que existe uma armadilha padrão, mas que é passível de ser alterada, ai vamos guardar isso numa lista. Ai também vamos informar pro R qual o nome da classe, lembre-se que nas funções genéricas, essa é a informação que que é passada pra frente, e se você digitar o nome do objeto, você vai receber de volta a lista de informações dele.
Agora a outra coisa legal, é que podemos validar as informações sendo entradas, veja que não permitimos armadilhas sem contagens, o número de contagens tem que ser igual ao número de distâncias para criarmos o objeto, poderíamos também não aceitar números negativos para as contagens, e muitas outras verificações, o que garantiria que sempre os dados vão estar consistentes e assim quando a gente usar algum método nele, esse método não vai falhar.

Agora que temos uma classe, guardando os dados de forma organizada e consistente, podemos criar os métodos para ela, ou melhor, fazer com que métodos básicos consigam entender os dados dessa classe.

Para tal vamos mexer em dois métodos básicos, o print, que é o que imprime o conteúdo da classe de forma “legível” ao usuário, e o método mean, pra calcular a média, que agora vai ser a dispersão média, que vamos calcular de forma diferente da média aritmética.

print.ArmadilhaDeSementes <- function(objeto, ...) {
    cat(paste("Distâncias: ",paste(objeto$distancias,collapse=", "),"\n"))
    cat(paste("Contagens : ",paste(objeto$contagem.sementes,collapse=", "),"\n"))
    cat(paste("Área da armadilha: ",objeto$area.armadilha,"\n"))
}
mean.ArmadilhaDeSementes <- function(objeto, ...) {
    return(weighted.mean(objeto$distancias, w = objeto$contagem.sementes))
}

A assim ficaram nossos métodos agora. Veja que no print, basicamente a gente escreve o conteúdo, os dados, de forma mais bonitinha. Talvez a gente poderia ja dar alguma medida derivada, por exemplo no print para sequência do pacote ape, a print não da a sequência, mas a porcentagem de cada base, porque dificilmente você quer ver a sequência inteira logo de cara, quando você abre uma sequência a primeira coisa que você quer ver é que conseguiu abrir, ai você da outra informação, e depois faz outro método para ver as letrinhas, pra não poluir o console com mil informações inúteis, mas aqui a gente tenta passar a informação de forma mais caprichada.

E no caso da média, a gente vai calcular a media da distância ponderada pela contagem de sementes, que vai ser nossa dispersão média, então temos que acessar o objeto em questão, extrair a informação que nos interessa para esse método e usa-la.
Agora não é porque não usamos a área da armadilha que ela não deveria estar la, pense que na classe agora, temos informação organizada, que podemos ir fazendo os métodos e sabemos como chamá-las.

Veja que se agora, a gente perguntar pro R se existem métodos da nossa classe, ele vai falar que existem dois, que acabamos de criar. Sendo que ele acha isso por usamos o .NomeDaClasse, isso que torna possível a associação entre função com classe.

> methods(class = "ArmadilhaDeSementes") [1] mean.ArmadilhaDeSementes print.ArmadilhaDeSementes

Veja que se olharmos os métodos do mean denovo, agora temos nosso método la.

> methods(mean) [1] mean.ArmadilhaDeSementes mean.Date mean.default [4] mean.difftime mean.POSIXct mean.POSIXlt

Assim como o print, existe nosso print.ArmadilhaDeSementes agora.

> methods(print) [1] print.acf* print.anova* [3] print.aov* print.aovlist* [5] print.ar* print.Arima* [7] print.arima0* print.ArmadilhaDeSementes [9] print.AsIs print.aspell* . . . [181] print.xngettext* print.xtabs* Non-visible functions are asterisked >

Agora podemos criar objetos do tipo ArmadilhaDeSementes e usar nossos métodos.

> s1 <- ArmadilhaDeSementes(distancias = 1:4, contagem.sementes = c(4, 3, 2, 0)) > s1 Distâncias: 1, 2, 3, 4 Contagens : 4, 3, 2, 0 Area da armadilha: 1e-04 > mean(s1) [1] 1.777778 >

Veja que a classe do objeto s1 que criamos é a nossa classe

> class(s1) [1] "ArmadilhaDeSementes"

Mas outra coisa, é que apesar de essa ser a nossa classe, veja que dentro da classe, a informação esta na forma de uma lista, então todas as características da classe list também são herdadas na nossa classe, ou seja podemos usar nossa classe como se fosse uma lista também.

> is.list(s1) [1] TRUE > s1[[1]] [1] 1 2 3 4

Criar toda essa infraestrutura pode parecer a princípio um monte de trabalho dispensável. A gente pode pensar, poxa porque simplesmente a gente não cria dois vetores e faz a média ponderada. Mas as vantagens aparecem quando a gente começa a ter que construir, manipular e analisar modelos de sistemas complexos. Orientação a objetos suporta uma prototipagem fácil e rápida, e um caminho para ir adicionando complexidade conforme essa se faz necessária. Por exemplo é fácil fazer uma nova função que calcula a média da nossa classe ArmadilhaDeSementes, e além disso, a gente faz isso sem correr risco de mexer em funções de outros pacotes, e cagar alguma coisa do R que pode ter impacto no R como um todo, como muitas funções dependem umas das outras, é perigoso mexer nas coisas, ainda mais funções básicas assim. A orientação a objeto oferece essa possibilidade de encapsulamento, você mexe nas suas coisas sem correr riscos desnecessários. Existe ainda uma consideração final. É possível que mais de um pacote tenha uma mesma função para lidar com um mesmo tipo de classe, so pensar que classes como sequência de dna são usadas por muitos pacotes, o que tornaria pacotes incompatíveis. Uma solução é usar os namespaces. Um pacote pode criar um namespace e se quisermos ter certeza que estamos acessando aquela função, é so usarmos o nome do pacote::função (help(“::”) para mais informações).

Mas se a gente digitar iris, a gente ve um dataset que esta no pacote datasets

> iris Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa . . . 149 6.2 3.4 5.4 2.3 virginica 150 5.9 3.0 5.1 1.8 virginica >

Isso porque esse pacote está aberto e o R sabe onde procurar.

> search() [1] ".GlobalEnv" "ESSR" "package:stats" "package:graphics" [5] "package:grDevices" "package:utils" "package:datasets" "package:methods" [9] "Autoloads" "package:base"

Veja que a função search() mostra o que o R está olhando, quais os namespace, o R ta olhando. Mas suponha que o R não tivesse aberto o pacote datasets, a gente pode fechar ele.

> detach("package:datasets", unload=TRUE) > search() [1] ".GlobalEnv" "ESSR" "package:stats" "package:graphics" [5] "package:grDevices" "package:utils" "package:methods" "Autoloads" [9] "package:base" >

Agora quando a gente olhar o search, o R não enxerga mais o datasets.

> iris Erro: objeto 'iris' não encontrado

Mas se você especificar o namespace da certo.

> datasets::iris Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa . . . 149 6.2 3.4 5.4 2.3 virginica 150 5.9 3.0 5.1 1.8 virginica >

Bem, acho que é isso ai, fiquei um tempo sem fazer posts, mas espero voltar a postar com regularidade denovo. Tem bastante coisa para continuar, e agora fiquei devendo um post sobre classes s4. Os scripts sempre no repositório recologia e é isso, até o próximo post, que espero que seja logo, se você leu até aqui e me viu escrevendo alguma abobrinha, por favor não deixe de me avisar do meu erro por e-mail ou fazendo um comentário.

Referência:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

mean
methods(mean)
 
objeto1<-c(30,1,2,3)
class(objeto1)
mean(objeto1)
 
 
datas<- c("30/01/14", "01/02/14", "02/02/14", "03/02/14")
objeto2<-as.Date(datas, "%d/%m/%y")
class(objeto2)
mean(objeto2)
 
methods(plot)
 
 
#Parametros para desenhar um circulo.
#a e b são o ponto central, r é o raio e t é a variavel parametrica que vai de 0 a 2pi
t<-seq(0, 2*pi, length = 100)
r<-1
a<-0
b<-0
 
plot(a+r*cos(t),b+r*sin(t),type="l",xlim=c(-2,5),ylim=c(-3,3))
points(0,0,pch=19,col="green")
points(c(0.5,1,1.5,2,2.5),rep(0,5),pch=22,cex=2)
 
methods(class = "ArmadilhaDeSementes")
 
ArmadilhaDeSementes <- function(distancias, contagem.sementes, area.armadilha = 0.0001) {
    if (length(distancias) != length(contagem.sementes))
        stop("Quantidade de distancias e contagens é diferente.")
    if (length(area.armadilha) != 1) stop("Area da armadilha é ambigua.")
    ArmadilhaDeSementes <- list(distancias = distancias,
                         contagem.sementes = contagem.sementes,
                         area.armadilha = area.armadilha)
    class(ArmadilhaDeSementes) <- "ArmadilhaDeSementes"
    return(ArmadilhaDeSementes)
}
 
print.ArmadilhaDeSementes <- function(objeto, ...) {
    cat(paste("Distâncias: ",paste(objeto$distancias,collapse=", "),"\n"))
    cat(paste("Contagens : ",paste(objeto$contagem.sementes,collapse=", "),"\n"))
    cat(paste("Area da armadilha: ",objeto$area.armadilha,"\n"))
}
 
mean.ArmadilhaDeSementes <- function(x, ...) {
    return(weighted.mean(x$distancias, w = x$contagem.sementes))
}
 
methods(class = "ArmadilhaDeSementes")
methods(mean)
methods(print)
 
s1 <- ArmadilhaDeSementes(distancias = 1:4, contagem.sementes = c(4, 3, 2, 0))
s1
mean(s1)
 
class(s1)
 
is.list(s1)
 
s1[[1]]
 
help("::")
envir
 
search()
 
iris
 
detach("package:datasets", unload=TRUE)
search()
 
iris
 
datasets::iris