<!-- mathjax: "/home/fpmayer/GitHub/MathJax/MathJax.js?config=TeX-AMS-MML_HTMLorMML" -->

```{r setup, include=FALSE, purl=FALSE, eval=TRUE}
opts_chunk$set(
    ## knitr options
    cache = TRUE,
    tidy = FALSE,
    comment = NA,
    fig.width = 7,
    fig.height = 5,
    fig.align = "center",
    #dpi = 60, ## aumentar o dpi para exibir maior.
    dev = "png"
    # fig.path = "figures/",
    )
```

# Experimentos fatoriais $2^k$

## Análise de um experimento $2^3$

Exemplo 14-2, Montgomery (EAPE), pg. 375.

Um engenheiro mecânico está estudando a rugosidade da superfície de
uma peça usinada. Três fatores são de interesse: taxa de alimentação
(A), profundidade do corte (B), e ângulo da ferramenta (C). Todos os
três fatores tem 2 níveis, e um planejamento fatorial foi feito com 2
réplicas.

A representação geométrica é

```
   nom A B C
 1   1 - - -
 2   a + - -
 3   b - + -
 4  ab + + -
 5   c - - +
 6  ac + - +
 7  bc - + +
 8 abc + + +
```

```
        bc ------- abc
        .|        .|
 [+]  c__|______ac |
      |  |      |  |
  C   |  b -----|- ab [+]
      | .       | .   B
 [-] (1)________a  [-]
     [-]   A   [+]
```

```{r}
## Importação
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-2.txt"
dados <- read.table(url, header = TRUE)
str(dados)

## Visualizacao
par(mfrow=c(2,2))
with(dados, interaction.plot(Alim, Prof, Rugosidade))
with(dados, interaction.plot(Alim, Angulo, Rugosidade))
with(dados, interaction.plot(Prof, Angulo, Rugosidade))
par(mfrow=c(1,1))

wireframe(Rugosidade ~ Alim * Prof, data = dados)
wireframe(Rugosidade ~ Alim * Angulo, data = dados)
wireframe(Rugosidade ~ Prof * Angulo, data = dados)

xyplot(Rugosidade ~ Alim | Prof, groups = Angulo, data = dados,
       type = c("p", "a"), auto.key=TRUE)
xyplot(Rugosidade ~ Alim, groups = interaction(Prof, Angulo, sep = ":"),
       data = dados, type = c("p", "a"), auto.key = TRUE)

##----------------------------------------------------------------------
## Analise usando as funcoes do R

## Transforma as variáveis em fator
dados <- transform(dados, Alim = factor(Alim),
                   Prof = factor(Prof),
                   Angulo = factor(Angulo))

## Ajusta o modelo, especificando o conraste soma zero
m0 <- lm(Rugosidade ~ Alim * Prof * Angulo, data = dados,
         contrasts = list(Alim = "contr.sum",
                          Prof = "contr.sum",
                          Angulo = "contr.sum"))
anova(m0)
summary(m0)
```

Repare que os sinais dos coeficientes estimados pela função `lm()` são
contrários aos sinais que calculamos a mão. Isso ocorre porque no
contraste soma zero do R, é assumido que o último nível é igual ao
negativo do primeiro, ou seja,
$$
\alpha_2 = -\alpha_1
$$

Como, por padrão, assumimos calcular os efeitos usando a diferença de
médias sempre no formato `efeito = nivel_alto - nivel_baixo` então
precisamos assumir que
$$
\alpha_1 = -\alpha_2
$$

Do ponto de vista prático, o inversão dos sinais não altera em nada a
interpretação e os resultados do modelo, mas podemos fazer a saída do R
ficar compatível com a nossa notação. Para isso, basta utilizarmos a
codificação *numérica* -1 e 1 para os níveis baixo e alto dos fatores,
ao invés de utilizarmos a codificação padrão do R.

```{r}
## Funcao para converter os niveis de um fator na codificacao -1 (baixo)
## e 1 (alto). Valido apenas para fatores com 2 niveis.
convert <- function(x){
    if(!is.factor(x)){
        x <- as.factor(x)
    }
    if(!identical(nlevels(x), 2L)){
        stop("O fator deve ter exatamente 2 niveis")
    }
    x <- as.numeric(x)
    cod <- (x - mean(x))/0.5
    return(cod)
}
```

```{r}
## Convertendo os níveis dos fatores para -1 e 1
dados <- transform(dados, Alim.n = convert(Alim),
                   Prof.n = convert(Prof),
                   Angulo.n = convert(Angulo))

## Ajusta o modelo com a nova codificação
m1 <- lm(Rugosidade ~ Alim.n * Prof.n * Angulo.n, data = dados)
anova(m1)
summary(m1)
```

### Suposições do modelo

Como todo modelo de ANOVA:

1. $\epsilon_{ijkl} \sim \text{N}(0, \sigma^2)$
2. $\epsilon_{ijkl}$ são independentes
3. Homogeneidade de variâncias

### Avaliação do modelo

Quando ajustamos um modelo a um conjunto de dados, podemos ter
basicamente dois objetivos:

1. Verificar quais variáveis influenciam na média da variável resposta
   $Y$ 
2. Fazer predições com base nos parâmetros estimados
3. Outros: estimar parâmetros e fazer inferência

Para verificar as variáves que infuenciam na rugosidade, vemos a ANOVA do
modelo, 

```{r}
anova(m1)
```

Se o nosso objetivo fosse simplesmente construir um modelo para avaliar
quais efeitos e interações influemciam no processo, então a análise
poderia parar aqui, pois já temos esse resultado. Então deveria se
proceder para a análise dos resíduos para verificar as suposições do
modelo.

Caso o nosso objetivo seja fazer uma predição para a variável resposta,
então devemos ajustar o modelo mais **parcimonioso** possível. Um modelo
parcimonioso é aquele que consegue explicar a maior parte da
variabilidade dos dados, com o menor número de parâmetros possível.

Asumindo uma avaliação não muito "fechada" (ou restrita) em um ponto
de corte de p-valor especifico, notamos que `Alim` e `Prof` são fatores
influentes. Como os dois são influentes, e a interação entre eles é
"marginalmente" significativa, vamos manter estes dois fatores
principais e a interação para construir um novo modelo.

```{r}
m2 <- update(m1, . ~ Alim.n * Prof.n)
```

Como, por construção, o modelo é ortogonal, então as somas de quadrados
e os coeficientes estimados também serão ortogonais. Dessa forma, a
retirada de termos do modelo não interfere nos coeficientes
estimados. No caso da ANOVA, as somas de quadrados permanecem as mesmas,
mas a $SQRes$ e os $GL$ associados aumentam, pois as $SQ$ dos fatores
retirados são incorporados aos resíduos. Veja:

```{r}
anova(m1)
anova(m2)
```

Porque se alteram os p-valores?

```{r}
## Exemplo para Alim.n
## Usando o modelo completo
pf(anova(m1)$"F value"[1], df1 = anova(m1)$"Df"[1],
   df2 = anova(m1)$"Df"[8], lower.tail = FALSE)
## Exemplo para o modelo sob avaliação
pf(anova(m2)$"F value"[1], df1 = anova(m2)$"Df"[1],
   df2 = anova(m2)$"Df"[4], lower.tail = FALSE)
```

```{r, eval=FALSE, include=FALSE}
## par(mfrow = c(1,2))
## plot(seq(0, 10, length=100), type = "l", xlab = "X", ylab = "f(x)",
##      y = df(x = seq(0, 10, length=100), df1 = 1, df2 = 8),
##      main = expression(list(nu[1] == 1, nu[2] == 8)))
## abline(v = qf(0.05, df1 = anova(m1)$"Df"[1],
##               df2 = anova(m1)$"Df"[8], lower.tail = FALSE))
## plot(seq(0, 10, length=100), type = "l", xlab = "X", ylab = "f(x)",
##      y = df(x = seq(0, 10, length=100), df1 = 1, df2 = 12),
##      main = expression(list(nu[1] == 1, nu[2] == 12)))
## abline(v = qf(0.05, df1 = anova(m2)$"Df"[1],
##               df2 = anova(m2)$"Df"[4], lower.tail = FALSE))

```

Porque mudam os graus de liberdade do numerador e do denominador, altera
a distribuição.

Para comparar se o modelo reduzido pode de fato ser utilizado no lugar
do modelo completo, podemos fazer o **teste da razão de
verossimilhança** (TRV) entre os dois modelos encaixados.

```{r}
anova(m2, m1)
```

Como a diferença das somas de quadrados residuais entre os dois modelos
não é significativa, então concluímos que o modelo com menos parâmetros
pode ser utilizado no lugar do modelo com mais parâmetros - este é o
**princípio da parcimônia**.

Podemos ainda estar em dúvida sobre a inclusão do termo de interação,
que é "marginalmente" significativo. Para avaliar se de fato esse termo
é importante, podemos ajustar um novo modelo sem esse termo

```{r}
m3 <- update(m2, . ~ Alim.n + Prof.n)
anova(m3)
```

Assim como antes, podemos usar o TRV para comparar os modelos encaixados

```{r}
anova(m3, m2, m1)
```

Com esse resultado, podemos então concluir que a interação não é
importante, e vamos dizer que o último modelo é o modelo escolhido.

### Análise dos resíduos

Para verificar as suposições feitas para o modelo, procedemos com a
análise dos resíduos do modelo final

```{r}
## Residuos de modelos
res <- residuals(m3)
## Quantis normais
qqnorm(res); qqline(res)
## Residuos vs fatores
library(lattice)
xyplot(res ~ Alim, data = dados,
       panel = function(...){
           panel.xyplot(...)
           panel.abline(h = 0, col = 2, lty = 2)
       })
xyplot(res ~ Prof, data = dados,
       panel = function(...){
           panel.xyplot(...)
           panel.abline(h = 0, col = 2, lty = 2)
       })
```

### Predição

A partir dessa avaliação, podemos fazer a predição usando os coeficientes
estimados do modelo

```{r}
summary(m3)
```

```{r}
## Predição para as combinações únicas dos fatores
pred <- expand.grid(Alim.n = unique(dados$Alim.n),
                    Prof.n = unique(dados$Prof.n))
pred$Rugosidade <- predict(m3, newdata = pred)
pred

## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(Alim.n = seq(-1, 1, length.out = 20),
                    Prof.n = seq(-1, 1 ,length.out = 20))
pred$Rugosidade <- predict(m3, newdata = pred)
## Vários formas de visualizar
wireframe(Rugosidade ~ Alim.n + Prof.n, data = pred)
wireframe(Rugosidade ~ Alim.n + Prof.n, data = pred, drape = TRUE)
levelplot(Rugosidade ~ Alim.n + Prof.n, data = pred)
levelplot(Rugosidade ~ Alim.n + Prof.n, data = pred, cuts = 90,
          col.regions = heat.colors)
```

## Fatoriais $2^k$ para $k > 3$

Á medida que o número de fatores cresce em um experimento fatorial,
o número de efeitos que podem ser estimados também cresce. Por exemplo,

- $2^4$ possui 4 efeitos principais, 6 interações de segunda ordem, 4
  interações de terceira ordem, e 1 interação de quarta ordem
- $2^6$ possui 6 efeitos principais, 15 interações de segunda ordem, 20
  interações de terceira ordem, 15 interações de quarta ordem, 6
  interações de quinta ordem, e 1 interação de sexta ordem

Para um experimento $2^4$, a representação geométrica fica

```
     nom A B C D
 1   (1) - - - -
 2     a + - - -
 3     b - + - -
 4    ab + + - -
 5     c - - + -
 6    ac + - + -
 7    bc - + + -
 8   abc + + + -
 9     d - - - +
 10   ad + - - +
 11   bd - + - +
 12  abd + + - +
 13   cd - - + +
 14  acd + - + +
 15  bcd - + + +
 16 abcd + + + +
```

```
           D [-]             |             D [+]
                             |
        bc ------- abc       |         bcd ------- abcd
        .|        .|         |          .|        .|
 [+]  c__|______ac |         |   [+] cd__|_____acd |
      |  |      |  |         |        |  |      |  |
  C   |  b -----|- ab [+]    |    C   | bd -----|- abd [+]
      | .       | .   B      |        | .       | .    B
 [-] (1)________a  [-]       |   [-]  d _______ ad  [-]
     [-]   A   [+]           |       [-]   A   [+]
```

Em muitas situações, o **princípio da esparsidade dos efeitos** se
aplica, ou seja, o sistema geralmente é dominado pelos efeitos
principais e pelas interações de ordens baixas. As interações de
terceira ordem ou superiores geralmente são negligenciadas.

Como consequência, quando o número de fatores for moderadamente grande,
como $k \geq 4$ ou $5$, uma prática comum é correr somente uma réplica
do planejamento, e então combinar as interações de ordem alta como uma
estimativa do erro. No entanto, se eventualmente alguma interação de
ordem alta for significativa, então esse procedimento não é adequado.

Um método simples para verificar o tamanho dos efeitos de um fator é
construindo um gráfico das estimativas dos efeitos em uma escala de
probabilidade normal. Por construção, a distribuição dos efeitos
estimados possui média zero e variância $\sigma^2$. Portanto, os efeitos
que forem desprezíveis estarão em cima de uma linha reta nesse gráfico,
enquanto que efeitos significativos não terão média zero e estarão mais
afastados dessa linha.

### Exemplo de um experimento $2^4$

Exemplo 14-5, Montgomery (EAPE), pg. 385: processo de ataque químico
localizado sobre nitreto, por meio de uma sonda de plasma de pastilha
única. Variável resposta: taxa de ataque do nitreto de
silício. Variáveis explicativas:

```
                Fatores do experimento
 Nível     Espaçamento   Pressão    Vazão de C2F6   Potência
              (cm)       (m Torr)   (cm³ pad/min)     (W)
 Baixo (-)    0.80         450           125          275
 Alto  (+)    1.20         550           200          325
```

```{r}
##----------------------------------------------------------------------
## Importação
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/mont_14-5.txt"
dados <- read.table(url, header = TRUE)
## Definições do experimento
k <- 4
n <- 1

##----------------------------------------------------------------------
## Visualização
xyplot(y ~ A, groups = interaction(B, C, D, sep = ":"), data = dados,
       type=c("p", "a"), auto.key = TRUE)
xyplot(y ~ A | D, groups = interaction(B, C, sep = ":"), data = dados,
       type = c("p", "a"), auto.key = TRUE)
xyplot(y ~ A | B + C, groups = D, data = dados,
       type = c("p", "a"), auto.key = TRUE)

##----------------------------------------------------------------------
## Montando a tabela de sinais
tab <- model.matrix(~A * B * C * D, data = dados)
tab

##----------------------------------------------------------------------
## Estimativa dos efeitos
## Podemos usar as colunas da tabela de sinais para calcular os
## contrastes
y <- dados$y
## Contrastes
(contr <- t(tab[, -1]) %*% y)
## Efeitos = contraste/(n2^{k-1})
(ef <- contr/(n * 2^(k - 1)))
```

Para testar quais efeitos são significativos, poderíamos pensar
inicialmente em fazer uma ANOVA. No entanto, este experimento não possui
repetição, e temos um parâmetro para cada observação, fazendo ocm que o
modelo seja saturado. Dessa forma, não temos como avaliar a
significância dos efeitos. Veja o resultado da ANOVA para o modelo
completo:

```{r}
anova(lm(y ~ (A * B * C * D), data = dados))
```

Uma forma de verificar quais efeitos são importantes é através de um
gráfico de probabilidade normal para a estimativa dos
esfeitos. Lembre-se que, sob a hipótses nula, todos os parâmetros do
modelo são iguais a zero, e pela definição do termo de erro do modelo,
se os parâmetros são iguais a zero, então sobra apenas o erro que possui
média 0 e variância contante $\sigma^2$. Dessa forma, se a hipótese nula
for verdadeira, esperamos que os efeitos tenham também média 0, e fiquem
em cima da linha em um gráfico de probabilidade normal. Os efeitos que
se afastarem muito da linha são aqueles que possuem média diferente de
zero, e portanto, são aqueles que temos interesse.

```{r}
##----------------------------------------------------------------------
## Gráfico de probabilidade normal dos efeitos estimados
qqaux <- qqnorm(ef, col = 2, pch = 19); qqline(ef)
text(qqaux$x, qqaux$y, rownames(qqaux$y), cex = 0.8, pos = 3)
```

Através do gráfico acima, vemos que termos de ordem maior que 2 ficaram
no centro, o que mostra que interações de ordem grande podem ser
consideradas como nulas. Portanto podemos ajustar agora um modelo
considerando apenas as interações de segunda ordem.

```{r}
##----------------------------------------------------------------------
## Ajuste do modelo
## Como vimos no gráfico anterior, as interações de ordem maior do que 2
## são desnecessárias, portanto vamos ajustar o modelo com interações de
## até segunda ordem apenas
m0 <- lm(y ~ (A + B + C + D)^2, data = dados)
anova(m0)
```

Note que a variabilidade das interações não consideradas foi integrada
na soma de quadrado dos resíduos, e assim podemos ter uma estimativa
para $\sigma^2$.

Com esse resultado, confirmamos a hipótese leventada no gráfico de
probabilidades normais de que apenas 2 fatores e uma interação são
importantes. Por isso, podemos atualizar o modelo apenas com estes
termos

```{r}
m1 <- update(m0, . ~ A * D, data = dados)
anova(m1)

## Fazendo o TRV para confirmar o modelo
anova(m1, m0)
```

#### Análise dos resíduos

```{r}
##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m1)
## Quantis normais
qqnorm(res); qqline(res)
## Residuos vs fatores
library(lattice)
xyplot(res ~ A, data = dados,
       panel = function(...){
           panel.xyplot(...)
           panel.abline(h = 0, col = 2, lty = 2)
       })
xyplot(res ~ D, data = dados,
       panel = function(...){
           panel.xyplot(...)
           panel.abline(h = 0, col = 2, lty = 2)
       })
xyplot(res ~ interaction(A, D, sep = ":"), data = dados,
       panel = function(...){
           panel.xyplot(...)
           panel.abline(h = 0, col = 2, lty = 2)
       })
```

#### Predições

```{r}
##----------------------------------------------------------------------
## Predições

## Predição para as combinações únicas dos fatores
pred <- expand.grid(A = unique(dados$A),
                    D = unique(dados$B))
pred$y <- predict(m1, newdata = pred)
pred
```

Note como fica a predição para o conjunto de dados

```{r}
cbind(dados, predict(m1))
```

As médias estimadas só levam em consideração a mudança nos níveis dos
fatores importantes para o processo (nesse caso A, D e A:D),
independente dos níveis dos outros fatores.

```{r}
## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(A = seq(-1, 1, length.out = 30),
                    D = seq(-1, 1 ,length.out = 30))
pred$y <- predict(m1, newdata = pred)
## Vários formas de visualizar
wireframe(y ~ A + D, data = pred, drape = TRUE)
levelplot(y ~ A + D, data = pred)
levelplot(y ~ A + D, data = pred, cuts = 90,
          col.regions = heat.colors)
```

#### EXTRA: visualização dos resultados

```{r}
##----------------------------------------------------------------------
## Carregando função para fazer contornos de nível
source("panel.3d.contour.R")
##----------------------------------------------------------------------
## Gráfico com contornos de nível
display.brewer.all()                        # palhetas de cores predefinidas
colr <- brewer.pal(11, "Spectral")          # escolhe uma palheta
colr <- colorRampPalette(colr, space="rgb") # interpola para ter mais níveis

wireframe(y ~ A + D, data = pred,
          scales = list(arrows = FALSE),
          zlim = c(300, 1100),
          col = "gray50", col.contour = 1,
          panel.3d.wireframe = "panel.3d.contour",
          type = c("on", "bottom"),
          col.regions = colr(100), drape = TRUE, alpha.regions = 0.5,
          screen = list(z = 40, x = -70))
```

```{r, eval=FALSE}
##----------------------------------------------------------------------
## Escolhendo a posição de observação do gráfico iterativamente
require(rpanel)
wire <- function(panel){
  print(
  wireframe(y ~ A + D, data = pred,
            scales = list(arrows = FALSE),
            zlim = c(300, 1100),
            col = "gray50", col.contour = 1,
            panel.3d.wireframe = "panel.3d.contour",
            type = c("on","bottom"),
            col.regions = colr(100), drape = TRUE, alpha.regions = 0.5,
            screen = list(x = panel$x, z = panel$z, y = panel$y))
    )
  panel
}

panel <- rp.control()
rp.slider(panel, x, -90, 90, initval=-70, showvalue=TRUE, action=wire)
rp.slider(panel, z, -90, 90, initval=40, showvalue=TRUE, action=wire)
rp.slider(panel, y, -90, 90, initval=0, showvalue=TRUE, action=wire)
```

A partir desse ponto, o experimentador possui duas opções:

1. Fazer um novo experimento, agora com apenas estes dois fatores e com
   um número maior de repetições: **projeção**
2. Fazer um novo experimento com estes dois fatores e, se desejado,
   outros fatores que não foram considerados no primeiro
   experimento. Nesse caso é preciso lembrar que as podem haver
   interações entre estes novos fatores e aqueles que foram descartados
   no primeiro experimento, e isto não está sendo considerado neste
   caso.

### Análise de um experimento fatorial $2^6$

```{r, eval=FALSE, echo=FALSE, include=FALSE}
## k <- 6
## input <- replicate(k, c(-1,1), simplify=FALSE)
## names(input) <- LETTERS[1:k]
## dados <- do.call(expand.grid, input)
## dados$y <- c(36,35,37,34,30,32,28,38,40,32,32,28,34,26,30,28,
##              40,30,35,35,32,35,36,32,20,35,35,28,27,36,36,35,
##              32,34,30,24,30,20,32,25,20,20,22,35,26,15,19,24,
##              22,23,22,18,20,20,20,36,20,11,20,35,16,32,20,20)
```

```{r}
##----------------------------------------------------------------------
## Dados
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/fatorial_2-6.txt"
dados <- read.table(url, header = TRUE)
str(dados)

## Definições
k <- 6
n <- 1

##----------------------------------------------------------------------
## Visualização
xyplot(y ~ A | B * C, groups = interaction(D, E, F), data = dados,
       type = c("p", "a"))
xyplot( y ~A | B * C * D, groups = interaction(E, F), data = dados,
       type = c("p", "a"))

##----------------------------------------------------------------------
## Montando a tabela de sinais
tab <- model.matrix(~A * B * C * D * E * F, data = dados)
tab

##----------------------------------------------------------------------
## Estimativa dos efeitos
## Podemos usar as colunas da tabela de sinais para calcular os
## contrastes
y <- dados$y
## Contrastes
(contr <- t(tab[, -1]) %*% y)
## Efeitos = contraste/(n2^{k-1})
(ef <- contr/(n * 2^(k - 1)))

##----------------------------------------------------------------------
## Gráfico de probabilidade normal dos efeitos estimados
qqaux <- qqnorm(ef, col = 2, pch = 19); qqline(ef)
text(qqaux$x, qqaux$y, rownames(qqaux$y), cex = 0.8, pos = 3)

##----------------------------------------------------------------------
## Ajuste do modelo
m0 <- lm(y ~ (A + B + C + D + E + F)^4, data = dados)
anova(m0)

m1 <- update(m0, . ~ (A + B + C + D + E + F)^3)
anova(m1)

m2 <- update(m1, . ~ (A + B + C + D + E + F)^2)
anova(m2)

m3 <- update(m2, . ~ D + F + (C + E)^2)
anova(m3)

m4 <- update(m3, . ~ D + F)
anova(m4)

## Fazendo o TRV para confirmar o modelo
anova(m4, m3, m2, m1, m0)

##----------------------------------------------------------------------
## Residuos do modelo
res <- residuals(m4)
## Quantis normais
qqnorm(res); qqline(res)

##----------------------------------------------------------------------
## Predições

## Predição para as combinações únicas dos fatores
pred <- expand.grid(D = unique(dados$D),
                    F = unique(dados$F))
pred$y <- predict(m4, newdata = pred)
pred

## Predição para um intervalo de valores entre os níveis baixo e alto
## dos fatores
pred <- expand.grid(D = seq(-1, 1, length.out = 30),
                    F = seq(-1, 1 ,length.out = 30))
pred$y <- predict(m4, newdata = pred)
## Vários formas de visualizar
wireframe(y ~ D + F, data = pred, drape = TRUE)
levelplot(y ~ D + F, data = pred)
levelplot(y ~ D + F, data = pred, cuts = 90,
          col.regions = heat.colors)
```
