R Under development (unstable) (2023-11-12 r85514 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(robustlmm)
Loading required package: robustlmm
Loading required package: lme4
Loading required package: Matrix
> 
> isPackageInstalled <- function(package) {
+     return(length(find.package(package, quiet = TRUE)) > 0)
+ }
> 
> 
> set.seed(1)
> oneWay <- generateAnovaDatasets(1, 1, 10, 4,
+                                 lmeFormula = y ~ 1,
+                                 heavyLmeRandom = ~ 1,
+                                 heavyLmeGroups = ~ Var2,
+                                 lqmmRandom = ~ 1,
+                                 lqmmGroup = "Var2",
+                                 groups = cbind(rep(1:4, each = 10), rep(1:10, 4)),
+                                 varcov = matrix(1, 4, 4),
+                                 lower = 0)
> test <- function(result, expectedLength = 1) {
+     stopifnot(length(result) == expectedLength,
+               !is(result[[1]], "try-error"))
+     print(result[[1]])
+ }
> test(expected <- fitDatasets_lmer(oneWay))
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | Var2)
   Data: data
REML criterion at convergence: 211.4088
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups:  Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> test(fitDatasets_lmer_bobyqa(oneWay))
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | Var2)
   Data: data
REML criterion at convergence: 211.4088
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups:  Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> test(fitDatasets_lmer_Nelder_Mead(oneWay))
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | Var2)
   Data: data
REML criterion at convergence: 211.4088
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups:  Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> fitDatasets_rlmer_custom <- function(datasets) {
+     return(fitDatasets_rlmer(datasets,
+                              method = "DASvar",
+                              tuningParameter = c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14),
+                              label = "fitDatasets_rlmer_custom"))
+ }
> test(fitDatasets_rlmer_custom(oneWay))
Robust linear mixed model fit by DASvar 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.516   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.046  
> test(fitDatasets_rlmer_DAStau(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.529   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.046  
> test(fitDatasets_rlmer_DAStau_lmerNoFit(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.529   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.046  
> test(fitDatasets_rlmer_DASvar(oneWay))
Robust linear mixed model fit by DASvar 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.516   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.046  
> test(fitDatasets_rlmer_DAStau_noAdj(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.404   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.047  
> test(fitDatasets_rlmer_DAStau_k_0_5(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.00    
 Residual             3.41    
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      2.109  
> test(fitDatasets_rlmer_DAStau_k_0_5_noAdj(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 0.000   
 Residual             3.003   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
       2.13  
> test(fitDatasets_rlmer_DAStau_k_2(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.499   
 Residual             3.247   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      1.919  
> test(fitDatasets_rlmer_DAStau_k_2_noAdj(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.664   
 Residual             3.203   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      1.914  
> test(fitDatasets_rlmer_DAStau_k_5(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> test(fitDatasets_rlmer_DAStau_k_5_noAdj(oneWay))
Robust linear mixed model fit by DAStau 
Formula: y ~ (1 | Var2) 
   Data: data 
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups: Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> if (isPackageInstalled("heavy")) {
+     test(fitDatasets_heavyLme(oneWay))
+ }
> if (isPackageInstalled("lqmm")) {
+     test(fitDatasets_lqmm(oneWay))
+ }
Call: lqmm::lqmm(fixed = y ~ 1, random = ~1, group = "Var2", covariance = "pdDiag", 
    data = data)

Quantile 0.5 

Fixed effects:
(Intercept)  
      2.154  

Covariance matrix of the random effects:
(Intercept) 
      2.189 

Residual scale parameter: 1.213 (standard deviation 3.432)
Log-likelihood: -106.4 

Number of observations: 40 
Number of groups: 10 
> ## if (isPackageInstalled("rlme")) {
> ##     test(fitDatasets_rlme(oneWay)) # won't work as dataset is balanced
> ## }
> if (isPackageInstalled("robustvarComp")) {
+     fitDatasets_varComprob_custom <- function(datasets, postFit) {
+         lcontrol <- robustvarComp::varComprob.control(lower = datasets[["lower"]])
+         return(
+             fitDatasets_varComprob(
+                 datasets,
+                 lcontrol,
+                 "fitDatasets_varComprob_custom",
+                 postFit = postFit
+             )
+         )
+     }
+     test(fitDatasets_varComprob_custom(oneWay))
+     test(fitDatasets_varComprob_compositeTau(oneWay))
+     test(fitDatasets_varComprob_compositeTau_OGK(oneWay))
+     test(fitDatasets_varComprob_compositeTau_2SGS(oneWay))
+     test(fitDatasets_varComprob_compositeS(oneWay))
+     test(fitDatasets_varComprob_compositeS_OGK(oneWay))
+     test(fitDatasets_varComprob_compositeS_2SGS(oneWay))
+     test(fitDatasets_varComprob_S(oneWay))
+     test(fitDatasets_varComprob_S_OGK(oneWay))
+     test(fitDatasets_varComprob_S_2SGS(oneWay))
+ }

 Method:  compositeTau 

Fixed effects: 
(Intercept) 
    1.90939 
Random effect variances: 
[1] 0.5533093
Residual variance: 
error variance 
      11.13874 
Value of the objective function at convergence: 
[1] 8.217236

 Method:  compositeTau 

Fixed effects: 
(Intercept) 
    1.90939 
Random effect variances: 
[1] 0.5533093
Residual variance: 
error variance 
      11.13874 
Value of the objective function at convergence: 
[1] 8.217236

 Method:  compositeTau 

Fixed effects: 
(Intercept) 
    1.90939 
Random effect variances: 
[1] 0.5533084
Residual variance: 
error variance 
      11.13874 
Value of the objective function at convergence: 
[1] 8.217236

 Method:  compositeTau 

Fixed effects: 
(Intercept) 
    1.90939 
Random effect variances: 
[1] 0.5533086
Residual variance: 
error variance 
      11.13874 
Value of the objective function at convergence: 
[1] 8.217236

 Method:  compositeS 

Fixed effects: 
(Intercept) 
   2.299836 
Random effect variances: 
[1] 1.472417
Residual variance: 
error variance 
      9.873617 
Value of the objective function at convergence: 
[1] 33.56952

 Method:  compositeS 

Fixed effects: 
(Intercept) 
   2.299836 
Random effect variances: 
[1] 1.472418
Residual variance: 
error variance 
      9.873616 
Value of the objective function at convergence: 
[1] 33.56952

 Method:  compositeS 

Fixed effects: 
(Intercept) 
   2.299836 
Random effect variances: 
[1] 1.472419
Residual variance: 
error variance 
      9.873615 
Value of the objective function at convergence: 
[1] 33.56952

 Method:  S 

Fixed effects: 
(Intercept) 
   1.848618 
Random effect variances: 
[1] 0
Residual variance: 
error variance 
      19.50105 
Value of the objective function at convergence: 
[1] 69.37515

 Method:  S 

Fixed effects: 
(Intercept) 
   1.848617 
Random effect variances: 
[1] 0
Residual variance: 
error variance 
      19.50105 
Value of the objective function at convergence: 
[1] 69.37515

 Method:  S 

Fixed effects: 
(Intercept) 
   1.848617 
Random effect variances: 
[1] 0
Residual variance: 
error variance 
      19.50105 
Value of the objective function at convergence: 
[1] 69.37515
> 
> ## test datasetIndices argument
> 
> set.seed(1)
> oneWay3 <- generateAnovaDatasets(3, 1, 10, 4,
+                                  lmeFormula = y ~ 1,
+                                  heavyLmeRandom = ~ 1,
+                                  heavyLmeGroups = ~ Var2,
+                                  lqmmRandom = ~ 1,
+                                  lqmmGroup = "Var2",
+                                  groups = cbind(rep(1:4, each = 10), rep(1:10, 4)),
+                                  varcov = matrix(1, 4, 4),
+                                  lower = 0)
> 
> test(actual <- fitDatasets_lmer(oneWay3, datasetIndices = 1))
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | Var2)
   Data: data
REML criterion at convergence: 211.4088
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 1.632   
 Residual             3.195   
Number of obs: 40, groups:  Var2, 10
Fixed Effects:
(Intercept)  
      1.899  
> stopifnot(all.equal(expected, actual, check.attributes = FALSE))
> 
> test(fitDatasets_lmer(oneWay3, datasetIndices = 2:3), expectedLength = 2)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | Var2)
   Data: data
REML criterion at convergence: 238.2448
Random effects:
 Groups   Name        Std.Dev.
 Var2     (Intercept) 4.455   
 Residual             3.980   
Number of obs: 40, groups:  Var2, 10
Fixed Effects:
(Intercept)  
      2.017  
> 
> ## test rlme implementation
> if (require(rlme)) {
+ 
+     data(schools)
+     datasets <- createDatasetsFromList(list(schools),
+                                        formula = y ~ 1 + sex + age + (1 | region) + (1 | region:school),
+                                        trueBeta = rep(NA_real_, 3),
+                                        trueSigma = NA_real_,
+                                        trueTheta = rep(NA_real_, 2))
+     fits <- fitDatasets_rlme(datasets)
+     result <- processFit(fits[[1]], all = TRUE)
+     stopifnot(length(result$coefficients) == 3,
+               length(result$standardErrors) == 3,
+               length(result$tValues) == 3,
+               length(result$sigma) == 1,
+               length(result$thetas) == 2,
+               length(result$residuals) == NROW(schools))
+ }
Loading required package: rlme
> 
> 
> proc.time()
   user  system elapsed 
  26.82    0.39   27.20