Jak používat bigstatsr R balíček pomocí dvou souborů dat pro odhad parametrů?

Mám nezávislé a závislé datové sady. Chci vyzkoušet všechny možné vztahy mezi závislé a nezávislé proměnné. V mém předchozím příspěvku (Jak replikovat funkci pomocí mapply s více argumenty pro výpočet výkonu metoda?), Chtěl jsem udělat power analýza pomocí data simulace. Nyní chci analyzovat reálná data pomocí stejné funkce. Problém je, že test_function potřeba více času jako můj datový soubor je velký (rozměr každého souboru dat větší než 10000 X 40000). Také chci použít paralelní výpočty pro urychlení výpočtu. Zjistil jsem, že bigstatsr paket (https://privefl.github.io/bigstatsr/index.html) zvládne matice, které jsou příliš velké, aby se vešly do paměti. Kromě toho, chci, aby se zabránilo rozšíření.grid , jak to je také výpočetně drahé pro big data. Nenašel jsem žádné místo, které může používat dva datové soubory současně pomocí bigstatsr balíček a odhad parametrů paralelně. Datové soubory a příklady kódu jsou uvedeny níže:

# dependent dataset
test_A <- data.frame(matrix(rnorm(100), nr=10, nc=10))
# independent dataset
test_B <- data.frame(matrix(sample(c(0,1,2), 500, replace = TRUE), nr=50, nc=10))
# Find all combination using dependent and independe datasets's variables
A_B_pair <- subset(expand.grid(c1=names(test_A), c2=names(test_B), 
                               stringsAsFactors = FALSE))
# Main function to estimate the parameter and p-values 
test_function <- function(test_A, test_B, x,y){
  c1 <- test_A [[x]]
  c2 <- test_B[[y]]
  Data <- data.frame(1, XX=c1, YY=c2)
  
  model_lm <- lm(YY ~ XX, Data)
  est_lm <- as.numeric(model_lm$coefficients)[2]
  pvalue_lm <- as.numeric(summary(model_lm)$coeffi[,4][2])
  
  return(unlist(data.frame(lm.estimator = est_lm, lm.pvalue =pvalue_lm)))
}
# Final output
output <- mapply(test_function, MoreArgs = list(test_A, test_B),
                 x = A_B_pair$c1, y = A_B_pair$c2)

Edit: Chci použít můj návrh metody pro odhad parametrů a porovnat výsledky s lm metoda. Moje navrhované metody je uveden níže:

library(pracma)
Proposed_method<- function(Data, Beta) 
{ 
  n = dim(Data)[1]
  Median <- t(apply(Data,2,median))
  Dist <- sqrt(rowSums((Data - as.matrix(rep(1,dim(Data)[1]))%*%Median)^2))
  Data0 <- as.matrix(Data[which(Dist <= as.numeric(quantile(Dist, p=.45, na.rm = TRUE))),])
  Yo <- as.matrix(Data0[,dim(Data0)[2]])
  Xo <- as.matrix(Data0[,-dim(Data0)[2]])
  Gama0 <- as.numeric(pinv(crossprod(Xo, Xo))%*%crossprod(Xo, Yo))
  Sigma2o <- var(Yo)
  Y <- as.matrix(Data[,dim(Data)[2]])
  X <- as.matrix(Data[,-dim(Data)[2]])
  
  DiffTol = 0.0001;
  DiffNorm = +10000;
  Iter = 0;
  ###########While loop################
  while (DiffNorm > DiffTol)
  {
    Const <- sqrt(2*pi*Sigma2o)
    devmat <- (Y-X%*%Gama0)
    Squaremat <- as.matrix(apply(devmat, c(1,2), function(x) x^2))
    Gauss <- exp(-Squaremat/(2*as.numeric(Sigma2o)))/as.numeric(Const)
    Wbeta <- exp(-(Beta*((Y-X%*%Gama0)*(Y-X%*%Gama0)))/(2*as.numeric(Sigma2o)))
    ONE1 <- rep(1,dim(X)[2]);
    Xb <- (X*(Wbeta%*%ONE1)) 
    Gama <- as.numeric(pinv(crossprod(X, Xb))%*%crossprod(Xb, Y)) 
    hedprod <- (Y-X%*%Gama)*(Y-X%*%Gama) 
    tWbeta <- as.matrix(t(Wbeta)) 
    One_1 <- as.matrix(rep(1,dim(X)[1])) 
    Sigma2 <- (tWbeta%*%hedprod)*pinv(tWbeta%*%One_1)
    
    LHb<-(sum(Gauss^Beta)/n-1)/Beta
    LH<-prod(Gauss)
    ##########
    Norm2 <- ((sum(Gama*Gama))^0.5 + abs(Sigma2))
    DiffNorm <-((sum((Gama-Gama0)*(Gama-Gama0)))^0.5 + abs(Sigma2 - Sigma2o))/Norm2
    ###
    Gama0 = Gama
    Sigma2o=Sigma2
    Iter = Iter + 1 
  }
  return(list(Gama=Gama,Sigma2=Sigma2,Wt=Wbeta,LHb=LHb,LH=LH))
}
# independent variable dataset
test_A <- data.frame(matrix(sample(c(0,1,2), 500, replace = TRUE), nr=10, nc=50))
# dependent variable dataset
test_B <- data.frame(matrix(rnorm(1000), nr=10, nc=100))
# Find all combination using dependent and independe datasets's variables
A_B_pair <- subset(expand.grid(c1=names(test_A), c2=names(test_B), 
                               stringsAsFactors = FALSE))
# Main function to estimate the parameter and p-values by proposed method and lm 
test_function <- function(x, y){
  c1 <- test_A[[x]]
  c2 <- test_B[[y]]
  Data <- data.frame(1, XX=c1, YY=c2)
  nn <- dim(Data)[1]
  Beta = 0.1
  Omit = 2
  ResL1 <- Proposed_method(Data, Beta)
  ResL0 <- Proposed_method(as.matrix(Data[,-Omit]), Beta)
  LR0 <- (-nn)*log(ResL1$Sigma2/ResL0$Sigma2)
  
  # Proposed estimator
  Proposed_estimator <- (ResL1$Gama)[2]
  Proposed_pvalue <- as.numeric(pchisq(q=LR0, df=1, lower.tail = FALSE))
  
  #lm model
  model_lm <- lm(YY ~ XX, Data)
  est_lm <- as.numeric(model_lm$coefficients)[2]
  pvalue_lm <- as.numeric(summary(model_lm)$coeffi[,4][2])
  
  return(unlist(data.frame(lm.estimator = est_lm, lm.pvalue =pvalue_lm, Proposed_estimator,Proposed_pvalue)))
}

# Output:
output <- mapply(test_function, x = A_B_pair$c1, y = A_B_pair$c2)
# transpose the output
output_t <- data.frame(t(output))

# Final output
output_final <- cbind(A_B_pair, output_t)
output_final <- structure(list(c1 = c("X1", "X2", "X3", "X4", "X5"), c2 = c("X1", 
"X1", "X1", "X1", "X1"), lm.estimator = c(-0.855708052636761, 
0.227250280548332, -0.128955946232531, 0.171650221327542, -0.701027831473379
), lm.pvalue = c(0.0361141129937136, 0.646905371365762, 0.816730073250761, 
0.780290676037238, 0.261013977519426), Proposed_estimator = c(-0.879232513006948, 
0.242368232504351, -0.110999951753211, 0.174574390311335, -0.76456493319124
), Proposed_pvalue = c(0.0131801103443272, 0.583155149115837, 
0.870570103632653, 0.783460676404866, 0.154142429946211)), row.names = c(NA, 
5L), class = "data.frame"))

Jak mohu použít bigstatsr a paralelně vypočítat tuto funkci získat výstupy? Děkuji moc za vaše úsilí a pomoc.