Почему предикат lm возвращает NA как прогноз в R? - PullRequest
Новая
       134

Почему предикат lm возвращает NA как прогноз в R?

0 голосов
/ 13 апреля 2020

Я пытаюсь воспроизвести симуляцию из статьи для проекта класса. Я создал симулированную матрицу х и симулированный вектор у. Я хочу сначала уменьшить размер модели с p = 1000 до p = n / log (n), используя SIS. Затем я хочу еще больше уменьшить векторное пространство до p ', используя селектор Primal Dual Dantzig.

В этом пакете нет функции прямого предсказания, поэтому я подключаю p 'выбранные предикторы и их коэффициенты в функцию lm для получения предсказаний. Я использую функцию Предсказание для получения предсказанных значений. Тем не менее, возвращаемый вектор - все NA.

Следующий код может дополнительно проиллюстрировать проблему. Я не получаю никаких ошибок после запуска следующего кода. Я должен следовать этой процедуре использования SIS, а затем Dantzig Selector, чтобы соответствовать симуляции в статье. 

library(SIS) #import the SIS library
library(fastclime)

errors_DS_small = c() #create empty array to store errors from each simulation
model_sizes_DS_small = c() #creat empty array to store selected model size from each simulation
start_time = Sys.time()

  set.seed(123*1) #re-create results at a later time but different seed for each data set
  n = 200 #sample size
  p = 1000 #variables
  x = matrix(rnorm(n*p, mean=0, sd=1), n, p) #creating IID standard Gaussian random predictors

  # gaussian response
  set.seed(456*1) #re-create results at a later time but different seed for each data set
  s = 8
  u = rbinom(s, 1, 0.4)
  z = rnorm(s, mean=0, sd=1)
  a = 4*log(n)/sqrt(n)
  b= ((-1)**(u))*(a + abs(z))
  y=x[, 1:s]%*%b 

  #creating SIS-DS model and gaussian response. iter=FALSE means not doing ISIS
  modelSIS_small = SIS(x, y, family='gaussian', iter = FALSE, nsis=(n/log(n)))

  modelDS = dantzig(x[,modelSIS_small$sis.ix0], y)
  selectDS = dantzig.selector(modelDS$lambdalist, modelDS$BETA0, lambda=min(modelDS$lambdalist))

  linearMod = lm(y~x[,modelSIS_small$sis.ix0])
  linearMod$coefficients = selectDS

  newx = x[,modelSIS_small$sis.ix0]

  predTest = predict(linearMod, data=newx) #create predictions using test data test
  a = modelDS$BETA0[, 9] != 0
  mse = mean((y - predTest)^2) #compare predictions to real values of test y
  rmse = sqrt(mse) #Square root of MSE (Square-loss function)
  errors_DS_small[1] = rmse #store the RMSE

  model_sizes_DS_small[1] = table(a)["TRUE"] #Store model size including intercept
  print(c(1, errors_DS_small[1], model_sizes_DS_small[1])) #print results

  end_time = Sys.time()
  DS_total_time_small = end_time - start_time
  DS_error_small_median = median(errors_DS_small)
  DS_model_Size_small_median = median(model_sizes_DS_small)

  print(c("DS", DS_model_Size_small_median, DS_error_small_median, DS_total_time_small))
...