reuse mixed effects structures when X and Z are constant

[andportnoy]

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]

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

Might be worth forking pymer4 and adding refit to it.

[andportnoy]

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.