Gráficos em R Base:

O mecanismo central de construção de gráficos no R estão contidos nos seguintes pacotes, carregadoa ao se incializar o R:

O foco será usar o sistema gráfico do R Base para criar plots no dispositivo da tela do R.

Salienta-se os pacotes do R Base mencionados acima oferecem um sistema muito poderoso para a criação de gráficos, sendo muito utilizados. Há duas fases na criação de gráficos em R Base:

1- Inicialização de novo gráfico, com abertura do dispositivo gráfico.

2- Anotação e/ou acréscimo e elementos gráfico no dispositivo em uso (aberto na fase 1).

Executar uma função gráfico de high-level [plot (), hist (), boxplot(), etc.] inicializará um dispositivo gráfico (se ainda não estiver aberto) e desenhará um novo gráfico no dispositivo.

O sistema gráfico base tem muitos parâmetros globais que podem ser definidos ou ajustados. Eles ão usados para controlar o comportamento global dos gráficos, tais como, margens, orientação do eixo, etc. Sempre consulte a documentação da função par() que controla os paramêtros gráficos mais relevantes.

Função plot():

Geralmente, a função plot () é utilizada para desenhar diagramas de dispersão e, nesse caso, tem como dados de entrada dois vetores de numéricos de mesmo tamanho, contendo as coordenadas dos pontos que serão plotados. : um para as coordenadas do eixo x e um para as coordenadas do eixo y.

Além disso, plot () é uma função genérica em R, o que significa que seu comportamento pode mudar dependendo de quais tipos de dados de entrada são utilizados como argumento para a função. Não entraremos em detalhes sobre esse comportamento por enquanto.

Comportamento padrão de plot ():

Seja o conjunto de dados iris. Ele contém 150 observações com quatro medidas morfológicas de representando três espécies da flor iris (Iris setosa, versicolor and virginica).

O conjunto de dados é interno ao R e é um data.frame, ou seja, é parecido com uma matriz, mas suas colunas pode ser de diferentes tipos, podendo ser acessadas pelo seu nome.

      # -------  Inspeção do data.frame

  # classe do objeto
class(iris)
## [1] "data.frame"
  # estrutura do objeto iris
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
  # nome das colunas
colnames(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
  # 1os. 6 elementos da variável Petal.Length
head(iris$Petal.Length)
## [1] 1.4 1.4 1.3 1.5 1.4 1.7
  # nome dos níveis do fator Species
levels(iris$Species)
## [1] "setosa"     "versicolor" "virginica"

Vamos construir um diagrama de dispersão simples das variáveis Petal.Length vs. Petal.Width:

plot(iris$Petal.Length, iris$Petal.Width)

O nomes dos eixos x e y foram criados automaticamente a partir dos nomes das variáveis (ou seja, iris$Petal.Length e iris$Petal.Width). Pode ser útil quando você está fazendo gráficos rapidamente, mas exige que você utilize nomes descritivos para suas variáveis e objetos.

Argumentos da função plot():

Muitas das funções gráficos compartilham um conjunto de parãmetros globais. Abaixo, estão apresentados alguns parãmetros chave:

  • pch: símbolo do gráfico (default é círculo).

  • cex: controla a escala dos símbolos e textos plotados, com valor especificado como um múltiplo da fonte básica (default é 1).

  • lty: tipo de linha (default é linha contínua), pode ser pontilhada, tracejada, etc.

  • lwd: largura da linha, especificada como um múltiplo da linha básica (default é 1).

  • col: cor da plotagem, especificado como um númro, string ou código hexadecimal. O comando colors() fornece um vetor dos nomes das cores padrão.

  • bg: cor de fundo dos símbolos (apenas de 21 a 25).

  • xlab: string de caracteres do rótulo do eixo x.

  • ylab: string de caracteres do rótulo do eixo x.

  • main: título do gráfico.

  • sub: sub-título do gráfico.

  • The background color of symbols (only 21 through 25) cex s

  • …: outros parâmetros gráficos.

No caso de nosso diagrama de dispersão, por exemplo, poderíamos acrescentar o título e os rótulos dos eixos x e y.

plot(iris$Petal.Length, iris$Petal.Width, 
     xlab = "Comprimento", ylab = "Largura",
     main = "Dimensões Morfológicas - Pétalas", 
     sub = "Conjunto de dados: iris"
     )

Há um fator no conjunto de dados (Species) e é interessante colorir ou diferenciar os caracteres dos pontos de acordo com os níveis do fator Species. Podemos usar o argumento pch (caractere do ponto) para mudar o tipo dos pontos. Podemos usar, por exemplo, pch = 21, para círculos preenchidos, pch = 22, para quadrados preenchidos e pch = 23, para losangos preenchidos, pch = 24 ou pch = 25 triângulos. Mudar os pontos por tipo de espécie é necessário listar as espécies de cada um das 150 flores do conjunto de dados, atribuindo a os valores 23, 24 e 25 a cada uma dos níveis do fator Species, usando esses 150 valores no argumento pch:

plot(iris$Petal.Length, iris$Petal.Width,
     pch = c(23, 24, 25)[unclass(iris$Species)], 
     xlab = "Comprimento", ylab = "Largura",
     main = "Dimensões Morfológicas - Pétalas", 
     sub = "Conjunto de dados: iris"
     )

O comando unclass(iris$Species) transforma os níveis do fator Species em números (1, 2 e 3). Observe como o comando c(23,24,25)[unclass(iris$Species)] atua.

c(23, 24, 25)[unclass(iris$Species)]
##   [1] 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
##  [26] 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
##  [51] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
##  [76] 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24
## [101] 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25
## [126] 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25 25

Verifique no gráfico abaixo a cofdificação utilizado no argumento pch para ajustar o tipo dos pontos.

plot(1:25, pch=1:25)
text(1:25, 1:25, pos = 3, cex = 0.90)

Podemos controlar o tamanho dos pontos com o argumento cex. Por exemplo, vamos aumentar em 50% o tamanho dos pontos do gráfico para verificar se melhoramos a visualização.

plot(iris$Petal.Length, iris$Petal.Width, 
     cex = 1.50,
     pch = c(23, 24, 25)[unclass(iris$Species)],
     xlab = "Comprimento", ylab = "Largura",
     main = "Dimensões Morfológicas - Pétalas", 
     sub = "Conjunto de dados: iris"
     )

Similarmente ao que foi efetuado com o tipo dos caracteres, podemos gerar uma lista de cores e usá-la para colorir cada ponto conforme sua espécie:

plot(iris$Petal.Length, iris$Petal.Width, 
     pch=21, bg = c("red","green3","blue")[unclass(iris$Species)],
     xlab = "Comprimento", ylab = "Largura",
     main = "Dimensões Morfológicas - Pétalas", 
     sub = "Conjunto de dados: iris"
     )

Todos os pontos são círculos preenchíveis (pch = 21) que serão preenchidos de acordo com a cor correspondente no vetor de cores. O argumento bgcolore o interior do ponto (se o tipo de caractere permitir). As cores usadas são internas do R. Há 657 cores disponíveis. Verifique colors():

  # Cores disponíveis no R Base
head(colors())
## [1] "white"         "aliceblue"     "antiquewhite"  "antiquewhite1"
## [5] "antiquewhite2" "antiquewhite3"

Há pacotes que oferecem uma ampla gama de cores para serem utilizadas nos gráficos, com paletas que podem ser adequadas para colorir gradações de acordo com temperatura, por exemplo. Um pacote bastante utilizado é o RColorBrewer. Ele usa o trabalho de http://colorbrewer2.org para ajudar a escolher esquemas sensíveis de cores para figuras em R. As cores são divididas em três grupos: sequencial, divergente e qualitativo.

As paletas sequenciais são adequadas para dados ordenados crescentes. As etapas gradativas de claridade dominam a aparência desses esquemas de cores, com cores claras para valores baixos dos dados e cores escuras, valores altos.

As paletas divergentes enfatizam de igual maneira os valores críticos do meio e em ambos os extremos do intervalo dos dados. A classe crítica ou quebra no meio da lengenda é enfatizada com cores claras e os extremos são enfatizados com cores escuras, possuindo matizes contrastantes.

As paletas qualitativas não implicam em diferenças de magnitude entre as classes da legenda, sendo que os matizes são usados para criar as principais diferenças visuais entre as classes. Os esquemas qualitativos são mais adequados para representar dados nominais ou categóricos.

As paletas do RColorBrewer podem ser mostradas com o comando display.brewer.all(),

Esse pacote oferece também cores adequadas a daltônicos (argumento colorblindFriendly = TRUE). Verifique as paletas disponíveis:

display.brewer.all(colorblindFriendly = TRUE)

Verifique o uso da paleta Set2 com os comandos abaixo:

plot(iris$Petal.Length, iris$Petal.Width, pch = 21, 
     bg = brewer.pal(n = 3, name = "Set2")[unclass(iris$Species)],
     xlab = "Comprimento", ylab = "Largura",
     main = "Dimensões Morfológicas - Pétalas", 
     sub = "Conjunto de dados: iris"
     )

Poderíamos ter um problema de overplotting, ou seja, a ocorr~encia de uma concentração de pontos em alguma região do diagrama de dispersão, embaraçando a visualização do padrão dos pontos. Vamos gerar duas amostras, xe y, formando 10.000 pontos:

set.seed(666)
  # amostra normal com              
  # média zero e desvio-padrão 2
x <- rnorm(10000, sd = 2)

  # amostra normal padrão
y <- rnorm(10000)

O plot default com pch = 19 será:

plot(x, y, pch = 19)

Percebe-se que as escalas dos eixos não concordam. Vamos então ajustar para que a escala dos eixos x e y coincidam:

plot(x, y, pch = 19, 
     xlim = c(-6, 6), ylim = c(-6, 6)
     )

Como era de se esperar, a massa de pontos tem o formato elípitico, mas continua um borrão contínua sem dar ideia das diferenças em conentração dos dados. Espera-se que os pontos estejam concentrado em torno de zero mas que pelo menos diminuam para os pontos que estejam além da região \([-2, 2] \times [-1, 1]\), Uma alternativa seria diminuir o símbolo dos pontos, para, por exemplo pch = ".":

plot(x, y, pch = ".", 
     xlim = c(-6, 6), ylim = c(-6, 6)
     )

As concentrações são mais visíveis. Podemos também alterar a transparência dos pontos com a função rgb().

plot(x, y, pch = 19, col = rgb(0, 0, 0, 0.15), 
     xlim = c(-6, 6), ylim = c(-6, 6)
     )

Percebe-se claramente uma diminuição de pontos ao nos afastarmos além de dois desvios-padrão em cada eixo. Os formatos das concentração aparentemente continuam elípticos. Como traçar três elipses, centradas na origem, e que, conforme a distribuições teóricas, contivessem, respectivamente, 68%, 95% e 99% dos dados?

Para isso, necessitamos de funções auxiliares, para desenhar em cima do gráfico construído, como se fossem camadas em uma tela.

Funções de anotação em gráficos (low-level):

A execução da função plot() constrói um plot na tela do dispositivo gráfico (abrirá o dispositivo gráfico se ele não estiver aberto). Depois disso, as funções de anotação (funções low-level) podem ser executadas para adicionar anotações no gráfico que estiver construído no dispositivo ativo. Abaixo, uma pequena lista de algumas dessas funções de anotação:

Vamos acompanhar alguns exemplos de construção de um gráfico base e adição de algumas anotações.

Funções points() e title():

São funções secundárias e apenas anotam em gráfico aberto (se não houver gráfico aberto elas não funcionam!) O objetivo aqui e construir um gráfico das dimensões morfológicas das pétalas e sépalas do conjunto de dados iris. O procedimento será criar o gráfico com as largura e o comprimento das pétalas, utilizando pch = 21 para os pontos e posteriormente acrescentar os pontos relativos à largura e comprimento das sépalas, com outro tipo de caractere (pch = 3). Vamos manter a distinção das espécies por cores. Primeiro temos de determinar os intervalos de variação da largura e do comprimento tanto de pétalas, quanto de sépalas. Para tanto usamos a função summary().

summary(iris[-5])
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500

Os comprimentos variam de 1,000 (Petal.Length) a 7,900 (Sepal.Length) e as larguras, de 0,100 (Petal.Width) a 4,400 (Sepal.Width). Vamos usar a função plot() para construir o gráfico principal com as medidas das pétalas, sem nomear os eixos e ajustando a escala do eixo x para o intervalo (1, 8) e do eixo y, para (0, 4.5). Posteriormente, será usada a função points() para plotar os pontos das dimensões das sépalas. Por fim, acrescentaremos os títulos e rótulos dos eixos com a função title().

   # cores das espécies
  cores <- c("coral1", "cyan4", "darkgoldenrod2")
  
          # -------  Gráfico primário  ----------

     # plot principal - medidas das pétalas
  plot(Petal.Width ~ Petal.Length, data = iris, 
       xlab = "", ylab = "",
       pch = 21, bg = cores[unclass(Species)],
       xlim = c(1, 8), ylim = c(0, 4.5) 
       )

       # -------  Anotações no gráfico primário  ----------
  
    # dimensões morfológicas das sépalas
  points(Sepal.Width ~ Sepal.Length, data = iris, pch = 3, cex = 1.2,
         col = cores[unclass(Species)])
  
  # títulos centralizados
  title(main = "Dimensões Morfológicas", 
        xlab = "Comprimento, em cm", ylab = "Largura, em cm",)
  
  # sub-título alinhado à esquerda
  title(sub = "Conjunto de dados: iris", cex.sub = 0.8, font.sub = 2, 
        adj = 0)

As escalas dos eixos x e y foram impostas pelos argumentos xlim = c(1, 8) e ylim = c(0, 4.5). O subtítulo teve sua fonte modificada para negrito (font.sub = 2), diminuída 80% do tamanho default (cex.sub = 0.8) e alinhado à esquerda (adj = 0). O gráfico possui agora todas as informações do conjunto de dados, pois as 5 variáveis estão representadas nele. Necessitamos de legenda para saber as cores que representam as espécies e os caracteres que representam as pétalas e sépalas.

Legenda:

Função usada nesta seção: legend():

Também é uma função secundária e atua apenas sobre o gráfico que estiver ativo. Temos de especificar suas coordenadas ou a posição em que queremos possicionar a legenda no gráfico. No gráfico que construímos há espaço para legenda no canto inferior direito e superior esquerdo. Vamos construir duas legendas: uma diferenciando as cores para as espécies (horizontal, no canto inferior direito) e outra, a morfologia da flor (vertical, no canto superior direito). Ambas terão título. No R Markdown, temos de construir o gráfico novamente, para a função legend() ter um dispositivo para anotas.

Diagrama de dispersão com reta de regressão:

Funções usadas nesta seção: abline(), polygon(), text() e expression():

Um dos interesses em Estatística é modelar a relação entre uma variável resposta (ou variável dependente) e uma ou mais variáveis explicativas (ou variáveis independentes), no sentido de estabelecer algum tipo de função matemática do tipo \(y = f(x_1, x_2, \dots, x_p)\), onde \(y\) é a variável resposta e \(x_1, x_2, \dots, x_p\), as variáveis explicativas. A regressão linear é uma abordagem para modelo essa relaçao por meio de uma função \(f\) linear. Caso haja uma única variável explicativa essa relação é denominada regressão linear simples e, para mais de uma variável explicativa, o procedimento é denominado regressão linear múltipla. Vamos apresentar rapidamente o modelo de regressão linear simples, no qual uma variável aleatória resposta contínua \(Y\) está relacionada linearmente a uma variável explicativa \(x\) e com uma variável aleatória \(\epsilon\), \(\epsilon \sim \mathcal{Normal}(0, \sigma^2)\). Essa variável aleatória \(\epsilon\) pode ser considerado um ruído do modelo, representando outras variáveis explicativas que porventura não foram incluídas no modelo (o ideal seria ter modelo com valores baixos de ruído, ou seja, \(\sigma^2\), baixos). O modelo probabilístico é então: \[ Y = \beta_0 + \beta_1 X + \epsilon, ~ \epsilon \sim\mathcal{Normal}(0,\sigma^2), \]

onde \(\beta_0\) é denominado intercepto, pois é o ponto de intersecção com o eixo y e \(\beta_1\) é a inclinação da reta. Dentre outras, o R oferece a função lm() que é bastante poderosa e, dentre outras funcionalidades, estima os valores doscoeficientes da regressão (intercepto e da inclinação), a partir de uma amostra de valores observados das variáveis resposta (\(y\)) e explicativa (\(x\)). vamos gerar uma amostra com 25 pontos para visualizar o funcionamento da função (lm()). Serão simulados 25 pontos do modelo \(Y = 1 + 1.5 x + \epsilon, ~ \epsilon \sim \mathcal{N}(0,1)\) e usada a semente 666, para replicação posterior dos resultados:

   # geração da amostra
  set.seed(666)
  ruido <- rnorm(25)
  explicativa <- runif(25, min = 0, max = 3)
  resposta <- 1 + 1.5 * explicativa + ruido

A função lm() estima os valores do intercepto e da inclinação dados 25 pontos de uma amostra com valores de \(x\) (explicativa) e \(y\) (resposta). Já a utilizamos anteriormente e vamos criar um objeto (obj.reg).

   # ajuste do modelo de regressão linear
  obj.reg <- lm(resposta ~ explicativa)
  
    # coeficentes do modelo
  coeficientes <- coef(obj.reg)
  coeficientes
## (Intercept) explicativa 
##   0.6905156   1.6937340

Os coeficientes estimados (obtidos pela função coef()) são intercepto (\(\hat{\beta}_0\)) com o valor de 0.691 e inclinação (\(\hat{\beta}_1\)), com valor de 1.694. Essas estimativas irão variar para cada amostra simulada de 25 valores. Perceba também que esses valores estimam os valores verdadeiros (que usamos na simulação), ou sejam, \(\beta_0 = 1\) e \(\beta_1 = 1,5\), ou seja, os erros amostrais serão variáveis também. O diagrama de dispersão dos pontos e a reta de regressão estimada estão plotadas abaixo:

      # -----  Gráficos dos pontos e da reta de reressão  ------
  
  plot(x = explicativa, y = resposta, pch = 21, bg = "cyan4", 
       xlim = c(0, 3.1), main = "Reta de Regressão")
  
    # reta de regressão
  abline(obj.reg, lty = 2, col = "gray48", lwd = 2)
  
  
    # expressão da reta de regressão estima
  text(0, 6, expression(y == 0.691 + 1.694 ~x), pos = 4)
  
    # vértices de triângulo com altura beta.1 e base unitária
  vert.x <- c(1, 2, 2)
  vert.y <- coeficientes[1] + coeficientes[2] * c(1, 1, 2)
 
     # triângulo com altura beta.1 e base unitária
  polygon(vert.x, vert.y, col = 'gray88', border = NA)
    # anotação do valor da inclinação
  text(2, 3.28, expression(hat(beta)[1] == 1.694), cex = 0.85, pos = 4)
  
  
    # ponto do intercepto
  points(0, coeficientes[1], pch = 21, bg = "gray48")
    # anotação do valor do intercepto
  text(0.176, 2.75, expression(hat(beta)[0] == 0.691), cex = 0.85)
  arrows(0.2, 2.4, 0, coeficientes[1], length = 0.15, angle = 15)

Verifique com atenção a documentação de todas as funções de low-level que utilizamos para elaborarmos anotações no gráfico ativo: abline(), text(), expression(), arrows(), polygon() e points() (já utilizada em gráfico anterior). O triângulo no gráfico foi desenhada com a função polygon(), que tem os vétices do polígono como valores de entrada. Ela pode ser usada, por exemplo, para desenhar a área entre duas curvas (descubra como!). Por sua vez função text() presta-se para anotar textos no área do gráfico, tendo como valores de entrada o texto a ser anotado (uma string) e as coordenadas em que essa string será posicionada no gráfico. Possui vários argumentos de ajuste do tamanho da fonte (cex), posição em torno da coordenada (pos), justificação em torno do ponto (adj), etc. Confira sua documentação, pois, ela apresenta muitos bons exemplos. Junto com ela foi usada a função expression() que cria strings com anotações matemáticas: símbolos de operadores, letras gregas, etc. Essa função permite utilizar uma ampla gama de anotações matemáticas e tem muita similaridade com o LaTeX. Confira um pouco mais sobre a função expression() aqui: http://magnusmetz.github.io/2013/04/mathematical-annotation-in-r. Por sua vez, a função abline() trabalha diretamente com o objeto obj.reg criado com o ajuste com a função lm(), o que facilita bastante seu uso em estudos de regressão. Ela desenha qualquer reta em gráficos ativos, bastando informar seus valores do intercepto e inclinação. Perceba o que acontece se desenharmos a reta verdadeira (\(y = 1 + 1,5x\)). Lembre-se que simulamos com esses parâmetros.

# -----  Gráficos dos pontos e da reta de reressão  ------
  
  plot(x = explicativa, y = resposta, pch = 21, bg = "cyan4", 
       xlim = c(0, 3.1), main = "Reta de Regressão")
  
    # reta de regressão
  abline(obj.reg, lty = 2, col = "gray48", lwd = 2)

  # reta de regressão verdadeira (populacional)
  abline(1, 1.5, col = "blue", lty = 3)
  
    # centróide da massa de pontos
  points(x = mean(explicativa), y = mean(resposta), pch = 21,
        bg = "blue")
  
    # coordenadas pelo centróide
  abline(h = mean(resposta), v = mean(explicativa), lty = 3, 
         col = "grey38")
    # anotação
  text(1.7, 1.3, "centróide", cex = 0.85)
  arrows(1.7, 1.5, mean(explicativa), mean(resposta), length = 0.10, 
         angle = 15)

Tivemos de criar o gráfico novamente, pois o comando abline(1, 1.5, col = "blue", lty = 3) necessita de gráfico ativo para atuar. Há uma pequena diferença entre a reta verdadeira (populacional) e a reta de regerssão estimada. Perceba que as retas se cruzam no centróide dos pontos (médias das duas variáveis).

A pergunta que fica é como a reta de regressão é estimada? Qual critério é utilizado para passar uma reta pela massa de pontos? A reta de regressão foi estimada pelo método dos mínimos quadrados ordinários. São consideradas os desvios quadráticos entre o valor ajustado (na reta) e cada valor observado da resposta. Observe no diagrama de dispersão de nossa amostra de 25 pontos.

# ----------  Desvios quadráticos  -----
  
  # -----  Gráficos dos pontos e da reta de reressão  ------
  
  plot(x = explicativa, y = resposta, pch = 21, bg = "cyan4", 
       xlim = c(0, 3.1), main = "Reta de Regressão")
  
  # reta de regressão
  abline(obj.reg, lty = 2, col = "gray48", lwd = 2)
  

  
    # --------- função para calcular ponto na reta ------
  reta <- function(x) coeficientes[1] + x * coeficientes[2]
  
    # ordenadas dos pontos sob a reta
  ajustado <- reta(explicativa)
   
    # ordenada dos pontos observados
  observado <- resposta
  
    # desvios da reta
  segments(x0 = explicativa, y0 = ajustado, 
           x1 = explicativa, y1 = resposta, 
           col = "gray58"
  )

A função reta() foi construída para calcular o valor da ordenada sob a reta (chamamos de valor ajustado) ou seja \(y= \hat{\beta}_0 + \hat{\beta}_1 x)\). Calculamos então um valor ajustado apar cada valor observado e desenhos os segmentos entre o ponto observado e a reta usando a função segments(). Essa função tem como valores de entrada os vetores das coordenadas dos extremos de cada segmento de reta [\((x_0, y_0)\) e \((x_1, y_1)\)]. Podem ser utilizados os mesmos argumentos utilizados com linhas (lty, lwd, col, etc.),

Diagrama de dispersão com linha suavizada:

Funções usadas nesta seção: lines() e lowess():

A função lowess() executa uma estimação de regressão não paarmétrica ponderada localmente (locally weighted scatter plot smoothing - LOWESS. Ela é denominada de estimação não paramétrica por não fornecer uma expressão matemática (com parâmetros ou coeficientes), já que a estimação de cada ponto é calculada por um algoritmo apropriado. Esses pontos estimados não parametricamente podem ser desenhados no gráfico pela função lines() que conecta todos os pontos com o tipo especificado de linha, formando assim uma linha suavizada. Vamos verificar o resultado em nosso último exemplo, em que é possível comparar também com a reta de regressão.

               # ------------------------------------
                #     Regressão suavizada - LOWESS
  
  
       # -----  Diagrama de dispersão dos pontos  ------
  
  plot(x = explicativa, y = resposta, pch = 21, bg = "cyan4", 
       xlim = c(0, 3.1), main = "Reta de Regressão")
  
  
    # suavização por LOWESS
  lines(lowess(explicativa, resposta), lty = 3, lwd = 3, col = "brown2")
  
    # reta de regressão
  abline(obj.reg, lty = 2, col = "gray48", lwd = 1)

Verifica-se que a reta de regressão e a linha suavizada são praticamente coincidente, com exceção do centro da massa de dados que, provavelmente por ter um menor quantidade de pontos, fez a linha suavizada infletir na direção dos dados existentes em sua vizinhança. Essas linhas suavizadas costumam ser muito úteis como auxilixar na exploração de relações entre variáveis, ajudando a perceber estruturas possíveis em modelos entre a variável resposta e a explicativa. Vamos observar na relação entre as variáveis price e carat do conjunto de dados diamonds, pertencente ao pacote ggplot2.

 library(ggplot2)
  
    # dimensões do conjunto de dados
  dim(diamonds)
## [1] 53940    10
    # diagrama de dispersão dos dados
  plot(price ~ carat, data = diamonds)
  
  # suavização por LOWESS
  lines(lowess(x = diamonds$carat, y = diamonds$price), lty = 3, 
        lwd = 3, col = "white")

O conjunto de dados diamonds é bastante grande (53940 observações) e o diagrama de dispersão apresenta o problema de overplotting devida a quantidade de pontos. A linha suavizada (em cor branca para poder ser visualizada) indica uma estrutura não linear que merece ser investigada (como é esse comportamento nos diversos estratos do conjunto de dados?). Uma outra possibilidade de visualização é tornar os pontos mais transparentes para permitir um gradiente cromático que auxilie na percepção das concentrações dos pontos.

 # diagrama de dispersão dos dados com transparência
  plot(price ~ carat, data = diamonds, pch = 19, 
       col = rgb(0, 0, 0, 0.15))
  
  # suavização por LOWESS
  lines(lowess(x = diamonds$carat, y = diamonds$price), lty = 3, 
        lwd = 3, col = "white")

Ao tornarmos os pontos mais transparentes, percebemos truncamentos na distribuição dos pontos que podem indicar outras relações, com variáveis categóricas do conjunto de dados ou mesmo com medidas discretas das medidas de carat. Mas essa pergunta, embora bastante interessante, foge um pouco de nosso objetivo atial.

Anotação matemática

Funções usadas nesta seção: curve(), axis(), expression(), substitute() e bquote():

Se, em alguma das funções de anotação de texto do R Base [text(), mtext(), axis(), legend(), title()], a entrada de texto estiver na forma de expressão, então esse argumento é interpretado como uma expressão matemática e a saída estará formatada de acordo com regras semelhantes ao TeX. Expressões não podem ser usadas nos rótulos dos eixos em plots construídos pela função persp(). Note que expressão matemática passa a fazer parte da string que será escrita no plot. Uma maneira de concatenar texto é usar a função expression com a forma "texto" ~ expressão como seu argumento ou, paste("texto" ~ expressão). No exemplo abaixo desenhamos as curvas da função seno e cosseno com a função curve().

    # gráfico da função seno
  curve(sin, from = -3*pi/2, to = 3*pi/2,
        xaxt = "n", ylab = "", xlab = "")

    # gráfico da função cosseno acrescentado ao gráfico ativo 
 curve(cos, from = -3*pi/2, to = 3*pi/2, lty = 2, col = "blue",
        add = TRUE)
  
   # eixo x com anotação matemática no eixo x
  axis(1, at = c(-pi, -pi/2, 0, pi/2, pi), 
       labels = expression(-pi, -pi/2, 0, pi/2, pi))
  
    # titulo e nomes dos eixos com anotação matemática
  title(main = expression(plain(sin)*phi~"e"~plain(cos)*phi), 
        col.main = "blue",
        xlab = expression("Ângulo"~ phi),
          # só funciona a primeira anotação 
        ylab = expression("sin"*phi, "cos"*phi)
        )

O gráfico da função seno foi facilmente desenhado com a função curve(), bastanto o nome da função e o intervalo sobre o qual ela será plotada. A curva da função cosseno foi acrescentada ao gráfico ativo com a mesma função, mas agora com o argumento add = TRUE. Se esse argumento estivesse em seu valor default, seria aberto um novo dispositivo gráfico agora com a função cosseno. No comando para abrir o dispositivo com o gráfico da função seno foi usado o argumento xaxt = "n" para nada ser anotado no eixo x. O objetivo foi para anotar o eixo x com anotação mateática, no caso o número \(\pi\). O título e o eixo x foram escritos com strings de texto e de anotação matemática, concantenadas pelo operador ~. Uma questão importante é utilizar o valor de uma variável numérica em uma anotação de expressão matemática. Para combinar expressão matemática e texto usamos os a função substitue() ou a função bquote().

      # ----  Combinando anotação matemática e variáveis numéricas ---
  
  
    # gráfico primário em branco para as anotações
  plot(1:10, type = "n", xlab = "", ylab = "", 
       main = "Anotação Matemática e Variáveis")
 
    # anotando as coordenadas de 4 pontos no gráfico (abcissas pares)
  for(i in (1:4)*2)
    text(i, i+1, 
       substitute(list("x"[x],"y"[y]) == 
                  group("(",list(x, y),"]"), list(x = i, y = i + 1)))
  
    # Símbolo matemático e valor de variável fora da região do gráfico
  media <- 1.23 
  mtext(bquote(hat(mu) == .(media)), line= .25)

A função substitute() está utilizando os valores de x e y de seu segundo argumento [list(x = i, y = i + 1)], os quais estão mudando a cada loop em i e está substituindo em todos os x’s e y’s de seu primeiro argumento [list("x"[x],"y"[y]) === group("(",list(x, y),"]")]. Note que tivemos de utilizar "x" e "y" para que eles não fossem substituídos também. Imagine que este valor pode variar em cada gráfico que você constrói com aquele código, você necessitaria de uma maneira para introduzi-lo no gráfico automaticamente (nos exemplos anteriores, inserimos os valores manualmente). A função bquote() concatena a anotação matemática referente a \(\hat{\mu}\) com o valor atual de media, o qual está estabelecido no código. Essa inserção deve se dar na forma .(media). Lembre-se que, caso necessite anotar uma igualdade, você necessita usar ==. O mais importante é que você vá se familiarizando com os comandos dos vários símbolos matemáticos, como, por exemplo, a finalidade do comando group, utilizado no código acima. Verifique os comandos disponíveis em R aqui: https://stat.ethz.ch/R-manual/R-devel/library/grDevices/html/plotmath.html.

Veja abaixo algumas anotações de fórmulas e seus comandos correspondentes;

        # ----  Algumas expressões matemáticas ---
  
     # criando quadro branco
  plot(1:10, type = "n", xlab = "", ylab = "", yaxt = "n", xaxt = "n",
       main = "Algumas Expressões Matemáticas")
    # expressão e comando para anotação de vetor de coeficientes
  text(2, 9, expression(hat(beta) == (X^t * X)^{-1} * X^t * y))
  text(4, 9, "expression(hat(beta) == (X^t * X)^{-1} * X^t * y)",
       cex = .8, pos = 4)
  
    # expressão e comando para anotação da média amstral
  text(2, 6.5, expression(bar(x) == sum(frac(x[i], n), i==1, n)))
  text(4, 6.5, "expression(bar(x) == sum(frac(x[i], n), i==1, n))",
       cex = 0.8, pos = 4)
  
    # expressão da função de densidade de probabilidade da normal
  text(2, 3, expression(paste(frac(1, sigma*sqrt(2*pi)), " ",
                              plain(e)^{frac(-(x-mu)^2, 2*sigma^2)})),
       cex = 1.2)
  text(4, 3, "expression(paste(frac(1, sigma*sqrt(2*pi)), \" \",
                   plain(e)^{frac(-(x-mu)^2, 2*sigma^2)}))",
       cex = 0.8, pos = 4)
  
    # cabeçalho da tabela
  mtext(c("Expressão", "Comando"), at = c(2.5, 7), line = 0.25, side = 3)

    # linhas de grade
  abline(v = 4, h = c(5.5, 8), lty = 3)

No comando da expressão da função de densidade de probabilidade da distribuição normal foi usada a função paste() para formar uma string de texto e de expressão matemática, de maneira a possibilitar sua anotação no gráfico. Temos também a função mtext() anota nas margens da figura especificada pelo argumento side (verifique os valores possíveis na documentação dessa função). seu argumento at posiciona o texto na coordenada dada pelo eixo paralelo à margem especificado e o argumento line especifica na direção perpendicular ao eixo da margem, em termos de múltiplos de uma linha default. A seguir alguns códigos úteis para se obter símbolos especiais:

  plot.new(); 
  plot.window(c(0,4), c(15,1))
  
  pos.d <- 1      # descrição
  pos.c <- 2.5    # comando
  pos.s <- 3.3    # símbolo
  
  title(main = "Símbolos Especiais")
  mtext(c("Descrição", "Comando", "Símbolo"), at = c(1, 2.2, 3), 
        line = 0.25, side = 3, adj = 0)

    # quantificação universal
  text(pos.d, 1, "universal", adj = 0)
  text(pos.c, 1,  "\\042")
  text(pos.s, 1, expression(symbol("\042")))
  
    # quantificação existencial
  text(pos.d, 2, "existencal", adj = 0)
  text(pos.c, 2,  "\\044")
  text(pos.s, 2, expression(symbol("\044")))
  
    # pertinência a conjunto
  text(pos.d, 3, "pertence a", adj = 0)
  text(pos.c, 3,  "\\047")
  text(pos.s, 3, expression(symbol("\047")))
  
    # portanto
  text(pos.d, 4, "portanto", adj = 0)
  text(pos.c, 4,  "\\134")
  text(pos.s, 4, expression(symbol("\134")))
  
    # perpendicular
  text(pos.d, 5, "perpendicular", adj = 0)
  text(pos.c, 5,  "\\136")
  text(pos.s, 5, expression(symbol("\136")))
  
    # círculo e vezes
  text(pos.d, 6, "circulo vezes", adj = 0)
  text(pos.c, 6,  "\\304")
  text(pos.s, 6, expression(symbol("\304")))
  
    # círculo e mais
  text(pos.d, 7, "círculo mais", adj = 0)
  text(pos.c, 7,  "\\305")
  text(pos.s, 7, expression(symbol("\305")))
  
    # conjunto vazio
  text(pos.d, 8, "conjunto vazio", adj = 0)
  text(pos.c, 8,  "\\306")
  text(pos.s, 8, expression(symbol("\306")))
  
    # ângulo
  text(pos.d, 9, "ângulo", adj = 0)
  text(pos.c, 9,  "\\320")
  text(pos.s, 9, expression(symbol("\320")))
  
    # ângulo esquerdo
  text(pos.d, 10, "ângulo esquerdo", adj = 0)
  text(pos.c, 10,  "\\341")
  text(pos.s, 10, expression(symbol("\341")))
  
    # ângulo direito
  text(pos.d, 11, "ângulo direito", adj = 0)
  text(pos.c, 11,  "\\361")
  text(pos.s, 11, expression(symbol("\361")))

Deseja-se: (i) simular 1.000 amostras exponencias padrão (média 1) de tamanho \(n=25\); (ii) calcular a média de cada amostra; (iii) fazer um histograma dessas 1.000 médias amostrais; (iv) colocar no mesmo gráfico o gráfico da função de densidade de probabilidade de uma normal com média 1 e desvio padrão \(1/\sqrt{n}\) e, por fim, (iv) anotar no gráfico o tamanho da amostra e o valor da média global (\(\bar{\bar{x}}\)), ou seja, a média de todas as médias amostrais. Verifique o código abaixo:

    # tamanho amostral
  n <- 25
  
  # geração de 1000 medias de amostras de tamanho n
  medias <- replicate(1000, mean(rexp(n)))
  
    # histograma das médias amostrais
  hist(medias, freq = F, #xlim = c(0.6, 1.4), 
       main = "Histograma de Médias Amostrais",
       xlab = "Médias amostrais", ylab = "Densidade")
  
    # curva da densidade da normal
  curve(dnorm(x, mean = 1, sd = 1/sqrt(n)), add = T)
  
    # cálculo da média global
  x.bar <- mean(medias)
 
    # tamanho das amostra - abaixo do título do gráfico
  mtext(paste("Amostras exponenciais (n =", n, ")"), cex = 0.9, 
        line = 0.25, font = 3)
  
    # valor da média global de todas as médias
  text(1.3, 1.5, bquote(bar(bar(x)) == .(x.bar)), 
       adj = 0, cex = 0.8)

Como estamos comparando histograma com densidade, temos de utilizar histograma de densidade (argumento freq = FALSE). verifique a documentação de algumas novas funções: dnorm(), que calcula a densidade de um normal, no caso com parâmetros mean = 1 e sd = 1/sqrt(n) (por que??); rexp(25) que sorteia uma amostra exponencialmente distribuída de tamanho 25 e a função replicate() que repete a quantidade de vezes especificado (1000) o sorteio desejado (rexp(25)). Cada vez que executarmos esses comandos o histograma e média global serão modificados (e as escalas dos eixos também). Vamos tentar automatizar mais o gráfico para podermos utilizar os comandos na geração de amostras de outros tamanhos, digamos \(n=100\). O melhor é montar uma função para poder repetir a operação mais facilmente:

amostra.exp <- function(n){
   
   # desvio padrão da média
   dp.xbar <- 1/sqrt(n)
   # limite eixo y - máxima densidade
   y.sup <- 1.1*dnorm(1, mean = 1, sd = dp.xbar)
   
   # posição média global - 
   x.glb <- 1 + 1.1 * dp.xbar
   y.glb <- 2/3 * y.sup
   
   
   # geração de 1000 medias de amostras de tamanho n
   medias <- replicate(1000, mean(rexp(n)))
   
   # histograma das médias amostrais
   hist(medias, freq = F, #xlim = c(0.6, 1.4), 
        ylim = c(0, y.sup),
        main = "Histograma de Médias Amostrais",
        xlab = "Médias amostrais", ylab = "Densidade")
   
   # curva da densidade da normal
   curve(dnorm(x, mean = 1, sd = dp.xbar), add = T)
   
   # cálculo da média global
   x.bar <- mean(medias)
   
   
   # tamanho das amostra - abaixo do título do gráfico
   mtext(paste("Amostras exponenciais (n =", n, ")"), cex = 0.9, 
         line = 0.25, font = 3)
   
   # valor da média global de todas as médias
   text(x.glb, y.glb, bquote(bar(bar(x)) == .(x.bar)), 
        adj = 0, cex = 0.8)
 } 
  # amostras de tamanho 25
amostra.exp(25)

  # amostras de tamanho 100
amostra.exp(100)

Tentem explorar em detalhes as funções gráficas do R Base. O importante é vocês se familiarizarem com a estrutura de objetos e funções do R. Não procurem gráfico prontos, aprendam a montá-los. O caminho pode parecer árduo no início, mas abrirá um universo de possibilidades gráficas! Há muito mais funções para apresentar.

Continua em breve!!