手动计算泊松回归饱和模型系数遇矩阵不匹配错误求助

Hey there! Let's work through this saturated Poisson model dimension mismatch error you're running into—super common pitfall with these manual implementations, so let's break it down.

First, Why the Error Happens #

The non-conformable arguments error tells us your matrix multiplication in the log-likelihood function is using incompatible dimensions. For saturated Poisson models, this almost always boils down to two issues:

1. Backwards matrix multiplication order in your linear prediction calculation
2. Mismatched length between your initial beta values and the saturated model's X matrix

Let's Fix It Step by Step #

1. Correct the Matrix Multiplication Order

Poisson regression's linear prediction for each observation is $\eta_i = X_i\beta$, where:

• $X$ is an $n \times p$ matrix (n observations, p predictors)
• $\beta$ is a $p \times 1$ vector of coefficients

In R, this translates to x %*% beta (not beta %*% t(x)). The latter will throw a dimension error because a $p \times 1$ vector can't multiply a $p \times n$ transposed matrix (the inner dimensions don't match).

Here's the right way to structure your negative log-likelihood function (we use negative because optim minimizes by default):

poisson_neg_loglik <- function(beta, x, y) {
  # Correct linear prediction: n x 1 vector matching y's length
  eta <- x %*% beta
  lambda <- exp(eta)
  # Calculate negative log-likelihood (ignore log(y!) since it's constant for optimization)
  -sum(y * log(lambda) - lambda)
}

2. Build the Saturated Model's X Matrix Properly

A saturated model requires one unique parameter per observation, so your X matrix should be $n \times n$. The easiest way to build this is with an identity matrix (each row has a 1 in its own column, 0 elsewhere):

n <- length(y)
x_saturated <- diag(n)
# Alternatively, using dummy variables:
# x_saturated <- model.matrix(~0 + factor(1:n))

3. Initialize Beta with the Right Length

When calling optim, your initial beta vector must have the same length as the number of columns in X (which is n for the saturated model). Don't reuse the length-1 initial value from your null model—use a vector of zeros (or any reasonable starting values) with length n:

beta0_saturated <- rep(0, n)

Full Working Example #

Let's put this all together with simulated data to confirm:

set.seed(123)
# Simulate data where each observation has its own lambda
n <- 5
y <- rpois(n, lambda = exp(rnorm(n)))

# Build saturated model matrix
x_saturated <- diag(n)

# Define negative log-likelihood
poisson_neg_loglik <- function(beta, x, y) {
  eta <- x %*% beta
  lambda <- exp(eta)
  -sum(y * log(lambda) - lambda)
}

# Optimize the saturated model
fit_saturated <- optim(
  par = rep(0, n),
  fn = poisson_neg_loglik,
  x = x_saturated,
  y = y,
  method = "BFGS"
)

# Check results: estimates should be close to log(y) (since saturated model fits lambda_i = y_i)
cat("Saturated model beta estimates:\n")
print(fit_saturated$par)
cat("\nlog(y) for comparison:\n")
print(log(y))

This should run without dimension errors, and you'll see the beta estimates align closely with log(y)—which makes sense, since the saturated model perfectly fits each observation's mean.

Quick Recap of the Fixes #

• Swap beta %*% t(x) to x %*% beta in your log-likelihood
• Ensure your saturated X matrix has n columns (one per observation)
• Initialize beta with a vector of length n, not 1

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