**[[.:00-inicio|Uso da Linguagem R para Análise de Dados em Ecologia]]** ====== 3. Funções Matemáticas e Estatísticas ====== ===== O R como uma Potente Calculadora ===== ==== Operações Aritméticas Básicas ==== A linha de comando do R funciona como uma potente calculadora. Todas operações aritméticas e funções matemáticas principais estão disponíveis. Exemplo: > 4 + 9 [1] 13 > 4 - 5 [1] -1 > 4 * 5 [1] 20 > 4 / 5 [1] 0.8 > 4^5 [1] 1024 > A notação básica de operações algébricas, como a aplicação hierárquica de parênteses, também pode ser utilizada: > (4 + 5 ) * 7 - (36/18)^3 [1] 55 > (2 * ( 2 * ( 2 * (3-4)))) [1] -8 > Note que somente os parênteses podem ser utilizados nas expressões matemáticas. As chaves ("{}") e os colchetes ("[]") têm outras funções no R: > (2 * { 2 * [ 2 * (3-4)]}) Error: syntax error in "(2 * { 2 * [" > Por que o R é uma **potente** calculadora? Experimente fazer a seguinte operação matemática na sua calculadora: > 1 - (1 + 10^(-15)) === Exercícios === Pollard (1971) propôs o seguinte estimador para estimar a densidade no método de quadrantes: \widehat{N}_p = \frac{4(4n-1)}{\pi \sum_{i=1}^{n}\sum_{j=1}^4 r_{ij}^2} onde, r_{ij} é a distância de árvore do quadrante //j// no ponto //i// ao centro do ponto quadrante e //n// é o número de pontos quadrantes. A variância desse estimador é: \textrm{Var}\{ \widehat{N}_p \} = \frac{\widehat{N}_p}{4n-2}. Se o valor da soma do quadrado das distâncias for: \sum_{i=1}^{30}\sum_{j=1}^4 r_{ij}^2} = 2531.794 qual a densidade estimada e sua variância? A área transversal de uma árvore é calculada assumindo que a secção transversal do tronco à altura do peito (1,3m) é perfeitamente circular. Se o diâmetro à altura do peito (DAP) de uma árvore for 13.5cm, qual a área transversal? Se uma árvore possui três fustes com DAPs de: 7cm, 9cm e 12cm, qual a sua área transversal? ==== Funções Matemáticas Comuns ==== As funções matemáticas comuns também estão disponíveis e podem ser aplicadas diretamente na linha de comando: > sqrt(9) # Raiz Quadrada [1] 3 > abs( - 1 ) # Módulo ou valor absoluto [1] 1 > abs( 1 ) [1] 1 > log( 10 ) # Logaritmo natural ou neperiano [1] 2.302585 > log( 10, base = 10) # Logaritmo base 10 [1] 1 > log10(10) # Também logaritmo de base 10 [1] 1 > log( 10, base = 3.4076) # Logaritmo base 3.4076 [1] 1.878116 > exp( 1 ) # Exponencial [1] 2.718282 > As funções trigonométricas: > sin(0.5*pi) # Seno [1] 1 > cos(2*pi) # Coseno [1] 1 > tan(pi) # Tangente [1] -1.224647e-16 > > asin(1) # Arco seno (em radianos) [1] 1.570796 > asin(1) / pi * 180 [1] 90 > > acos(0) # Arco coseno (em radianos) [1] 1.570796 > acos(0) / pi * 180 [1] 90 > atan(0) # Arco tangente (em radianos) [1] 0 > atan(0) / pi * 180 [1] 0 > Funções para arredondamento: > ceiling( 4.3478 ) [1] 5 > floor( 4.3478 ) [1] 4 > round( 4.3478 ) [1] 4 > round( 4.3478 , digits=3) [1] 4.348 > round( 4.3478 , digits=2) [1] 4.35 > Funções matemáticas de especial interesse estatístico: > factorial( 4 ) # Fatorial de 4 [1] 24 > choose(10, 3) # Coeficientes binomiais: combinação de 10 3-a-3 [1] 120 > === Exercícios === Se uma árvore possui três fustes com DAPs de: 7cm, 9cm e 12cm, qual o diâmetro (único) que é equivalente à sua área transversal? O modelo alométrico de biomassa ajustado para árvores do Cerradão estabele que a biomassa é dada pela expressão: \widehat{b_i} = \exp(-1.7953) \, d_i^{2.2974} onde b_i é a biomassa em //kg// e d_i é o DAP em //cm//. Já um outro modelo para biomassa das árvores na mesma situação tem a forma: \widehat{\ln(b_i)} = -2.6464 + 1.9960 \ln(d_i) } + 0.7558 \ln(h_i) onde h_i é a altura das árvores em //m//. Para uma árvore com DAP de 15//cm// e altura de 12//m//, os modelos resultarão em estimativas muito distintas? ==== Mantendo a Coerência Lógica-Matemática ==== O R também lida com operações matemáticas que envolvem **elementos infinitos** e **elementos indeterminados**: > 1/0 [1] Inf > -5/0 [1] -Inf > 500000000000000000/Inf [1] 0 > 0/0 [1] NaN > Inf/Inf [1] NaN > log(0) [1] -Inf > exp(-Inf) [1] 0 > sqrt(Inf) [1] Inf > sqrt( - 1 ) [1] NaN Warning message: NaNs produced in: sqrt(-1) > 2 * NA [1] NA > 2 * NaN [1] NaN > NA / 10 [1] NA > NaN / -1 [1] NaN > Note que determinadas **palavras** (além do nome das funções) estão reservadas no R, pois são utilizadas com significado especial: * ''pi'' - constante \pi = 3.141593 ; * ''Inf'' - infinito; * ''NaN'' - indeterminado (Not a Number), normalmente resultado de uma operação matemática indeterminada; * ''NA'' - indeterminado (Not Avaiable), normalmente caracterizando uma observação perdida (//missing value//). Na operações matemáticas, ''NaN'' e ''NA'' atuam sempre como **indeterminado**. === Exercícios === Como se caracteriza uma **observação perdida**? Quando o diâmetro de uma árvore deve ter o valor **zero** ou o valor **NA**? E o peso de um animal? E a biomassa de uma florestal? E a espécie de uma ave? ==== Criando Variáveis ==== Mais do que simples operações aritméticas, o R permite que executemos operações **algébricas** operando sobre variáveis prédefinidas. Para definir uma variável, basta escolher um nome (//lembre-se das regras de nomes no R//) e atribuir a ela um valor: > > a = 3.6 > b = sqrt( 35 ) > c = -2.1 > a [1] 3.6 > b [1] 5.91608 > c [1] -2.1 > > a * b / c [1] -10.14185 > b^c [1] 0.02391820 > a + exp(c) - log(b) [1] 1.944782 > > a - b * c / d Error: object "d" not found > Não esqueça de definir as variáveis previamente!! === Exercícios === O que acontece se você criar uma variável com o nome ''pi''? Por exemplo, > pi = 10 O que acontece com a constante \pi? E se for criada uma constante de nome ''sqrt''? O que acontece com a função raíz quadrada (''sqrt()'')? **DICA:** O que faz a função ''search'', no comando: > search() ==== Sinais de Atribuição ==== Um **detalhe importante** em linguagens de programação é diferenciar um **sinal de igualdade** de um **sinal de atribuição**. O sinal de igualdade faz uma comparação entre dois elementos. No R o sinal de igualdade é **"=="** e será visto em detalhes quando a questão de comparações for tratada. O R possui vários sinais de atribuição. No tópico acima utilizamos o sinal '''=''', mas o mesmo efeito pode ser obtido por dois outros meios: * o sinal **"<****-"** (composto pelo sinal de menor **"<"** seguido do hífen **"-"**) atribui o valor **à direita** ao elemento **à esquerda**: > k <- 2.46 > k [1] 2.46 *o sinal **"-****>"** (composto pelo hífen **"-"** seguido do sinal de maior **">"**) atribui o valor **à esquerda** ao elemento **à direita**: > 2.46 -> k2 > k2 [1] 2.46 > Essas três formas de atribuir valores a variáveis são equivalentes e cada usuário trabalha com a forma que julgar mais conveniente. ===== O R como uma Calculadora Vetorial ===== ==== Criando Vetores ==== O R, e a linguagem S, foram criados para operar não apenas //número-a-número// como uma calculadora convencional. O R é um ambiente **vetorial**, isto é, quase todas suas operações atua sobre uma //seqüência de números//, que genericamente chamaremos de vetores. Para criar um vetor, podemos usar a função ''c'' (c = colar, concatenar). Essa função simplesmente cola todos os argumentos dados a ela, formando um vetor: > a = c(1, 10, 3.4, pi, pi/4, exp(-1), log( 2.23 ), sin(pi/7) ) > a [1] 1.0000000 10.0000000 3.4000000 3.1415927 0.7853982 0.3678794 0.8020016 0.4338837 > Para criar vetores de números com intervalo fixo unitário (intervalo de 1) se utiliza o //operador seqüencial// ('':''): > b = 1:8 > b [1] 1 2 3 4 5 6 7 8 > c = 20:32 > c [1] 20 21 22 23 24 25 26 27 28 29 30 31 32 > d = 2.5:10 > d [1] 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 Uma forma mais flexível de criar seqüências de números (inteiros ou reais) é usando a função '''seq''': > seq(10, 30) [1] 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 > seq(10, 30, by=2) [1] 10 12 14 16 18 20 22 24 26 28 30 > seq(1.5, 7.9, length=20) [1] 1.500000 1.836842 2.173684 2.510526 2.847368 3.184211 3.521053 3.857895 [9] 4.194737 4.531579 4.868421 5.205263 5.542105 5.878947 6.215789 6.552632 [17] 6.889474 7.226316 7.563158 7.900000 Também é fácil criar uma seqüência de números repeditos utilizando a função '''rep''': > rep(5, 3) [1] 5 5 5 > rep(1:5, 3) [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 > rep(1:5,each=3) [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 > === Exercícios === Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7. Crie um vetor '''dap''' com os diâmetros acimas e uma seqüência que enumera os fustes. ==== Vetores: Operações Matemáticas ==== Todas operações matemáticas aplicadas sobre um vetor, serão aplicadas sobre cada elemento desse vetor: > 2 * a [1] 2.0000000 20.0000000 6.8000000 6.2831853 1.5707963 0.7357589 1.6040032 [8] 0.8677675 > sqrt( a ) [1] 1.0000000 3.1622777 1.8439089 1.7724539 0.8862269 0.6065307 0.8955454 [8] 0.6586985 > > log( a ) [1] 0.0000000 2.3025851 1.2237754 1.1447299 -0.2415645 -1.0000000 -0.2206447 [8] -0.8349787 > Se as variáveis que trabalhamos são vetores, operações matemáticas entre variáveis serão realizadas pareando os elementos dos vetores: > a* b [1] 1.000000 20.000000 10.200000 12.566371 3.926991 2.207277 5.614011 [8] 3.471070 > a - b [1] 0.0000000 8.0000000 0.4000000 -0.8584073 -4.2146018 -5.6321206 -6.1979984 [8] -7.5661163 > a^(1/b) [1] 1.0000000 3.1622777 1.5036946 1.3313354 0.9528356 0.8464817 0.9689709 [8] 0.9008898 > > sqrt( a ) [1] 1.0000000 3.1622777 1.8439089 1.7724539 0.8862269 0.6065307 0.8955454 [8] 0.6586985 > log( b ) [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101 [8] 2.0794415 > Nesse tipo de operação é importante termos clareza sobre o que pretendemos quando fazemos uma operação matemática entre dois vetores de comprimentos (''length'') diferentes: > length( a ) [1] 8 > length( b ) [1] 8 > length( c ) [1] 13 > a * c [1] 20.000000 210.000000 74.800000 72.256631 18.849556 9.196986 [7] 20.852041 11.714861 28.000000 290.000000 102.000000 97.389372 [13] 25.132741 Warning message: longer object length is not a multiple of shorter object length in: a * c > === Exercícios === Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7. Calcule a área transversal de cada fuste dessa palmeira. Como construir uma seqüência que representa o aumento do número de bits por byte de computador, quando se dobra o tamanho dos bytes? Essa seqüência numérica parte do 2 e dobra os valores a cada passo. ==== Vetores: Operações Estatísticas ==== Operações matemática sobre vetores, opera //elemento-a-elemento// do vetor. Já as funções estatísticas operam no vetor **como um todo**: > mean( a ) [1] 2.491344 > var( b ) [1] 6 > max( c ) [1] 32 > sd( a ) [1] 3.259248 > sum( c ) [1] 338 > min( b ) [1] 1 > range( c ) [1] 20 32 > Algumas funções úteis que não são estatísticas, mas operam no vetor são: > a [1] 1.0000000 10.0000000 3.4000000 3.1415927 0.7853982 0.3678794 0.8020016 [8] 0.4338837 > sort(a) [1] 0.3678794 0.4338837 0.7853982 0.8020016 1.0000000 3.1415927 3.4000000 [8] 10.0000000 > rev(sort(a)) [1] 10.0000000 3.4000000 3.1415927 1.0000000 0.8020016 0.7853982 0.4338837 [8] 0.3678794 > cumsum(sort(a)) [1] 0.3678794 0.8017632 1.5871613 2.3891629 3.3891629 6.5307556 9.9307556 [8] 19.9307556 > cumsum(a) [1] 1.00000 11.00000 14.40000 17.54159 18.32699 18.69487 19.49687 19.93076 > > diff(a) [1] 9.0000000 -6.6000000 -0.2584073 -2.3561945 -0.4175187 0.4341221 -0.3681178 > diff( seq(10, 34, length=15) ) [1] 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 [9] 1.714286 1.714286 1.714286 1.714286 1.714286 1.714286 > === Exercícios === Uma palmeira perfilhada possui 10 fustes com os seguintes diâmetros: 5, 6, 7, 5, 10, 11, 6, 8, 9 e 7. Encontre para essa palmeira as seguintes medidas: - a área transversal total (área basal) da palmeira (//g//); - o diâmetro médio dos fustes (\overline{d}); - o desvio padrão dos diâmetros (//s//); Como podemos representar essa palmeira com vários fustes por uma única medida de diâmetro? Calculemos a área transversal a partir de três opções: A) A partir do diâmetro total, ou soma dos diâmetros: g_A = \frac{\pi}{4} \left( 10\,\overline{d} \right)^2 B) A partir do diâmetro médio e o número de fustes: g_B = 10\, \frac{\pi}{4} \left( \overline{d} \right)^2 C) A partir do diâmetro médio e do desvio padrão: g_C = \frac{\pi}{4} \left( (10)\overline{d}^2 + (10-1)s^2 \right) Quais dessas opções se aproxima do valor correto de área basal (//g//)? ==== As Funções no R ==== Já foi visto que ao se digitar o nome de uma função na linha de comando, o R retorna o **código** da função. Veja a diferença de: > ls() para: > ls A maioria das funções precisa de certas **informações** para orientar o seu procedimento, tais informações são chamados de **argumentos**. Os argumentos de qualquer função são detalhadamente explicados nas páginas de ajuda sobre a função. Mas para uma rápida consulta dos argumentos de uma função podemos usar a função '''args''': > args(ls) function (name, pos = -1, envir = as.environment(pos), all.names = FALSE, pattern) NULL > args(q) function (save = "default", status = 0, runLast = TRUE) NULL > args(save.image) function (file = ".RData", version = NULL, ascii = FALSE, compress = !ascii, safe = TRUE) NULL > Algumas funções, entretanto, são primitivas ou internas e seus argumentos não são apresentados. Geralmente, nesses casos os argumentos são bastante óbvios: > args(sin) NULL > sin .Primitive("sin") > Outras funções simplesmente não possuem argumentos: > args(getwd) function () NULL > getwd function () .Internal(getwd()) > Ao observar o resultado da função '''args''', você notará que alguns argumentos são seguindos de uma expressão que se inicia com o sinal de igualdade ('''='''). A expressão após o sinal de igualdade é chamada de **valor default** do argumento. Se o usuário não informar o valor para um dado argumento, a função usa o valor default. Como exemplo veja a função '''save.image''': > args(save.image) function (file = ".RData", version = NULL, ascii = FALSE, compress = !ascii, safe = TRUE) NULL > Se o usuário simplesmente evocar a função '''save.image()''', sem informar o nome do arquivo onde á área de trabalho deve ser gravada, o R gravará as informações num arquivo com nome '''.RData'''. === Exercícios === Quais são os argumentos (e seus valores default) das seguintes funções: * **mean**, * **sd**, * **range**, e * **cumsum**. Quais são os argumentos (e seus valores default) das funções: * **sort**, * **log**, e * **seq**. O que é o argumento **". . ."**? ===== Distribuições Estatísticas: Funções no R ===== Sendo um ambiente para análise de dados, o R dispõe de um grande conjunto de funções para trabalhar com //Distribuições Estatísticas//. Essas funções ajudam não só na análise de dados, como também permitem a //simulação// de dados. ==== Distribuição Normal ==== A distribuição Normal é a distribuição central da teoria estatística. Para gerar uma amostra de observações de uma distribuição normal utilizamos a função '''rnorm''': > args( rnorm ) function (n, mean = 0, sd = 1) NULL > vn1 = rnorm( 1000, mean = 40, sd = 9 ) > mean( vn1 ) [1] 39.47248 > sd( vn1 ) [1] 8.523735 > range( vn1 ) [1] 14.93126 62.11959 > > vn2 = rnorm( 100000, mean = 40, sd = 9 ) > mean( vn2 ) [1] 40.02547 > sd( vn2 ) [1] 9.025218 > range( vn2 ) [1] 3.40680 78.25496 > Se quisermos saber a //probabilidade acumulada// até um certo valor de uma variável com distribuição normal utilizamos a função '''pnorm''': > args(pnorm ) function (q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) NULL > > pnorm( 1.96, mean = 0 , sd = 1 ) [1] 0.9750021 > pnorm( 1.96 ) [1] 0.9750021 > > pnorm( 27, mean = 20, sd = 7 ) [1] 0.8413447 > pnorm( 13, mean = 20, sd = 7 ) [1] 0.1586553 > Se quisermos obter o valor de um //quantil// da distribuição normal utilizamos a função '''qnorm''': > args( qnorm ) function (p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) NULL > qnorm( 0.90 ) [1] 1.281552 > qnorm( 0.30 ) [1] -0.5244005 > > qnorm( 0.90, 20, 7) [1] 28.97086 > qnorm( 0.30, 20, 7) [1] 16.32920 > A função '''dnorm''' fornece a //densidade probabilística// para cada valor de uma variável Normal: > args( dnorm ) function (x, mean = 0, sd = 1, log = FALSE) NULL > x = seq(-4, 4, length=10000) # Seqüência de -4 a 4 com 10.000 valores > > plot(x, dnorm(x)) # Curva da Dist. Normal com média 0 e desvio padrão 1 > points(x, dnorm(x, sd=2)) # Curva da Dist. Normal com média 0 e desvio padrão 2 (adicionada ao gráfico) > === Exercícios === Tomando uma variável que segue a Distribuição Normal, o que acontence com a //amplitude de variação// dos dados à medida que o tamanho da amostra cresce (por exemplo n= 100, 1000, 10000)? Qual o intervalo da Distribuição Normal Padronizada que têm a média no centro e contem 50% das observações? Qual a probabilidade de uma observação da variável Normal Padronizada estar no intervalo [-1.96 , 1.96]? ==== As Funções que Operam em Distribuições Estatísticas ==== O que foi apresentado para Distribuição Normal pode ser generalizado para todas as distribuições que o R trabalha. Há quatro funções para se trabalhar com distribuições estatísticas: * **d**//distrib// - retorna a //densidade probabilística// para um dado valor da variável; * **p**//distrib// - retorna a //probabilidade acumulada// para um dado valor da variável; * **q**//distrib// - retorna o //quantil// para um dado valor de probabilidade acumulada; * **r**//distrib// - retorna //valores// (números aleatórios) gerados a partir da distribuição; No caso da Distribuição Normal: //distrib// = ''norm''. Para outras distribuições temos: ^ Distribuição ^ Nome no R ^ Parâmetros((os argumentos de cada função incluem estes parãmetros, entre outras coisas)) ^ | beta | beta | shape1, shape2, ncp | | binomial | binom | size, prob | | Cauchy | cauchy | location, scale | | qui-quadrado | chisq | df, ncp | | exponential | exp | rate | | F | f | df1, df2, ncp | | gamma | gamma | shape, scale | | geométrica | geom | prob | | hypergeométrica | hyper | m, n, k | | log-normal | lnorm | meanlog, sdlog | | logística | logis | location, scale | | binomial negativa | nbinom | size, prob | | normal | norm | mean, sd | | Poisson | pois | lambda | | t de Student | t | df, ncp | | uniforme | unif | min, max | | Weibull | weibull | shape, scale | | Wilcoxon | wilcox | m, n | === Exercícios === Você realizou um teste //t// de Student bilateral e obteve o valor //t = 2.2// com 19 graus de liberdade. O teste é significativo ao nível de probabilidade de 5%? E se o valor observado fosse //t = 1.9//? Você realizou um teste //F// e obteve o valor //F = 2.2// com 19 graus de liberdade no numerador e 24 graus de liberdade no denominador. O teste é significativo ao nível de probabilidade de 5%? E se o valor observado fosse //F = 2.5//? Gere duas amostras (p.ex.: //x// e //y//) de tamanho 1000 (n=1000) de números da distribuição Uniforme. Faça um gráfico plotando uma amostra contra a outra (''plot(x,y)''). Qual o padrão espacial observado? Você consegue explicá-lo? Gere duas amostras (p.ex.: //xp// e //yp//) de tamanho 10 (n=10) de números da distribuição Uniforme. Gere duas amostras (p.ex.: //xf// e //yf//) de tamanho 1000 (n=1000) de números da distribuição Normal com média zero e desvio padrão 0.2) Faça um gráfico plotando a soma das amostras X (//xp+xf//) contra a soma das amostras Y (//yp+yf//) (''plot(xp+xf,yp+yf)''). Qual o padrão espacial observado? Você consegue explicá-lo? Construa uma seqüência **ordenada** de 1000 números entre 0 e 1: > p = seq(0, 1, length=1000) O vetor '''p''' representa um vetor de probabilidades acumuladas. Gere 1000 números aleatórios da distribuição Normal Padronizada e coloque os números em ordem: > x = sort( rnorm(1000) ) Faça um gráfico dos quantis da distribuição Normal, tomando o vetor '''p''' de probabilidades, contra os valores de '''x''': > plot( qnorm(p), x ) Como é o gráfico resultante? Repita o mesmo processo para as distribuições Exponencial ( '''rexp''' ) e Lognormal ( '''rlnorm''' ). Como são os gráficos resultantes? Por que? Suponha que terremotos ocorram de forma completamente aleatória no decorrer do tempo. Qual a melhor distribuição estatística para representar o **intervalo de tempo** entre dois terremotos **sucessivos**?