set.seed(123)
dt <- 
  data.frame(
    x1 = (1:10)+rnorm(10),
    x2 = (11:20)+rnorm(10),
    x3 = (10:1)+rnorm(10),
    x1r = (10:1)+rnorm(10),
    x2r = (20:11)+rnorm(10),
    x3r = (1:10)+rnorm(10),
    y1 = (1:10)+rnorm(10),
    y2 = (1:10)+rnorm(10),
    y3 = (1:10)+rnorm(10)
  )
pairs(dt)
cat("\n\n=== Correlation coefficients ===\n")
cor(dt)

dt.x <- dt[, 1:3]
dt.xr <- dt[, 4:6]
dt.y <- dt[, 7:9]
# remind that dim(dt.x) == dim(dt.xr)

betas <- matrix(NA_real_, 3, 3, dimnames = list(paste0("x", 1:ncol(dt.x)), paste0("y", 1:ncol(dt.y))))
is.reverse <- matrix(NA, 3, 3, dimnames = list(paste0("x", 1:ncol(dt.x)), paste0("y", 1:ncol(dt.y))))
for(i in 1:ncol(dt.y)) {
  for(j in 1:ncol(dt.x)) {
    beta <- coef(lm(dt.y[, i] ~ dt.x[, j]))[2]
    if(beta < 0) {
      betas[i, j] <- coef(lm(dt.y[, i] ~ dt.xr[, j]))[2]
      is.reverse[i, j] <- T
    } else {
      betas[i, j] <- beta
      is.reverse[i, j] <- F
    }
  }
}


cat("\n\n=== 'Adjusted' coefficients ===\n")
betas
cat("\n\n=== Sum of 'adjusted' coefficients by y ===\n")
colSums(betas)
cat("\n\n=== Is 'adjusted'? ===\n")
is.reverse