Introdução

Este documento apresenta uma metodologia de Superfície de Resposta (MSR) utilizando o R. A análise inclui ajustes de modelos de primeira e segunda ordem, análise de resíduos e gráficos que representam a relação entre variáveis preditoras (x1, x2) e a variável resposta (y).

Instalação e Carregamento de Pacotes

Função para verificar se os pacotes necessários estão instalados e carregá-los:

# Função para verificar e instalar pacotes
check_and_install <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
  library(pkg, character.only = TRUE)
}

# Verificar e instalar pacotes necessários
check_and_install("rsm")
check_and_install("car") # Alternativa ao alr3

Ajuste de Modelos

Ajuste de Modelo de Primeira Ordem

O modelo de primeira ordem ajusta uma relação linear simples entre as variáveis preditoras (x1, x2) e a variável resposta (y).

# Ajuste de modelo de primeira ordem e análise inicial
ajuste_modelo_primeira_ordem <- function(dados) {
  modelo <- lm(y ~ x1 + x2, data = dados)
  print(summary(modelo))
  print(anova(modelo))
  return(modelo)
}

Ajuste de Modelo com Termos Quadráticos e Interação

Este modelo considera termos mais complexos, incluindo interações (x1:x2) e efeitos quadráticos (x1^2). Ele é usado para capturar relações não lineares na superfície de resposta.

# Ajuste de modelo com termos quadráticos e interação
ajuste_modelo_segunda_ordem <- function(dados) {
  modelo <- lm(y ~ x1 + x2 + x1:x2 + I(x1^2), data = dados)
  print(anova(modelo, test = "F"))
  return(modelo)
}

Construção de Gráficos

Superfície de Resposta

Os gráficos de superfície de resposta mostram como as variáveis x1 (tempo) e x2 (temperatura) afetam a variável resposta y (rendimento).

# Construção de gráficos
construir_graficos <- function(modelo, x_range, y_range, z_function, title_3d, title_2d) {
  z <- outer(x_range, y_range, z_function)

  # Gráfico tridimensional
  persp(x_range, y_range, z, theta = -35, phi = 5, expand = 0.5, col = "green",
        xlab = "Tempo", ylab = "Temperatura", zlab = "Rendimento", main = title_3d)

  # Gráfico de contorno
  contour(x_range, y_range, z, xlab = "Tempo", ylab = "Temperatura", main = title_2d)

  # Representações adicionais
  image(x_range, y_range, z, col = topo.colors(12))
  contour(x_range, y_range, z, add = TRUE)
  image(x_range, y_range, z, col = terrain.colors(12), main = "Gráfico com terrain.colors")
  image(x_range, y_range, z, col = heat.colors(12), main = "Gráfico com heat.colors")
}

Gráfico com Ponto Estacionário

Este gráfico destaca o ponto estacionário calculado na superfície de resposta.

# Gráfico com ponto estacionário destacado
grafico_ponto_estacionario <- function(modelo, x_range, y_range, z_function, x0) {
  z <- outer(x_range, y_range, z_function)
  image(x_range, y_range, z, col = topo.colors(12), main = "Superfície com Ponto Estacionário")
  contour(x_range, y_range, z, add = TRUE)
  points(x0[1], x0[2], col = "red", pch = 19)
}

Determinação do Ponto Estacionário

Esta função calcula o ponto estacionário do modelo ajustado, verificando se a matriz é válida.

# Determinação do ponto estacionário
determinar_ponto_estacionario <- function(modelo) {
  b <- matrix(c(modelo$coef[2], modelo$coef[3]))
  B <- matrix(c(modelo$coef[4], modelo$coef[6] / 2, modelo$coef[6] / 2, modelo$coef[5]), ncol = 2)

  if (det(B) == 0 || any(is.na(B))) {
    cat("Erro: A matriz B é singular ou contém valores ausentes. Não é possível determinar o ponto estacionário.\n")
    return(NULL)
  }

  x0 <- -0.5 * solve(B) %*% b

  tempo <- x0[1] * 5 + 85
  temperatura <- x0[2] * 5 + 175

  cat("Coordenadas do ponto estacionário:", x0, "\n")
  cat("Tempo:", tempo, "\n")
  cat("Temperatura:", temperatura, "\n")

  return(list(x0 = x0, tempo = tempo, temperatura = temperatura))
}

Análise Canônica

Verifica o tipo de ponto estacionário (máximo, mínimo ou sela).

analise_canonica <- function(modelo) {
  B <- matrix(c(modelo$coef[4], modelo$coef[6] / 2, modelo$coef[6] / 2, modelo$coef[5]), ncol = 2)

  if (any(is.na(B)) || any(is.infinite(B)) || det(B) == 0) {
    cat("Erro: A matriz B é inválida (valores ausentes, infinitos ou singular). Não é possível realizar a análise canônica.\n")
    return(NULL)
  }

  raizes <- eigen(B)
  cat("Raízes do modelo:", raizes$values, "\n")

  if (all(raizes$values < 0)) {
    cat("O ponto estacionário é um ponto de máximo.\n")
  } else if (all(raizes$values > 0)) {
    cat("O ponto estacionário é um ponto de mínimo.\n")
  } else {
    cat("O ponto estacionário é um ponto de sela.\n")
  }
}

Análise de Resíduos

Os gráficos abaixo avaliam se o modelo ajustado é adequado, verificando normalidade e homocedasticidade dos resíduos.

# Análise de resíduos e verificação de suposições
analise_residuos <- function(modelo) {
  residuos <- residuals(modelo)

  # Gráfico de probabilidade normal
  qqnorm(residuos, main = "Gráfico de Probabilidade Normal")
  qqline(residuos, col = "red")

  # Gráfico de preditos x resíduos
  preditos <- predict(modelo)
  plot(preditos, residuos, main = "Valores Preditos x Resíduos", pch = 16, xlab = "Valores Preditos", ylab = "Resíduos")
  abline(h = 0, col = "red")

  # Gráfico de resíduos x observações
  plot(residuos, main = "Resíduos x Observações", pch = 16, xlab = "Observações", ylab = "Resíduos")
  abline(h = 0, col = "red")
}

Execução do Código

A função principal executa todos os passos: ajuste dos modelos, gráficos de superfície, ponto estacionário e análise de resíduos.

# Executa todo o fluxo de análise
main <- function() {
  # Dados iniciais
  dados1 <- data.frame(
    x1 = c(-1, -1, 1, 1, 0, 0, 0, 0, 0),
    x2 = c(-1, 1, -1, 1, 0, 0, 0, 0, 0),
    y = c(39.3, 40, 40.9, 41.5, 40.3, 40.5, 40.7, 40.2, 40.6)
  )
  
  modelo1 <- ajuste_modelo_primeira_ordem(dados1)
  modelo2 <- ajuste_modelo_segunda_ordem(dados1)
  
  x_range <- seq(-1, 1, length = 10)
  y_range <- seq(-1, 1, length = 10)
  z_function <- function(x, y) modelo1$coef[1] + modelo1$coef[2] * x + modelo1$coef[3] * y
  
  construir_graficos(modelo1, x_range, y_range, z_function, "Superfície de Resposta 3D", "Contornos")
  
  ponto_estacionario <- determinar_ponto_estacionario(modelo2)
  if (!is.null(ponto_estacionario)) {
    grafico_ponto_estacionario(modelo2, x_range, y_range, z_function, ponto_estacionario$x0)
  }
  
  analise_canonica(modelo2)
  analise_residuos(modelo2)
}

# Executa o fluxo principal
main()
## 
## Call:
## lm(formula = y ~ x1 + x2, data = dados)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.244444 -0.044444  0.005556  0.055556  0.255556 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 40.44444    0.05729 705.987 5.45e-16 ***
## x1           0.77500    0.08593   9.019 0.000104 ***
## x2           0.32500    0.08593   3.782 0.009158 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1719 on 6 degrees of freedom
## Multiple R-squared:  0.941,  Adjusted R-squared:  0.9213 
## F-statistic: 47.82 on 2 and 6 DF,  p-value: 0.0002057
## 
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value   Pr(>F)    
## x1         1 2.40250 2.40250  81.339 0.000104 ***
## x2         1 0.42250 0.42250  14.304 0.009158 ** 
## Residuals  6 0.17722 0.02954                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## x1         1 2.40250 2.40250 55.8721 0.001713 **
## x2         1 0.42250 0.42250  9.8256 0.035030 * 
## I(x1^2)    1 0.00272 0.00272  0.0633 0.813741   
## x1:x2      1 0.00250 0.00250  0.0581 0.821316   
## Residuals  4 0.17200 0.04300                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Erro: A matriz B é singular ou contém valores ausentes. Não é possível determinar o ponto estacionário.
## Erro: A matriz B é inválida (valores ausentes, infinitos ou singular). Não é possível realizar a análise canônica.

Conclusão

Esta versão do código apresenta diversas melhorias em relação ao original, tornando-o mais organizado, robusto e compreensível. Primeiramente, a modularização do código em funções claramente definidas facilita a manutenção e o reaproveitamento em outros projetos. Cada função foi projetada para desempenhar uma tarefa específica, como ajuste de modelos, construção de gráficos ou análise de resíduos, o que melhora a legibilidade e reduz redundâncias.

Além disso, o novo código introduz verificações automáticas de pacotes necessários, garantindo que o ambiente esteja configurado corretamente antes da execução, algo ausente no original. Esse controle não apenas evita erros, mas também facilita a execução por outros usuários que possam não estar familiarizados com os requisitos do código.

O tratamento de erros foi significativamente aprimorado. Por exemplo, ao calcular o ponto estacionário ou realizar a análise canônica, são feitas verificações para garantir que os dados e as matrizes sejam válidos, evitando falhas inesperadas, como matrizes singulares ou coeficientes ausentes. Isso aumenta a confiabilidade dos resultados e reduz os riscos de interpretações equivocadas.