# É er(R)ando que se aprende - 3
Uma introdução a linguagem R para análise de dados - 3.

Parte 3: http://bit.ly/aprendeerrando3
Parte 4: http://bit.ly/aprendeerrando4
Parte 1: http://bit.ly/aprendeerrando
Parte 2: http://bit.ly/aprendeerrando2

## Obtenção de dados

> Link para download dos conjuntos de dados: https://drive.google.com/drive/folders/1Up7qIFrxjaAlCkdOQZpNrEM_JOYIQvsZ?usp=sharing

##### paises:
> Area, tamanho populacional, taxas de nascimento e óbito, expectativa de vida e forma de governo de 249 países

##### dados_politicos:
> Conjunto de dados políticos e institucionais de 36 países

#Carregando texto separado por tabulação (.tsv ou .txt)
t <- read.table(f, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
head(t)
tail(t)

#Carregando texto separado por vírgula (.csv)
v <- read.table(f, header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(v)
tail(v)

#Carregando arquivo do excel (.csv)
install.packages("readxl")
library("readxl")
e <- read_excel(f, sheet = 1, col_names = TRUE)
head(e)
tail(e)
str(e)

# OU 
install.packages("XLConnect")
install.packages("XLConnectJars")
library("XLConnect")
library("XLConnectJars")
e1 <- readWorksheetFromFile(f, sheet = 1, header = TRUE)
head(d)

#Criando um arquivo excel no R
writeWorksheetToFile(e2, e, sheet = "myData", clearSheets = TRUE)

### Explorando os dados

d <- read.csv(f, header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d)
summary(d)
names(d)

#Como saber a densidade populacional dos países?
d$density <- d$population/d$area

#Quais são os 10 países com maior densidade populacional?
d <- d[order(-d$density), ]
d[1:10, ]

#E os países com menor densidade populacional?
d <- d[order(d$density), ]
d[1:10, ]

#Extraindo os países que começam com "A" até "F" e calculando a média e o tamanho populacional desses países
new <- d[grep("^[A-F]", d$country), ]
summary(new)
mean(new$population, na.rm = TRUE)
mean(new$area, na.rm = TRUE)

#Criando um histograma com os dados em log() de tamanho populacional e area
par(mfrow = c(1, 2))  # Nos dá dois paineis
attach(d) # Com essa função, não precisamos usar $ para acessar as colunas do conjunto de dados
hist(log(population), freq = FALSE, col = "red", main = "População", xlab = "log(Tamanho populacional)", 
    ylab = "density", ylim = c(0, 0.2))
    
hist(log(area), freq = FALSE, col = "red", main = "Tamanho da area", xlab = "log(Area)", 
    ylab = "density", ylim = c(0, 0.2))
    
#Inserindo linha de média e estimativa de distribuição
par(mfrow = c(1, 1))  # Estabelece um painel e redesenha o histograma acima de log(Tamanho populacional)
hist(log(population), freq = FALSE, col = "white", main = "Histograma com média e densidade", 
    xlab = "log(Tamanho populacional)", ylab = "Densidade", ylim = c(0, 0.2))
abline(v = mean(log(population), na.rm = TRUE), col = "blue")
lines(density(log(population), na.rm = TRUE), col = "green")
detach(d)

### Medidas de tendência e variância
> São essenciais para analises descritivas e teste de hipótese
##### Conceitos importantes:
* População: Conjunto com todos os dados do conjunto de dados;
* Amostra: Uma ou mais observações da população;
* Parâmetro: Uma característica mensurável na população;
* Estatística: Uma característica mensurável na amostra.

#### Medidas de tendência
* Moda: Valor mais comum observado;
* Mediana: Valor medio em uma seria ordenado de valores;
* Média: Soma de todos os valores, divido pela quantidade de valores

#### Medidas de disperção
* Amplitude: Diferença entre o máximo e mínimo;
* Interquartil: 25th até 75th quartil;

summary()

* Soma dos quadrados: Soma dos desvios quadrados da média


* Variância populacional: Média da soma dos quadrados de toda a população;


* Variância amostral: Média da soma dos quadrados da amostra. Aqui perde-se um grau de liberdade;


* Desvio padrão populacional: Raiz quadrada da variância populacional;


* Desvio padrão amostral: Raiz quadrada da variância amostral;



> Descrição de incerteza na estimação de parametros
* Erro padrão: Verifica a confiabilidade da média. É a medida de variação da média amostral em relação à media da população.
* Através do erro padrão pode-se inferir um intervalo de confiança.
* Enquanto o desvio padrão mostra o indice de disperção da amostra em relação à media, o erro padrão avalia a confiabilidade da média.


> Intervalo de confiança: Provável faixa de valores que uma estimativa cairia caso a amostragem fosse repetida.
> Ex: Um intervalo de confiança de 95% descreve os valores em que uma estatística, calculada a partir de uma amostra repetida, deverá cair 95% das vezes.

#Analise de uma distribuição normal de números aleatórios
set.seed(1)
x <- rnorm(10000, 0, 1)
hist(x)

#Plotar a densidade e a distribuição de probabilidade
x <- seq(from = -4, to = 4, by = 0.01)
plot(x, dnorm(x), cex = 0.4)

plot(x, pnorm(x), cex = 0.4)

x <- seq(from = 0, to = 1, by = 0.01)
plot(qnorm(x), x, cex = 0.4)

#Cálculo do intervalo de confiança
x <- c(1:15)
m <- mean(x)
n <- length(x)
v <- var(x)
s <- sd(x)
e <- sqrt(v/n)
upper <- mean(x) + qnorm(0.975, mean = 0, sd = 1) * se(x)
lower <- mean(x) + qnorm(0.025, mean = 0, sd = 1) * se(x)  # or lower <- mean(x) - qnorm(0.975)*se(x)
ci <- c(lower, upper)
ci

#Interpretação do IC
    #Baseados no conjunto de dados, temos 95% probabilidade que a verdadeira média da população esteja entre o máximo e o mínimo;
    #É esperado que em uma reamostragem, os valores da distribuição caiam no mesmo intervalo 95% das vezes

### Teste de hipótese
* Hipótese nula (H0): Afirma que um parametro da população é um valor hipotético
* Hipótese alternativa(H1): É a hipotese acreditada ser verdadeira por quem faz a análise
* valor *p*: A probabilidade de que a hipótese nula seja verdadeira

#Peso de macacos em 2016. Estima-se que em 2016, os macacos estiveram mais pesados em comparação com os macacos de outros anos com média de 4,9kg. A média de 2016 é de 5.324kg. Ela é significantemente maior?

#H0?

#H1?

###### Teste T
> Usado para comparar duas amostras idependentes

library(curl)
f <-curl("https://raw.githubusercontent.com/fuzzyatelin/fuzzyatelin.github.io/master/AN597_Fall17/vervet-weights.csv")
d <- read.csv(f, header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d)
mean(d$weight)

#Calculo de media, desvio padrão e o erro padrão da média
pesoinicial <- 4.9 
pesos <- d$weight
m <- mean(x)
s <- sd(x)
n <- length(x)
sem <- s/sqrt(n)

t <- t.test(x = pesos, mu = pesoinicial, alternative = "greater")
# x é o vetor de dados a ser analisado
# mu é a media verdadeira, a ser comparada com os novos dados
# alternative diz como a hipótese alternativa deve ser em relação à hipotese nula, nesse caso h1 deve ser maior que h0

## Análises filogenéticas
> Aprendendo sobre arvores filogenéticas e os diferentes métodos de construção no R. Usaremos um conjunto de dados de amostras de influenza dos eua de 1993-2008.

install.packages(c("adegenet", "phangorn", "stats", "ape", "ade4"))

### Tipos de árvores filogenéticas
>Árvore enraizada

> Árvore não enraizada

>Diferentes tipos de representação


### Principais métodos de reconstrução de árvores
* Baseada em distância
Encontra a "distância genética" entre taxons e agrupa espécies usando a distância
* Máxima parsimônia
* Máxima verossimilhança

## Iniciando as análises

#Importando os dados
dna <- fasta2DNAbin(file="http://adegenet.r-forge.r-project.org/files/usflu.fasta")
dna

#Lendo as anotações dos dados
annot <- read.csv("http://adegenet.r-forge.r-project.org/files/usflu.annot.csv", header=TRUE, row.names=1)
annot

* Análise baseada em distância

>UPGMA: Método mais simples, assume a mesma velocidade evolutiva em todas as linhagens. Todos os braços tem a mesma distância da raiz.

> Neighbor-joining: Seleciona os dois nóss mais próximos da arvore e os define como vizinhos. Faz isso até todos os nós terem sido pareados.


#Encontrar as distâncias genéticas
D <- dist.dna(dna, model = "TN93")
length(D) #numero de distancias pareadas

temp <- as.data.frame(as.matrix(D))
table.paint(temp, cleg=0, clabel.row=.5, clabel.col=.5) #Quanto mais escuro, maior a distância

#Construindo a árvore
tre <- nj(D)
class(tre)

tre <- ladderize(tre)
tre 

plot(tre, cex = 0.6)
title("Árvore neighbor-joining")
 
#ou

h_cluster <- hclust(D, method = "average", members = NULL) 
plot(h_cluster, cex = 0.6)

#Melhorando a visualização
plot(tre, show.tip=FALSE)
title("Árvore NJ não enraizada")
myPal <- colorRampPalette(c("red","yellow","green","blue"))
tiplabels(annot$year, bg=num2col(annot$year, col.pal=myPal), cex=.5) 
temp <- pretty(1993:2008, 5)
legend("bottomleft", fill=num2col(temp, col.pal=myPal), leg=temp, ncol=2)

#Árvore não enraizada
plot(tre, type = "Não enraizado", show.tip = FALSE)
title("Árvore NJ não enraizada")
tiplabels(tre$tip.label, bg = num2col(annot$year, col.pal = myPal), cex = 0.5)

# Criando uma árvore enraizada
tre2 <- root(tre, out = 1)
tre2 <- ladderize(tre2)
plot(tre2, show.tip=FALSE, edge.width=2)
title("Rooted NJ tree")
tiplabels(tre$tip.label, bg=transp(num2col(annot$year, col.pal=myPal),.7), cex=.5, fg="transparent")
axisPhylo()
temp <- pretty(1993:2008, 5)
legend("topright", fill=transp(num2col(temp, col.pal=myPal),.7), leg=temp, ncol=2)

#Verificando se o método foi apropriado para o conjunto de dados

#Método NJ
x <- as.vector(D)
y <- as.vector(as.dist(cophenetic(tre2)))
plot(x, y, xlab="original pairwise distances", ylab="pairwise distances on the tree", main="Is NJ appropriate?", pch=20, col=transp("black",.1), cex=3)
abline(lm(y~x), col="red")
cor(x,y)^2

#Método UPGMA
tre3 <- as.phylo(hclust(D,method="average"))
y <- as.vector(as.dist(cophenetic(tre3)))
plot(x, y, xlab="original pairwise distances", ylab="pairwise distances on the tree", main="Is UPGMA appropriate?", pch=20, col=transp("black",.1), cex=3)
abline(lm(y~x), col="red")

cor(x,y)^2

#Bootstraping
myBoots <- boot.phylo(tre2, dna, function(e) root(nj(dist.dna(e, model = "TN93")),1))
myBoots

plot(tre2, show.tip=FALSE, edge.width=2)
title("NJ tree + bootstrap values")
tiplabels(frame="none", pch=20, col=transp(num2col(annot$year, col.pal=myPal),.7), cex=3, fg="transparent")
16

axisPhylo()
temp <- pretty(1993:2008, 5)
legend("topright", fill=transp(num2col(temp, col.pal=myPal),.7), leg=temp, ncol=2)
nodelabels(myBoots, cex=.6)

temp <- tre2
N <- length(tre2$tip.label)
toCollapse <- match(which(myBoots<70)+N, temp$edge[,2])
temp$edge.length[toCollapse] <- 0
tre3 <- di2multi(temp, tol=0.00001)
plot(tre3, show.tip=FALSE, edge.width=2)
title("NJ tree after collapsing weak nodes")
tiplabels(tre3$tip.label, bg=transp(num2col(annot$year, col.pal=myPal),.7), cex=.5, fg="transparent")
axisPhylo()
temp <- pretty(1993:2008, 5)
legend("topright", fill=transp(num2col(temp, col.pal=myPal),.7), leg=temp, ncol=2)