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

У меня проблема с использованием функции 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, z, x_endo, y): число строк матриц должно совпадать (см. аргумент 2)

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

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

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

Ошибка в R_Z[c(n_G, n_y), c(n_G, n_y)] : нижний индекс выходит за пределы

Я также не понимаю, почему он не опускает NA сам по себе, поскольку глобальная опция для na.action — na.omit. Хотя бывает еще страннее. Если вы удалите «data=df» из функции, а затем повторно запустите модель, сообщение об ошибке вернется к

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

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


r nas cbind
person avocado1    schedule 04.10.2015    source источник
comment
попробуйте запустить с в 10 раз больше данных. Я думаю, что вы нарушаете предположения алгоритма.   -  person pcantalupo    schedule 04.10.2015

Ответы (1)


arrow_upward
0
arrow_downward

Похоже, в вашем примере есть пара ошибок. Вы должны определить переменные в 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 и предлагает вариант LIML. Однако он находится в постоянном движении и не подвергался серьезным испытаниям, поэтому было бы лучше подтвердить результаты.

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
person potterzot    schedule 04.10.2015
comment
Да, я сделал из него фрейм данных, не знаю, почему его нет в первой части кода Вот как это было: df <- data.frame(a,b,c,d,e,f,g,h). Я тоже так думал, он опускает NA из h, но не из других переменных. Есть идеи, почему это так? Я также попробую ваш пакет и сообщу, как только получу. - person avocado1; 04.10.2015