R: Ошибки в cbind() с функцией kclass() и отсутствующие значения в кадре данных

Question

У меня проблема с использованием функции kclass() пакета RCompAngrist, когда у меня отсутствуют значения в моем df. Это функция, которая должна вычислять оценку «ограниченного объема информации максимального правдоподобия» с зависимой переменной с левой стороны и двумя частями в правой части уравнения. Первый для эндогенных переменных и второй для инструментов. Он основан на функции ivreg() пакета AER. Ниже приведен минимальный рабочий пример, который будет воспроизводить ошибку.

library(magrittr) 
library(devtools) 
install_github(repo = "RCompAngrist", 
         username = "MatthieuStigler", 
         subdir = "RcompAngrist") 
library(RcompAngrist) 
a <- runif(10, 5, 90) 
b <- runif(10, 4, 10) 
c <- runif(10, 0, 1) 
d <- runif(10, 5, 65) 
e <- runif(10, 1, 2) 
f <- runif(10, 1, 100) 
g <- runif(10, 80, 90) 
h <- c(1,12,3,5,NA,16,17,NA,9,10) 
dummy <- kclass(a ~ b + c + d | d + e + f + g + h, 
model = T, 
data=df) 

Если запустить этот код, вы должны получить сообщение об ошибке от R:

Ошибка в cbind (x_exo, г, x_endo, у): число строк матриц должны совпадать (см arg 2)

Это связано с NA в кадре данных, но я не могу понять, что именно происходит неправильно. Он работает, если вы создаете переменную «h» без NA. Однако, если опустить Nas через

df <- data.frame(a,b,c,d,e,f,g,h) %>% na.omit() 

прежде чем переоценку модели, R дает мне сообщение об ошибке:

Ошибка в R_Z [с (n_G, n_y), C (n_G, n_y)]: индекс за пределами

Я также не понимаю, почему он не опускает НС самостоятельно, поскольку глобальная опция для na.action - это na.omit. Это становится еще более странным. При удалении «данных = ДФ» из функции, а затем повторно запустить модель сообщение об ошибке переключается на

Ошибка в cbind [...]

Почему это делает никакой разницы здесь, если «data = df» находится в коде или нет? Есть ли у кого-нибудь идеи, в которых может возникнуть проблема? Я не понимаю, что здесь происходит.

Answer 1

Похоже, в вашем примере есть несколько ошибок. Вы должны определить переменные в data.frame для аргумента data=df для работы. Я сделал некоторые корректировки ниже:

df <- data.frame(
    a = runif(10, 5, 90), 
    b = runif(10, 4, 10), 
    c = runif(10, 0, 1), 
    d = runif(10, 5, 65), 
    e = runif(10, 1, 2), 
    f = runif(10, 1, 100), 
    g = runif(10, 80, 90), 
    h = c(1,12,3,5,NA,16,17,NA,9,10) 
) 

Вызов RcompAngrist дает мне:

m_rca <- RcompAngrist::kclass(a ~ b + c + d | d + e + f + g + h, 
           data = df) 
#>Error in cbind(x_exo, z, x_endo, y) : 
#>number of rows of matrices must match (see arg 2) 

m_rca <- RcompAngrist::kclass(a ~ b + c + d | d + e + f + g + h, 
          data = na.omit(df)) 
#> Error in R_Z[c(n_G, n_y), c(n_G, n_y)]: subscript out of bounds 

В первом случае, это выглядит как NA, исчезают из h, но не от другого строки, которые дадут вам ошибку number of rows must match. Во-вторых, я не уверен, что происходит.

Есть несколько альтернатив. Лучше всего, если вы не требуете LIML, было бы использовать AER::ivreg, что даст вам двухэтапную оценку наименьших квадратов.

Вы также можете использовать rkclass пакет (отказ от ответственности: я написал, и репо Стиглера была огромная помощь в то время как я делал это), что согласуется с результатами TSLS из AER::ivreg, и предлагает вариант ММПОИ. Он находится в состоянии потока, хотя и не подвергся серьезной проверке, поэтому лучше было бы подтвердить результаты.

m_aer <- AER::ivreg(a ~ b + c + d | d + e + f + g + h, 
        data = df) 

m_rkc <- rkclass::kclass(a ~ b + c + d | d + e + f + g + h, 
         model.type = "TSLS", 
         data = df) 

summary(m_aer) 
#> 
#> Call: 
#> AER::ivreg(formula = a ~ b + c + d | d + e + f + g + h, data = df) 
#> 
#> Residuals: 
#>  Min  1Q Median  3Q  Max 
#> -30.296 -16.990 5.157 11.238 28.988 
#> 
#> Coefficients: 
#>    Estimate Std. Error t value Pr(>|t|) 
#> (Intercept) 64.7730 87.3657 0.741 0.500 
#> b    1.3600  9.0175 0.151 0.887 
#> c   -7.8303 52.9975 -0.148 0.890 
#> d   -0.5200  0.5537 -0.939 0.401 
#> 
#> Residual standard error: 27.63 on 4 degrees of freedom 
#> Multiple R-Squared: 0.2457, Adjusted R-squared: -0.3201 
#> Wald test: 0.43 on 3 and 4 DF, p-value: 0.743 


summary(m_rkc) 
#> Call: 
#> rkclass::kclass(formula = a ~ b + c + d | d + e + f + g + h, 
#>  data = df, model.type = "TSLS") 
#> 
#> Residuals: 
#> Min  1Q Median  3Q Max 
#> -30.30 -16.99 5.16 11.24 28.99 
#> 
#> Coefficients: 
#>    Estimate Std. Error t value Pr(>|t|) 
#> (Intercept) 64.773  87.366 0.74  0.50 
#> b    1.360  9.018 0.15  0.89 
#> c    -7.830  52.997 -0.15  0.89 
#> d    -0.520  0.554 -0.94  0.40 

Результаты rkclass::kclass с LIML указанной дико разница, так что это, безусловно, стоит проверить, используя что-то вроде Stata:

m_liml <- rkclass::kclass(a ~ b + c + d | d + e + f + g + h, 
        model.type = "LIML", 
        data = df) 
summary(m_liml) 
#>Call: 
#>rkclass::kclass(formula = a ~ b + c + d | d + e + f + g + h, 
#>    data = df, model.type = "LIML") 
#> 
#>Residuals: 
#> Min  1Q Median  3Q Max 
#>-19.45 -13.84 -1.77 10.82 25.89 
#> 
#>Coefficients: 
#>   Estimate Std. Error t value Pr(>|t|) 
#>(Intercept) 18.784  36.142 0.52  0.63 
#>b    10.678  11.677 0.91  0.41 
#>c   -13.035  52.957 -0.25  0.82 
#>d    -1.101  0.911 -1.21  0.29