kutaslab / fitgrid

Multichannel event-related time-series regression modeling for EEG, MEG, and sensor array data
https://kutaslab.github.io/fitgrid
BSD 3-Clause "New" or "Revised" License
8 stars 4 forks source link

reuse mixed effects structures when X and Z are constant #108

Open andportnoy opened 5 years ago

andportnoy commented 5 years ago

Fitting mixed effects models is slow, so it's worth investigating whether we can reuse the output of lFormula when the predictors are constant across time. We are also guaranteed that they stay constant across channels.

andportnoy commented 5 years ago

lme4 has refit: https://cran.r-project.org/web/packages/lme4/lme4.pdf

Might be worth forking pymer4 and adding refit to it.

andportnoy commented 5 years ago
library(lme4)                                                                                                         

set.seed(42)                                                                                                        

m = 3000                                                                                                            
n = 300                                                                                                             

# predictors                                                                                                        
x <- rnorm(m)                                                                                                       
a <- sample(LETTERS[1:5], m, replace=TRUE)                                                                          

# responses                                                                                                         
Y <- matrix(rnorm(m * n), ncol=n)                                                                                   

# fit normally                                                                                                      
fit_measurements <- numeric(n)                                                                                      
for (i in 1:n) {                                                                                                    
    y <- Y[, i]                                                                                                     
    start = Sys.time()                                                                                              
    fit <- lmer(y ~ x + (x | a))                                                                                    
    fit_measurements[i] <- as.numeric(Sys.time() - start)                                                           
}                                                                                                                   

# fit the same data using fit object from before                                                                    
refit_measurements <- numeric(n)                                                                                    
for (i in 1:n) {                                                                                                    
    y <- Y[, i]                                                                                                     
    start = Sys.time()                                                                                              
    fit2 <- refit(fit, y)                                                                                           
    refit_measurements[i] <- as.numeric(Sys.time() - start)                                                         
}                                                                                                                   

fit_time <- sum(fit_measurements)                                                                                   
refit_time <- sum(refit_measurements)                                                                               

print(paste("Fit time: ", fit_time))                                                                                
print(paste("Refit time: ", refit_time))                                                                            
print(paste("Refit/fit ratio: ", refit_time / fit_time))

prints:

Loading required package: Matrix
[1] "Fit time:  21.2563896179199"
[1] "Refit time:  9.30353307723999"
[1] "Refit/fit ratio:  0.437681715684999"

So might actually be worth looking into.