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.

> ## test the availability of generic functions.
> ae <- function(target, current, ...) {
+     ret <- all.equal(target, current, ...)
+     if (isTRUE(ret)) return(ret)
+     print(ret)
+     stop("Objects not equal")
+ }
> chgClass <- function(object) {
+     class(object) <- sub("(merMod|mer)", "rlmerMod", class(object))
+     object
+ }
> 
> require(robustlmm)
Loading required package: robustlmm
Loading required package: lme4
Loading required package: Matrix
> 
> 
> set.seed(3)
> sleepstudy2 <- within(sleepstudy, {
+     Group <- letters[1:4]
+     Covar <- rnorm(180)
+ })
> rfm <- rlmer(Reaction ~ Days + (Days|Subject) + (1|Group), sleepstudy2,
+              rho.e = cPsi, rho.b = cPsi, doFit=FALSE)
boundary (singular) fit: see help('isSingular')
> fm <- lmer(Reaction ~ Days + (Days|Subject) + (1|Group), sleepstudy2)
boundary (singular) fit: see help('isSingular')
> ## all three print the same:
> print(rfm)
Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.
> show(rfm)
Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.
> summary(rfm)
Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.
> 
> ## object information
> ## ae(df.residual(fm), df.residual(rfm))
> ae(formula(fm), formula(rfm))
[1] TRUE
> stopifnot(isLMM(rfm))
> stopifnot(isREML(rfm))
> ae(model.frame(fm), model.frame(rfm))
[1] TRUE
> ae(model.matrix(fm), model.matrix(rfm), check.attributes=FALSE)
[1] TRUE
> nobs(rfm)
[1] 180
> ## ae(getInitial(fm), getInitial(rfm))
> ae(terms(fm), terms(rfm))
[1] TRUE
> weights(rfm)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> 
> ## basic accessors for the results
> ae(chgClass(coef(fm)), coef(rfm))
[1] TRUE
> ## dummy.coef(rfm)
> stopifnot(inherits(try(deviance(rfm), silent=TRUE), "try-error"))
> stopifnot(inherits(try(extractAIC(rfm), silent=TRUE), "try-error"))
> family(rfm)

Family: gaussian 
Link function: identity 

> ## after version 1.1 fitted values are named
> if (packageVersion("lme4") > "1.1") ae(fitted(fm), fitted(rfm))
[1] TRUE
> ae(fixef(fm), fixef(rfm))
[1] TRUE
> stopifnot(inherits(try(logLik(rfm), silent=TRUE), "try-error"))
> ae(chgClass(ranef(fm, condVar=FALSE)), ranef(rfm))
[1] TRUE
> ## after version 1.1 fitted values are named
> if (packageVersion("lme4") > "1.1") ae(resid(fm), resid(rfm))
[1] TRUE
> ae(sigma(fm), sigma(rfm))
[1] TRUE
> ## weighted.residuals(rfm)
> 
> ## var-covar methods
> ## VarCorr(rfm)
> ae(chgClass(VarCorr(fm)), VarCorr(rfm))
[1] TRUE
> ae(vcov(fm), vcov(rfm), tolerance=1e-4)
[1] TRUE
> ## vcov(rfm)
> 
> ## confidence intervals
> ## confint(rfm)
> 
> ## other (deprecated?)
> theta(rfm)
     Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
               0.9667371                0.0151686                0.2309095 
       Group.(Intercept) 
               0.0000000 
> 
> ## other methods
> getInfo(rfm)
$data
sleepstudy2

$coef
(Intercept)        Days 
  251.40510    10.46729 

$stderr
(Intercept)        Days 
   6.824576    1.545788 

$varcomp
(Intercept) | Subject        Days | Subject   (Intercept) | Group 
            24.740552              5.922129              0.000000 

$sigma
[1] 25.59181

$correlations
(Intercept) x Days | Subject 
                  0.06554938 

$rho.e
[1] "classic"

$rho.sigma.e
[1] "classic"

$rho.b_1
[1] "classic"

$rho.sigma.b_1
[1] "classic"

$rho.b_2
[1] "classic"

$rho.sigma.b_2
[1] "classic"

> compare(fm, rfm)
                             fm           rfm         
Coef                                                  
(Intercept)                  251.4 (6.82) 251.4 (6.82)
Days                          10.5 (1.55)  10.5 (1.55)
                                                      
VarComp                                               
(Intercept) | Subject        24.74        24.74       
Days | Subject                5.92         5.92       
(Intercept) | Group           0.00         0.00       
                                                      
Correlations                                          
(Intercept) x Days | Subject 0.0655       0.0655      
                                                      
sigma                        25.6         25.6        
                                                      
REML                         1744                     
rho.e                                     classic     
rho.sigma.e                               classic     
rho.b_1                                   classic     
rho.sigma.b_1                             classic     
rho.b_2                                   classic     
rho.sigma.b_2                             classic     
> update(rfm, ~ . + Covar)
boundary (singular) fit: see help('isSingular')
Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.
> 
> ## predict method (various examples)
> ae(predict(fm), predict(rfm))
[1] TRUE
> if (packageVersion("lme4") > "1.1") {
+     ae(predict(fm,re.form=NA), predict(rfm,re.form=NA))
+     newdata <- with(sleepstudy, expand.grid(Subject=unique(Subject),
+                                             Days=3:5, Group=letters[1:2]))
+     ae(predict(fm,newdata), predict(rfm,newdata))
+     ae(predict(fm,newdata,re.form=NA), predict(rfm,newdata,re.form=NA))
+     ae(predict( fm,newdata, re.form=~(1|Subject)),
+        predict(rfm,newdata, re.form=~(1|Subject)))
+ }
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   2.79    0.28    3.03