Как сгенерировать данные о заданном значении коэффициента множественного определения в R?

Question

Мне нужно сгенерировать данные о заданном значении коэффициента множественного определения. Например, если я указал R^2 = 0,77, я хочу сгенерировать данные, которые создают модель регрессии с R^2 = 0,77

но эти данные должны быть в определенном диапазоне. Например, sample= 100 и мне нужно 4 переменные (x1 - зависимая переменная), где значения находятся в диапазоне от 5-15. Как это сделать? Я использую Optim

optim(0.77, fn, gr = NULL,
      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",
                 "Brent"),
      lower = 5, upper = 15,
      control = list(), hessian = FALSE)

но я не знаю, как создать функцию fn для моей цели. Пожалуйста, помогите написать эту функцию

Answer 1

Сначала вот решение:

library(mvtnorm)

get.r <-  function(x) c((x+sqrt(x**2+3*x))/(3),(x-sqrt(x**2+3*x))/(3))

set.seed(123)
cv <- get.r(0.77)[1]
out <- rmvnorm(100,sigma=matrix(c(1,cv,cv,cv,cv,1,cv,cv,cv,cv,1,cv,cv,cv,cv,1),ncol=4))
out1 <- as.data.frame(10*(out-min(out))/diff(range(out))+5)

range(out1)
# [1]  5 15

lm1 <- lm(V1~V2+V3+V4,data=out1)

summary(lm1)
# Call:
#   lm(formula = V1 ~ V2 + V3 + V4, data = out1)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -1.75179 -0.64323 -0.03397  0.64770  2.23142 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.36180    0.50940   0.710 0.479265    
# V2           0.29557    0.09311   3.175 0.002017 ** 
# V3           0.31433    0.08814   3.566 0.000567 ***
# V4           0.35438    0.07581   4.674 9.62e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.927 on 96 degrees of freedom
# Multiple R-squared:  0.7695,  Adjusted R-squared:  0.7623 
# F-statistic: 106.8 on 3 and 96 DF,  p-value: < 2.2e-16

Теперь позвольте мне объяснить, как я туда попал. Мы можем построить это статистически. Для начала нам нужно немного понять корреляцию и ковариацию. Одна формула для корреляции

Corr(X, Y) = Cov(X,Y)/sqrt(Var(X)Var(Y))

И одна формула для ковариации это:

Cov (X, Y) = E (XY) - E (X) E (Y)

В своем вопросе вы хотите получить множественную корреляцию регрессионной модели:

Y = X1 + X2 + X3

Давайте сделаем это максимально простым и сделаем так, чтобы дисперсия всех переменных была равна 1, и давайте сделаем попарную корреляцию между любыми двумя переменными равной и назовем ее r.

Теперь мы ищем квадрат корреляции между Y и X1 + X2 + X3, который:

R^2 = [Cov (Y, X1 + X2 + X3)] ^ 2 / [Var (Y) Var (X1 + X2 + X3)]

Обратите внимание, что

Cov(Y,X1 + X2 + X3) = Cov(Y,X1) + Cov(Y,X2) + Cov(Y,X3)

Также отметим, что дисперсия каждой переменной равна 1, а попарная корреляция равна r, поэтому приведенный выше результат эквивалентен 3r.

Также обратите внимание, что

Var(X1 + X2 + X3) = Var(X1) + Var(X2) + Var(X3) + Cov(X1,X2) + Cov(X1,X3) + Cov(X2,X3).

Поскольку дисперсия каждого равна 1, это эквивалентно 3 + 6r, поэтому

R^2 = 9r^2/(3 + 6r) = 3r^2/(1 + 2r)

Мы можем использовать квадратное уравнение, чтобы решить для г и получить

r = (R^2 +/- sqrt((R^2)^2+3R^2))/3

Если подставить R^2 = 0,77, то r = -0,3112633 или 0,8245966. Мы можем использовать любой, чтобы получить то, что вам нужно, используя rmvnorm() в пределах mvtnorm пакет. А поскольку R^2 инвариантен к линейным преобразованиям, мы можем преобразовать результирующие переменные так, чтобы они попадали между 5 и 15.

Обновить:

Если мы хотим смоделировать с n Предикторы, мы можем использовать следующее (обратите внимание, что я не преобразую диапазон каждого предиктора, но это можно сделать после факта, не изменяя кратное число R^2):

get.r <- function(x,n) c(((n-1)*x+sqrt(((n-1)*x)**2+4*n*x))/(2*n),
                         ((n-1)*x-sqrt(((n-1)*x)**2+4*n*x))/(2*n))

sim.data <- function(R2, n) {
  sig.mat <- matrix(get.r(R2,n+1)[1],n+1,n+1)
  diag(sig.mat) <- 1

  out <- as.data.frame(rmvnorm(100,sigma=sig.mat))

  return(out)
}

Answer 2

Это не ответ, но я хотел поделиться тем, что я сделал. Я не верю optim можно использовать так, как вы хотите. Я попытался использовать метод "грубой силы", чтобы найти набор данных, который мог бы работать, но самый высокий r-квадрат, который я "получил", был 0,23:

# Initializing our boolean and counter.
rm(list = ls())
Done <- FALSE
count <- 1
maxr2 <- .000001

# I set y ahead of time.
y <- sample(5:15, 100, replace = TRUE)

# Running until an appropriate r-squared is found.
while(!Done) {

  # Generating a sample data set to optimize y on.
  a <- sample(5:15, 100, replace = TRUE)
  b <- sample(5:15, 100, replace = TRUE)
  c <- sample(5:15, 100, replace = TRUE)
  data <- data.frame(y = y, a = a, b = b, c = c)

  # Making our equation and making a linear model.
  EQ <- "y ~ a + b + c" # Creating the equation.
  model <- lm(EQ, data) # Running the model.
  if (count != 1) { if (summary(model)$r.squared > maxr2) { maxr2 <- summary(model)$r.squared } }
  r2 <- summary(model)$r.squared # Grabbing the r-squared.
  print(r2) # Printing r-squared out to see what is popping out.
  if (r2 <= 0.78 & r2 >= 0.76) { Done <- TRUE } # If the r-squared is satfisfactory, pop it out.
  count <- count + 1 # Incrementing our counter.
  if (count >= 1000000) { Done <- TRUE ; print("A satisfactory r-squared was not found.") } # Setting this to run at most 1,000,000 times.

}

# Data will be your model that has an r-squared of 0.77 if you found one.

Вопрос с optim заключается в том, что он оптимизирует отдельные параметры, отдельные значения. Первый аргумент в optim это par аргумент, который должен быть списком значений, которые вы хотите оптимизировать. Это может быть использовано при оптимизации r-квадрата с помощью некоторой функции затухания, которая зависит от нескольких значений (это будет ваш par ценности). Однако в этом случае вы просите оптимизировать целые столбцы для максимизации r-квадрата, что не имеет смысла (насколько я знаю) с optim,