Laboratórios de Estatistica - Aulas
Rafael Oliveira
10/14/2021
- Inicio
- 1 - Estrutura de dados
- Vetores numéricos
- Operaçoes aritmeticas
- Vetores logicos
- Comandos logicos
- Vetores de Caracteres
- Fatores | levels | tapply
- Manipulaçao de vetores (1d) | length, max, min, sum, prod, sort…
- Matrizes (2d)
- Arrays (3d)
- Indexação de arrays e matrizes
- Listas
- Data frames
- Concatenando arrays | cbind, rbind
- Modos dos objetos
- Importação e Exportação
- Estruturas de Controlo e Funções
- 2 - Analise Exploratoria de Dados
- 3 - Distribuiçoes e simulaçao, reamostragem
Inicio
#https://github.com/Toniiiio/imageclipr
test1 = c(1,2,3)
test = as.data.frame(test1, row.names = c('Test', 'test1','test3'))
#library(cluster, lib.loc = "C:\Users\Rafael Oliveira\Documents\R\win-library\4.0")
# library(boot)
# library(dplyr)
# detach("package::boot", unload =TRUE)
# install.packages("")
# ls()1 - Estrutura de dados
Vetores numéricos
A=c(1,2,3,4,5,6,7)
1:20 #cria a sequencia por passos de 1 unidade## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 2020:1 # cria a sequencia por passos de 1 unidade de reversa## [1] 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1seq(1,20, by=2) #slicing 2 de cada vez## [1] 1 3 5 7 9 11 13 15 17 19seq(1,20, length.out=50) #cria 50 numeros (com decimais) de 1 até 20## [1] 1.000000 1.387755 1.775510 2.163265 2.551020 2.938776 3.326531
## [8] 3.714286 4.102041 4.489796 4.877551 5.265306 5.653061 6.040816
## [15] 6.428571 6.816327 7.204082 7.591837 7.979592 8.367347 8.755102
## [22] 9.142857 9.530612 9.918367 10.306122 10.693878 11.081633 11.469388
## [29] 11.857143 12.244898 12.632653 13.020408 13.408163 13.795918 14.183673
## [36] 14.571429 14.959184 15.346939 15.734694 16.122449 16.510204 16.897959
## [43] 17.285714 17.673469 18.061224 18.448980 18.836735 19.224490 19.612245
## [50] 20.000000Operaçoes aritmeticas
#; para colocar varios parametros. Fazemos operaçoes elemento a elemento
a=c(1,2)
b=c(5,10)
c=c(-1,-5,0)
d=c(1,2,3,4)
a+b;a-b;a*b;a/b## [1] 6 12## [1] -4 -8## [1] 5 20## [1] 0.2 0.2a+c;b+d## Warning in a + c: longer object length is not a multiple of shorter object
## length## [1] 0 -3 1## [1] 6 12 8 14#R faz reciclagem do vetor mais pequeno - Numa operação entre dois vetores de tamanhos diferentes, o mais pequeno é “reciclado” de forma a atingir o tamanho do maior e assim ser executada a operação
#Quando vetor nao tem tamanho multiplo do mais pequeno - longer object is nota multiple of shorter object lenght
#Quando vetor maior tem tamanho multiplo do mais pequeno (tamanho 4 e 2, respetivamente, por exemplo) nao dá aviso no termial (exemplo do a+d)
a*b## [1] 5 20d*d## [1] 1 4 9 16d^2## [1] 1 4 9 16c*d## Warning in c * d: longer object length is not a multiple of shorter object
## length## [1] -1 -10 0 -4#produto interno soma vetores, produto de elemento guarda em vetores
a%*%b #produto interno de vetores## [,1]
## [1,] 25#a = [a1,a2]
#b = [b1,b2]
#a*b = a1b1 + a2b2Vetores logicos
a = c(1,2)
b = c(1,5)
c=c(-1,3)
a>b #comparaçoes entre vetores sao feitas entre elemento a elemento consoante os seus respetivos indices. Neste caso vai dar [FALSE, FALSE]## [1] FALSE FALSEa>=b## [1] TRUE FALSEa==b## [1] TRUE FALSEa!=b## [1] FALSE TRUE###Disjunçao, Conjunçao, Negaçao
#TRUE = 1
#FALSE = 0
a>b & a>c #FALSE + FALSE## [1] FALSE FALSEa>b | a>c ## [1] TRUE FALSE!(a>b) #would be the same as saying "a<b"## [1] TRUE TRUE#E
#TT da T, TF da F, FT da F e FF da F
#OU
#TT da T, TF da T, FT da T FF da F
sum(a>b & a>c) #R interpreta logica de forma numerica onde T é 1 e F é 0; ## [1] 0#Neste caso da 0, logo, nao ha nenhum elemento de a que seja superior a b ou a c. Diz-nos se há pelo menos 1 elemento (ou nenhum)
sum(a>b | a>c) #O a consegue ser superior a primeira componente do b, logo, dá 1. Como o operador é "ou" basta que 1 seja verdade para validar o pressuposto## [1] 1Comandos logicos
#which #quais os indices (posiçoes) do vetor x que são TRUE
which(a>b);which(a>c)## integer(0)## [1] 1any(a)# avaliar se ha algum elemento de X que é True ou nao. Para fazer o contrario fazemos !any(x) (se ha algum elemento de x que é false ou nao)## Warning in any(a): coercing argument of type 'double' to logical## [1] TRUEany(a>b);any(a>c) ## [1] FALSE## [1] TRUEall(a) #retorna falso quando pelo menos 1 dos elementos do vetor é falso## Warning in all(a): coercing argument of type 'double' to logical## [1] TRUEall(a>b);all(a>c) #FALSO e FALSE - O 1º retorna false pois temos 2 falsos e no 2º retorna false pois temos pelo menos um false## [1] FALSE## [1] FALSE!all(a>b)## [1] TRUEall(!(a>b))## [1] TRUEtable(a) #contagem frequencias## a
## 1 2
## 1 1Vetores de Caracteres
# São criados por concatenação de expressões entre aspas:
a = c("expressão1","expressão2")
a## [1] "expressão1" "expressão2"# A função paste permite concatenar ﴾reciclando se necessário﴿ um número arbitrário de argumentos ﴾caracteres ou não﴿ um por um definindo o separador entre eles #como se fosse o .format no python ou as template strings no javascript
#paste(arg1,arg2,sep=",")
paste("Teste", a[2])## [1] "Teste expressão2"paste("Neste momento são",date())## [1] "Neste momento são Thu Dec 30 22:34:24 2021"Fatores | levels | tapply
Definem grupos de valores, previamente especificados
Criados pela função factor() e os grupos por si definidos podem ser obtidos com a função levels()
x = factor(c("Suficiente","Bom", "Bom", "Suficiente"))
x## [1] Suficiente Bom Bom Suficiente
## Levels: Bom Suficientelevels(x) #Vai mostrar apenas os valores unicos (nao repete os valores de x). Funçao que vai ao x para ir buscar os levels## [1] "Bom" "Suficiente"Podemos ainda definir a ordem desejada dos niveis:
y=factor(c("Bom","Suficiente", "Bom", "Suficiente"),
levels = c("Suficiente","Bom")) # este levels é uma opçao dentro do factor
y## [1] Bom Suficiente Bom Suficiente
## Levels: Suficiente BomOs valores pre-especificados podem tambem nao estar no vetor em causa:
z=factor(c("Suficiente", "Bom", "Bom", "Suficiente"),
levels = c("Insuficiente","Suficiente","Bom")) #ao fazer grafico podemos ver que nao teriamos a barra do insuficiente
#levels = as.character(1:100)
#sort - ordem alfabetica
#sort(decreasing=True); -> ao contrario
?sort## starting httpd help server ... donez## [1] Suficiente Bom Bom Suficiente
## Levels: Insuficiente Suficiente BomA função table() permite obter um resumo das frequencias absolurtas de cada grupo de um fator. Para usar as frequencias relativas devemos usar prop.table(table(factor))
table(z)## z
## Insuficiente Suficiente Bom
## 0 2 2prop.table(table(z))## z
## Insuficiente Suficiente Bom
## 0.0 0.5 0.5A função tapply permite aplica uma funçao a cada grupo de elementos de um objecto (grupos definidos por um fator). Vai no fundo aplicar uma determinada transformaçao definida por uma funçao tapply(objeto,fator,função). O resultado corresponde à soma dos elementos de x agrupados pelos niveis de z. Aplica uma funçao a um objeto demonstrando o resultado consoante os fatores existentes.
Ex: Grupo etario, peso e queremos calcular media do peso. Podemos calcular media global, ou media dos jovens, adultos e seniors (fatores) de forma separada com o tapply
x=c(10,8,5,4)
tapply(x,z,sum) #soma elementos do vetor x que fazem parte de cada um dos grupos que estao definidos em z -> 10 + 4 para a mesma posiçao em indice do vetor z ([1] e [4] sao suficiente e [2] e [3] sao bom). O resultado da soma dá 14 para suficiente (10+4 das posiçoes 1 e 4) e 13 para bom (8 e 5 nas posiçoes 2 e 3). O tamanho dos vetores tem de ser iguais para comparar os indices## Insuficiente Suficiente Bom
## NA 14 13Manipulaçao de vetores (1d) | length, max, min, sum, prod, sort…
length(x) #para obter o comprimento do vetor x## [1] 4max(x) #ou min(x) para obter o valor máximo e mínimo, respetivamente, do vetor x## [1] 10sum(x) #para obter a soma de todos os elementos de x## [1] 27prod(x) #para obter o produto de todos os elementos de x## [1] 1600sort(x) #para ordenar os elementos de x por ordem crescente## [1] 4 5 8 10sort(x,decreasing=TRUE) #para ordem decrescente﴿## [1] 10 8 5 4Podemos também introduzir missing values usando a expressão NA; o resultado das operações anteriores(exceto length) em vetores com NA é NA
Concatenação e replicação (rep()…)
Criar novos vetores mais compridos, concatenando vetores já existentes usando a função c()
a=c(1,2)
B=3:10
c(a,B)## [1] 1 2 3 4 5 6 7 8 9 10Podemos replicar um determinado vetor usando rep()
n = 2
m = 3
rep(x,n) #ou: ## [1] 10 8 5 4 10 8 5 4rep(x,times=n) #concatena n replicações de x## [1] 10 8 5 4 10 8 5 4rep(x,each=m) #concatena m replicações de cada elemento de x - pega no primeiro elemento do vetor x e replica 4 vezes, pega no segundo elemento e faz o mesmo... "Por cada elemento no vetor X, replicar N vezes)## [1] 10 10 10 8 8 8 5 5 5 4 4 4rep(x,times=n,each=m) #mistura das anteriores## [1] 10 10 10 8 8 8 5 5 5 4 4 4 10 10 10 8 8 8 5 5 5 4 4 4rep(x,length.out=10) #replica x com eventual reciclagem até um comprimento de n. Vetor tem comprimento 3 (original) e faz 3 replicaçoes inteiras e depois so falta uma posiçao ate termos 10 (pega no primeiro elemento do X e completa o que falta)## [1] 10 8 5 4 10 8 5 4 10 8Indexação de vetores: Seleçao | x[indice]
Para selecionar um subconjunto de elementos de um determinado vetor, usamos o comando x[vetor de indices]
A definição do vetor de índices pode ser realizada usando:
- um vetor lógico, que é reciclado se necessário até ao comprimento de x; os elementos de x correspondentes ao valor TRUE do vetor lógico são selecionados e aqueles correspondentes a FALSE são omitidos
x=c(10,-32,40)
tf=c(F,F,T)
x[tf]## [1] 40- um vetor de inteiros positivos, que representam os índices dos elementos a selecionar
x[c(1,3)]## [1] 10 40- um vetor de inteiros negativos, que representam os índices dos elementos a omitir
x[-c(1,3)] #diz que nao queremos escolher os elementos no indice 1 nem no indice 3## [1] -32- um vetor de strings, quando um objeto tem associado um atributo
names
x=c(3,1,4,2,5,6,8,7,9)
x## [1] 3 1 4 2 5 6 8 7 9x[3]## [1] 4x[c(3,5,8)]## [1] 4 5 7x[-c(3,5,8)] #todos menos os elementos do indice 3,5,8## [1] 3 1 2 6 8 9x[x>3] #Vai ao X e retorna quais os que sao superiores a 3## [1] 4 5 6 8 7 9x>3 # Verfica quais deles sao superiores a 3## [1] FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUEx[(x%%2)==0] #retorna os pares cuja divisao da resto 0## [1] 4 2 6 8x[!(x%%2)==0] #retorna os impares - Negaçao## [1] 3 1 5 7 9(x%%2)==0## [1] FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSEModificação
Tambem podemos usar a indexação para modificar elementos de um vetor com x[vetor de indices] = novos valores
Notas: * se novos valores é escalar, então a todos os elementos selecionados de x é atribuído esse valor
x## [1] 3 1 4 2 5 6 8 7 9aux = x[x>3]
x[x>3]=10
x## [1] 3 1 10 2 10 10 10 10 10x[x>3]= aux
x## [1] 3 1 4 2 5 6 8 7 9x[(x%%2)==0]=c(0,-1) #substitui todos os pares por 0 e -1, respetivamente
x## [1] 3 1 0 -1 5 0 -1 7 9se
novos valoresé um vetor com comprimento superior a x, apenas os seus primeiroslength(x[vetor de indices])valores são utilizadosse
novos valoresé um vetor com comprimento inferior a x, então é reciclado atélength(x[vetor de indices])
Matrizes (2d)
Matriz
Podem ser criadas com o comando matrix(x,nrow,ncol) organizando os valores de x por coluna em nrow linhas e ncol colunas (com a opção byrow=TRUE preenche por linha)
A=matrix(1:12,nrow = 4,ncol = 3)
A## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12B=matrix(1:12,nrow=4,ncol=3,byrow=TRUE) #preencher por linha
B## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9
## [4,] 10 11 12A[2:3] ## [1] 2 3B[,2] #todos os elementos da 2º coluna## [1] 2 5 8 11B[4,] #todos os elementos da 4ª linha## [1] 10 11 12B[2:4,2:3]#da linha 2 â 4 e da coluna 2 a 3## [,1] [,2]
## [1,] 5 6
## [2,] 8 9
## [3,] 11 12B[c(1,3),1] #Linha 1 e linha 3 da coluna 1## [1] 1 7dim(A) #obtém a dimensão ﴾nº linhas e nº de colunas﴿ de A## [1] 4 3t(A) #obtém a transposta de A## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 5 6 7 8
## [3,] 9 10 11 12ncol(A); nrow(A) # o número de colunas e linhas, respetivamente, de A## [1] 3## [1] 4A * B #obtém o produto elemento a elemento entre as matrizes. Multiplicaçao elemento a elemento## [,1] [,2] [,3]
## [1,] 1 10 27
## [2,] 8 30 60
## [3,] 21 56 99
## [4,] 40 88 144#A%*%B #obtém o produto matricial usual entre as matrizes. Numero de colunas de A tem que ser igual a numero de linhas da B, por isso é que executando isto dá erro. Temos que fazer a transposta de uma delas para conseguir
#Produto matricial -> 1*1 + 5*2 + 9 * 3 (elemento das linhas de A * elemento das colunas de B e a sua soma)
A%*% t(B)## [,1] [,2] [,3] [,4]
## [1,] 38 83 128 173
## [2,] 44 98 152 206
## [3,] 50 113 176 239
## [4,] 56 128 200 272diag(A) #para obter as entradas da diagonal de A (quando sao quadradas)## [1] 1 6 11k = 4
diag(4) #com k inteiro positivo, cria a matriz identidade de dimensão k. Matriz identidade que só tem 0, com exceçao da diagonal. Mtriz identidade é a matriz neutra (1 é o elemento neutro da multiplicaçao)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1Arrays (3d)
Forma de matriz com varias dimensoes 
Para criar uma array usamos o comando array(x,dim) onde x é um vetor contendo as entidades e dim é um vetor especificando as dimensoes pretendidas (o preenchimento ocorre por coluna)
É efetuada reciclagem no caso de x não conter elementos suficientes para preencher todas as posições do array
x=1:12
array(x,dim=c(2,3,3))## , , 1
##
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 2 4 6
##
## , , 2
##
## [,1] [,2] [,3]
## [1,] 7 9 11
## [2,] 8 10 12
##
## , , 3
##
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 2 4 6Podemos também criar um array a partir de x, forçando a alteração das suas dimensões: dim(x)=c(a,b,c,...) ﴾o vetor x transforma‐se no array﴿
- Neste caso o número de elementos de x tem de ser igual ao produto das dimensões
a×b×c×· · ·﴾não é efetuada reciclagem﴿
#dim(x)=c(2,3,3) #a array tem que ter as mesmas dimensoes de x, neste caso, 12 (na array estamos a tentar passar para dimensao 18)
dim(x)=c(1,6,2)
x## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 2 3 4 5 6
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 7 8 9 10 11 12Indexação de arrays e matrizes
- Para aceder a um elemento de um array A especificamos a sua posição:
A[dim1,dim2,dim3,...]onde dim1, dim2, etc são vetores índice para cada uma das dimensões do array - Se algum dos vetores índice for deixado em branco, então todas as entradas na respetiva dimensão são selecionadas; podem também ser utilizados vetores índice lógicos ou de inteiros negativos com o mesmo objetivo visto para vetores
- Podemos também usar matrizes índice com tantas colunas quanto a dimensão do array a indexar e um número de linhas igual ao número de elementos a selecionar ﴾cada linha na matriz índice contém a posição de um elemento a selecionar﴿
X=1:18
dim(X)=c(3,2,3)
X## , , 1
##
## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6
##
## , , 2
##
## [,1] [,2]
## [1,] 7 10
## [2,] 8 11
## [3,] 9 12
##
## , , 3
##
## [,1] [,2]
## [1,] 13 16
## [2,] 14 17
## [3,] 15 18X[2,1,3]## [1] 14X[c(2,3),1,2]## [1] 8 9X[,,1]## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6X[-1,,3]## [,1] [,2]
## [1,] 14 17
## [2,] 15 18X[(X%%2)==0]## [1] 2 4 6 8 10 12 14 16 18Listas
Coleções ordenadas de objetos, nao necessariamente do mesmo tipo, criadas assim: A = list(nome1=objeto1, nome2=objeto2...)
#Sem atribuir nomes:
L1 = list(1:3,"a",c(TRUE,FALSE,TRUE),c(2.3, 5.9))
str(L1)## List of 4
## $ : int [1:3] 1 2 3
## $ : chr "a"
## $ : logi [1:3] TRUE FALSE TRUE
## $ : num [1:2] 2.3 5.9names(L1)## NULLIndexação de listas
- A cada objeto pode ser associado um nome, pelo qual pode ser acedido usando
A$nomedoobjeto - A indexação pode tambem ser realizada na forma
A[[indice]], onde indice é um inteiro positivo (diferente deA[indice]) - Os elementos de um array dentro da lista podem também ser acedidos:
A[[indice do array]][indice], onde em index se aplicam as mesmas regras de indexação de arrays. [] retorna a lista enquanto que [[]] retorna o conteudo.
L2=list(Obj1=1:3,Obj2="a",Obj3=c(TRUE,FALSE,TRUE),Obj4=c(2.3,5.9))
names(L2)## [1] "Obj1" "Obj2" "Obj3" "Obj4"L2$Obj3## [1] TRUE FALSE TRUEL2[[3]] #podemos chamar pelo nome do elemento ou o indice do elemento## [1] TRUE FALSE TRUEL2[3]## $Obj3
## [1] TRUE FALSE TRUEL2[[3]][2] # entra no objeto 3 e vai buscar o 2º elemento## [1] FALSEL2[3][2]## $<NA>
## NULLModificação e concatenaçao de listas
- Para substituir um determinado elemento de uma lista A
A[[indice]] = novo objetoouA[indice] = list(novo objeto) - Para alterar o nome de um elemento da lista:
names(A)[indice] = novo nome - Para adicionar novas componentes num determinado indice:
A[indice novo] = list(nome=objeto) - Para concatenar listas A,B,C…:
c(A,B,C...)
L2[[2]] = c("a", "b")
L2## $Obj1
## [1] 1 2 3
##
## $Obj2
## [1] "a" "b"
##
## $Obj3
## [1] TRUE FALSE TRUE
##
## $Obj4
## [1] 2.3 5.9names(L2) #nome dos objetos por ordem## [1] "Obj1" "Obj2" "Obj3" "Obj4"names(L2)[3] = "Objeto 3"
L2[5] = list(pi)
names(L2)[5] = "Objeto 5"
str(L2)## List of 5
## $ Obj1 : int [1:3] 1 2 3
## $ Obj2 : chr [1:2] "a" "b"
## $ Objeto 3: logi [1:3] TRUE FALSE TRUE
## $ Obj4 : num [1:2] 2.3 5.9
## $ Objeto 5: num 3.14Data frames
Linhas representam os casos e colunas as variaveis recolhidas
Uma data frame é uma lista em que cada objeto vai representar uma variavel. Nas data frames cada coluna deve ter o mesmo numero de elementos (valores para 150 pacientes, por exemplo)
Colunas devem ter nomes (das variaveis, nas listas chamam-se objetos).
data.frame(nomevar1=valores, nomevar2=valores, ...)
df1 = data.frame(var1 = 1:3, var2 = c("a", "b", "c"), row.names= c("Identificador1", "Identificador2", "Identificador3")) #temos 3 linhas e 2 variaveis. Podemos adicionar como factor em vez de chr
str(df1)## 'data.frame': 3 obs. of 2 variables:
## $ var1: int 1 2 3
## $ var2: chr "a" "b" "c"df1## var1 var2
## Identificador1 1 a
## Identificador2 2 b
## Identificador3 3 c#identificadores de pacientes, por exemplo, podem ser usados no "indice" de cada linha (que nao é uma variavel). Os row names pode ser nomes de genes, identificadores de pacientes, etc###Indexação de data frames Data frame é como se fosse matriz misturada com caracterisitcas de uma lista
df1$var1; df1[,1]## [1] 1 2 3## [1] 1 2 3df1$var2[2]; df1[2,2] #vai a variavel 2 e dentro desta variavel vai buscar o elemento na posiçao 2## [1] "b"## [1] "b"Podemos referenciar diretamente o nome das colunas
#var2[2] #nao resulta
attach(df1) #para entrar no df1 e buscar o var2 podemos usar o $, mas para evitar esta parte fazemos o attach da data frame. Este attach vai abrir o df1 - colunas de df1 estao diretamente acessiveis
var2[2]## [1] "b"#Este metodo nao é o melhor principalmente se tivermos data frames com as mesmas variaveis
detach(df1)Outros comandos de interesse
names(df1) #obtém os nomes das colunas ﴾variáveis﴿ da data frame A## [1] "var1" "var2"colnames(df1) #é semelhante ao anterior## [1] "var1" "var2"rownames(df1) #obtém os nomes das linhas de A## [1] "Identificador1" "Identificador2" "Identificador3"cbind(df1,df1) #concatena as colunas de A com as colunas de B ﴾A e B têm de ter o mesmo número de linhas﴿## var1 var2 var1 var2
## Identificador1 1 a 1 a
## Identificador2 2 b 2 b
## Identificador3 3 c 3 crbind(df1,df1) #concatena as linhas de A com as linhas de B ﴾A e B têm de ter o mesmo número e nome de colunas﴿## var1 var2
## Identificador1 1 a
## Identificador2 2 b
## Identificador3 3 c
## Identificador11 1 a
## Identificador21 2 b
## Identificador31 3 cdf2 = df1 #vai criar indices novos (1...)
# names(df2)[2] = "var3" #vai dar erro pois o numero de variaveis em df1 e df2 vai ser diferente
# rbind(df1,df2)Concatenando arrays | cbind, rbind
cbind(arg1,arg2,arg3,…)- concatena horizontalmente os objetos ﴾do mesmo tipo﴿ arg1, arg2, arg3, etc; para tal, deverão ter o mesmo número de linhasrbind(arg1,arg2,arg3,…)- concatena verticalmente os objetos ﴾do mesmo tipo﴿ arg1, arg2, arg3, etc; para tal, deverão ter o mesmo número de colunas
No caso de alguns argumentos serem vetores com dimensão desadequada, são automaticamente reciclados
Modos dos objetos
Int, logic, char…
Se precisarmos de saber qual o modo podemos perguntar: is.??() e diz se é true ou false (onde ?? = numeric, character…)
Podemos ainda alterar (coerce) o tipo usando as.??() (onde ?? = numeric, character, data.frame…)
is.numeric(df1$var1)## [1] TRUEImportação e Exportação
Podemos usar datasets para testes com data(package = .packages(all.available = TRUE)) e data()
x = data("EuStockMarkets")
x## [1] "EuStockMarkets"Importar Excel
- Requer packages como o
readxle/ouxlsx range(A1:C11)vai escolher as 3 primeiras variaveis e as primeiras 10 colunas (colocamos 11 pois a primeira linha tambem conta que é o header)col_types=c("numeric")para importar as colunas como numericas
?read_excel## No documentation for 'read_excel' in specified packages and libraries:
## you could try '??read_excel'Importar CSV
#swimming_pools = read.csv("D:/swimming_pools.csv")
#swimming_pools = read.csv("D:/swimming_pools.txt", sep=" ") #separa por espaçosImportar from Text
- Com a importação via From Text (readr) é necessario instalar package
readr
# sp_text = read.csv(...)
?read.csv #first column as namesImportar de Bases de Dados | header, sep, quote, dec, read.table…
Podemos usar read.table("ficheiro") ou read.table(file, opção1=valor, opção2=valor,...) destacando-se as seguintes opções:
header=valor lógicoindicando se a primeira linha do ficheiro contém os nomes das variáveissep="separador"indicando qual o separador entre colunas ﴾para separação por espaço, deixar em branco﴿quote="caractere"indicando o caractere envolvente dos valores das variáveisdec="."indicando qual o separador decimal utilizado no ficheirorow.names=vetorpara especificar um vetor com os nomes das linhas ou o número da coluna contendo esses nomescol.names=vetorpara especificar um vetor com os nomes das variáveisnrows=valorpara especificar o número de linhas a importarskip=valorpara especificar o número de linhas a ignorar antes de importar
Exportação | .RData
Para exportar uma base de dados (matriz ou data frame) usamos comando write.table - write.table(objeto,"ficheiro",opções)
Destacam-se as seguintes opções: * append valor lógico para indicar se o objeto é adicionado ao ficheiro * quote valor lógico ou vetor numérico para indicar, no primeiro caso, se as colunas de modo character ou factor são ”quoted”; no segundo caso, o vetor indica as colunas a aplicar quotes * sep para definir o separador entre colunas * dec para definir a separação decimal
#write.table(S1,"LE.txt", sep=":")
#write.table(S1,"LE1.txt", sep=":", row.names=F) #vai retirar os indices de cada linha
#S1
#save.image guarda tudo
#save(S1,file="S1.RData")Estruturas de Controlo e Funções
if statement
Executa instruções dependendo do valor logico de uma expressao
if (expr1) { expr2 }
onde expr1 representa uma condição ﴾lógica﴿ a ser verificada e expr2 é uma instrução ou conjunto de instruções a serem executadas caso expr1 for verdadeira
a=5
if (a>2) {
a+10 #temos que colocar print() para mostrar resultados intermedios
a+50
}## [1] 55a = 1
if (a>2) {
a+10
a+50
}Na execução anterior, verifica‐se que apenas o último resultado é apresentado; caso expr2 contenha várias instruções que importa visualizar o resultado, podemos usar a função print()
Podemos também ter interesse em executar um conjunto de instruções caso expr1 seja falsa
if (expr1) { expr2 } else { expr3 }
onde expr3 é executada se expr1 for falsa
Disjunções, conjunçoes e negaçoes
- expr1 pode conter disjunções, conjunções, negações, etc; em particular, os operadores
&&﴾AND﴿ e||﴾OR﴿ são normalmente usados em vez de&e|, respetivamente, por uma questão de eficiência ﴾evitam a avaliação da segunda expressão caso não seja necessário﴿
a=1
if (a>2 && a<=5) { #por questao de eficiencia como o python se for F no && nao lemos o 2º
print(a+10)
print(a+50)
} else {
print("O valor de a não está entre 2 e 5")
}## [1] "O valor de a não está entre 2 e 5"- Estrutura if permite ainda uma estrutura em escada para testar varias condiçoes
if (expr1) { executar1 } else if (expr2) { executar2 } else if (expr3) { executar3 } else { executar4 }
x = 0
if (x < 0) {
print("Número negativo")
} else if (x > 0) {
print("Número positivo")
} else {
print("Zero")
}## [1] "Zero"
.
x= -2
if (x<=-2) {
print(x**2)
} else if (x > -2 && x <5 ) {
print(exp(x))
} else {
print(x/2)
}## [1] 4for loop
Usadas para executar repetições ﴾ciclos﴿ de um determinado número de instruções: for, while, repeat
Para o ciclo for
for (nome in expr1) { expr2 }
onde nome representa uma variável incrementadora ﴾índice﴿, expr1 é tipicamente uma sequência numérica e expr2 as instruções a executar;
for (i in 1:6) {
print(i)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6#for (i in aux)...while
Para o ciclo while while (expr1) { expr2 } onde expr1 representa uma condição que, enquanto verdadeira, implica a execução das instruções em expr2;
i=0
while (i<=5) {
i=i+1
print(i)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6i = 0
while (i<=3) {
i=runif(1,0,5) # 1 numero aleatorio entre 0 e 5. Semelhante ao sample?
print(i)
}## [1] 1.684002
## [1] 2.883149
## [1] 2.884999
## [1] 4.015532repeat
Para o ciclo repeat:
repeat { expr1 } onde expr1 é o conjunto de instruções a executar
- Não existe qualquer condição de paragem neste ciclo! O ciclo tornar‐se‐á infinito a menos que se utilize o comando
break
i=0
repeat {
i=i+1
if (i==7) {
break
}
print(i)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6Criar funções
A sua estrutura básica é da forma nome = function(arg1,arg2,…) {expr} onde nome é, como o próprio indica, o nome que pretendemos para a função, arg1, arg2 são argumentos ﴾objetos por exemplo﴿ passados para o “interior” da função e expr é o conjunto de instruções a executar
Numa função, os argumentos podem ser dados sob qualquer ordem desde que nomeados, caso contrário devem ser dados pela ordem definida na função
Exemplo: dada a função
funLE = function(x,y,z,w) {expr}
onde x e y são números, z é um valor lógico e w é uma string qualquer, podemos executá‐la das seguintes formas:
- Pela ordem:
funLE(1,2,TRUE,"blue") - Nomeando os argumentos:
funLE(w="blue",x=1,z=TRUE,y=2) - Mistura dos anteriores ﴾os argumentos não nomeados devem vir primeiro, por ordem﴿:
funLE(1,2,w="blue",z=TRUE)
Valores default
É comum ocorrerem argumentos que o utilizador utiliza com maior frequência com um determinado valor
Nesse sentido, é possível estabelecer valores default que poupam o utilizador de os definir sempre que usa a função
Exemplo: suponhamos que na função anterior se pretende estabelecer como default que z=TRUE e w="blue"; então definimos
funLE = function(x,y,z=TRUE,w="blue") {expr}
x e y não têm valores default pelo que têm sempre de ser fornecidos; z e w se não forem fornecidos vão ser considerados como TRUE e “blue”, respetivamente
Assim, são equivalentes as seguintes execuções: funLE(1,2), funLE(1,2,TRUE,"blue") ou funLE(1,2,w="blue"); já funLE(1,2,z=FALSE) altera um dos defaults
Função break
A função break pode ser usada para terminar uma função, mas desde que dentro de um ciclo for, while ou repeat; para outras situações, usar a função stop
funlog = function(x,y,out=TRUE, pl=FALSE) {
if (x+y<=0) {
stop("x+y nao é positivo")
}
z=log(x+y) #so vamos para aqui se o numero for positivo, caso contrario, paramos o if statement pois nao podemos usar numeros negativos no log
if (out) {# mostrar valor de z se for true
print(z)
}
if (pl) { # so faz o plot se for true
plot(x+y,z)
}
}funlog(1,2)## [1] 1.098612funlog(1,2,FALSE)
#funlog(-2,1) - Error in funlog(-2, 1) : x+y nao é positivo
funlog(1,2,pl=TRUE)## [1] 1.098612funlog(1,2,FALSE, TRUE)
#Exercício: altere a função de modo a que x e y possam ser vetores numéricos.Argumento …
Muitas vezes as funções criadas contêm outras funções às quais necessitamos de passar parâmetros
Isso é possível utilizando o argumento … ﴾dot‐dot‐dot ou ellipsis ﴿
fun = function(x,y,z,...) {expr1 outrafunc(arg1,arg2,...) expr2}
onde expr1 e expr2 são instruções a executar e outrafunc é uma outra função a executar com argumentos arg1 e arg2 ﴾eventualmente criados em expr1 ou dados como entrada explícita na função principal﴿
Desta forma o R sabe que quaisquer parâmetros ou objetos passados para além de x,y,z serão destinados a outra função chamada internamente que também deverá conter o argumento …
Semelhante ao spread operator no javascript ou o *args ou **kwargs no python Exemplo:
funLE1 = function(x,y,...) {
z = x+y
plot(x,y)
mean(x,...)
}
funLE1(c(1,2,3,4,5,60),c(3,1,5,4,7,3))## [1] 12.5funLE1(c(1,2,3,4,5,60),c(3,1,5,4,7,3),trim=0.4)
## [1] 3.5O argumento … permite passar argumentos para outras funções que sejam executadas na função funLE1
funLE1(c(1,2,3,4,5,60),c(3,1,5,4,7,3))
funLE1(c(1,2,3,4,5,60),c(3,1,5,4,7,3),trim=0.4)
O argumento trim é assim passado para a função mean
Saidas de uma funçao
As funções retornam apenas o último objeto criado; para retornar vários objetos, constrói‐se ﴾por exemplo﴿ uma lista na última instrução da função
funlog3 = function(x,y,out=TRUE){
aux=c(x+y,x*y,x-y)
if (any(aux<=0)) { # garantir que qualquer uma das operaçoes no vetor é positivo
stop("argumento não positivo")
}
z1=log(x+y);z2=log(x*y);z3=log(x-y)
if (out) {
print(z1)
}
list(soma=z1,produto=z2,diferença=z3) # para exportar as variaveis de forma local (dentro da funçao) para modo global (Acesso global)
}funlog3(2:10,1:9)## [1] 1.098612 1.609438 1.945910 2.197225 2.397895 2.564949 2.708050 2.833213
## [9] 2.944439## $soma
## [1] 1.098612 1.609438 1.945910 2.197225 2.397895 2.564949 2.708050 2.833213
## [9] 2.944439
##
## $produto
## [1] 0.6931472 1.7917595 2.4849066 2.9957323 3.4011974 3.7376696 4.0253517
## [8] 4.2766661 4.4998097
##
## $diferença
## [1] 0 0 0 0 0 0 0 0 0A =funlog3(2:10,1:9) #a vantagem de ter uma lista no final é que podemos aceder as variaveis com o $## [1] 1.098612 1.609438 1.945910 2.197225 2.397895 2.564949 2.708050 2.833213
## [9] 2.944439A$soma## [1] 1.098612 1.609438 1.945910 2.197225 2.397895 2.564949 2.708050 2.833213
## [9] 2.944439A$produto## [1] 0.6931472 1.7917595 2.4849066 2.9957323 3.4011974 3.7376696 4.0253517
## [8] 4.2766661 4.4998097A$diferença## [1] 0 0 0 0 0 0 0 0 02 - Analise Exploratoria de Dados
Censo e analise de dados amostrais
- Quando a recolha e análise se centra em todos os elementos da população de estudo, realiza‐se um Censo
- Quando a recolha e análise se centra numa parte dos elementos da população de estudo, uma amostra, realiza‐se um Análise Estatística dos Dados Amostrais
Parametros vs Medidas amostrais
- Parâmetro ‐ é uma medida numérica que descreve alguma característica de uma população ﴾habitualmente representado por letras gregas﴿:
- Ex: μ ﴾média﴿, σ ﴾desvio padrão﴿, ρ ﴾coeficiente de correlação﴿
- Medida amostral ‐ é uma medida numérica que descreve alguma característica de uma amostra ﴾habitualmente representadas por letras latinas﴿:
- Ex: x¯ ﴾média﴿, s ﴾desvio padrão﴿, r ﴾coeficiente de correlação﴿
Analise exploratoria vs Inferencia
- A Inferência Estatística centra o seu estudo na estimação dos parâmetros desconhecidos de uma população *Para tal, pressupõe a representatividade da amostra face à população
- A inferência é sempre realizada com um certo grau de incerteza que não pode ser eliminado mas pode ser medido e controlado recorrendo à Teoria das Probabilidades
- Algumas técnicas:
- Testes t para amostras independentes e dependentes
- Análises de variância ﴾ANOVA﴿
- Análise fatorial
- Regressão simples e múltipla
Variaveis e Natureza e Escala de uma variavel
Variável ‐ característica mensurável ou atributo de interesse para estudo na população *Ex: presença da doença; idade em que foi diagnosticada a doença
Uma amostra pode conter mais de uma característica de interesse para cada uma das unidades observadas. Com duas características temos uma amostra bivariada; com mais do que duas temos uma amostra multivariada.
- Ex: na população de alunos de 1a matrícula da U. Aveiro em 2012/13 pode interessar a análise de várias características em simultâneo como, por exemplo,“Nota de ingresso” e “Região de proveniência”
Quantitativa | Discretas ou Continuas
- Variável Quantitativa ‐ assume valores intrinsecamente numéricos; associadas a medições ou contagens; distância quantificável entre valores da variável
- Ex. A Maria pesa 50 Kg e o António 100 Kg. O António pesa duas vezes mais do que a Maria.
Podemos definir dois grupos de variáveis quantitativas:
- Discretas ‐ assumem um conjunto finito ou infinito numerável ﴾que se consegue contar﴿ de possibilidades; probabilidade alta de haver valores repetidos; nas discretas ha possivel repetiçao dos mesmos valores, enquanto que as continuas ja nao ha essa mesma probabilidade
- Ex: Número de peças defeituosas por lote; número de chamadas telefónicas; número de crianças
- Contínuas ‐ assumem valores em intervalos reais; probabilidade baixa de haver valores repetidos
- Ex: Comprimentos, pesos, tempos, pressão arterial
Qualitativa
- Podemos ter variaveis dummy aqui
- Variável Qualitativa ﴾ou categórica﴿ ‐ assume valores intrinsecamente não‐numéricos; associadas à classificação em níveis, estados ou classes; não existência de distância ou “distância” não quantificável entre valores
- Ex: Considere‐se a variável “cor preferida”, cujos valores possíveis são azul, preto e branco. Não existe qualquer distância ou relação entre azul e preto
- Outros exemplos: Região de proveniência de um produto agrícola, grupo etário, hipertensão ﴾presente/ausente﴿, tipo de tratamento
Escala de razoes e escala de intervalos
- Escala de razões ‐ os valores numéricos possuem ordem e as diferenças têm significado. O valor zero representa ausência da característica ﴾zero absoluto﴿ e múltiplos de valores possuem significado. É sempre constante independentemente da escala.
- Ex: Medidas de comprimento, pesos, intervalos de tempo, nº de filhos de um casal
- Escala de intervalos ‐ os valores numéricos possuem ordem e diferenças têm significado. O valor zero não representa ausência da característica ﴾zero arbitrário﴿
- Ex: Temperatura medida em graus C ou F
Escala ordinal e escala nominal
- Escala Ordinal ‐ os valores pertencem a uma de várias categorias existentes e possuem uma ordem intrínseca. Já não é sempre constante, depende da escala.
- Ex: Classificações obtidas no 2º e 3º ciclos do ensino básico ﴾1 a 5﴿, grupos etários ﴾crianças, jovens, adultos e idosos﴿, estado de saúde do doente ﴾muito bom, bom, razoável, fraco e péssimo﴿
- Escala Nominal ‐ os valores pertencem a uma de várias categorias distintas e não possuem uma ordem intrínseca
- Ex: Grupo sanguíneo ﴾O, A, B e AB﴿, Género ﴾0 ‐ Feminino, 1 ‐ Masculino﴿
Tabela de frequencias
A frequência de um determinado valor é uma medida do número de vezes que esse valor ocorre
Se uma variável tem um grande número de diferentes valores ﴾quantitativa contínua﴿ então convém agrupar os dados em classes
Em geral, são calculadas as:
- Frequências absolutas ﴾acumuladas ou não acumuladas﴿
- Frequências relativas ﴾acumuladas ou não acumuladas﴿
O cálculo de frequências acumuladas ﴾Ni e Fi﴿ pressupõe a existência de ordem nos valores, pelo que o seu cálculo não se aplica a variáveis nominais
Nao faz sentido calcular as frequencias acumuladas numa nominal (preto, azul, amarelo.. ordem nao interessa) pois estas acumuladas pressupoe que ha uma ordem (ordinais).
Para variaveis nao continuas
Para variáveis não contínuas ﴾discretas ou qualitativas, portanto﴿ podemos usar uma combinação dos comandos
table,cumsumeprop.tableprop.tablepega nas frequencias absolutas e vai calcular em proporçao %cumsumfaz a soma acumuladaFreq_rel_acum=cumsum(prop.table(tbl))soma das frequencias relativas acumuladasFreq_acum=cumsum(tbl)guardar na variavel como a somaFreq_rel=prop.table(tbl)calcular frequencias relativas
load("C:/Users/Rafael/Desktop/main/University/BioinformaticaClinica/1Semestre/LE/Material 2021/Scripts/LE-Fichas/DQ.RData")
tbl = table(dados$Escolaridade)
tbl##
## Primaria Escolaridade obrigatoria Liceu e Curso Superior
## 17 60 23cbind(Freq=tbl,
Freq_acum=cumsum(tbl),
Freq_rel=prop.table(tbl),
Freq_rel_acum=cumsum(prop.table(tbl)))## Freq Freq_acum Freq_rel Freq_rel_acum
## Primaria 17 17 0.17 0.17
## Escolaridade obrigatoria 60 77 0.60 0.77
## Liceu e Curso Superior 23 100 0.23 1.00Para variaveis quantitativas continuas
- Aqui definimos intervalos de valor e contabilizamos quantos caem nesse intervalo de valores
.
Para variáveis contínuas começamos por definir as classes e tudo o resto é semelhante
#Para as continuas primeiro criamos as classes
classes=table(cut(dados$PAS,
breaks=nclass.Sturges(dados$PAS))) #numero de classes baseado na variavel dados$PAS
classes##
## (110,114] (114,118] (118,122] (122,126] (126,129] (129,133] (133,137] (137,141]
## 1 0 1 3 7 25 41 22cbind(Freq=classes,Freq.acum=cumsum(classes), Freq.rel=prop.table(classes), Freq.rel.acum=cumsum(prop.table(classes)))## Freq Freq.acum Freq.rel Freq.rel.acum
## (110,114] 1 1 0.01 0.01
## (114,118] 0 1 0.00 0.01
## (118,122] 1 2 0.01 0.02
## (122,126] 3 5 0.03 0.05
## (126,129] 7 12 0.07 0.12
## (129,133] 25 37 0.25 0.37
## (133,137] 41 78 0.41 0.78
## (137,141] 22 100 0.22 1.00#freq.rel.acum -> 12% tem uma PAS menor ou igual a 129- Podemos definir manualmente os intervalos usando
breaks=c(a,b,c,…)onde a,b,c,… são os extremos dos intervalos
Para 2+ variaveis qualitativas | Tabelas de dupla entrada
Também podemos fazer o cruzamento de duas ou mais variáveis qualitativas, numa tabela de dupla entrada
table(dados$Gr_etario,dados$Escolaridade)##
## Primaria Escolaridade obrigatoria Liceu e Curso Superior
## Ate aos 34 anos 3 9 4
## 35 - 64 anos 12 40 16
## 65 e mais anos 2 11 3prop.table(table(dados$Gr_etario,dados$Escolaridade))##
## Primaria Escolaridade obrigatoria Liceu e Curso Superior
## Ate aos 34 anos 0.03 0.09 0.04
## 35 - 64 anos 0.12 0.40 0.16
## 65 e mais anos 0.02 0.11 0.03ftable(xtabs(~Gr_etario+Escolaridade,data=dados))## Escolaridade Primaria Escolaridade obrigatoria Liceu e Curso Superior
## Gr_etario
## Ate aos 34 anos 3 9 4
## 35 - 64 anos 12 40 16
## 65 e mais anos 2 11 3+## Medidas Amostrais ### Media vs Mediana vs Media aparada
Podemos agrupar as medidas amostrais em três grupos * Medidas de Localização Central * Central: Média, Mediana, Moda, Média Aparada
- Média ﴾x¯﴿ ‐ A média pode ser pensada como o centro de massa dos valores das observações, i.e., o ponto de equilíbrio após dispormos as observações sobre uma régua
- Mediana ﴾Me﴿ ‐ é a observação central da amostra ordenada
- pelo menos 50% dos dados são inferiores ou iguais a Me
- pelo menos 50% dos dados são superiores ou iguais a Me
- Média aparada ‐ Uma média aparada a α% é a média calculada quando se exclui α% das observações mais extremas ou seja, α/2 % em cada extremo da amostra ﴾geralmente α é um valor pequeno﴿
- Moda ﴾Mo﴿ ‐ é o valor mais frequente de uma amostra; ao contrário do que acontece com a mediana e a média, uma amostra pode possuir mais do que uma moda
Medidas de localizaçao relativa
Relativa: Mínimo, Máximo, Quantis
- Mínimo ‐ é o valor mais reduzido da amostra
- Máximo ‐ é o valor mais elevado da amostra
- Quantil de ordem p ﴾0 ≤ p ≤ 1﴿ ‐ é um valor, xp, que divide a amostra ordenada em duas partes, tal que pelo menos p × 100% das observações da amostra são menores ou iguais a xp e pelo menos (1 − p) × 100% das observações são maiores ou iguais a xp
- Consoante o valor de p alguns quantis recebem designações próprias. Temos o
- Quartis ﴾p=1/4, 2/4, 3/4﴿
- Decis ﴾p=1/10, 2/10,…, 9/10﴿
- Percentis ﴾p=1/100, 2/100,…, 9/100﴿
- Os quartis são habitualmente denotados por Q1, Q2 e Q3. Note que Q2 coincide com a mediana.
- Consoante o valor de p alguns quantis recebem designações próprias. Temos o
O R possui vários comandos de base para obter as medidas de localização de uma variável x
mean(x)oumean(x,trim=valor)para obter a média ou média aparada respetivamente ﴾valor varia entre 0 e 0.5﴿median(x)` para obter a medianamax(x)emin(x)para obter os valores máximo e mínimo, respetivamentequantile(x)ouquantile(x,probs=vetor)para obter os quantis ﴾vetor contém valores em [0,1]﴿summary(x)para obter um resumo ﴾simplificado﴿ de todas estas medidasnames(which(x == max(x)))para obter a moda de uma tabela usandorep()- exercicio 3,e,i) ficha 2
Não existe um comando imediato para obter a moda; pode ser identificada analisando a tabela de frequências
Muitas vezes interessa calcular determinadas medidas por grupos de indivíduos; para tal podemos usar a função tapply ou by:
quantile(dados$Peso, probs=c(0.1,0.3,0.8,0.95))## 10% 30% 80% 95%
## 64.003 66.788 73.810 77.483tapply(dados$PAS,dados$Gr_etario, mean) #no tapply podemos definir as nossas proprias funçoes em vez de ter as pre-estabelecidas como mean, summary, etc ## Ate aos 34 anos 35 - 64 anos 65 e mais anos
## 127.1875 134.7500 138.7500by(dados$PAS,dados$Gr_etario, mean)## dados$Gr_etario: Ate aos 34 anos
## [1] 127.1875
## ------------------------------------------------------------
## dados$Gr_etario: 35 - 64 anos
## [1] 134.75
## ------------------------------------------------------------
## dados$Gr_etario: 65 e mais anos
## [1] 138.75by(dados$PAS,dados$Gr_etario, median)## dados$Gr_etario: Ate aos 34 anos
## [1] 129
## ------------------------------------------------------------
## dados$Gr_etario: 35 - 64 anos
## [1] 135
## ------------------------------------------------------------
## dados$Gr_etario: 65 e mais anos
## [1] 139by(dados$PAS,dados$Gr_etario, summary)## dados$Gr_etario: Ate aos 34 anos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 110.0 126.5 129.0 127.2 130.0 134.0
## ------------------------------------------------------------
## dados$Gr_etario: 35 - 64 anos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 125.0 133.0 135.0 134.8 136.0 140.0
## ------------------------------------------------------------
## dados$Gr_etario: 65 e mais anos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 135.0 138.0 139.0 138.8 139.2 141.0Medidas de dispersao
- Medidas de Dispersão
- Amplitude, Distância Inter‐Quartil, Variância, Desvio padrão, Desvio absoluto mediano
- Amplitude ‐ é a diferença entre o valor máximo e o valor mínimo da amostra.
- Distância inter‐quartil ‐ é a diferença entre o 3º e o 1º quartis ﴾H = Q3 − Q1﴿
- Variância ‐ é a média dos quadrados dos desvios das observações em relação à média da amostra
- A variância não vem representada na mesma unidade das observações. Se tomarmos a raiz quadrada da variância obtemos o desvio padrão que também é uma medida de dispersão e vem na mesma unidade das observações
- Desvio absoluto mediano (mad - median absolute…) ‐ constitui uma alternativa ao desvio padrão, mais robusta a outliers
- Amplitude, Distância Inter‐Quartil, Variância, Desvio padrão, Desvio absoluto mediano
- As medidas de dispersão podem ser obtidas com os seguintes comandos
max(x)-min(x)para obter a amplitudeIQR(x)para obter a distância inter‐quartilsd(x)para obter o desvio‐padrãovar(x)ousd(x)^2para obter a variânciamad(x)para obter o desvio absoluto mediano - Comparar o sd com o mad pode ajudar para ver se temos valores extremos (outliers)
Medidas de forma
- A forma da distribuição de valores da amostra pode ser avaliada pela sua simetria ou assimetria. Diz‐se que uma distribuição possui assimetria positiva ﴾assimetria negativa﴿ quando existe uma concentração de valores na zona de valores mais reduzidos ﴾elevados﴿ da amostra.
Coeficiente de assimetria
- O coeficiente de assimetria ﴾skewness﴿ é uma medida amostral que assume o valor zero quando a distribuição de frequências da amostra é completamente simétrica e assume valores diferentes de zero ﴾positivos ou negativos, resp.﴿ quando a distribuição não é simétrica ﴾assimétrica positiva ou assimétrica negativa, resp.﴿
- O coeficiente de assimetria pode ser obtido com a função
skewnessdo pacotemoments
Relaçao entre 2 variaveis: correlaçao
Medida que permite analisar o grau de relaçao (associaçao) linear entre 2 variaveis to coeficiente de correlaçao linear (coeficientes de correlaçao variam de -1 e 1)
Um valor negativo ﴾positivo﴿ indica que a relação é negativa ﴾positiva﴿, ou seja, com o aumento dos valores de uma variável verifica‐se uma ﴾um﴿ diminuição ﴾aumento﴿ dos valores da outra
Corr.=1 ou Corr.=‐1 significa correlação perfeita positiva ou negativa, respetivamente
O coeficiente de correlação de Pearson é calculado para duas variáveis medidas numa escala de intervalos ou razões ﴾Quantitativa vs Quantitativa﴿
O coeficiente linear de Spearman, RS, aplica‐se no caso de duas variáveis medidas pelo menos numa escala ordinal ﴾adequado para o caso Ordinal vs Ordinal﴿
O coeficiente de Kendall é uma alternativa ao coeficiente de Spearman quando existem empates entre as ordens das observações
No caso de uma amostra multivariada, podemos calcular uma matriz de correlações, constituída por todos os pares de correlações entre as várias variáveis
cor(dados$PAS,dados$PAD,method="pearson")## [1] 0.9691827cor(dados$PAS,dados$PAD,method="spearman")## [1] 0.946051cor(dados$PAS,dados$PAD,method="kendall")## [1] 0.855748# Matriz de correlações (apenas as variáveis numéricas)
View(cor(dados[,sapply(dados,is.numeric)])) #so faz correlaçao de variaveis numericas logo neste codigo estamos a fazer a correlaçao de todas as colunas de dados cujas colunas sao numericsa (sapply vai a data frame dados e verfica se é numerica). Como estamos a fazer indexaçao so vai considerar aquelas que tem true (numericas)
#Processo devia estar como row.names (index) e nao como variavel
#Correlaçao negativa- quanto maior a altura, menor o IMC Correlaçao vs Causalidade
- Encontrar uma correlação forte entre as variáveis não significa necessariamente um efeito de causalidade ﴾e a existir, não indica qual é a direção dessa causalidade﴿
- A correlaçao indica apenas o grau de associaçao linear entre as duas variaveis
- Duas variaveis podem ter correlaçao forte devido a uma 3º variavel (confounding variable)
- Note-se tambem que a correlaçao nula nao significa que nao haja relaçao entre as variaveis (apenas nao se encontra associaçao linear)
Representaçoes graficas
Grafico de barras
- Representadas proporçoes, contagens ou percentagens
- Uteis para variaveis discretas e variaveis qualitativas (ordinais e nominais)
- podem ser obtidos com o comando
barplot(), geralmente associado ao comandotable: barplot(table(x));também complot(f)onde f é um fator
barplot(table(dados$Escolaridade))
plot(dados$Escolaridade)
Histograma
- Usado para variaveis continuas
- Divisao da amplitude das observaçoes em k subintervalos (classes) iguais (usando, por exemplo, regra de sturges)
- As barras representam a frequencia (asboluta ou relativa) em cada classe
- Usamos função
hist(x)- Destaca-se a opção
breakspara fornecer informação sobre a divisão em classes ﴾pode ser um vetor com os pontos de divisão, uma função para calcular os pontos de divisão, uma estimativa para o número de divisões, etc﴿
- Destaca-se a opção
hist(dados$PAS)
hist(dados$PAS,breaks = nclass.Sturges(dados$PAS))
hist(dados$PAS,breaks = 15)
O que procurar?
Outliers ‐ Observações discrepantes ou extremas * devidas a erros de leitura ou medição ﴾corrigir se possível﴿ * valores reais ﴾nada a fazer?﴿ * Estas observações causam problemas! Enviesamento de futuras análises
Multimodalidade ‐ “picos” mais salientes; várias modas ou classes modais são indicadores de possíveis grupos distintos nos dados 
Caixa de bigodes
Usado principalmente nas quantitativas (continuas e discretas).
- Ex: idade em funçao do grupo etario, etc As caixas de bigodes ﴾boxplot﴿ correspondem a uma representação gráfica que dá informação sobre
Localização central: mediana
Localizações relativas: 1º e 3º quartis e mínimo e máximo.
Dispersão: amplitude e distância inter‐quartil
Assimetria: posição relativa da mediana na caixa e o comprimento dos bigodes
Uteis para comparar varias amostras (mesma variavel medida em diferentes grupos)
No R colocamos boxplot(x) Para caixa de bigodes comparativas boxplot(x,y,z...) Se alguma das variaveis é um fator (definindo por isso, grupos) entao boxplot(x ~ grp) ou plot(grp,x) tambem gera um boxplot comparativo (x é a variavel a analisar e grp é o fator)
boxplot(dados$PAS)
boxplot(dados$PAS,dados$PAD)
boxplot(dados$PAS~dados$Gr_etario)
plot(dados$Gr_etario, dados$PAS)
Diagramas de dispersão
O diagrama de dispersão é uma representação gráfica, dos pares de pontos (x, y) de uma amostra bivariada. Este tipo de gráficos permitem relacionar duas variáveis entre si.
Diagrama de dispersao entre as variaveis x e y
plot(x,y)NOTA sao este tipo de graficos que fazemos para ver se ha uma relaçao linear entre as 2 variaveis quantitativas (ver se faz sentido fazer uma regressao, etc)
plot(dados$PAS,dados$PAD)
Grafico para matriz de correlaçao
Sendo dificil a analise de uma matriz de correlaçao, podemos optar por analise grafica
library(corrplot) #apenas para testar correlaçoes lineares## corrplot 0.90 loadeddados1=dados[,-1] # Para excluir a variável Processo
corrplot(cor(dados1[,sapply(dados1,is.numeric)])) # passamos a matriz de correlaçao e depois temos uma funçao sapply que aplica is.numeric coluna a coluna. Escolhemos apenas as colunas de dados1 que sao numericas para fazer as correlaçoes dados1[row, column]. em que queremos dados[todas as linhas, das colunas que sao numericas] - é por isso que deixamos as linhas vazias
Analise de dados longitudinais | Medidas historicas
Quando os dados são recolhidos sobre os mesmos sujeitos em diferentes momentos temporais ﴾ao longo do tempo﴿ dizemos que estamos perante um problema de dados longitudinais ﴾ou medidas repetidas﴿
Geralmente neste tipo de dados pretende‐se analisar a evolução de uma qualquer característica com o tempo, eventualmente sujeita a um qualquer “tratamento”
Desta forma, os dados podem ser vistos como séries temporais para cada sujeito e a visualização pode ser feita por sujeito, por grupo, por valor médio, entre outras
Em geral os dados podem ser apresentados em dois formatos: formato “long” e formato “wide”
- No formato “long” as medidas repetidas constituem diferentes linhas da matriz de dados, enquanto que no formato “wide” cada sujeito é representado por uma linha e as medidas repetidas são diferentes colunas
#install.packages("HSAUR3")
library(HSAUR3)## Warning: package 'HSAUR3' was built under R version 4.0.5## Loading required package: toolsdata("epilepsy")
head(epilepsy,10) # Estes dados estão no formato long - pacote lattice pressupoe que os dados estao em formato long## treatment base age seizure.rate period subject
## 1 placebo 11 31 5 1 1
## 110 placebo 11 31 3 2 1
## 112 placebo 11 31 3 3 1
## 114 placebo 11 31 3 4 1
## 2 placebo 11 30 3 1 2
## 210 placebo 11 30 5 2 2
## 212 placebo 11 30 3 3 2
## 214 placebo 11 30 3 4 2
## 3 placebo 6 25 2 1 3
## 310 placebo 6 25 4 2 3- O pacote
latticedispõe de ferramentas para visualização destes dados ﴾assumindo o formato “long”﴿ nos chamados Trellis plots
#install.packages("lattice")
library(lattice)## Warning: package 'lattice' was built under R version 4.0.5xyplot(seizure.rate/base ~ period, data = epilepsy) #seizure.rate/base para todos os sujeitos para cada um dos periodos (59 sujeitos para periodo 1, 59 para periodo 2... mas temos mistura com os de placebo e tratamento e tambem nao sabemos a ligaçao entre sujeitos por cada um dos periodos. A bolinha no topo do periodo 1 nao é a mesma que a bolinha no topo do periodo 2.
xyplot(seizure.rate/base ~ period | treatment, data = epilepsy) #Este permite distinguir por tratamento. Os casos extremos (pessoas com muitas seizures) parecem ter reduzido. A concentraçao das bolinhas no progabide tambem estao mais concentradas em baixo (mais cor azul).
xyplot(seizure.rate/base ~ period | treatment, group = subject, data = epilepsy, type = "b") # separa por sujeito para associar os sujeitos durante os periodos. A tendencia das linhas tambem parece reduzir (linhas mais compactas)
epilepsy$gretario=cut(epilepsy$age,breaks =c(0,20,30,40,50) ) #vamos pegar na idade e dividmos em niveis/grupo etario. Fazemos o cut com os niveis que definimos
xyplot(seizure.rate/base ~ period | treatment+gretario, group = subject, data = epilepsy, type = "b") #separa todos os periodos por tratamento e grupo etario
Funçao plot()
- É a função ﴾para gráficos﴿ mais utilizada e talvez mais versátil do R; é designada de função genérica pois adequa o seu comportamento ao tipo ou classe dos seus argumentos
- Vários tipos de gráficos podem ser obtidos com esta função
plot(x)ouplot(x,type="l")para um gráfico do tipo “índice vs valores de x” ﴾estilo série temporal﴿plot(x,y)para um diagrama de dispersão entre as variáveis numericas x e yplot(f)gera um gráfico de barras para o fator fplot(f,x)gera caixas de bigodes dos valores de x para cada grupo do fator f
- A função possui vários argumentos para definir o resultado final do gráfico, podendo ainda passar parâmetros gráficos ﴾ver no seguimento
par﴿type="p"﴾pode ser “l”,“b”,etc, ver ?plot﴿main="Titulo"sub="Subtítulo"xlab="Nome var. x"ylab="Nome var. y"asp=número﴾relação de escala y/x, aspect ratio﴿ X de 10 em 10 e y de 5 em 5
plot(dados$Colesterol,dados$Peso,main="Relação Peso - Colesterol", xlab="Colesterol",ylab="Peso",asp=2)
Graficos multivariados
- Se
Aé uma data frame entãoplot(A)gera uma matriz de diagramas de dispersão entre cada par de variáveis da data frame Aplot(∼ a + b + …, data=A)gera uma matriz de diagramas de dispersão entre as variáveis de A com nome a, b, …plot(a ∼ b + c + …,data=A)gera gráficos de dispersão individualizados entre a variável a e as variáveis b, c, …
plot(~ Idade + Altura + Peso + IMC, data=dados)
#Graficos de dispersao entre uma variavel e varias outras de forma individualizada
plot(Idade ~ Peso + IMC, data=dados)
- Se A é uma matriz, então plot(A) produz o diagrama de dispersão das duas primeiras colunas de A
- Para obter uma matriz de diagramas de dispersão de uma matriz ou data frame A também podemos usar
pairs(A)
Considerando os dados DQ.RData
#load("DQ.RData")
aux=sapply(dados,is.numeric) # Verificar que vars. são numéricas
plot(dados[,aux]) #plot todas as linhas e apenas variaveis numericas
plot(~ Colesterol+FreqCard+Peso,data=dados)
plot(Colesterol ~ FreqCard+Peso,data=dados)
pairs(dados[,aux])
- A função
coplotpode também ser usada para gerar gráficos multivariados, neste caso condicionais; permitem relacionar 3 ou 4 variáveis da seguinte forma ﴾particularmente útil quando c e d são fatores﴿. Semelhantes ao lattice, mas nao pressupoe evoluçao temporal. Aqui é uma so medidacoplot(a∼b| c, data=A)para uma matriz de diagramas de dispersão entre as variáveis a e b para cada grupo do fator c. a em funçao de b condicionada em c. Diagrama de dispersao de a para b por cada grupo do fator ccoplot(a∼b| c+d, data=A)para uma matriz de diagramas de dispersão entre as variáveis a e b para cada combinação de grupos dos fatores c e d
coplot(FreqCard ~ Colesterol|Tabaco,data=dados) # analisar da esquerda para a direita, de baixo para cima (o grafico no topo sozinho corresponde ao fumador) - Fumadores tem freq cariacas mais elevadas e colesterol mais elevado. 
coplot(FreqCard ~ Colesterol|Tabaco+Sexo,data=dados) 
- Para uma seleção mais completa de gráficos multivariados, consultar, por exemplo, o pacote
lattice
Comandos de baixo nivel
- São vários os comandos de baixo nível que permitem adicionar informação a um gráfico previamente criado
points(x,y)para adicionar pontos cujas coordenadas são especificadas por x e ylines(x,y)para adicionar pontos ligados por linhas com coordenadas especificadas por x e ytext(x,y,labels)para adicionar texto contido emlabelsnas coordenadas (x,y)abline(a,b)para adicionar a reta de equação y = a + bxpolygon(x,y)para adicionar um polígono cujos vértices ordenados estão em (x,y)legend(x,y,legend=leg,col=vetor1,pch=vetor2,lty=vetor3,…)para adicionar uma legenda, onde(x,y)é a posição do canto superior esquerdo,legé um vetor, tipicamente de caracteres com as etiquetas,vetor1é um vetor que especifica as cores,vetor2é um vetor para especificar o símbolo a marcar cada ponto,vetor3especifica o tipo de linha, etctitle("titulo principal","subtitulo")para adicionar titulo ﴾no topo﴿ e subtítulo ﴾em baixo﴿. Podemos usarexpressionepastepara concatenar strings, variaveis e notaçao matematica.
- Nos comandos que exigem localizações, o par
(x,y)pode ser substituído porlocator(n)ondené o número de inserções a realizar; esta opção permite definir as localizações usando o rato
plot(dados$Colesterol,dados$Peso)
title("Relação Peso - Colesterol")
Col.medio=mean(dados$Colesterol)
Peso.medio=mean(dados$Peso)
points(x=Col.medio,y=Peso.medio,col="blue",pch="+") #fez o ponto medio entre colesterol e peso. Se passasemos uma reta de regressao ela passava por este ponto medio
#text(locator(1),"Sem relação aparente")
x=seq(-5,5,by=0.1)
f1=x^2; f2=exp(x)
plot(x,f1,col="blue",type="l") #plot é uma funçao de alto nivel logo cria um novo grafico sempre que executado. É por isso que usamos o "lines" na proxima linha pois queremos adicionar essa parte do grafico ao grafico pre-existente feito pelo "plot"
lines(x,f2,col="red",type="l")
legend("bottomright",legend=c(expression(paste(x^2)),expression(paste(exp(x)))), col=c("blue","red"),lty=c(1,1)) #lty (line type) define as linhas que vamos definir
Parametros graficos
Cada graphics device tem os seus defaults que podem ser acedidos com o comando
par()﴾veja‐se ?par﴿Os parâmetros gráficos pode ser alterados de forma temporária ou de forma permanente
- permanentemente se alterados com a função
par - temporariamente e apenas para o gráfico em questão quando passados dentro de um comando de alto ou baixo nível ﴾como o
plot﴿
- permanentemente se alterados com a função
par(bg="darkolivegreen1")
hist(dados$Peso,main="Distribuição de Peso",xlab="Peso", ylab="Frequência",col="aquamarine",border="chocolate1")
plot(dados$Colesterol,dados$Peso)
par(bg="white")
plot(dados$Colesterol,dados$Peso, col=dados$Escolaridade,cex=3,pch="+")
plot(dados$Colesterol,dados$Peso, col=dados$Escolaridade, pch=as.numeric(dados$Tabaco)) #em cada nivel, coloca uma cor diferente. Escolaridade temos diferentes cores e tabaco temos diferentes simbolos. Este "as.numeric" so faz sentido se colocarmos como fator. Se tivessemos as.character teriamos de passar para fator primeiro - pch= as.character(as.factor(c("a","b","c")))
#cor por escolaridade
plot(dados$Colesterol,dados$Peso, col=dados$Escolaridade, pch=as.character(dados$Tabaco)) #neste caso vamos ter diferentes cores para escolaridade e os dados do tabaco vao ser atribuidos a 1ª letra dos fatores 
Figuras multiplas
- Para criar criar a matriz n × m onde “cairão” os gráficos alteramos o parâmetro gráfico
mfrowoumfcolcom a funçãopar(mfrow=c(n,m)); a diferença reside no preenchimento, que no primeiro caso é por linha e no segundo por coluna- A janela de gráficos será dividida ﴾invisivelmente﴿ em n × m células bastando, para preenchê‐las, gerar n × m gráficos
par(mfrow=c(1,3))
plot(dados$Colesterol,dados$Peso)
boxplot(Colesterol~Tabaco,data=dados)
hist(dados$Idade)
par(mfrow=c(1,1))ggplot2
Visão geral
- Camadas ou comandos podem categorizar-se da seguinte forma:
- Esteticos: comandos que regulam a forma como as variaveis sao mapeadas para objetos visuais no grafico
aes()-> o que queremos desenhar - Geometricos: comandos que controlam o tipo de grafico e os objetos geometricos no grafico
geom-> como queremos desenhar
- Esteticos: comandos que regulam a forma como as variaveis sao mapeadas para objetos visuais no grafico
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.0.5data("iris")
p1 = ggplot(data = iris, aes(x = Petal.Length, y = Petal.Width)) #gera grafico
p1
#Diagrama de dispersao entre as 2 variaveis
p2 = p1 + geom_point()
p2
# Nestes dados as observaçoes correspondem a 1 de 3 especies: iris versicolor, iris virginica ou iris setosa
p3 = p1 +geom_point(aes(color=Species))
p3
p4 = p3 + xlab("Comprimento da petala (cm)") + ylab("largura da petala (cm)") + ggtitle("Comprimento da petala vs largura da petala")
p4
#Alterar o fundo, eixos, etc
a1 = ggplot(data=iris, aes(x=Petal.Length, y = Petal.Width))
a1+theme_classic()
a1+theme_dark()
#install.packages("ggthemes")
library(ggthemes)## Warning: package 'ggthemes' was built under R version 4.0.5a1+theme_economist()
a1+theme_excel_new()
a1+theme_solarized()
Grafico de barras
Usados para analisar frequencias absolutas ou relativas de uma ou mais variaveis qualitativas
Apresentar sumario (media) de uma variavel quantitativa por grupos definidos por uma variavel qualitativa
Funçao usada é
geom_bar
library(ggplot2)
#Com frequencias absolutas
#labs permite estabelecer titulos e legendas
ggplot(data=dados, aes(x=Gr_etario))+geom_bar()+labs(x="Grupo etario", y="Frequencia",title="Distribuiçao do grupo etario")
#Com frequencias relativas. geom_bar por defeito faz ..count..; Se fizessemos y = count ia procurar variavel count, mas isso nao existe. Os ".." é uma funçao built in de geom_bar. Se queremos a proporçao fazemos o count pela soma do count (contagem de cada barra a dividir pelo total.)
ggplot(data=dados, aes(x=Gr_etario, y=..count../sum(..count..)))+geom_bar()+labs(x="Grupo etario", y = "Percentagem", title="Distribuiçao do grupo etario")+scale_y_continuous(labels=scales::percent) #scale_y_continous para colocar o eixo do Y como percentagem. Na escala do y queremos que os labels venham em percentagem
#Barras podem ser ordenadas por ordem crescente, decresente e apresentarem uma orientaçao horizontal
tbl=as.data.frame(table(dados$Gr_etario))
names(tbl)[1]="GrupoEtario"
#Tivemos que criar uma contagem dos casos (dataframe) e fazer uma reordenaçao dos valores
ggplot(data=tbl, aes(x=reorder(GrupoEtario,Freq), y=Freq)) + geom_bar(stat="identity")+labs(x="Grupo etario", y="Frequencia", title="Distribuiça do grupo etario")+coord_flip() #alterar os eixos para ficar a horizontal
#reorder vai organizar o grupo etario a partir de um sort pelo valor de frequencia por ordem decrescente. No x vai colocar o grupo etario
#aes(x=reorder(GrupoEtario.sort(Freq.decreasing=T))
#Stat = identity -> contagem ja foi feita pelo table. A estatistica a usar deve ser identity (para usar os valores que estao em data=tbl)
#Depois de fazer coord_flip com barras horizontais, colocavamos Freq.decreasing = F e a barra maior passava para a primeira
#ggplot(data=tbl, aes(x=reorder(GrupoEtario,sort(Freq.decreasing=F)), y=sort(Freq.decreasing=F))) + geom_bar(stat="identity")+labs(x="Grupo etario", y="Frequencia", title="Distribuiça do grupo etario")+coord_flip()- Podemos ainda comparar duas variaveis simultaneamente; no caso de duas variaveis qualitativas, existem diferentes abordagens
ggplot(data=dados, aes(x=Gr_etario, fill=Tratamento)) + geom_bar(position="stack")
ggplot(data=dados, aes(x=Gr_etario, fill=Tratamento)) + geom_bar(position="dodge")
ggplot(data=dados, aes(x=Gr_etario, fill=Tratamento)) + geom_bar(position="fill") #coloca cada barra para 100% e permite analisar proporçoes 
- No caso de uma quantitativa e qualitativa, é comum descrever-se uma determinada medida amostral (como a media) da primeira comparando nos diferentes niveis da segunda
ggplot(dados) + geom_bar(aes(x=Gr_etario, y = PAS), stat="summary", fun="mean") + labs(x="Grupo etario", y="PAS") #para y=PAS queremos fazer um summary em que a funçao a fazer é a media. Eixo do y vai ter o PAS mas com o valor medio. Cada barra representa valor medio 
Treemap
geom_treemappodem ser vistos como uma modificaçao dos graficos de barra, permitindo mais facilmente efetuar uma comparaçao face ao total
#install.packages("treemapify")
library(treemapify)## Warning: package 'treemapify' was built under R version 4.0.5ggplot(data=tbl, aes(fill=GrupoEtario, area=Freq)) + geom_treemap() + labs(title="Distribuiçao do grupo etario")
tbl## GrupoEtario Freq
## 1 Ate aos 34 anos 16
## 2 35 - 64 anos 68
## 3 65 e mais anos 16ggplot(data=tbl, aes(fill=GrupoEtario, area= Freq, label=GrupoEtario))+geom_treemap()+geom_treemap_text(colour="white", place="centre") + labs(title="Distribuiçao do grupo etario") + theme(legend.position = "none")
Histograma
geom_hist()uteis para distribuiçao de uma amostra relativa a uma variavel quantitativa
ggplot(data=dados,aes(x=Altura)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=dados, aes(x=Altura)) + geom_histogram(fill="dodgerblue", color="black")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#O numero de barras (sub-intervalos) a construir pode ser redefinido de varias formas (numero, largura ou regra)
ggplot(data=dados, aes(x=Altura)) + geom_histogram(fill="dodgerblue", color="black", bins=20)
ggplot(data=dados, aes(x=Altura)) + geom_histogram(fill="dodgerblue", color="black", binwidth=5)
ggplot(data=dados, aes(x=Altura)) + geom_histogram(fill="dodgerblue", color="black", bins=nclass.Sturges(dados$Altura
))
Density plot
- Os histogramas podem ser vistos como uma aproximação descontínua à função densidade da variável aleatória subjacente à variável em estudo
- Para obter uma aproximação contínua, é possível usar métodos de estimação de densidades baseados em funções kernel ﴾kernel density estimation﴿ contínuas
- O resultado permite obter uma curva contínua, mais condizente com o esperado
ggplot(data=dados, aes(x=Altura)) + geom_density()
ggplot(data=dados, aes(x=Altura)) + geom_density(fill="dodgerblue", color="firebrick", size=1)
ggplot(data=dados, aes(x=Altura)) + geom_density(fill="dodgerblue", color="firebrick", size=1, bw=1) #densidade estimada depende do kernel bandwidth, bw. Se bw = 0.1, estimativa esta a ser muito local. Quanto mais alargarmos o bw, mais "smooth" fica
ggplot(data=dados,aes(x=Altura,fill=Gr_etario)) + geom_density(alpha=0.5) #density plot para cada grupo etario
#OU
#install.packages("ggridges")
library(ggridges)## Warning: package 'ggridges' was built under R version 4.0.5ggplot(data=dados, aes(x=Altura, y = Gr_etario, fill=Gr_etario)) + geom_density_ridges()## Picking joint bandwidth of 2.01
Diagramas de dispersao
Pretende-se analisar a existencia de relaçao entre duas variaveis quantitativas
- estimando um modelo linear
grafico=ggplot(dados,aes(x = Altura,y = Peso)) + geom_point(color="cornflowerblue",size = 2,alpha=.8) + scale_y_continuous(limits = c(50,90)) + scale_x_continuous(breaks = seq(155, 180, 5), limits=c(155, 180))
grafico + geom_smooth(method="lm")## `geom_smooth()` using formula 'y ~ x'
- estimando um modelo polinomial
grafico + geom_smooth(method = "lm", formula = y ~ poly(x, 2), color = "indianred3")
- usando metodo de estimaçao local
grafico + geom_smooth(method = "loess", color = "indianred3")## `geom_smooth()` using formula 'y ~ x'
Caixa de bigodes
- Obtidos com o comando
geom_boxplotpodendo-se escolher a orientaçao horizontal ou vertical comcoord_flip
#orientaçao horizontal
ggplot(data=dados, aes(x=Status)) + geom_boxplot()
#Vertical
ggplot(data=dados, aes(x=Status)) + geom_boxplot()+coord_flip()
ggplot(data=dados, aes(y=Status)) + geom_boxplot()
#- É tambem uma boa opçao quando se pretende comparar as distribuiçoes de uma variavel quantitativa separadamente por niveis de uma qualitativa
ggplot(data=dados,aes(x=Escolaridade,y=Status)) + geom_boxplot(color="black",fill="cadetblue1", outlier.colour = "red", outlier.shape = 2)
* A opção notch permite obter uma aproximaçao a um IC a 95% para as medianas e assim efetuar uma comparaçao entre as mesmas (a nao sobreposiçao dos respetivos notches indica diferenças significativas)
ggplot(data=dados,aes(x=Escolaridade,y=Status))+ geom_boxplot(color="black",fill="cadetblue1", outlier.colour = "red", outlier.shape = 2,notch = TRUE) #intervalos de confiança sao sobrepostos, logo nao ha diferença significativa entre medianadas de status entre gurpos etarios## notch went outside hinges. Try setting notch=FALSE.
- Os violin plots sao uma combinaçao entre diagramas de caixa e os graficos de densidade
ggplot(data=dados,aes(x=Escolaridade,y=Status))+ geom_violin(fill="cornflowerblue")+ geom_boxplot(color="black",fill="cadetblue1", outlier.colour = "red", outlier.shape = 2, width=0.3) #conjugar info de boxplot e density plot, o width controla o tamanho da caixa
Error plots
- Permitem analisar o comportamento medio e respetiva variabilidade de uma variavel quantitativa por grupos definidos por variaveis qualitativas
- É necessario preparar os dados na forma de uma matriz (com o pacote
psych)
#install.packages("psych")
library(psych)## Warning: package 'psych' was built under R version 4.0.5##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphaaux=describeBy(dados$PAS,group=dados$Gr_etario,mat = T) #describeBy calcyla sd, media, etc
View(aux) #mat=T passa a formato matriz para ser usado no ggplot
ggplot(aux,aes(x = group1, y = mean,group = 1)) + geom_point(size = 3) + geom_line() + geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width = .1) + labs(x="Grupo Etário",y="PAS médio")
# com duas qualitativas
aux=describeBy(dados$PAS,group=list(dados$Gr_etario,dados$Tratamento),mat = T)
ggplot(aux,aes(x = group1, y = mean,group = group2,color=group2)) + geom_point(size = 3) + geom_line() + geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width = .1) + labs(x="Grupo Etário",y="PAS médio")
#cruzamento de linhas pode significar que variaveis tem interaçao entre si. PAS pode ter resultados diferentes consoante os niveis das variaveis. Ha interaçao entre os niveis dos grupo etario com os tratamentos e o PAS Graficos de correlaçoes
#install.packages("ggcorrplot")
library(ggcorrplot)## Warning: package 'ggcorrplot' was built under R version 4.0.5dados1=dados[,2:ncol(dados)]
MatCorr=cor(dados1[,sapply(dados1,is.numeric)])
ggcorrplot(MatCorr)
- Podemos reorganizar as variaveis, juntando aquelas com padroes semelhantes de correlaçao ou alterar o aspeto do grafico mostrando apenas uma das partes triangulares (Recordar que a matriz é simetrica)
ggcorrplot(MatCorr,hc.order = TRUE,type = "lower",lab = TRUE) #hc.order - clustering das correlaçoes
Dados longitudinais - Dados em que pacientes sao medidos ao longo do tempo para perceber se ha alguma evoluçao face a um tratamento, p.exemplo
library(HSAUR3)
data("epilepsy")
ggplot(epilepsy, aes(x = period, y = seizure.rate/base)) + geom_line(aes(color = treatment, group = subject)) + geom_point() #group = subject - da os pontos do grafico 
ggplot(epilepsy, aes(x = period, y = seizure.rate/base)) + geom_line(aes(group = subject)) + facet_wrap(~treatment,labeller = label_both) #facet wrap separa por tratamento
ggplot(epilepsy, aes(x = period, y = seizure.rate/base)) + geom_line(aes(color = treatment, group = subject)) + facet_wrap(~subject,labeller=label_both) #linha individualizada, por sujeito
ggplot(epilepsy, aes(x = period, y = seizure.rate/base)) + stat_summary(aes(group = treatment, color = treatment), geom = "line", fun = mean, size = 3) #valores medios das seizure rates de todos do progabide vs todos do placebo
ggplot(epilepsy, aes(x = period, y = seizure.rate/base)) + geom_line(aes(color = treatment, group = subject)) + stat_summary(aes(group = treatment, color = paste("mean", treatment)), geom = "line", fun = mean, size = 3)
3 - Distribuiçoes e simulaçao, reamostragem
Distribuições de probabilidade
Tipos | dnorm, pnorm, qnorm, rnorm…
- R possui varias distribuiçoes de probaiblidade que permitem calcular valores da função distribuiçao, da funçao densidade, gerar aleatorios…
- Primeira letra significa o tipo de valor a calcular (probabilidade/densidade, probabilidade acumulada, numero aleatorio) e as restantes especificam a distribuiçao
dnormobtém valores da densidade de uma distribuição normal - f minusculopnormobtém valores de probabilidade acumulada ﴾função distribuição﴿ de uma distribuição normal - F maiusculoqnormobtém os quantis de uma distribuição normal - damos probabilidade acumulada e retorna o X onde atingimos essa probabilidade (inversa da pnorm)rnormgera números aleatórios ﴾random﴿ de uma distribuição normal
Exemplos
- Tipo discreto | Binomial
- Numero de sucessos em n provas independentes (nº de faces saidas em 20 lançamentos de uma moeda)
- Depende de 2 parametros: n e p (probabilidade de ocorrer um sucesso ao acaso)
- Representa-se por X ~ B(n,p). Assim:
- P(X = x) =
dbinom(x,n,p)- probabilidade de ocorrer X sucessos - P(X ≤ x) =
pbinom(x,n,p)- probabilidade acumulada até ao X - P(X ≤ x) = k ⇔ x =q
qbinom(k,n,p)﴾k é um valor de probabilidade dado﴿ rbinom(N,n,p)gera N números aleatórios com distribuição B(n, p)
- P(X = x) =
#B(20, 0.5) ﴾gráfico conjunto da amostra e da distribuição exata﴿
p.exato=dbinom(0:20,20,0.5)
x=rbinom(1000,20,0.5)
B=barplot(prop.table(table(factor(x,levels=0:20))))
lines(B,p.exato,type="h",col="blue",lwd=2)
- Tipo continuo | Normal
- Definia por 2 parametros, media e desvio padrao, sigma (ou variancia). A area abaixo da curva tem que ser sempre 1, por isso, se as medidas de forma (desvio padrao) forem pequenos (mais compactos), o tamanho do sino vai subir.
- Representa‐se por X ∼ N(μ, σ2). Assim,
- f(x) =
dnorm(x,mu,sigma)﴾numa distribuição contínua f(x) diz‐se a função densidade e não representa uma probabilidade﴿. dnorm como nao estamos nas discretas, ja nao representa uma probabilidade, mas sim cada um dos pontos. - P(X ≤ x) =
pnorm(x,mu,sigma) - P(X ≤ x) = k ⇔ x =
qnorm(k,mu,sigma)﴾k é um valor de probabilidade dado﴿ retorna o ponto até a area acumulada rbinom(N,mu,sigma)gera N números aleatórios com distribuição N(mu, sigma2)
- f(x) =
x=rnorm(1000,mean=0,sd=1)
y=rnorm(1000,mean=5,sd=2)
hist(x,freq=F) #freq = F, nao vai calcular frequencias. Temos de colocar freq = F pois se tirassemos o Freq, iamos criar um histograma com escala da nossa amostra (0 a 200) e a curva iria ter escala de 0 a 1 (vai estar escondida). Freq = F faz frequencias relativas e depois normaliza os valores para que a soma das barras dê 1
lines(sort(x),dnorm(sort(x),mean=0,sd=1))
hist(y,freq=F)
lines(sort(y),dnorm(sort(y),mean=5,sd=2))
* Os QQ-plots tambem sao uteis para avaliar o ajustamento de uma distribuiçao. Funções qqnorm e qqline podem ser usadas para analisar a proveniencia de uma dist. normal
x=rnorm(1000,0,1)
qqnorm(x)
qqline(x,col="blue")
y=rchisq(1000,3)
qqnorm(y)
qqline(y,col="red")
Introdução à simulação e reamostragem
Estudos Simulados
- Realizar numero elevado de experiencias controladas
- Permite estabelecer propriedades de estimadores dos quais é dificil obter por via analitica
- Perceber diferenças de modelos com diferentes amostras
- Ver as consequencias da falha de pressupostos, testes de hipoteses, ICs…
- Comparar performance de modelos
Simulaçao vs Reamostragem
Simulaçao | Simulação de monte carlo
- Nas simulaçoes de monte carlo, temos o completo controlo sobre todo o processo que leva à geraçao da amostra (distribuiçao…); podemos assim gerar varias amostras que simulam um determinad oprocesso a partir de um modelo teorico
- Envolve geraçao de varias amostras aleatorias por um processo gerador pre-definido (estabelece o modelo dos dados)
Exemplo 1
Consideremos a estatística T, usada na construção de IC e em TH 
- Na hipótese de X ∼ N(μ, σ2) então T ∼ tn−1 (por nao conhecermos o sigma, é por isso que temos graus de liberdade n-1)
- Suponhamos que era analíticamente difícil estabelecer a distribuição de T; podemos simular instâncias de T e testar sobre algumas distribuições
N=1000; # número de repetições
media=0
desviop=1
n=10
T=NULL
for (i in 1:N){
x=rnorm(n,media,desviop)
T[i]=(mean(x)-media)/(sd(x)/sqrt(n))
}
#Perceber qual a dist. que tem o T:
qqnorm(T);qqline(T)
shapiro.test(T) #vai ver se t obs esta ou nao na regiao critica. Rejeita a normalidade pois se é uma t student rejeita a normalidade. Temos que usar a distribuiçao correta para confirmar a distribuiçao de T, que neste caso será uma KS .##
## Shapiro-Wilk normality test
##
## data: T
## W = 0.97988, p-value = 1.571e-10#NOTA: se aumentarmos muito o n, o shapiro.teste ja vai confirmar a normalidade apesar de ser a distribuiçao errada uma vez que aumentando o n, o T vai tender para a normalidade
#Como rejeitamos com o shapiro.test, vamos agora ajustar cm os graus de liberdade n-1, escolhendo a disribuiçao "pt" (t-student)
ks.test(T,"pt",df=(n-1))##
## One-sample Kolmogorov-Smirnov test
##
## data: T
## D = 0.033222, p-value = 0.2197
## alternative hypothesis: two-sidedExemplo 2
- Comparar três estimadores para a média μ de uma distribuição baseada em amostras i.i.d.:
- Média amostral ¯X
- Média amostral aparada a 20% ¯Xap
- Mediana mostral Med (poderiamos usar a mediana em vez da media se a distribuiçao fosse simetria)
- Se a distribuição dos dados for simétrica, os três estimadores são plausíveis para estimar a média, o que não acontece se for assimétrica
set.seed(3)
N=1000; n=15; media.ex=0
media=NULL;mediaAp=NULL; mediana=NULL
for (i in 1:N){
x=rnorm(n,media.ex,1)
media[i]=mean(x)
mediaAp[i]=mean(x,trim=0.2)
mediana[i]=median(x)
}
resultados=data.frame(media=media, media.aparada=mediaAp, mediana=mediana)
# viés
colMeans(resultados)-media.ex # media dos valores esperados - media exata (Teorica) = vies## media media.aparada mediana
## -0.011501470 -0.010150869 -0.006407224# Desvio padrão
sqrt(diag(var(resultados))) #variabilidade de um estimador - o quao diferente o estimador pode ser consoante diferentes amostras. Um sd pequeno significa que por muito que a amostra varie, o estimador vai ser consistente## media media.aparada mediana
## 0.2563014 0.2695634 0.3080002# Erro quadrático médio - quanto menor o erro quadratico melhor. É a variancia da diferença
diag(var(resultados-media.ex))## media media.aparada mediana
## 0.06569039 0.07266444 0.09486415# Eficiência relativa face ao estimador standard da média (media amostral)
var(resultados$media)/diag(var(resultados)) #mediana seria um estimador menos eficiente face à media tradicional## media media.aparada mediana
## 1.0000000 0.9040239 0.6924680#Perante todos os resultados, a media da amostra é o melhor estimador para estimar a media populacionalReamostragem
- Sao geradas multiplas amostras mas agora nao a partir de um modelo teorico pre-estabelecido, mas sim a partir de uma amostra de dados disponivel. Originamos estas amostras multiplas a partir da amostra orginal por metodos como:
- Testes de permutaçao ou testes de aleatorizaçao (randomizaçao)
- Jackknife
- bootstrap
Testes de permutaçao
- Gerar varias amostras baralhando aleatoriamente a amostra original, quebrando qualquer relaçao que exista na amostra original
- Compara-se amostra original com as amostras baralhadas tentando-se verificar se uma potencial relaçao existente na amostra original nao tera ocorrido ao acaso
- Tipicamente usado para testar a H0 de que um dado tratamento nao tem efeito (quando temos n = 10 para t student, por exemplo)
Exemplo
Exemplo: comparar as médias de uma variável x entre dois grupos de indivíduos dos quais 3 receberam um tratamento e 3 não receberam esse tratamento
library(combinat)##
## Attaching package: 'combinat'## The following object is masked from 'package:utils':
##
## combnset.seed(123)
individuos=letters[1:6]
tratamento=factor(c("Sim","Sim","Sim","Não","Não","Não"),levels = c("Sim","Não"))
x=rnorm(6,mean=c(5,5,5,0,0,0),sd=1) #se alterarmos a meida para 1, vai ser quase indistinguivel pois vao estar sobrepostas. Nesse caso, a diferença ja é irrelevante
dados=data.frame(individuos,tratamento,x)
dados## individuos tratamento x
## 1 a Sim 4.43952435
## 2 b Sim 4.76982251
## 3 c Sim 6.55870831
## 4 d Não 0.07050839
## 5 e Não 0.12928774
## 6 f Não 1.71506499aux=by(dados$x,dados$tratamento,mean); aux## dados$tratamento: Sim
## [1] 5.256018
## ------------------------------------------------------------
## dados$tratamento: Não
## [1] 0.638287dif.obs=aux[1]-aux[2]Por inspeção visual ﴾e por construção/simulação do próprio problema﴿ sabemos que existe uma diferença nítida entre as médias dos dois grupos
Uma primeira abordagem seria utilizar um teste de hipóteses para a diferença de médias ﴾teríamos de assumir a Normalidade dos dados﴿ ou um teste não paramétrico
Como nesta abordagem são geradas todas as possíveis permutações da amostra original, é possível calcular a probabilidade exata da diferença observada ﴾na amostra original﴿ ocorrer por acaso
Neste baralhamento, estamos a quebrar as ligaçoes com a amostra original. Faz todas as trocas possiveis e calcula em cada uma delas a media por grupo e faz a diferença.
# Gerar todas as permutações possíveis de tratamento
combin=unique(permn(dados$tratamento)) #todas as permutaçoes para a variavel tratamento - 720 permutaçoes, mas trocar um sim com outro sim da a mesma permutaçao, é por isso que colocamos o unique para ter apenas as diferentes
baralha.tratamento = matrix(unlist(combin), nrow = length(combin), ncol = length(dados$tratamento), byrow = TRUE)
baralha.tratamento## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] "Sim" "Sim" "Sim" "Não" "Não" "Não"
## [2,] "Sim" "Sim" "Não" "Sim" "Não" "Não"
## [3,] "Sim" "Não" "Sim" "Sim" "Não" "Não"
## [4,] "Não" "Sim" "Sim" "Sim" "Não" "Não"
## [5,] "Sim" "Sim" "Não" "Não" "Sim" "Não"
## [6,] "Sim" "Não" "Sim" "Não" "Sim" "Não"
## [7,] "Não" "Sim" "Sim" "Não" "Sim" "Não"
## [8,] "Não" "Sim" "Não" "Sim" "Sim" "Não"
## [9,] "Sim" "Não" "Não" "Sim" "Sim" "Não"
## [10,] "Não" "Não" "Sim" "Sim" "Sim" "Não"
## [11,] "Não" "Não" "Sim" "Sim" "Não" "Sim"
## [12,] "Não" "Sim" "Não" "Sim" "Não" "Sim"
## [13,] "Não" "Sim" "Sim" "Não" "Não" "Sim"
## [14,] "Sim" "Não" "Sim" "Não" "Não" "Sim"
## [15,] "Sim" "Não" "Não" "Sim" "Não" "Sim"
## [16,] "Sim" "Sim" "Não" "Não" "Não" "Sim"
## [17,] "Sim" "Não" "Não" "Não" "Sim" "Sim"
## [18,] "Não" "Sim" "Não" "Não" "Sim" "Sim"
## [19,] "Não" "Não" "Sim" "Não" "Sim" "Sim"
## [20,] "Não" "Não" "Não" "Sim" "Sim" "Sim"# Calcular a diferença de médias para cada combinação
baralha.dif=NULL
for (i in 1:nrow(baralha.tratamento)){ #percorre a matriz
aux=by(dados$x,baralha.tratamento[i,],mean) #pegamos no x e fazemos as medias por grupo (sim e nao) - calculamos a distribuiçao das diferenças ao acaso. Depois vemos onde é que a nossa observaçao inicial (a primeira com sim sim sim nao nao nao), onde é que se localiza na distribuiçao - se ela se posicionar no meio, entao nao é mt diferente do acaso, se ficar nas extermidades, entao esta mt disante do acaso. A diferença que encontramos, estando nas extermidades, é estatisticamente significativa
baralha.dif[i]=aux[1]-aux[2]
}
hist(baralha.dif)
lines(c(dif.obs,dif.obs),c(0,4),col="red")
quantile(baralha.dif,probs = seq(0,1,0.05))## 0% 5% 10% 15% 20% 25% 30%
## -4.6177314 -2.8922404 -2.6032462 -1.8697882 -1.7128913 -1.5692946 -1.4966111
## 35% 40% 45% 50% 55% 60% 65%
## -1.4223126 -0.7543249 -0.3098985 0.0000000 0.3098985 0.7543249 1.4223126
## 70% 75% 80% 85% 90% 95% 100%
## 1.4966111 1.5692946 1.7128913 1.8697882 2.6032462 2.8922404 4.6177314- O histograma apresenta‐nos a distribuição das diferenças obtidas por acaso, isto é, na hipótese ﴾nula﴿ de não existência de diferença ﴾veja‐se que o valor médio é próximo de zero﴿ entre tratamento e controlo
- Verificamos assim que a probabilidade da diferença observada ocorrer por acaso é muito baixa:
- o valor dif.obs=4.6 está na cauda do histograma
- o valor dif.obs=4.6 está acima do quantil 95%
- Conclui‐se assim pela rejeição da hipótese nula, pelo que existe diferença significativa entre controlo e tratamento
Testes de aleatorizaçao
- Metodo de permutaçoes nao é bom para amostras muito grande devido ao crescimento exponencial de combinaçoes possiveis. Nos testes de aleatorizaçao, apenas se analisam um numero elevado de possiveis permutaçoes, mas nao toda.
tratamento=factor(c(rep("Sim",30),rep("Não",40)), levels = c("Sim","Não"))
x=rnorm(70,mean=c(rep(5,30),rep(0,40)),sd=1)
dados=data.frame(tratamento,x)
aux=by(dados$x,dados$tratamento,mean)
dif.obs=aux[1]-aux[2];dif.obs## Sim
## 4.952294- Esta amostra dava algo muito grande de permutaçoes, por isso, vamos usar apenas
repspermutaçoes
reps=1000
combin=unique(t(replicate(reps,sample(dados$tratamento,70)))) #vamos aos dados tratamento fazer uma amostragem aleatoria - Em vez de termos sim todos seguidos e nao todos seguidos, vamos ter reamostragem -replicate repete 1000 vezes a amostragem
baralha.dif=NULL
for (i in 1:nrow(combin)){
aux=by(dados$x,combin[i,],mean)
baralha.dif[i]=aux[1]-aux[2]
}
hist(baralha.dif)
lines(c(dif.obs,dif.obs),c(0,100),col="red") #nao aparece porque o valor obs está nos 4 e a escala so vai até 2 - é uma diferença significativa
quantile(baralha.dif,probs = seq(0,1,0.05))## 0% 5% 10% 15% 20% 25%
## -1.98673863 -0.98688545 -0.76627098 -0.62990576 -0.52213776 -0.41765093
## 30% 35% 40% 45% 50% 55%
## -0.32449068 -0.24499685 -0.17668453 -0.10327525 -0.01961262 0.07284918
## 60% 65% 70% 75% 80% 85%
## 0.14873653 0.22780426 0.30354989 0.40457066 0.47581681 0.61243908
## 90% 95% 100%
## 0.77724018 1.00235838 1.97366720dif.obs## Sim
## 4.952294!!! Duvida porque usamos quantis como teste nao parametrico | Metodo de bootstrap
- Tecnica de amostragem com reposiçao. Gera amostras a partir da amostra disponivel.
- 3 caracteristicas:
- Amostragem deve ser independente ou seja qualquer observaçao da amostra original deve ter a mesma prob. de ser escolhida
- Amostra bootstrap devera ter o mesmo tamanho da amostra original
- Amostragem é feita com reposiçao
- Suponhamos que pretendíamos estimar uma propriedade de um estimador ˆθ de um parâmetro θ de uma qualquer população F ﴾desconhecida﴿ tendo apenas disponível uma amostra S de tamanho n
- Com o método bootstrap
- Obter uma amostra aleatória, Sboot, de tamanho n por amostragem com reposição de S
- Calcular ˆθ a partir de Sboot e guardar o valor num vetor V
- Repetir N vezes os passos 1 e 2.
- Se N for suficientemente grande, o vetor V ﴾que contém uma amostra de valores de θˆ﴿ pode ser agora usado para estimar o erro padrão, um intervalo de confiança, etc., para θ
- É fácil perceber que cada amostra Sboot pode conter repetições de uma ou mais observações da amostra original ou mesmo não conter algumas das observações originais
- É esta característica que permite obter variabilidade de amostras e reproduzir a população F desconhecida, permitindo obter diferentes estimativas para θ
paises = c("EUA", "Canadá", "México",
"Inglaterra", "França","Espanha", "Alemanha", "Itália",
"China", "Japão","Portugal")
sample(paises, replace = FALSE) #sample faz reamostragem com mesmo tamanho da amostra inicial quando nao colocamos reps, e com replace = F, vai buscar todos os paises por ordem diferente de forma aleatoria ## [1] "França" "México" "China" "Inglaterra" "Canadá"
## [6] "Japão" "Alemanha" "Portugal" "Itália" "EUA"
## [11] "Espanha"sample(paises, replace = TRUE) #com replace = T, pode ter repetiçoes do mesmo pais## [1] "Alemanha" "EUA" "França" "Canadá" "Espanha"
## [6] "França" "Inglaterra" "Canadá" "Inglaterra" "Itália"
## [11] "Canadá"sample(paises, replace = TRUE)## [1] "Japão" "Espanha" "Portugal" "Japão" "Japão"
## [6] "Espanha" "Inglaterra" "EUA" "China" "Canadá"
## [11] "Portugal"- Neste exemplo pretendemos comparar o valor do erro padrao da media amostral (dado por (desvio padrao / raiz de n) com estimativa obtida por simulaçao
set.seed(123)
n = 500
N = 1000
amostra = rnorm(n, 4, 5)
media = mean(amostra)
erro.pad = 5/sqrt(n) #sigma / raiz de n
erro.pad.est = sd(amostra)/sqrt(n) #desvio padrao da amostra / raiz de n - podemos ver que o sd da amostra é um estimador do sd populacional pois sao proximos
erro.pad;erro.pad.est## [1] 0.2236068## [1] 0.2175178- Por definiçao, o erro padrao da media é o desvio padrao da sua distribuiçao amostral
- Podemos usar bootstrap para construir a distribuiçao amostral de medias
media.boot = NULL
for (i in 1:N) {
amostra.boot = sample(amostra,size=n,replace = TRUE)
media.boot[i] = mean(amostra.boot)
}
hist(media.boot,main="Distribuição amostral de médias")
sd(media.boot) #erro padrao é o desvio padrao da dist. amostral das medias. ## [1] 0.2151349# erro padrao da media amostral = sigma / raiz de n (sigma = desvo padrao populaciona) que é aproximado pelo s / raiz de n em que o s é o desvio padrao da amostra original.
# O erro padrao é o desvio padrao da dist. das medias amostrais
#valor obtido é proximo do real (erro.pad e erro.pad.est)- Suponhamos que agora queremos determinar um IC para media com (1-alfa)*100% de confiança:
- Podemos fazelo de diferentes maneiras:
- ]¯x − z1−α/2EPboot, ¯x + z1−α/2EPboot[, onde EPboot é uma estimativa para o erro padrao (podemos obter uma estimativa deste EP)
- Usando uma alternativa nao parametrica (note-se que a anterior pressupoe normalidade) a partir dos quantis da amostra de medias bootstrap
- Permite replicar muitas amostras (até quando temos, por exemplo, apenas uma media, um ponto), criar medias e com isto criar uma distribuiçao de medias com a nossa original. É neste dist. de medias que podemos fazer o corte e perceber onde é que a nossa media teorica cai no grafico
- Podemos fazelo de diferentes maneiras:
#IC 9% parametrico
media - qnorm(0.975)*sd(media.boot) # Limite inf.## [1] 3.751295media + qnorm(0.975)*sd(media.boot) # Limite sup.## [1] 4.594609# IC 95% (não-paramétrico) - como temos dist. de medias, vamos ver onde ocorre 2.5 e 97.5% - os valores sao quase iguais com o IC parametrico. O parametrico é mais proximo do teorico pois assume normalidade dos dados. Nesta nao parametrica, baseia-se na replicaçao bootstrap para simular a distribuiçao dos nossos dados originais
quantile(media.boot, .025) # Limite inf.## 2.5%
## 3.757187quantile(media.boot, .975) # Limite sup.## 97.5%
## 4.606669# Como os dados gerados sao normais, o parametrico e nao parametrico vao dar parecidos, podendo tambem fazer um t.test(amostra) para confirmar os intervalos.- Valores sao proximos pois a dist. verdadeira é normal
Reamostragem para seleçao de modelo | Nao sai no teste?
- A existência de uma grande quantidade de modelos estatísticos diferentes implica a realização de um processo de seleção
- Na verdade, os no free lunch thoerems dizem‐nos que não existe um método que seja o melhor de todos em todos os problemas
- Além disso, a generalidade dos modelos têm parâmetros que é necessário especificar/estimar ﴾por exemplo parâmetros de controlo de flexibilidade do modelo﴿
- Designamos de conjunto de treino ao conjunto ﴾amostra﴿ de dados utilizado para estimar um determinado modelo e de conjunto de teste ao conjunto de dados ﴾da mesma população do conjunto de treino﴿ utilizado para avaliar a performance do modelo criado ﴾dados não utilizados na construção do modelo﴿
O problema de over-fit
- O erro de treino ou resubstituição é uma estimativa de performance obtida com o próprio conjunto de treino sendo por isso demasiado otimista ﴾sujeita a over‐fit﴿
- De facto e em geral, aumentando a flexibilidade do modelo leva a um sobre‐ajuste aos dados de treino resultando num erro de treino cada vez menor; no entanto, a capacidade de generalização ﴾medida pelo erro de teste﴿ fica diminuída (modelo fica muito preso às caracteristicas dos dados de treino nao conseguindo “inferir” em dados reais)
O compromisso vies-variancia
A forma em U para o erro de teste que se pode observar na figura está ligada ao compromisso vies-variancia (bias-variance trade-off)
Por exemplo, prova-se para o caso de modelos de regressao que, para uma dada observaçao de teste x0, se tem:
onde t é a variavel a prever e y(x) é o modelo de regressaoVerifica-se que a reduçao do erro quadratico medio de teste deve ser realizada por reduçao simultanea da variancia e do vies do modelo
- Problema -> Em geral, para reduzir a variancia, devemos diminuir a flexibilidade o que faz aumentar o vies. Reduzir o vies implica aumentar a flexibilidade e, portanto, aumenta a variancia.
- É importante encontrar um compromisso entre ambos (ponto minimo da curva de teste)
- Problema -> Em geral, para reduzir a variancia, devemos diminuir a flexibilidade o que faz aumentar o vies. Reduzir o vies implica aumentar a flexibilidade e, portanto, aumenta a variancia.
Seleçao de modelos
- Qualquer modelo deve ser avaliado num conjunto de dados independente daquele que é usado para ajustar o modelo -> conjunto de teste
- Normalmente, é possivel reduzir o conjnto de treino da nossa amostra inicial (30% dos dados) e deixar o restante para testar o modelo
- Nestas situaçoes devem usar-se preferencialmente estrategias de validaçao para obter a estimativa do erro de teste
Conjunto de validaçao
- Dado um conjunto de dados deveriamos dividi-lo da seguinte forma:
- O conjunto de validaçao, que funciona como uma replica de um conjunto de teste, nao é usado para construir o modelo; é usado para escolher o melhor modelo que se ajusta ao conjunto de treino
Cross-validation | k-folds
- Quando o conjunto de dados original é pequeno, a divisao anterior pode nao ser a melhor opçao (poucos dados para treino)
- Uma forma de aliviar este problema é realizar uma amostragem do tipo
k-fold cross-validation
- Assim, para cada modelo, aplica-se o procedimento anterior para calcular o erro medio das folds a rosa
- Aleatoriza-se o conjunto de dados e repete-se o processo para outro modelo
- Escolhe-se o modelo com menor erro medio de cross-validation
- Avalia-se o modelo escolhido no conjnto de teste (caso esteja disponivel)
Leave one out CV
- É um caso particular do anterior em que k = N (N é o tamanho do conjunto de dados)
- Cada observaçao de treino é deixada de fora e o modelo construido com as outras N-1 observaçoes
- A observaçao deixada de parte é usada como conjunto de validaçao
- Mais exigente a nivel computacional
Como escolher o numero de folds
- Depende de fatores como tamanho do conjunto de dados, distribuiçao dos dados…
- O bias-variance tradeoff ajuda um pouco na escolha: valores maiores de k implicam maior variancia nas estimativas de erro; valores menores de k implicam maior vies da estimativa do erro
- K = 5 ou k=10 permitem obter estimativas que controlam vies e variancia
Exemplo
- Dataset iris, 150 obs e 3 especies da flor iris
- Objetivo é estimar a performance de um modelo que a partir das mediçoes das petalas e sepalas, seja capaz de reconhecer a especie
- Usaremos aqui o metodo dos 10-vizinhos mais proximos
- Uma vez que a amostra é pequena para extrair m conjunto de teste independente, usamos 10-CV e LOOCV
data("iris")
library(caret)## Warning: package 'caret' was built under R version 4.0.5library(class)## Warning: package 'class' was built under R version 4.0.5ind=createFolds(iris[,5],k=10,list=TRUE)
# 10-fold CV
erro.cv10=rep(0,10)
for (i in 1:10){
pred=knn(iris[-ind[[i]],-5],iris[ind[[i]],-5],iris[-ind[[i]],5],k=10)
erro.cv10[i]=mean(iris[ind[[i]],5]!=pred)
}
mean(erro.cv10)## [1] 0.03333333erro.loocv=rep(0,nrow(iris))
# LOOCV
for (i in 1:nrow(iris)){
pred=knn(iris[-i,-5],iris[i,-5],iris[-i,5],k=10)
erro.loocv[i]=(iris[i,5]!=pred)
}
mean(erro.loocv)## [1] 0.03333333- No exemplo anterior estimatmos a perfromance do modelo usando os metodos de 10-fold cross validation (10-fold CV) e leave one out CV (LOOCV)
- Suponhamos agora a necesidade de escolher entre diferentes modelos: sera que 10 vizinhos mais proximos é o melhor modelo?
- Nesta situaçao construimos um modelo para cada valor de k que se pretende testar e avalia-se com CV; no final, escolhe-se o modelo com melhor performance CV
set.seed(123)
ind.CV=createFolds(iris$Species,k=10,list=TRUE)
# 10-fold CV ; testar k=1,...,30
erro.k=rep(0,30)
for (k in 1:30){
erro.cv10=rep(0,10)
for (i in 1:10){
pred=knn(iris[-ind.CV[[i]],-5],iris[ind.CV[[i]],-5],
iris[-ind.CV[[i]],5],k=k)
erro.cv10[i]=mean(iris[ind.CV[[i]],5]!=pred)
}
erro.k[k]=mean(erro.cv10)
}
plot(erro.k,type="l",xlab="nº de vizinhos",ylab = "Erro")