Как провести реверс-инжиниринг набора данных?
Мой учитель дал мне результат регрессии, и упражнение состоит в том, чтобы перепроектировать набор данных, который привел к этой регрессии. Затем нам нужно провести регрессию и получить точно такие же результаты.
Коэффициент мне удалось довести до -19,93, но тот же SE у меня совсем не получается. Я не знаю, как это сделать. Я не знаю, следует ли мне использовать какие-то формулы, связывающие SE оценщика и стандартную ошибку регрессии (они у меня есть, но я не вижу способа реализовать их в R)... Заранее спасибо за ваша помощь!
Мой вывод R:
## Given values
n <- 1592
se_β1 <- 1.47
β1hat <- -19.93
## Create a dummy variable for control vs treatment condition
set.seed(123)
Low_anchor <- rbinom(n,1,0.5)
## Formula of standard error of beta 1 (assuming homoskedasticity)
calculate_standard_error <- function(u, Low_anchor) {
sqrt((1/(n - 2))*sum(u^2)/(n*sd(Low_anchor)^2))
}
## Define initial values of u
u <- rnorm(n)
## Tolerance for convergence
tolerance <- 0.1
## Iteratively adjust u until the standard error matches the target
while (abs(calculate_standard_error(u, Low_anchor) - se_β1) > tolerance) {
## Generate new set of values for u from a normal distribution
u <- rnorm(n)
}
print(u)
## regression
Yc <- -19.93*Low_anchor + u
model1 <- lm(Yc ~ Low_anchor - 1)
## Print the summary of the model
summary(model1)
Что такое Nc выше? Он нигде не появляется в наборе?
Ответы 1
Похоже, Nc не определен. Я думаю, вы вызываете независимую переменную. Здесь я буду использовать x. Обратите внимание, что этот вопрос, похоже, просит вас использовать свойства стандартной ошибки и коэффициента в случае (множественной) регрессии в форме y = bx + u. Ты должен это знать
- Если вы умножите x на некоторое число a, расчетный коэффициент будет b/a.
- Если вы умножите u на некоторое число a, предполагаемая стандартная ошибка составит a*s.
Благодаря этому вы можете написать простой взгляд, который сначала настраивает u, а затем настраивает x. Сначала мы определим некоторые праймериз:
n<- 1592
se_b1 <- 1.47
b1hat <- -19.93
set.seed(2)
x <- rnorm(n)
y_mean <- -19.93*x
# we are going to create a random variable to be the residuals
u <- rnorm(n,0,1)
error <- 1
tol <- 0.01
Обратите внимание, что их распределение не имеет отношения к ответу. Вы можете убедиться, что среднее значение u далеко от всего, что нам хотелось бы. Вы также можете проверить, что произойдет после запуска y <- y_mean + u резюме (lm (y ~ x)) Коэффициент и стандартная ошибка будут отличаться от желаемых. Как это исправить? Используя два свойства, которые мы упомянули выше.
error <- 1
tol <- 0.01
while (error > tol) {
# remember that the se_b1 is constructed as the rood of the diagonal of sigma^2 * (X'X)^-1
# determine the matrix of X (assuming there is an intercept here)
X <- matrix(c(rep(1, n), x), ncol = 2)
XX_minus_one <- solve(t(X) %*% X)
# so far, we would get a standard deviation os
present_se <- sqrt(var(u) * XX_minus_one[2,2])
# this is different. Let's adjust the residuals to have the desired variance
u_fitting <- u * se_b1 / present_se
y <- y_mean + u
reg <- lm(y ~ x)
estimated_b1 <- reg$coefficients[2]
estimated_se_b1 <- summary(reg)$coefficients[2,2]
error <- max(abs(estimated_se_b1 - se_b1), abs(estimated_b1 - b1hat))
# but now we need to refit the x
x <- x * estimated_b1/b1hat
u <- u_fitting
}
Вы можете проверить, работает ли он:
summary(lm(y ~ x))
Возможно, это сложнее, чем вам нужно, но по теме и интересно: Research.autodesk.com/app/uploads/2023/03/…