请求协助编写R语言DOE总体合意度3D绘图优化代码

Got it, let's tackle this step by step. You've got a DOE-generated model and want to visualize overall desirability in 3D to find optimal points—totally makes sense. I'll walk you through R code that aligns with the desirability framework you referenced, with clear explanations for each part.

First, let's recap the core idea: overall desirability (D) is a weighted geometric average of individual response desirabilities (d_i), where each d_i maps a response value to a 0-1 scale based on your optimization goals (maximize, minimize, target a specific value).

Here's the complete, runnable code with comments:

Step 1: Load Required Packages #

We'll use rsm for response surface model fitting, desirability for calculating desirability scores, and either rgl or plotly for 3D visualization.

# Install packages if you haven't already
# install.packages(c("rsm", "desirability", "rgl", "plotly"))

# Load libraries
library(rsm)
library(desirability)
library(rgl)
library(plotly)

Step 2: Prepare/Fit Your DOE Model #

I'll simulate a response surface dataset (matching typical DOE output) but you can replace this with your actual data. We'll fit quadratic response surface models for two responses (adjust this to match your number of responses).

# Set seed for reproducibility
set.seed(123)

# Create factor grids (replace with your actual factor ranges)
x1 <- seq(-2, 2, length.out = 20)  # Factor 1 (coded units)
x2 <- seq(-2, 2, length.out = 20)  # Factor 2 (coded units)
grid <- expand.grid(x1 = x1, x2 = x2)

# Simulate response data (replace with your actual response values)
grid$y1 <- 10 + 2*grid$x1 + 1.5*grid$x2 - 0.8*grid$x1^2 - 0.6*grid$x2^2 + rnorm(nrow(grid), 0, 0.5)
grid$y2 <- 50 - 3*grid$x1 + 2*grid$x2 - 1.2*grid$x1^2 - 0.9*grid$x2^2 + rnorm(nrow(grid), 0, 0.8)

# Fit quadratic response surface models (SO = second-order)
model_y1 <- rsm(y1 ~ SO(x1, x2), data = grid)
model_y2 <- rsm(y2 ~ SO(x1, x2), data = grid)

Step 3: Define Desirability Functions #

Customize these based on your optimization goals:

• dMax(): Maximize the response (set lower/upper bounds based on your acceptable range)
• dMin(): Minimize the response
• dTarget(): Target a specific response value

# Define individual desirability functions
# Example: Maximize y1 (acceptable range 10 to 15)
d_y1 <- dMax(model_y1, lower = 10, upper = 15)

# Example: Maximize y2 (acceptable range 40 to 55)
d_y2 <- dMax(model_y2, lower = 40, upper = 55)

# Calculate overall desirability (adjust weights if some responses are more important)
overall_desirability <- desirability(d_y1, d_y2, weights = c(1, 1))

Step 4: Generate 3D Desirability Surface #

We'll create a fine grid of factor values, predict overall desirability for each point, then plot it.

Option 1: Interactive 3D Plot with rgl (great for exploratory analysis)

# Create a dense grid for smooth visualization
x1_grid <- seq(-2, 2, length.out = 50)
x2_grid <- seq(-2, 2, length.out = 50)
plot_grid <- expand.grid(x1 = x1_grid, x2 = x2_grid)

# Predict overall desirability for each grid point
plot_grid$D <- predict(overall_desirability, newdata = plot_grid)

# Plot 3D surface
open3d()
surface3d(
  x = x1_grid, y = x2_grid, z = matrix(plot_grid$D, nrow = 50),
  xlab = "Factor X1", ylab = "Factor X2", zlab = "Overall Desirability",
  col = "lightblue", alpha = 0.7
)
# Add contour lines for better readability
contour3d(
  x = x1_grid, y = x2_grid, z = matrix(plot_grid$D, nrow = 50),
  levels = seq(0, 1, 0.1), col = "gray", add = TRUE
)
title3d(main = "3D Surface of Overall Desirability")

Option 2: Publication-Friendly Interactive Plot with plotly

# Use plotly for a sharper, web-friendly plot
plot_ly(plot_grid, x = ~x1, y = ~x2, z = ~D) %>%
  add_surface(colorscale = "Viridis") %>%
  layout(
    scene = list(
      xaxis = list(title = "Factor X1"),
      yaxis = list(title = "Factor X2"),
      zaxis = list(title = "Overall Desirability")
    ),
    title = "3D Surface of Overall Desirability"
  )

Key Notes to Adapt to Your Data: #

• Replace the simulated grid$y1 and grid$y2 with your actual response data.
• Adjust the desirability functions (dMax, dMin, dTarget) and their bounds to match your optimization objectives.
• If your model isn't a quadratic response surface, modify the rsm formula (e.g., use FO for first-order) or fit a standard lm model instead.
• Tweak the factor ranges in x1/x2 to match your actual experimental factors.

Once you run this, you'll get an interactive 3D plot where you can rotate/zoom to find the regions with the highest overall desirability (closest to 1)—those are your optimal factor combinations!

内容的提问来源于stack exchange，提问作者John