1 Introdução

Olá, prezado aluno! Este documento apresenta os blocos de código R necessários para realizar as análises do experimento de geometria fractal com as bolinhas de papel. Antes de mais nada, precisamos garantir que temos os dados no formato correto para importação. Para o exercício, precisaremos construir duas tabelas, descritas a seguir.

1.1 Tabela de tratamentos

A tabela de tratamentos deve conter três colunas, a primeira (Tamanho) em formato texto, que serve de identificador para cada um dos tratamentos (dimensões de papel), a segunda em formato numérico, com os valores referentes à área (m²) dos papeis, e a terceira em formato numérico, com os valores de massa (g) de cada bolinha de papel, conforme o exemplo:

Tamanho Area Massa
256 (A4) 0.0624 4.6800
128 (A5) 0.0312 2.3400
64 (A6) 0.0156 1.1700
32 (A7) 0.0078 0.5850
16 (A8) 0.0039 0.2920
8 (A9) 0.0019 0.1420
4 (A10) 0.0010 0.0750
2 (A11) 0.0005 0.0375
1 (A12) 0.0002 0.0150

Na tabela apresentada a massa de cada bolinha foi calculada pela multiplicação da área (m²) e da gramatura do papel (g/m²), que neste caso era de 75 (g/m²).

1.2 Tabela de medições

A segunda tabela deve conter 10 colunas, sendo a primeira referente a um identificador da medição, e as demais referentes a cada uma das medições realizadas para cada tratamento. O número de medições de diâmetro para cada tamanho de bolinha de papel deve ser o número de linhas da tabela (sem contar o cabeçalho), ou seja, dez linhas referentes às dez medições de cada tratamento.

Medicao 256 (A4) 128 (A5) 64 (A6) 32 (A7) 16 (A8) 8 (A9) 4 (A10) 2 (A11) 1 (A12)
Med1 3.88 3.06 2.60 1.52 1.32 0.96 0.64 0.52 0.41
Med2 3.87 2.96 2.20 1.45 1.21 1.07 0.75 0.44 0.45
Med3 3.89 2.89 2.68 1.59 1.13 0.92 0.68 0.62 0.46
Med4 3.77 2.94 2.44 1.51 1.55 0.84 0.72 0.60 0.51
Med5 3.90 2.82 2.42 1.46 1.21 1.01 0.74 0.54 0.49
Med6 3.78 2.77 2.16 1.42 1.16 1.00 0.83 0.60 0.44
Med7 3.89 2.94 2.56 1.55 1.28 0.80 0.84 0.69 0.44
Med8 3.88 2.77 2.29 1.82 1.19 0.99 0.71 0.61 0.43
Med9 3.94 2.70 2.31 1.59 1.26 0.88 0.67 0.68 0.45
Med10 3.98 2.83 2.53 1.52 1.19 0.80 0.70 0.61 0.43

2 Entrada de dados no R

Existem diferentes maneiras de importar dados para o ambiente do R, sendo os mais comuns a importação de arquivos .csv ou .txt, e o input manual.
Para este experimento, vamos construir as tabelas manualmente, utilizando as funções data.frame e c. Lembre-se que o marcador decimal no R é o ponto, e que a vírgula é o separador de argumentos das funções. Começaremos com a tabela de tratamentos, em que devemos informar o identificador de cada tratamento (tamanho do papel), e seus respectivos valores de área (m²) e massa..

# criar um dataframe com os dados dos tratamentos
dados_tratamentos <- data.frame(
  Tamanho = c('256 (A4)','128 (A5)','64 (A6)','32 (A7)','16 (A8)','8 (A9)','4 (A10)','2 (A11)','1 (A12)'),
  
  Area = c(0.0624, 0.0312, 0.0156, 0.0078, 0.0039, 0.0019, 0.001, 0.0005, 0.0002),
  
  Massa = c(4.68, 2.34, 1.17, 0.585, 0.292, 0.142, 0.075, 0.0375, 0.015))

# retornar o dataframe criado
dados_tratamentos
##    Tamanho   Area  Massa
## 1 256 (A4) 0.0624 4.6800
## 2 128 (A5) 0.0312 2.3400
## 3  64 (A6) 0.0156 1.1700
## 4  32 (A7) 0.0078 0.5850
## 5  16 (A8) 0.0039 0.2920
## 6   8 (A9) 0.0019 0.1420
## 7  4 (A10) 0.0010 0.0750
## 8  2 (A11) 0.0005 0.0375
## 9  1 (A12) 0.0002 0.0150

O próximo passo será a criação da tabela com as medições dos diâmetros das bolinhas de papel. Observe que ao declarar o nome da coluna, se há um ou mais espaços no rótulo, devemos colocá-lo entre crases. No final, podemos declarar o argumento check.names = FALSE, para que o R salve os nomes das colunas exatamente como digitamos. Do contrário, todos os caracteres especiais seriam tranformados em pontos.

# criar um dataframe com os dados coletados
dados_medicoes <- data.frame(
  Medicao = c('Med1', 'Med2', 'Med3', 'Med4', 'Med5', 'Med6', 'Med7', 'Med8', 'Med9', 'Med10'),
  `256 (A4)` = c(3.88, 3.87, 3.89, 3.77, 3.9, 3.78, 3.89, 3.88, 3.94, 3.98),
  `128 (A5)` = c(3.06, 2.96, 2.89, 2.94, 2.82, 2.77, 2.94, 2.77, 2.7, 2.83),
  `64 (A6)` = c(2.6, 2.2, 2.68, 2.44, 2.42, 2.16, 2.56, 2.29, 2.31, 2.53),
  `32 (A7)` = c(1.52, 1.45, 1.59, 1.51, 1.46, 1.42, 1.55, 1.82, 1.59, 1.52),
  `16 (A8)` = c(1.32, 1.21, 1.13, 1.55, 1.21, 1.16, 1.28, 1.19, 1.26, 1.19),
  `8 (A9)` = c(0.96, 1.07, 0.92, 0.84, 1.01, 1, 0.8, 0.99, 0.88, 0.8),
  `4 (A10)` = c(0.64, 0.75, 0.68, 0.72, 0.74, 0.83, 0.84, 0.71, 0.67, 0.7),
  `2 (A11)` = c(0.52, 0.44, 0.62, 0.6, 0.54, 0.6, 0.69, 0.61, 0.68, 0.61),
  `1 (A12)` = c(0.41, 0.45, 0.46, 0.51, 0.49, 0.44, 0.44, 0.43, 0.45, 0.43),
  
  check.names = FALSE
)

# retornar o dataframe criado
dados_medicoes
##    Medicao 256 (A4) 128 (A5) 64 (A6) 32 (A7) 16 (A8) 8 (A9) 4 (A10) 2 (A11)
## 1     Med1     3.88     3.06    2.60    1.52    1.32   0.96    0.64    0.52
## 2     Med2     3.87     2.96    2.20    1.45    1.21   1.07    0.75    0.44
## 3     Med3     3.89     2.89    2.68    1.59    1.13   0.92    0.68    0.62
## 4     Med4     3.77     2.94    2.44    1.51    1.55   0.84    0.72    0.60
## 5     Med5     3.90     2.82    2.42    1.46    1.21   1.01    0.74    0.54
## 6     Med6     3.78     2.77    2.16    1.42    1.16   1.00    0.83    0.60
## 7     Med7     3.89     2.94    2.56    1.55    1.28   0.80    0.84    0.69
## 8     Med8     3.88     2.77    2.29    1.82    1.19   0.99    0.71    0.61
## 9     Med9     3.94     2.70    2.31    1.59    1.26   0.88    0.67    0.68
## 10   Med10     3.98     2.83    2.53    1.52    1.19   0.80    0.70    0.61
##    1 (A12)
## 1     0.41
## 2     0.45
## 3     0.46
## 4     0.51
## 5     0.49
## 6     0.44
## 7     0.44
## 8     0.43
## 9     0.45
## 10    0.43

3 Trabalhando com os dados

Bom, agora que já temos os dados do experimento no R, começaremos a analisá-los. Na seção introdutória foi mencionado que o número de linhas da tabela de medições deve ser o número de medições realizadas para cada tamanho de folha de papel. Para extraírmos esse valor na nossa própria base de dados, basta usarmos a função nrow, que conta o número de linhas excluindo o cabeçalho.

# calcula numero de medicoes
n <- nrow(dados_medicoes)
# retorna o numero de medicoes
n
## [1] 10

A função retornou o valor 10, que está correto, já que temos 10 medições na nossa tabela de medições.
Agora vamos gerar estatísticas básicas para cada um dos tratamentos do experimento. Utilizaremos as funções mean e sd para calcular a média e o desvio padrão de cada tratamento. Conforme o exemplo a seguir.

mean(dados_medicoes[,2])
## [1] 3.878
sd(dados_medicoes[,2])
## [1] 0.06356099

Acabamos de calcular a média e o desvio padrão da coluna 2 da nossa tabela de medições, que se refere às medições do tratamento “256 (A4)”. No entanto, precisamos realizar esse cálculo para todas as colunas, exceto a primeira que é o campo identificador das medições. Para isso, usaremos a função sapply. Esta função recebe dois argumentos obrigatórios: o primeiro deve ser um data.frame ou uma matriz, e o segundo é a função que deverá ser aplicada em cada uma das colunas da base declarada no primeiro argumento. Vejamos.

# calcular media
medias <- sapply(dados_medicoes[,-1], mean)
medias
## 256 (A4) 128 (A5)  64 (A6)  32 (A7)  16 (A8)   8 (A9)  4 (A10)  2 (A11) 
##    3.878    2.868    2.419    1.543    1.250    0.927    0.728    0.591 
##  1 (A12) 
##    0.451
# calcular desvio padrao
desvpad <- sapply(dados_medicoes[,-1], sd)
desvpad
##   256 (A4)   128 (A5)    64 (A6)    32 (A7)    16 (A8)     8 (A9)    4 (A10) 
## 0.06356099 0.10921945 0.17546446 0.11255122 0.11962906 0.09416888 0.06511528 
##    2 (A11)    1 (A12) 
## 0.07445356 0.02960856

A partir das médias e desvios padrão, podemos calcular o coeficiente de variação de cada tratamento.

# calcular coeficiente de variacao
cvar <- desvpad / medias * 100
cvar
##  256 (A4)  128 (A5)   64 (A6)   32 (A7)   16 (A8)    8 (A9)   4 (A10)   2 (A11) 
##  1.639015  3.808210  7.253595  7.294311  9.570325 10.158455  8.944407 12.597896 
##   1 (A12) 
##  6.565090

Por fim, podemos completar a tabela de parâmetros dos tratamentos com as estatísticas calculadas.

# adicionar os vetores calculados como colunas na tabela
dados_tratamentos$Media_Diametros <- medias
dados_tratamentos$DesvPad_Diametros <- desvpad
dados_tratamentos$CV_Diametros <- cvar

dados_tratamentos
##    Tamanho   Area  Massa Media_Diametros DesvPad_Diametros CV_Diametros
## 1 256 (A4) 0.0624 4.6800           3.878        0.06356099     1.639015
## 2 128 (A5) 0.0312 2.3400           2.868        0.10921945     3.808210
## 3  64 (A6) 0.0156 1.1700           2.419        0.17546446     7.253595
## 4  32 (A7) 0.0078 0.5850           1.543        0.11255122     7.294311
## 5  16 (A8) 0.0039 0.2920           1.250        0.11962906     9.570325
## 6   8 (A9) 0.0019 0.1420           0.927        0.09416888    10.158455
## 7  4 (A10) 0.0010 0.0750           0.728        0.06511528     8.944407
## 8  2 (A11) 0.0005 0.0375           0.591        0.07445356    12.597896
## 9  1 (A12) 0.0002 0.0150           0.451        0.02960856     6.565090

Agora que completamos a tabela de parâmetros com médias e medidas de dispersão, vamos visualizar como a massa das bolinhas de papel varia em função das médias dos diâmetros das mesmas, utilizando a função plot.

plot(Massa ~ Media_Diametros, 
     data = dados_tratamentos, 
     type = 'b', # tipo de grafico (b = linha + ponto, p = ponto, l = linha)
     
     xlab = 'Diâmetro (cm)', # texto do eixo x
     ylab = 'Massa (g)', # texto do eixo y
     main = 'Massa em função dos diâmetros médios das bolas de papel', # titulo
     
     col = 'black', # cor da linha/pontos (red, brown, gray, green, darkgreen, darkred, etc.)
     lty = 2, # tipo de linha
     pch = 20, # tipo de marcador
     cex = 2 # tamanho do marcador
     )

O gráfico acima demonstra que a massa cresce exponencialmente a medida que o diâmetro da bolinha de papel cresce.
Outro gráfico interessante que podemos gerar é o boxplot das medições de diâmetro das bolinhas de papel. Por meio deste, podemos enxergar a dispersão das medidas dos diâmetros para cada um dos tratamentos.

# visualizar dispersao das medicoes com o grafico boxplot
boxplot(dados_medicoes[, -1], 
        ylab = 'Diâmetro (cm)', # texto do eixo y
        xlab = 'Tamanho do papel', # texto do eixo x
        main = 'Boxplot das medições dos diâmetros' # titulo
        )
# adicionar marcador para a media
points(colMeans(dados_medicoes[-1]), 
       pch = 20, # tipo de marcador
       cex = 2, # tamanho do marcador
       col = 'red' # cor do marcador
       )
# adicionar legenda
legend(8.5, # posicao da legenda no eixo x
       4, # posicao da legenda no eixo y
       legend = c('Média', 'Mediana'), # texto dos itens
       col = c('red','black'), # cor dos itens
       pch = c(20,NA), # tipos de marcadores
       lty = c(NA,1), # tipos de linhas
       lwd = c(NA, 3) # espessura de linhas
       )

4 Modelagem

Na seção anterior foi possível observar que a relação entre massa e diâmetro da bolinha de papel é exponencial. Esta constatação permite gerarmos um modelo matemático que nos retorne a massa esperada para uma bolinha de diâmetro médio qualquer, desde que a gramatura do papel seja a mesma da utilizada para ajuste do modelo.
Para aplicarmos uma regressão linear, vamos linearizar a relação entre a variável dependente e a independente. Podemos fazer isso transformando ambas as variáveis ao calcularmos o logaritmo natural das mesmas (função log). Vejamos como elas se relacionam após a transformação.

# visualizar relacao linearizada entre massa e diametros
plot(log(Massa) ~ log(Media_Diametros), 
     data = dados_tratamentos, 
     type = 'b', # tipo de grafico (b = linha + ponto, p = ponto, l = linha)
     
     xlab = 'ln do Diâmetro', # texto do eixo x
     ylab = 'ln da Massa', # texto do eixo y
     main = 'Linearização da reação Massa x Diâmetro', # titulo
     
     col = 'black', # cor da linha/pontos (red, brown, gray, green, darkgreen, darkred, etc.)
     lty = 2, # tipo de linha
     pch = 20, # tipo de marcador
     cex = 2 # tamanho do marcador
)

Observa-se que a transformação realizada foi suficiente para linearizarmos a relação entre as mesmas. Deste modo, podemos aplicar uma regressão linear em que o logaritmo natural da massa é função do logaritmo natural do diâmetro. Mas antes, vamos calcular a correlação entre as duas variáveis transformadas. Utilizaremos a função cor.test para calcular o coeficiente de correlação de Pearson e avaliar a significância do valor resultante.

cor.test(log(dados_tratamentos$Massa), log(dados_tratamentos$Media_Diametros),
         method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  log(dados_tratamentos$Massa) and log(dados_tratamentos$Media_Diametros)
## t = 30.782, df = 7, p-value = 9.858e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9819304 0.9992575
## sample estimates:
##       cor 
## 0.9963264

O valor próximo de 1 indica que as variáveis transformadas são altamente correlacionadas. O p-valor próximo de zero permite concluir que a correlação é estatisticamente significativa e não fruto de mero acaso.
Seguiremos com a regressão linear (função lm) e na sequência usaremos a função summary para retornarmos as estatísticas de ajuste.

# ajustar modelo de regressao
modelo <- lm(log(Massa)~log(Media_Diametros),
             data = dados_tratamentos)

# estatisticas de ajuste
summary(modelo)
## 
## Call:
## lm(formula = log(Massa) ~ log(Media_Diametros), data = dados_tratamentos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.24975 -0.07046  0.01655  0.12150  0.23319 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.89088    0.06250  -30.25 1.11e-08 ***
## log(Media_Diametros)  2.58584    0.08401   30.78 9.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1768 on 7 degrees of freedom
## Multiple R-squared:  0.9927, Adjusted R-squared:  0.9916 
## F-statistic: 947.5 on 1 and 7 DF,  p-value: 9.858e-09

As estatísticas demonstram que os dois coeficientes estimados para a equação ajustada são significativos. No entanto, para estimarmos o erro padrão obtido precisaremos transformar os valores estimados pelo modelo para a sua unidade original. Para tal, basta aplicarmos a função exp nos valores estimados pelo modelo (estimativas são geradas pela função predict).

# predicoes
dados_tratamentos$Massa_estimada <- exp(predict(modelo))

# retornar os valores estimados
dados_tratamentos$Massa_estimada
## [1] 5.02164441 2.30158381 1.48190818 0.46332362 0.26877822 0.12407198 0.06641926
## [8] 0.03874021 0.01925556

Após estimarmos a massa das bolinhas com o modelo ajustado, podemos calcular os desvios entre os dados estimados e observados, calcular a soma dos quadrados dos resíduos e, por fim, o erro padrão das estimativas.

# calculo do erro padrao das estimativas
dados_tratamentos$Desvio <- dados_tratamentos$Massa - dados_tratamentos$Massa_estimada
sqr <- sum(dados_tratamentos$Desvio^2) # soma dos quadrados dos residuos
n <- nrow(dados_tratamentos) # numero de observacoes
syx <- sqrt( sqr / (n - 2) ) # erro padrao absoluto

# retornar o erro padrao absoluto
syx
## [1] 0.1817543
# retornar o erro padrao relativo
syx / mean(dados_tratamentos$Massa) * 100
## [1] 17.52037

Os cálculos realizados, indicam que o erro padrão a ser considerado quando aplicamos o modelo ajustado para estimar a massa de uma bolinha é de 0,2 gramas, ou 19,9%, para mais ou para menos. Agora vamos observar os desvios relativos dos dados observados em relação à linha de ajuste.

# calcular o desvio relativo
dados_tratamentos$Desvio_relativo <- dados_tratamentos$Desvio/dados_tratamentos$Massa_estimada

# estimados vs observados
plot(Desvio_relativo ~ Massa_estimada, 
     data = dados_tratamentos,
     
     xlab = 'Massa estimada (g)', # texto do eixo x
     ylab = 'Desvio da estimativa (%)', # texto do eixo y
     main = 'Dispersão de Resíduos do Modelo Ajustado', # titulo
     ylim = c(-.5,.5), # limites do eixo x
     
     col = 'black', # cor dos pontos (red, brown, gray, green, darkgreen, darkred, etc.)
     pch = 20, # tipo de marcador
     cex = 2 # tamanho dos marcadores
     )
# adicionar linha horizontal
abline(h = 0, # altura do eixo y onde a linha horizontal passa
       lty = 2, # tipo de linha (1 = solida, 2 = tracejada, 3 = pontilhada)
       lwd = 2, # espessura da linha
       col = 'black' # cor da linha
       )

5 Dimensão Fractal

A dimensão fractal é equivalente ao coeficiente angular do modelo que ajustamos na seção anterior.
Para extrairmos o coeficiente angular da equação gerada, basta aplicarmos a função coef no modelo salvo. Esta função retornará dois valores, sendo o primeiro o intercepto, e o segundo o coeficiente angular.

# retorna os coeficientes do modelo
coef(modelo)
##          (Intercept) log(Media_Diametros) 
##            -1.890882             2.585840
# armazena o coeficiente angular do modelo
coef_angular <- coef(modelo)[[2]]

# retorna a dimensão fractal
coef_angular
## [1] 2.58584
# arredondar o valor da dimensao fractal (uma casa decimal)
round(coef_angular, 1)
## [1] 2.6

Fim do exercício.