Метод извлечения коэффициента с использованием пакета lme4 в R: большой набор данных - PullRequest
Новая
       19

Метод извлечения коэффициента с использованием пакета lme4 в R: большой набор данных

1 голос
/ 30 марта 2020

Я использую пакет lme4 для функциональной линейной модели со смешанным эффектом. Процедура состоит из создания базы сплайнов и последующего использования их в качестве переменных в модели.

Набор данных, который я использую, состоит из наблюдения 3654090, а моя модель выглядит следующим образом:

Параметр ~ 0 + X1 + X2 + X3 + ZONE + (0 + X1 | ANNI) + (0 + X2 | ANNI) + (0 + X3 | ANNI) + (0 + X1 | ID) + (0 + X2 | ID) + (0 + X3 | ID)

, где X1, X2, X3 - 3 сплайна, а ZONE, ANNI, ID - категориальные переменные. ZONE и ANNI имеют 5 уровней, ID вместо 89942 уровней. Модель сходится, но когда я пытаюсь экстраполировать коэффициенты, у меня появляется следующее сообщение об ошибке:

Ошибка в .local (a, b, ...): ошибка Cholmod «проблема слишком большая» в файле ../ Core / cholmod_memory. c, строка 334

Я понял, что проблема в матрице случайных эффектов (269841 x 3654090).

Мои вопросы:

  1. Есть ли способ извлечь коэффициенты из очень большого набора данных?
  2. Если нет, то какой максимальный размер матрицы я могу использовать, чтобы избежать этой проблемы?

I ' m с использованием R 3.6.3 и RStudio 1.2.5033 на компьютере Windows с 8 ГБ ОЗУ. Я также пытался запустить свой код на Linux виртуальной машине с 64 ГБ ОЗУ, но это не решило проблему.

Я создал для вас также пример: 

library(lme4)
library(splines)
library(tidyverse)

memory.limit(100000)

parabola<-function(a,b,c,x){
  a*x^2+b*x+c
}

#time interval
x <- seq(2,28,length.out = 100)

#I extract only 5 value for every ID (15000)
t<-matrix(0,nrow = 5,ncol = 15000)
set.seed(123)
for(i in 1:15000){
t[,i]=sample(x,5)
}

#creation of first case
y1<-matrix(0,nrow = 5,ncol = 15000)
set.seed(123)
for(i in 1:8000){
  a<-sample(seq(-0.01,0.01,length.out = 50),1)
  b<-sample(seq(-0.50,0.50,length.out = 50),1)
  c<-sample(seq(0,1,length.out=50),1)
  for(j in 1:5){
    y1[j,i]=parabola(a,b,c,t[j,i])
  }
}

#creation of second case
y2<-matrix(0,nrow = 5,ncol = 15000)
set.seed(123)
for(i in 1:8000){
  a<-sample(seq(-0.01,0.01,length.out = 50),1)
  b<-sample(seq(-0.50,0.50,length.out = 50),1)
  c<-sample(seq(1.5,2.5,length.out=50),1)
  for(j in 1:5){
    y2[j,i]=parabola(a,b,c,t[j,i])
  }
}

#creation of third case
y3<-matrix(0,nrow = 5,ncol = 15000)
set.seed(123)
for(i in 1:8000){
  a<-sample(seq(-0.01,0.01,length.out = 50),1)
  b<-sample(seq(-0.50,0.50,length.out = 50),1)
  c<-sample(seq(3,4,length.out=50),1)
  for(j in 1:5){
    y3[j,i]=parabola(a,b,c,t[j,i])
  }
}

#creation of forth case
y4<-matrix(0,nrow = 5,ncol = 15000)
set.seed(123)
for(i in 1:8000){
  a<-sample(seq(-0.01,0.01,length.out = 50),1)
  b<-sample(seq(-0.50,0.50,length.out = 50),1)
  c<-sample(seq(3.5,4.5,length.out=50),1)
  for(j in 1:5){
    y4[j,i]=parabola(a,b,c,t[j,i])
  }
}

# creation of the dataset
t<- t %>% as.data.frame() %>%  gather(value="time")
y1 <- y1 %>%
  as.data.frame() %>%
  gather(value="y") %>%
  add_column(case=rep('A',75000))

y2 <- y2 %>%
  as.data.frame() %>%
  gather(value="y") %>% 
  add_column(case=rep('B',75000))

y3 <- y3 %>%
  as.data.frame() %>%
  gather(value="y") %>% 
  add_column(case=rep('C',75000))

y4 <- y4 %>%
  as.data.frame() %>%
  gather(value="y") %>% 
  add_column(case=rep('D',75000))


# final dataset
X<-rbind(cbind(t,y1)[,-3],cbind(t,y2)[,-3],cbind(t,y3)[,-3],cbind(t,y4)[,-3])

# model
k=3
d=2

mean_t <- X$time #times
sort_t_list  <- sort(mean_t ,index.return=TRUE)
sort_t<- sort_t_list [[1]] # ordered times
sort_t_ind  <- sort_t_list [[2]] # indeces for permuting the original variable accordingly



spline <- bs(sort_t ,df=k,degree=d,intercept = T)

data <- data.frame(spline,Parameter=X$y[sort_t_ind ],case=as.factor(X$case[sort_t_ind]),ID=X$key[sort_t_ind]) # create a new dataframe for the regression model of the first variable

set.seed(1993)
fit_k3d2 <- lmer(Parameter~0+X1+X2+X3+(0+X1|case)+(0+X2|case)+(0+X3|case)+(0+X1|ID)+(0+X2|ID)+(0+X3|ID),
                 data=data,
                 control=lmerControl(optCtrl=list(xtol_abs=1e-10, ftol_abs=1e-10)),
                 na.action=na.exclude,
                 verbose = T)

# this command gives the error
coef_id<-coef(fit_k3d2)$ID
...