How do I add a conditional column based on an inequality with matrix multiplication (in R)?-CodePudding

I have a data frame of values and want to add a column based on an inequality condition that involves matrix multiplication.

The data frame looks like this

# Set possible values for variables
b0 <- b1 <- b2 <- seq(0, 2, by=0.1)
# Create data frame of all the different combos of these variables
df <- setNames(expand.grid(b0, b1, b2), c("b0", "b1", "b2"))

There are a lot of precursor objects I have to define before adding this column:

##### Set n
n = 100

#### Generate (x1i, x2i)
# Install and load the 'MASS' package
#install.packages("MASS")
library("MASS")
# Input univariate parameters
rho <- 0.5
mu1 <- 0; s1 <- 1
mu2 <- 0; s2 <- 1
# Generate parameters for bivariate normal distribution
mu <- c(mu1, mu2) 
sigma <- matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2), nrow=2, ncol=2) 
# Generate draws from bivariate normal distribution
bvn <- mvrnorm(n, mu=mu, Sigma=sigma ) # from MASS package
x1 <- bvn[, 1]
x2 <- bvn[, 2]
##### Generate error
error <- rnorm(n)

##### Generate dependent variable 
y <- 0.5   x1   x2   error

##### Create the model
lm <- lm(y ~ x1   x2)


# Setup parameters
n <- 100
K <- 3
c <- qf(.95, K, n - K)

# Define necessary objects
sigma_hat_sq <- 1
b0_hat <- summary(lm)$coefficients[1, 1]
b1_hat <- summary(lm)$coefficients[2, 1]
b2_hat <- summary(lm)$coefficients[3, 1]
x <- cbind(1, x1, x2)

I am trying to add this conditional column like this:

# Add a column to the data frame that says whether the condition holds
df <- transform(df, ueq = (
  (1/(K*sigma_hat_sq))*
    t(matrix(c(b0_hat-b0, b1_hat-b1, b2_hat-b2)))%*%
    t(x)%*%x%*%
    matrix(c(b0_hat-b0, b1_hat-b1, b2_hat-b2))
    <= c 
))

...but doing so generates the error message

Error in t(matrix(c(b0_hat - b0, b1_hat - b1, b2_hat - b2))) %*% t(x) : 
  non-conformable arguments

Mathematically, the condition is [1/(Ksigmahat^2)](Bhat-B)'X'X(Bhat-B) <= c, where, for each triple (b0,b1,b2), (Bhat-B) is a 3x1 matrix with elements {B0hat, B1hat, B2hat}. I'm just not sure how to write this condition in R.

Any help would be greatly appreciated!

CodePudding user response：

In order to only work with one row of df at a time (and get a separate answer for each 1 x 3 matrix, you need a loop. A simple way to do this in R is mapply.

df <- transform(df, ueq = mapply(\(b0, b1, b2)
  (1/(K*sigma_hat_sq)) *
    t(c(b0_hat-b0, b1_hat-b1, b2_hat-b2)) %*%
    t(x) %*% x %*%
    c(b0_hat-b0, b1_hat-b1, b2_hat-b2)
  <= c,
  b0 = b0, b1 = b1, b2 = b2
))

This leads to 91 TRUE rows.

sum(df$ueq)
[1] 91