在R语言中绘制基于x、z的倾向得分三维曲面图的最优方法

Got it, let's sort out this 3D plot problem! The issue with scatter3d() is that 10,000 data points are way too dense—they overlap so much you can't make out the actual trend of propensity scores relative to x and z. Instead of plotting every single point, we'll focus on visualizing the underlying surface of your data, either by fitting a smooth model or binning the data to show averages. Here are two robust methods:

Method 1: Fit a Smooth Surface with GAM + Interactive Plotly #

This approach uses a Generalized Additive Model (GAM) to capture non-linear relationships between x, z, and your propensity score, then plots the predicted surface interactively with Plotly (great for exploring the data by rotating/zooming).

First, load the necessary packages and your data:

library(tidyverse)
library(mgcv)
library(plotly)

# Load your data (make sure mydata2.csv is in your working directory)
mydata <- read_csv("mydata2.csv")

Next, fit the GAM model to create a smooth surface:

# Fit a GAM with smooth terms for x, z, and their interaction
gam_model <- gam(propensity_score ~ s(x) + s(z) + ti(x, z), data = mydata)

# Create a grid of x and z values to predict the surface
x_grid <- seq(min(mydata$x), max(mydata$x), length.out = 100)
z_grid <- seq(min(mydata$z), max(mydata$z), length.out = 100)
surface_data <- expand.grid(x = x_grid, z = z_grid)

# Predict propensity scores across the grid
surface_data$pred_propensity <- predict(gam_model, newdata = surface_data)

Now plot the interactive 3D surface:

plot_ly(data = surface_data, x = ~x, y = ~z, z = ~pred_propensity, type = "surface") %>%
  layout(
    scene = list(
      xaxis = list(title = "Variable x"),
      yaxis = list(title = "Variable z"),
      zaxis = list(title = "Propensity Score")
    ),
    title = "Smooth Surface of Propensity Scores vs x & z"
  )

You can rotate, zoom, and pan this plot to explore how propensity scores change with x and z. If you want to add a subset of raw points to ground the surface in real data, add this line after the plot:

# Add 1,000 random raw points (avoids overcrowding)
sample_points <- slice_sample(mydata, n = 1000)
add_markers(p = ., data = sample_points, x = ~x, y = ~z, z = ~propensity_score, 
            color = I("red"), size = I(2))

Method 2: Binned Averages with RGL (Static/Interactive 3D) #

If you prefer not to fit a model, you can bin x and z into intervals, calculate average propensity scores for each bin, then plot the binned surface with RGL.

library(tidyverse)
library(rgl)

# Load data
mydata <- read_csv("mydata2.csv")

# Bin x and z into 50 intervals each, calculate average propensity score per bin
binned_data <- mydata %>%
  mutate(
    x_bin = cut(x, breaks = 50),
    z_bin = cut(z, breaks = 50)
  ) %>%
  group_by(x_bin, z_bin) %>%
  summarise(
    avg_propensity = mean(propensity_score),
    x_center = mean(x),
    z_center = mean(z)
  ) %>%
  ungroup()

# Convert binned data to a matrix for RGL's persp3d
propensity_matrix <- matrix(binned_data$avg_propensity, nrow = 50)
x_centers <- unique(binned_data$x_center)
z_centers <- unique(binned_data$z_center)

# Plot the binned surface
persp3d(
  x_centers, z_centers, propensity_matrix,
  xlab = "x", ylab = "z", zlab = "Average Propensity Score",
  col = "lightblue", alpha = 0.7, main = "Binned Propensity Scores vs x & z"
)

# Optional: Add a subset of raw points
sample_points <- slice_sample(mydata, n = 1000)
points3d(sample_points$x, sample_points$z, sample_points$propensity_score, 
         col = "darkred", size = 1)

Why This Works Better Than scatter3d() #

10,000 points create a "cloud" that hides the actual trend. By focusing on a smooth surface or binned averages, you're highlighting the relationship between x, z, and propensity scores instead of just showing every data point. Both methods let you explore the data clearly, whether you want an interactive Plotly plot or a rotatable RGL surface.

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