Mini - Curso

Estatística em Hidrologia

Software aplicado

a recursos hídricos

Setembro de 2011

Fernanda Valente

2

O que é o R?

É um conjunto integrado de ferramentas computacionais que permitem a manipulação e análise de dados, o cálculo numérico e a produção de gráficos.

É uma linguagem de programação simples (uma variante da linguagem S).

É uma linguagem interpretada: – os comandos são imediatamente executados.

É uma linguagem orientada por objectos: – os dados manipulados são armazenados na memória activa

do computador na forma de objectos, que têm um nome e aos quais se podem aplicar acções.

3

Instalar o R

O R é uma aplicação de distribuição gratuita e de código público (http://cran.r-project.org/), existindo versões já compiladas para execução nos principais sistemas operativos (Windows, Linux e Macintosh).

Depois de fazer o download da versão adequada ao sistema operativo do computador (por ex. R-2.12.2-win32.exe, para o Windows), para instalar o R basta executar esse ficheiro.

4

Iniciar uma sessão de R

Criar na área de trabalho uma nova pasta (por exemplo, CursoR) onde irão ser guardados os ficheiros de dados – pasta de trabalho.

Iniciar o R:Start → All Programs → R → R 2.12.1

(em alternativa, pode-se criar um atalho no desktop e alterar as suas propriedades de modo a iniciar a aplicação na pasta de trabalho)

5

RGui (Graphical user interface)

O R funciona fundamentalmente no modo “pergunta-resposta”: – Os comandos (frases que se escrevem numa certa sintaxe)

introduzem-se a seguir à prompt ( > ) e são executados após pressionar Enter ( )

Edição de comandos:– Seleccionar comandos executados anteriormente: ↑ ou ↓– Percorrer a linha de comandos: ← ou → – Colocar o cursor no início / fim da linha de comandos:

Home / End

Para terminar o R:– Executar o comando q() ou fechar a janela da aplicação.

6

Comandos elementares

Expressão - o resultado é apresentado no visor e não é registado em memória.

> 2+3/4*7^2[1] 38.75> exp(-2)/log(sqrt(2))[1] 0.3904951 Atribuição – o resultado da expressão é atribuído à variável e

não é apresentado no visor.> x <- 2 # <- é o operador de atribuição> x[1] 2> x <- 2+3/4*7^2> x[1] 38.75

7

Designações das variáveis

Os nomes das variáveis podem ser uma combinação de letras (maiúsculas são diferentes de minúsculas), algarismos e pontos. O 1º caracter deve ser uma letra.

Alguns nomes são utilizados pelo sistema, pelo que devem ser evitados (por ex., c, q, t, C, D, F, I, T, diff, df, pt).

Exemplo: aula1.ex

8

Ajuda Para obter informação sobre uma função específica, por exemplo a

função mean(), pode utilizar-se qualquer uma das opções:

> help(mean)> help("mean")> ?mean> ?"+" # Nalguns casos é obrigatório utilizar aspas Para pesquisar uma sequência de caracteres no help do R o comando é

help.search(), por exemplo

> help.search("standard desviation")> ??"standard desviation" Para navegar nas páginas de ajuda do R utilizando um browser de

internet o comando é

> help.start()Nota: qualquer um dos comandos anteriores pode também ser executado

através do menu Help.

9

Objectos

As entidades que o R cria e manipula e que ficam guardadas na memória do computador. – Para listar os objectos disponíveis

> ls()– Para mostrar alguma informação sobre os objectos:

> ls.str()– Para remover objectos:

> rm(objecto1, objecto2,...)– Para remover todos os objectos:

> rm(list=ls())

10

Workspace

O conjunto de todos os objectos disponíveis numa sessão.

Pode ser guardado num ficheiro com vista à sua utilização em futuras sessões de R.

menu File → Save Workspace...

– No final de cada sessão é perguntado ao utilizador se quer guardar o workspace. Caso a resposta seja afirmativa, todos os objectos disponíveis em memória são guardados no ficheiro “.RData”, na pasta de trabalho em uso, podendo ser carregados na próxima sessão de R.

Podem ser guardados só alguns dos objectos> save(x,file="Exemplo.RData”)

11

Pasta de trabalho (Working Directory)

Localização da pasta de trabalho> getwd()

Especificação da pasta de trabalho> setwd("C:/Users/fernanda/Documents/CursoR")ou

menu File → Change dir...

Importar um ficheiro de objectos> load("Exemplo.RData”)ou

menu File → Load Workspace...

12

Exercício 1

a) Importar os ficheiros de objectos DadosMeteo.RData e Precipitacao.RData

b) Ver todos os objectos em memória (com e sem informação adicional)

Nota: 1 - Os objectos precip, temp e vento contêm dados diários para o ano de 2006 de precipitação (mm), temperatura média do ar (ºC) e velocidade média do vento (ms-1), respectivamente, registados em Évora (Portugal).2 – O objecto pg.Jan contém valores totais de precipitação (mm) registados em Lisboa e no Porto no mês de Janeiro, no período 1980 - 1999.3 – O objecto pg.PGarcia contém valores totais anuais de precipitação (mm) registados em Penha Garcia, Castelo Branco (Portugal), entre 1910 e 1994.

13

Scripts em R

Um script em R é um ficheiro de texto com a extensão “.R”, que contém uma sequência de comandos e que pode ser executado com o comando source("ficheiro") ou através do menu File → Source R code...

Tipicamente, um script é útil quando se pretende repetir as mesmas tarefas várias vezes ou quando se quer divulgar a outros os comandos necessários para resolver um dado problema.

Um modo simples de criar um script em R é abrir uma janela do editor de R (menu File → New script), ir escrevendo os diversos comandos e, no final, guardar o ficheiro de texto (menu File → Save). Caso se pretenda ver ou alterar um programa já existente

menu File → Open script...

14

Alguns objectos

– vector– factor– matrix– list– data Frame– ts– function

15

Vector

É uma colecção ordenada de elementos do mesmo tipo (valores numéricos, lógicos, alfanuméricos, ...).

Uma das maneiras de criar um vector no R é através da função c().

> x <- c(5.4, -3.7, 11.2, 0.78, 21.6) > x[1] 5.40 -3.70 11.20 0.78 21.60> y <- c(FALSE, TRUE, TRUE, TRUE, FALSE) > y[1] FALSE TRUE TRUE TRUE FALSE> nomes <- c("Ana","Paulo","Zé")> nomes[1] "Ana" "Paulo" "Zé"

16

Vector

Um vector pode conter os símbolos especiais NA e NaN, que designam um valor desconhecido e um valor não definido, respectivamente.

> y <- c(NA, 53, 31, 15, 62)> sqrt(c(-1,1,2))[1] NaN 1.000000 1.414214Warning message: NaNs produced in: sqrt(c(-1, 1, 2))

17

Vector

c() é a função concatenação. Os argumentos desta função podem ser eles próprios vectores, pelo que c(v1,v2) devolve a concatenação dos vectores v1 e v2. Esta função pode ter um qualquer número de argumentos.

> c(x,0,x)[1] 5.40 -3.70 11.20 0.78 21.60 0.00 5.40 -3.70 11.20 0.78 21.60

> nomes <- c(nomes,"Maria") > nomes[1] "Ana" "Paulo" "Zé" "Maria"

18

Geração de sequências regulares

O R permite gerar sequências de valores numéricos > 1:15 # define o vector c(1,2,...,14,15)> 15:1 # define o vector c(15,14,...,2,1)> 2*1:15 # define o vector c(2,4,...,28,30)> (2*1):15 # define o vector c(2,3,...,14,15)> seq(1,15)> seq(to=15,from=1)> seq(-5,-1,by=0.5)[1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 > seq(length=9, from=-5, by=0.5)> rep(1:3,2) # define o vector c(1,2,3,1,2,3)> rep(1:3,each=2) # é o vector c(1,1,2,2,3,3)> rep(1:3,1:3) # é o vector c(1,2,2,3,3,3)

19

Vector numérico (operações aritméticas)

As operações entre vectores são realizadas componente a componente.

Alguns operadores aritméticos: +, -, *, /, ^> v1 <- c(5,7)> v2 <- 10:13> 1/v1[1] 0.2000000 0.1428571> v1+v2 # o vector de menor dimensão é repetido> v1*v2 # (total ou parcialmente) até ficar com o

mesmo número de elementos do maior

> v1*4 # cada componente de v1 é multiplicada por 4

20

Vector numérico (algumas funções)

log,exp,sin,cos,tan,atan,sqrt,abs, ... length(x)devolve o número de elementos do vector x sum(x)devolve a soma dos elementos do vector x prod(x)devolve o produto dos elementos do vector x cumsum(x)devolve um vector cujos elementos são a soma

acumulada dos elementos do vector x cumprod(x)devolve um vector cujos elementos são o

produto acumulado dos elementos do vector x sort(x)devolve um vector com os elementos do vector x

ordenados por ordem crescente (decrescente com o argumento decreasing=TRUE)

choose(n,k)devolve o número de combinações nCk

21

Vector lógicoOs resultados das operações que envolvem os

operadores relacionais: <, <=, >, >=, ==, !=

são valores lógicos.> 2==sqrt(4)[1] TRUE> 2!=sqrt(4)[1] FALSEQuando aplicados a vectores, produzem vectores

lógicos.> x[1] 5.40 -3.70 11.20 0.78 21.60> x>11[1] FALSE FALSE TRUE FALSE TRUE

22

Vector lógico (operações)Os operadores lógicos conjunção, disjunção e

negação são, respectivamente: & | !

> x[1] 5.40 -3.70 11.20 0.78 21.60> x>13[1] FALSE FALSE FALSE FALSE TRUE> x<5[1] FALSE TRUE FALSE TRUE FALSE> x>13 | x<5[1] FALSE TRUE FALSE TRUE TRUE

23

Vector lógico (função is.na)

A função is.na(x) devolve um vector lógico em que cada elemento é TRUE se o elemento homólogo em x é do tipo NA ou NaN, e ou FALSE caso contrário.

> z[1] 5.4 -3.7 11.2 NA 21.6> is.na(sqrt(z))[1] FALSE TRUE FALSE TRUE FALSEWarning message:NaNs produced in: sqrt(z)> !is.na(z)[1] TRUE TRUE TRUE FALSE TRUE

24

Selecção de componentes de vectores

Cada elemento de um vector pode ser seleccionado indicando a sua posição entre parêntesis rectos a seguir ao nome do vector:

> letters> letters[3][1] "c" Pode também criar-se um novo vector constituído por alguns

dos elementos do vector indicando as suas posições num vector de índices colocado entre parêntesis rectos:

> letters[c(1,2,4)][1] "a" "b" "d"> letters[-(1:20)] # devolve os elementos de letters

que não estão nas posições indicadas

[1] "u" "v" "w" "x" "y" "z"

25

Selecção de componentes de vectores (cont.)

Pode também seleccionar-se elementos de um vector escrevendo uma condição entre parêntesis rectos a seguir ao nome do vector.

> x[1] 5.40 -3.70 11.20 0.78 21.60> x[x>0][1] 5.40 11.20 0.78 21.60> x[x>0 & x<10][1] 5.40 0.78

> doy <- 1:365> doy[precip>5] # equivalente a which(precip>5)[1] 80 85 87 92 131 136 ...> temp[precip>0 & temp<5][1] 38 60

26

Exercício 2

O ficheiro “serrasPortugal.RData” contém os objectos serras.altitude e serras.nome com dados da altitude (em m) e do nome (pela mesma ordem), respectivamente, das principais serras de Portugal Continental (Fonte: Instituto Nacional de Estatística, Portugal).

a) Importe este ficheiro de objectos.

b) Indique o número de serras de Portugal Continental existentes nestes dados.

c) Quais as serras que têm uma altitude inferior a 500 m?

d) Determine o nome e o número das serras que têm uma altitude entre os 800 e os 1000 m, inclusivé.

e) Determine a altitude da Serra da Estrela (Estrela).

27

Revisões de Estatística

Estatística - ciência que se ocupa da recolha, organização e análise de informação, com a finalidade de inferir de um conjunto limitado de informação para o todo e eventualmente prever a evolução futura de um fenómeno.

A estatística deve ser considerada um instrumento, uma ajuda à compreensão dos fenómenos, e nunca deve substituir a necessária e importante reflexão sobre o problema em estudo.

28

Dados estatísticos

Os dados (e respectiva variável) podem ser de natureza– qualitativa

• nominal (a ordem das categorias não tem significado)• ordinal (há uma ordenação natural das categorias)

– quantitativa• discreta (provêm de contagens)• contínua (provêm de medições)

Exemplos:– precipitação diária numa estação meteorológica– tipo de ocupação do solo numa bacia hidrográfica

29

Análise exploratória

A análise inicial de dados univariados tem como principais objectivos: – a exploração dos dados para descobrir/identificar

aspectos ou padrões de maior interesse;– a representação dos dados de modo a destacar esses

aspectos ou padrões:• condensar os dados observados sob a forma de quadros;• representar graficamente os dados;• calcular indicadores numéricos de localização e de dispersão.

Ao explorar dados bivariados o objectivo é:– o estudo da variação conjunta das variáveis, procurando

pôr em evidencia “relações” eventualmente existentes entre elas, utilizando

• indicadores numéricos;• análise gráfica.

30

Indicadores numéricos no R

mean(x)devolve a média dos elementos do vector x median(x)devolve a mediana dos elementos do vector x

quantile(x,probs=p)devolve o quantil de ordem p dos elementos do vector x; por omissão p=seq(0, 1, 0.25), isto é, devolve os extremos e os quartis de x

max(x)/min(x)devolve o máximo/mínimo dos elementos do vector x

range(x) devolve o vector c(min(x),max(x)) IQR(x) devolve a amplitude inter-quartil dos elementos

do vector x

31

Indicadores numéricos no R

var(x) devolve a variância dos elementos do vector x; é igual a sum((x-mean(x))^2)/(length(x)-1)

sd(x)devolve o desvio padrão dos elementos do vector x; é igual a sqrt(var(x))

summary(x) devolve os extremos, os quartis e a média do vector x

cov(x,y)devolve a covariância entre x e y cor(x,y)devolve o coeficiente de correlação (de

Pearson) de x e y cor(x,y, method="spearman")devolve o

coeficiente de correlação de Spearman de x e y

32

Exercício 3

1. Considere os vectores precip, temp e vento.a) Determine para cada vector o número total de elementos.

b) Determine indicadores numéricos de localização e de dispersão para estes dados meteorológicos.

c) Calcule o coeficiente de variação para cada um dos conjuntos de dados.

d) Calcular o coeficiente de assimetria para cada um dos conjuntos de dados.

e) Determine a precipitação anual medida nesta estação meteorológica.

Nota: coeficiente de variação:

coeficiente de assimetria:

C v=sx

C s=n∑ (x i− x)3

(n−1)(n−2) s3

33

Nota: Funções e valores NA Muitas das funções estatísticas anteriores podem ser

usadas com um argumento adicional, na seguinte formafunção(x, na.rm=FALSE)

em que na.rm é uma variável lógica que indica se os valores desconhecidos devem ser ou não ignorados.Caso o vector x inclua elementos NA,• se na.rm=FALSE (o que ocorre por omissão) a função

devolve NA ou uma mensagem de erro, • se na.rm=TRUE os valores NA são ignorados.

> y[1] NA 53 31 15 62> mean(y)> mean(y,na.rm=TRUE)

34

Factor

Um factor é um vector de categorias. É utilizado para dados de natureza qualitativa.

> nomes[1] "Ana" "Paulo" "Zé" "Maria" > sexo <- c("M","H","H","M")> sexo <- factor(sexo)> sexo[1] M H H MLevels: H M Os níveis (categorias) de um factor podem ser visualizados e

modificados através do atributo levels.> levels(sexo) <- c("Male","Female")> sexo[1] Female Male Male FemaleLevels: Male Female

35

Matrix

Uma matrix é uma colecção de dados, todos do mesmo tipo, referenciados por dois índices. É uma generalização para duas dimensões de um vector.

Quando a dimensão é superior a dois, a colecção é designada array. Uma matrix é um array de dimensão 2.

Uma matriz é definida pelo número de linhas (n), número de colunas (m) e um conjunto de (n×m) valores.

36

Matrix (cont.)

A função matrix()serve para criar matrizes. > M <- matrix(1:12,nrow=3,ncol=4)> M [,1] [,2] [,3] [,4][1,] 1 4 7 10[2,] 2 5 8 11[3,] 3 6 9 12 A opção byrow=TRUE faz com que os elementos sejam

dispostos por linha.> M <- matrix(1:12,3,4,byrow=TRUE)

37

Matrix (cont.)

Um objecto do tipo matrix tem associado o atributo dim que indica a sua dimensão, i.e., nº de linhas e nº de colunas.

> dim(M) # ver a dimensão de M [1] 3 4Um vector pode ser transformado numa matriz.> A <- c(3,2,-4,0,-1,8) # A é um vector > dim(A) <- c(3,2) # transforma A numa

matriz com 3 linhas e 2 colunas

> A [,1] [,2][1,] 3 0[2,] 2 -1[3,] -4 8

38

Matrix (cont.)

Uma matriz pode ser transformada num vector.> A [,1] [,2][1,] 3 0[2,] 2 -1[3,] -4 8> as.vector(A) # devolve um vector cujos

elementos são os da matriz

[1] 3 2 -4 0 -1 8

39

Matrix (cont.)

Uma matriz pode ser criada por concatenação de vectores e/ou matrizes:

• A função cbind faz a concatenação dos objectos por coluna> u <- 1:3 ; v <- 4:6> U <- cbind(u,v)

• A função rbind faz a concatenação dos objectos por linha> V <- rbind(u,v)

40

Matrix (indexação e selecção de componentes)

M[i,j] é o elemento de M que está na linha i e na coluna j > M[3,1] M[i,] é a linha i e M[,j] é a coluna j da matriz M. Estes

objectos são do tipo vector> M[2,] > M[,1] Para seleccionar parte de uma linha/coluna indicam-se os índices

correspondentes

> M[c(1,3),1] # igual a M[-2,1]

Se a selecção incluir mais do que uma linha e mais do que uma coluna, o resultado é uma submatriz

> M[2:3,c(1,3,4)]

41

List

Uma lista é uma colecção ordenada de objectos, que podem ser de tipos diferentes (vectores numéricos, vectores lógicos, matrizes, listas, funções, ...).

Chama-se componente a cada um dos objectos da lista. As componentes são numeradas e podem ter um nome associado.

> uma.lista <-list(um.vector=1:10, + uma.palavra="Olá",+ uma.matriz=matrix(1:20,ncol=5),+ outra.lista=list(a=5,b=factor(c("a","b")))) Muitas funções no R devolvem uma lista, ou seja, um conjunto

de resultados de tipo e dimensão distintos.

42

List (selecção de componentes)

As componentes de uma lista são indexadas e podem ser referidas pelo seu índice, utilizando parêntesis rectos duplos; se tiverem nome, também podem ser referidas pelo nome.

> uma.lista[[2]]> uma.lista$uma.palavra

43

Data frame

Uma data frame é um caso especial de lista em que as componentes têm todas o mesmo nº de elementos.

É um quadro de dados em que as colunas podem ser de tipos diferentes.

> meteo.Evora <-data.frame(ano=rep(2006, + 365), dia=doy, vento,temperatura=temp,+ precipitacao=precip)

> pb.Jan

44

Data frame (selecção de componentes)

Os elementos de uma data frame podem ser referidos indicando a linha e a coluna:

> meteo.Evora[3,1]> meteo.Evora[,2]

A notação nomeDataFrame$nomeComponente (como numa list) também pode ser usada:

> meteo.Evora$dia > meteo.Evora[[2]]

45

Data frame (função attach)

As consultas a data frames podem ser simplificadas utilizando a função attach()

A função attach() coloca o nome da data frame no search path (o que permite aceder directamente às colunas de uma data frame, sem necessidade de referir o seu nome).

> diaError: object 'dia' not found> attach(meteo.Evora)> dia[1] [1] 1 2 3 4 5 …> search() # o search path> detach(meteo.Evora)> diaError: object 'dia' not found

46

Leitura de ficheiros

Uma tabela de valores, registada num ficheiro de texto (ASCII), pode ser atribuída a uma data frame, através da função

read.table(ficheiro,header=FALSE,sep="", dec = ".",row.names, col.names, as.is = FALSE, na.strings = "NA")

em que,

ficheiro é o nome do ficheiro de dados (entre aspas); se este não se encontrar na pasta de trabalho, é necessário indicar o endereço (utilizando / ou \\ como separador de pastas);

header é uma variável lógica que indica se o ficheiro tem ou não os nomes das variáveis na primeira linha. Por omissão o valor é FALSE;

sep é o caracter que separa os valores em cada linha; por omissão é "" que corresponde a um ou mais espaços em branco;

47

Leitura de ficheiros (cont.)

dec é o caracter utilizado como separador decimal. Por omissão é o ponto;

row.names é um vector com os nomes das linhas. Por omissão 1,2,3,...;

col.names é um vector com os nomes das colunas. Por omissão V1,V2,V3,...;

as.is é uma variável lógica que controla a conversão das variáveis alfanuméricas em factores. Se o seu valor é TRUE a leitura mantém o tipo original dos dados;

na.strings é o valor usado para os dados desconhecidos no ficheiro e que será convertido em NA.

48

Leitura de ficheiros (exemplo)

O ficheiro de texto (com valores separados por vírgulas e com cabeçalho) “EstacMeteo.csv” contem dados referentes às estações meteorológicas automáticas do Instituto de Meteorologia de Portugal.Criar a data frame estac.meteo com o conteúdo deste ficheiro.

> estac.meteo <- read.table("EstacMeteo.csv",+ header=TRUE, as.is=TRUE,sep=",")

> estac.meteo <- read.csv("EstacMeteo.csv",+ as.is=TRUE)

49

Importar dados de uma folha de cálculo

Há packages que têm funções para importar e exportar dados de folhas de cálculo (por exemplo, xlsReadWrite, ROpenOffice).

Em alternativa, pode-se:– copiar os dados a importar para uma nova folha de

cálculo,deste modo evitam-se problemas com a importação de colunas e linhas desnecessárias

– salvar essa folha como ficheiro de texto separado por tabs (.txt) ou por vírgulas (.csv), com ou sem cabeçalho,

note que o separador utilizado pode ser a vírgula ou o ponto e vírgula, dependendo da configuração regional do Windows

– ler o ficheiro no R com a função read.table (ou read.csv)

50

Escrita em ficheiros

Para escrever num ficheiro o conteúdo de uma data frame ou de uma matriz usa-se a função

write.table(x,file="",quote=TRUE,sep=" ",na="NA",dec=".",row.names=TRUE,col.names=TRUE)em que,

x é o nome do objecto a guardar;

file é o nome do ficheiro, incluindo o caminho de pastas;

quote é uma variável lógica que indica se as variáveis alfanuméricas são ou não escritas entre aspas;

sep é o caracter que separa os valores em cada linha; por omissão é um espaço em branco;

51

Escrita em ficheiros (cont.)

na é a representação no ficheiro dos NA’s de x;

dec é o caracter utilizado como separador decimal. Por omissão é o ponto;

row.names , col.names ou são variáveis lógicas que indicam se o ficheiro vai ou não conter os nomes das linhas/colunas de x, ou são vectores alfanuméricos com os nomes a atribuir às linhas/colunas no ficheiro.

52

Escrita em ficheiros (exemplo)

Guardar no ficheiro “meteoEvora.txt” a data frame meteoEvora, utilizando como separador de campos o espaço em branco.

> write.table(meteo.Evora,"meteoEvora.txt",+ row.names=FALSE) Guardar no ficheiro “meteoEvora.csv” a data frame meteoEvora, utilizando como separador de campos a vírgula

> write.table(meteo.Evora,"meteoEvora.csv",+ row.names=FALSE, sep=",")> write.csv(meteo.Evora,"meteoEvora.csv",+ row.names=FALSE)

53

Ler/Escrever dados de/em outros softwares

No package foreign, há funções para importar e exportar dados em formato binário de Minitab, S, SAS, SPSS, Stata,Systat, ….

É também possível ler/escrever em Bases de Dados. Para mais informações consultar o manual “R Data Import/Export”:

menu Help → Manuals (in PDF) → R Data Import/Export

54

Gráficos

O R permite construir uma grande variedade de gráficos:> demo(graphics)> demo(persp) Os comandos para criar gráficos dividem-se

fundamentalmente em dois grupos:– funções gráficas de alto nível (que permitem criar um

novo gráfico);– funções gráficas de baixo nível (que permitem

acrescentar informação a um gráfico existente).

A função par permite ver e modificar uma lista de parâmetros gráficos.

> par() # devolve a lista dos valores dos parâmetros > ?par # informação sobre a função par

55

Função gráfica: plot plot(x), com x um vector numérico, produz um diagrama dos

valores de x versus o seu índice:> plot(temp)

plot(x,y), com x e y vectores numéricos, produz um diagrama de pontos de y versus x:

> attach(pg.Jan)> plot(Lisboa,Porto)> attach(estac.meteo)> plot(-Longitude,Latitude)

plot(f,from=a,to=b), com f uma função, traça o gráfico de f entre a e b (por omissão, a=0 e b=1):

> plot(sin)> plot(sin,-2*pi,2*pi)

56

Função gráfica: plot (argumentos)

Com a função plot podem ainda ser utilizados outros argumentos. O argumento type controla o tipo de gráfico produzido. O seu valor pode ser "p" (pontos) (valor por omissão), "l" (linhas), "b" (pontos ligados por linhas), "h" (linhas verticais dos pontos ao eixo das abcissas), "s" e "S" (função em escada).

> plot(precip,type="h")> plot(temp,type="l")> plot(Ano, Lisboa,type="h")

57

Função gráfica: plot (argumentos)

O título, sub-título e as legendas dos eixos num gráfico podem ser criados/alterados através dos argumentos main, sub, xlab e ylab, respectivamente. Os extremos dos eixos podem ser definidos com os argumentos xlim e ylim, cujos valores são vectores de duas componentes.

> plot(temp,type="l",main=+ "Temperatura média do ar em 2006",+ xlab="dia do ano",ylab="Temp (ºC)") Qualquer um dos argumento da função par pode ser

também utilizado como argumento da função plot para modificar as características do gráfico produzido.

58

Função gráfica: barplot

Cria um gráfico de barras> barplot(Lisboa,names.arg=Ano)

Um exemplo mais elaborado> LxPr.transp <- t(as.matrix(pb.Jan[,2:3]))> barplot(LxPr.transp,names.arg=Ano,+ beside=T, col=c("red","blue"),+ legend=c("Lisboa","Porto"),+ main="Precipitação mês de Janeiro (mm)")

Nota: a função colors() devolve o nome das cores que o R reconhece. Ver também palette()

59

Visualização de vários gráficos

As funções gráficas de alto nível criam, por omissão, um novo gráfico, apagando, se necessário, aquele que se encontra na janela gráfica. – Em alguns casos, a utilização do argumento add=TRUE

permite sobrepor diferentes gráficos na mesma janela gráfica.

> plot(sin,-pi,pi)> plot(cos,-pi,pi,add=TRUE,col="red",lty=2)

Pode-se utilizar a função par para visualizar simultaneamente vários gráficos.

> op<-par(mfrow=c(2,1))> plot(dia,temp,type="l");plot(dia,precip,type="l")> par(op) # recupera os valores iniciais dos parâmetros

60

Funções gráficas de baixo nível

Permitem acrescentar informação a um gráfico existente na janela gráfica:– as funções points(x,y) e lines(x,y) permitem

acrescentar, respectivamente, pontos e pontos ligados por linhas;

– a função abline(a,b) acrescenta uma recta de declive b e ordenada na origem a;

– as funções abline(v=x)e abline(h=y) permitem adicionar rectas verticais (de abcissa x) e horizontais (de ordenada y), respectivamente;

– a função legend permite acrescentar uma legenda ao gráfico.

61

Funções gráficas de baixo nível (exemplos) Um exemplo> plot(Ano,Lisboa,ylab="precipitação", type="b")> points(Ano,Porto,col="red",type="b")> legend(1980,350,c("Lisboa","Porto"),+ col=c("black","red"),lty=1,pch=1)

Outro exemplo mais complexo> plot(dia,temp,type="l",ylab="Temperatura média+ diária (T, mm)")> abline(h=mean(temp))> abline(h=c(mean(temp)-sd(temp),+ mean(temp+sd(temp)), lty=2,col="red")> legend(0,30,c(expression(bar(T)),+ expression(bar(T) %+-% s[T])), + col=c("black","red"),lty=1:2)

62

Guardar gráficos

Os gráficos produzidos podem ser guardados em ficheiros de vários tipos (pdf, jpeg, ps, …):– num ficheiro de format pdf> pdf("exemploG1.pdf")> plot(Ano,Lisboa,type="h")> dev.off()

– em formato ps segue-se um procedimento idêntico > postscript("exemploG1.ps")...

ou utiliza-se o botão direito do rato e escolhe-se a opção “Save as postscript ...”

63

Histograma

A função hist(x), cria um histograma com os elementos do vector x.

> par(mfrow=c(1,2)) > hist(Porto); hist(Lisboa) O parâmetro breaks=x da função hist() permite definir

o número de rectângulos do histograma (por omissão, calculado pela regra de Sturges, 1 + log2(n)) em que :

– x é o número de classes em que se pretende agrupar os dados, ou

– x é um vector numérico com os limites das várias classes.> par(mfrow=c(2,2))> hist(Pg)> hist(Pg,breaks=seq(300,1700,by=200))> hist(Pg,breaks=seq(300,1600,by=100))> hist(Pg,breaks=seq(300,1600,by=100),prob=T)

64

A função hist O parâmetro prob permite controlar se a área total do

histograma é um (prob=TRUE) ou se o histograma é construído com base nas frequências absolutas (prob=FALSE)– por omissão, prob=FAlSE sse as classes têm igual amplitude

Quando utilizada com o parâmetro plot=FALSE, a função hist devolve uma lista

> hist(Pg,breaks=seq(300,1600,by=100),plot=F)breaks – vector com os limites das classescounts – vector com as frequência absoluta das classesintensities – vector com (frequência relativa/amplitude) das classesdensity – idemmids – vector com os pontos médios das classesequidist – valor lógico que indica se as classes têm ou não amplitude

constante

65

Caixa-de-bigodes (box plot)

É um gráfico que permite representar em simultâneo medidas de localização e dispersão de uma amostra:– transmite de modo imediato uma ideia da localização,

dispersão e forma da população de que foi extraída a amostra;

– é um instrumento adequado para comparar várias amostras.

Estes gráficos podem ser obtidos em R utilizando a função boxplot().

> boxplot(Pg)> boxplot(pg.Jan[,2:3])

66

Exercício 4

Considere os dados do vector serras.altitude.a) Construa um diagrama de extremos e quartis e um

histograma.b) Construa três novos histogramas: um com classes de

amplitude de 100m, outro com classes de amplitude de 250m e outro com classes de amplitude 200m mas cujo o primeiro intervalo de classe é ]100,300].

c) Apresente os quatro histogramas construidos nas alíneas anteriores na mesma janela gráfica (sugestão: utilize a função par).

d) Dos quatro histogramas anteriores, qual permite ter uma melhor representação gráfica da distribuição das altitudes das serras de Portugal Continental? Justifique.

67

Packages

Todas as funções e comandos do R estão armazenados em módulos (packages).

Numa sessão de R, o conteúdo de um package só fica disponível quando este é carregado em memória.– Para ver quais os packages disponíveis: >(.packages())– Para ver quais os packages instalados localmente: >library()– Para carregar em memória um package já instalado: >library(nome_package)

ou menu Packages → Load Package ... O R inclui um conjunto de packages desenvolvidos pelos

utilizadores para implementar métodos específicos, quer na área da estatística, quer noutras áreas (por ex. R Commander-Rcmdr). – Para instalar um package fazer:

menu Packages → Install Package(s) ...

68

Package: fBasics

fBasics: estatísticas básicas e de mercado (ambiente criado para o ensino de “Engenharia Financeira” e “Finanças Computacionais”)

Indicadores de forma: – coeficiente de assimetria

> skewness(precip); skewness(pg.PGarcia$Pg)

– coeficiente de achatamento> kurtosis(precip); kurtosis(pg.PGarcia$Pg)

Resumos de indicadores:> basicStats(precip)> basicStats(pg.PGarcia$Pg)

69

Package: Rcmdr

R commander: uma interface gráfica do R para análises estatísticas básicas.

Para um correcto funcionamento em Windows, o R commander deve ser utilizado quando o R está a ser executado numa janela gráfica de um único documento (SDI). – Pode-se alterar as propriedades de um atalho do R para iniciar

esta aplicação em modo SDI: acrescentar --sdi ao caminho indicado em Target.

A função Commander() inicia a interface gráfica do R commander

70

Modelos probabilísticos no R

O software R disponibiliza “tabelas estatísticas” para um conjunto de distribuições.

Para cada distribuição estatística é possível calcular1. função massa de probabilidade ou função densidade 2. função distribuição cumulativa3. quantis4. números pseudo-aleatórios

Nota: o quantil de ordem p é o menor valor de x tal que F(x) ≥ p, sendo F a função distribuição cumulativa

71

Modelos probabilísticos no R

Distribution R name additional argumentsbeta beta shape1, shape2, ncpbinomial binom size, probCauchy cauchy location, scalechi-squared chisq df, ncpexponential exp rateF f df1, df2, ncpgamma gamma shape, scalegeometric geom probhypergeometric hyper m, n, klog-normal lnorm meanlog, sdloglogistic logis location, scaleneg. binomial nbinom size, probnormal norm mean, sdPoisson pois lambdaStudent's t t df, ncpuniform unif min, maxWeibull weibull shape, scaleWilcoxon wilcox m, n

72

Modelos probabilísticos no R

Ao nome R da distribuição acrescentam-se os prefixos:– d para função massa de probabilidade ou função densidade – p para função distribuição cumulativa– q para quantis– r para números pseudo-aleatórios

As funções pxxx e qxxx têm o argumentos lógico: lower.tail – se for TRUE (por omissão), as probabilidades são P[X ≤ x], caso contrário são P[X > x].

73

Distribuição binomial

n inteiro positivo e 0 < p < 1

– pode tomar os valores

– com probabilidade

Exemplo: Seja

– > dbinom(2,4,0.3)– > pbinom(2,4,0.3)

Alguns gráficos da função massa de probabilidade: > x <- 0:4> plot(x,dbinom(x,4,0.3),type="h")> x <- 0:50> plot(x,dbinom(x,50,0.3),type="h")

X ∩B n , p

P [ X =x ]=nx px 1− pn−x

x=0, 1,2,⋯ , n

X ∩B 4, 0.3

P [ X =2]=0.2646P [ X ≤2]=0.9163

74

Exercício 5

a) Qual a probabilidade do caudal instantâneo de um dado curso de água atingir ou ultrapassar (pelo menos uma vez) o valor associado a um período de retorno de T=50 anos durante os n=10 anos de construção de uma barragem nesse curso de água?

b) Construa um gráfico do risco R (a probabilidade de um determinado acontecimento hidrológico ser igualado ou excedido pelo menos uma vez em n anos sucessivos) em função do período de retorno T.Nota: utilize escalas logarítmicas nos dois eixos com a função par(xlog=TRUE,ylog=TRUE) e plot(...,log="xy")

75

Exercício 5 b)

76

Distribuição normal

– pode tomar qualquer valor real– função densidade– valor médio E[ X ] = μ– variância Var[ X ] = σ2

Exemplo: Seja – > pnorm(2,2,4)– > 1-pnorm(2.5,2,4)

Alguns gráficos da função densidade:> curve(dnorm(x), from=-6, to=13)> curve(dnorm(x,0,2), from=-6, to=6, add=T, col=2)> curve(dnorm(x,8,1), from=-6, to=13, add=T, col=3)

X ∩N ,

f ( x)= 1√2πσ

e−1

2 ( x−μσ )2

, x∈ℝ

∈ℝ , 0

P [ X ≤2]=0.5P [ X >2.5]=0.4502618

X ∩N (2, 4)

77

Outras distribuições com interesse emhidrologia Log-normal, Gumbel, Pearson tipo III, Weibull ...

– aparecem em diferentes packages: stats, evd, lmom, lmomco, …

Exemplo: Seja – > plnorm(2,0,1)– > 1-plnorm(2.5,0,1)

Alguns gráficos da função densidade:> curve(dlnorm(x),from=0,to=6,ylim=c(0,2))> curve(dlnorm(x,0,0.5),from=0,to=6,add=T,col=2)> curve(dlnorm(x,0,0.25),from=0,to=6,add=T,col=3)> curve(dlnorm(x,1,1),from=0,to=6,add=T,col=4)

P [ X ≤1]=0.5P [ X >2.5]=0.1797572

X ∩ln N (1, 2)⇔ ln ( X )∩N (1, 2)

78

Qual a distribuição teórica que se ajusta?

Métodos gráficos:– Histograma e função densidade (estimada e teórica)

> attach(pg.PGarcia)> hist(Pg)> hist(Pg,breaks=seq(200,1600,100),prob=TRUE)> lines(density(Pg)) # adiciona curva densidade estimada por

kernel> curve(dnorm(x,mean(Pg),sd(Pg)),from=200,to=1600,+ add=TRUE,lwd=2,col=2)

– Função distribuição cumulativa empírica> plot(ecdf(Pg))> plot(ecdf(Pg), verticals=TRUE, do.points=FALSE)> curve(pnorm(x,mean(Pg),sd(Pg)),from=200,to=1600,+ add=TRUE,lwd=2,col=2)

79

Qual a distribuição teórica que se ajusta?

Métodos gráficos (cont.):– Diagrama quantil-quantil (q-q plot)

• é um diagrama de pontos dos quantis de uma amostra versus os quantis, da mesma ordem, de outra amostra (ou distribuição teórica).

– Se as amostras têm a mesma dimensão, o q-q plot não é mais do que um diagrama de pontos das amostras ordenadas.

• Objectivo: avaliar se as duas amostras provêm de populações com uma distribuição comum ou se a amostra se ajusta à distribuição teórica considerada.

– Neste caso, os pontos devem-se situar aproximadamente ao longo da recta bissectriz dos quadrantes ímpares; quanto mais se afastarem desta linha de referência, mais provável é as duas populações terem diferentes distribuições.

80

Qual a distribuição teórica que se ajusta?

Métodos gráficos (cont.):– Diagrama quantil-quantil (q-q plot)

Exemplo: com a precipitação total anual de Penha Garcia

> par(pty="s") # gráfico numa região quadrada> qqnorm(Pg) # qqplot com a distribuição normal

idêntico ao comando seguinte> qqplot(qnorm(ppoints(length(Pg)),Pg)> qqline(Pg) # acrescenta uma linha recta ao

gráfico anterior que passa no 1º e 3º quartis

Nota: a função ppoints (também utilizada por qqnorm) gera uma sequência de n pontos de probabilidade utilizando a fórmula de Cunnane

(1:n – a)/(n + 1 – 2a) com a = 3/8 se n ≤ 10 ou a = 1/2 caso contrário.

81

Qual a distribuição teórica que se ajusta?

Outro exemplo: Pretende analisar-se se a distribuição de Gumbel se ajusta adequadamente à precipitação diária máxima anual registada em Penha Garcia, Portugal (dados no ficheiro “Pmad.PGarcia.csv”)

Relembrando a distribuição Gumbel:

– função densidade– valor médio– Variância

Nota: as funções associadas à distribuição de Gumbel (dgumbel, pgumbel, qgumbel, rgumbel) encontram-se no package evd.

X ∩Gum (α ,β) ,α∈ℝ , β>0

f (x)= 1α exp [− x−β

α −exp(− x−βα )] , x∈ℝ

E [ X ]=α+0.5772β

Var [ X ]=π2

6β2

82

Qual a distribuição teórica que se ajusta?

Outro exemplo (cont):

> Pmad.PGarcia<-read.csv("Pmad.PGarcia.csv")> lybrary(evd)> ?qgumbel# Cálculo das estimativas dos parâmetros de# escala (β) e localização (α)> escala<-sqrt(6)*sd(Pmad.PGarcia$Pmad)/pi> local<-mean(Pmad.PGarcia$Pmad)-+ 0.5772*escala> qqplot(qgumbel(ppoints(66),loc=local,+ scale=escala),Pmad.PGarcia$Pmad)> abline(a=0,b=1)

83

Qual a distribuição teórica que se ajusta?

Outro exemplo: simulação de 66 valores provenientes de uma Dist. de Gumbel com os mesmos parâmetros de localização e de escala do exemplo anterior

# fixar a semente do gerador de nº aleatórios> set.seed(1234) > x <- rgumbel(66,loc=local,scale=escala)> qqplot(qgumbel(ppoints(66),loc=local,+ scale=escala),x)> abline(a=0,b=1)

84

Qual a distribuição teórica que se ajusta?

Testes de ajustamento:– Teste de Shapiro-Wilk (teste de normalidade)shapiro.test(...)

– Teste do qui-quadrado de Pearsonchisq.test(...)

– Teste de Kolmogorov-Smirnovks.test(...)

Exemplo:> shapiro.test(Pg)> ks.test(Pmad.PGarcia$Pmad, + "pgumbel", loc = 45, scale=16.3)

85

Testes de hipóteses:

Procedimento:1º Formular duas hipóteses:

Hipótese nula (H0) onde se especifica o valor do parâmetro ou distribuição a verificar (a que se dá o beneficio da dúvida)Hipótese alternativa (H1) onde se específica a hipótese a considerar se rejeitarmos H0 (a que se dá o ónus da prova)

2º Identificar a estatística do teste (e as condições em que é válida).

3º Definir o nível de significância (α)

4º Calcular a estatística do teste5º Conclusão: rejeitar H0 se o valor p (p-value) ≤ α

o p-value é a probabilidade da estatística do teste tomar valores mais extremos que o valor calculado com base na amostra, sob H0 verdadeira.

α=P [Rejeitar H 0∣ H 0 verdade ]

86

Alguns testes paramétricos no R

Testes paramétricos – supõe-se conhecida a distribuição teórica da população em estudo e as hipóteses a formular dizem respeito ao(s) parâmetro(s).

t.test(…) Teste t-Student para o valor da média de uma população normal ou para a comparação das médias de duas populações (utilizando amostras independentes ou amostras emparelhadas)

Exemplo: Comparação da precipitação média no mês de Janeiro nas cidades do Porto e Lisboa (amostras emparelhadas – factor de emparelhamento: ano)> attach(pg.Jan)> shapiro.test(Lisboa-Porto) #teste à normalidade

da população das diferenças> t.test(Lisboa, Porto, pair=TRUE)var.test(…) Teste F à razão das variâncias de duas populações normais.

87

Exercício 6

É desencadeado um programa de controlo da poluição de um rio em que foram efectuadas medições, antes de lançar a campanha anti-poluição e um ano após. As medições são combinações de vários índices: quanto maior for o valor resultante maior é a poluição. Obtiveram-se os seguintes resultados:

Será que a campanha anti-poluição reduziu de facto a poluição? Explicite e verifique todas as hipóteses necessárias à resolução do problema

Ponto de controlo

1 2 3 4 5 6 7 8 9 10

Antes da campanha

68 88 101 82 96 74 65 74 52 99

Um ano após 67 87 90 76 98 69 68 65 59 70

88

O modelo linear

Designação genérica que engloba vários modelos estatísticos como:– a regressão linear (simples e múltipla)– a análise de variância– a análise de covariância

Relaciona, Y - a variável resposta (dependente, modelada), com x1, x2, …, xp - as variáveis preditoras (independentes, explicativas).

Pressupostos,Y i=β0+β1 x1(i)+β2 x2(i)+...+βp x p(i)+ϵi , i=1,... , nϵi∩N (0,σ2) ,∀i

{ϵ}i são v.a. independentes

89

O modelo linear no R

lm(…) função básica no estudo do modelo linear (linear model). Tem como primeiro argumento (obrigatório) um objecto do tipo formula que indica qual a variável resposta (e.g., y) e quais as variáveis preditoras (e.g., x1, x2, x3)

y ~ x1 + x2 + x3 Exemplos:> lm(y ~ x) # RLS de y sobre x > lm(y ~ -1+x) # RLS de y sobre x com ordenada na

origem nula (recta passa na origem)> lm(y ~ x1+x2+x3)# RLM de y sobre x1, x2, x3> lm(log(y) ~ x1+x2)# RLM da variável transformada,

log(y), sobre x1, x2

90

O modelo linear no R (cont.).

A função lm produz um objecto de classe “lm”, uma list com várias componentes. A maioria dessas componentes não é mostrada automaticamente.

Informação sobre o modelo ajustado pode ser obtida utilizando várias funções como, por exemplo,summary devolve um resumo alargado do resultado da

análise de regressão;coef devolve o vector das estimativas dos fitted devolve os valores ajustados da variável respostaresiduals devolve os resíduosconfint devolve intervalos de confiança para cadapredict permite ajustar novas observações à regressão

ajustada ou obter intervalos de confiança para o valor esperado da variável resposta ou intervalos de predição.

(e i= yi− yi)( yi)

βi

βi

91

Um exemplo

Um dos métodos de estimar a capacidade de armazenamento de um copado arbóreo (parâmetro utilizado na modelação da intercepção) utiliza a relação linear que se pode estabelecer entre a precipitação (bruta) e precipitação sob-coberto para chuvadas independentes com completa saturação do copado.

> pinhal<- read.csv("Pinhal.csv")> pinhal[1:3,]> summary(pinhal)> boxplot(pinhal)> cor(pinhal)> plot(pinhal) # idêntico a plot(pinhal$tf ~ pinhal$pg)> lm(tf ~ pg, data=pinhal)> pinhal.lm<-lm(tf ~ pg, data=pinhal)# recta ajustada

Tf = -0.3798 + 0.86642 pg> abline(pinhal.lm)> summary(pinhal.lm)

92

Os resultados de summary (lm(...)) Coefficientes: Estimate – estimativa de cada parâmetro do modelo

Std.Error – estimativa do erro padrão (desvio padrão associado a cada estimador dos parâmetros)t-value – valor da estatística num teste-t às hipótesesPr(>|t|) – p-value associado ao teste anterior

Residual standard error – estimativa do desvio padrão erros aleatórios (ε) e respectivos graus de liberdade.

Multiple R-squared – coeficiente de determinação (R2) Adjusted R-squared – R2 modificado

F-statistic – estatística do teste F de ajustamento global com as hipótesesrespectivos graus de liberdade e p-value

H 0 : β j=0 vs. H 1 : β j≠0

H 0 : β j=0, ∀ j=1,... , p vs. H 1: ∃ j t.q. β j≠0

=1−(1−R2)× n−1n−( p+1)

93

Análise de resíduos e outros diagonósticos

Para identificar eventuais violações grosseiras aos pressupostos do modelo é conveniente analisar os resíduos.

Outras ferramentas de diagnóstico visam identificar observações individuais com características especiais (por exemplo, distância de Cook e “observações alavanca”- leverage points)

O R disponibiliza uma forma gráfica rápida para fazer esta análise: aplicar o comando plot a um objecto da classe “lm”> op <- par(mfrow=c(3,2))> plot(pinhal.lm,which=1:6)> par(op)

94

Exercício 7

O ficheiro “Homogeneidade.csv” contém dados de precipitação anual (cm). A série AY corresponde a valores médios das precipitações anuais medidas em 25 estações homogéneas e a série Z às precipitações anuais registadas na estação z, situada na mesma região das anteriores. Sabe-se que na década de 1930 a estação z foi deslocada de lugar.

a) Utilizando o método da dupla-massa, verifique se a série Z é homogénea.b) Caso não haja homogeneidade, efectue as correcções necessárias para que a série Z fique homogénea.

95

Mais algumas informações

O package topmodel faz a implementação deste modelo hidrológico em Rhttp://cran.r-project.org/web/packages/topmodel/index.html

Outro package com interesse EcoHydRologyhttp://cran.r-project.org/web/packages/EcoHydRology/index.html

Para mais informações sobre análise de dados hidrológicos e modelação http://rwiki.sciviews.org/doku.php?id=guides:tutorials:hydrological_data_analysis

96

Bibliografia geral para o R

Documentação disponível no CRAN (Comprehensive R Archive Network - http://www.r-project.org/):

The R Manuals, editados por, R Development Core Team, Current Version: 2.13.1 (Julho 2011)

(http://cran.r-project.org/manuals.html).

Contributed Package: Available Bundles and Packages(http://cran.r-project.org/web/packages/)

Contributed Documentation, escritos por utilizadores de R. Por exemplo, “Introdução à Programação em R” por Luis Torgo. (http://cran.r-project.org/other-docs.html)

Fox, J., Getting Started With the R Commander (August 2006) (http://socserv.socsci.mcmaster.ca/jfox/Misc/Rcmdr/Getting-Started-with-the-Rcmdr.pdf)