R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> library(coxme)
Loading required package: survival
Loading required package: bdsmatrix

Attaching package: 'bdsmatrix'

The following object is masked from 'package:base':

    backsolve

> require(nlme)
Loading required package: nlme
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
> tdata <- eortc
> tdata$center2 <- factor(tdata$center)
> tdata$trt2 <- factor(tdata$trt)
> 
> fit1 <- lme(y ~ trt, random= ~ 1|center2, data=tdata,
+             method="ML")
> fit2 <- lmekin(y ~ trt + (1|center), tdata)
> 
> #I get the same loglik but slightly different coefficients.
> # (The loglik is flat on top, and we are using different starting
> #  estimates.)  So we need to use a toloerance for the test.
> vcomp <- function(m1, m2, ...) {
+     temp1 <- as.numeric(VarCorr(m1)[,2])  #lme gives a matrix of text
+ #    temp2 <- m2$sigma^2 * c(unlist(VarCorr(m2)), 1) #lmekin the variance struct
+     temp2 <- c(unlist(VarCorr(m2)), m2$sigma^2) #lmekin the variance struct
+     aeq(temp1, sqrt(temp2), ...)
+ }
> 
> aeq(fit1$logLik, fit2$loglik)
[1] TRUE
> aeq(fit1$sigma, fit2$sigma, tol=1e-4)
[1] TRUE
> aeq(fixef(fit1), fixef(fit2), tol=1e-4)
[1] TRUE
> vcomp(fit1, fit2, tol=1e-3)
[1] TRUE
> 
> # A second fit, using the matrix form
> ucen <- sort(unique(tdata$center))
> kmat <- diag(length(ucen))
> dimnames(kmat) <- list(ucen, ucen)
> fit3 <- lmekin(y ~ trt + (1|center), tdata, varlist=kmat)
> aeq(fit3$coefficients, fit2$coefficients)
[1] TRUE
> aeq(unlist(VarCorr(fit3)), unlist(VarCorr(fit2)))
[1] TRUE
> 
> # Force the same coefs. 
> temp <- as.numeric(unclass(fit1$modelStruct$reStruct)[[1]]) 
> fit3 <- lmekin(y~trt + (1|center), tdata,
+                vfixed= exp(2*temp))
> aeq(fit1$logLik, fit3$loglik)
[1] TRUE
> aeq(fit1$sigma, fit3$sigma)
[1] TRUE
> aeq(fixef(fit1), fixef(fit3))
[1] TRUE
> aeq(ranef(fit1), ranef(fit3), check.attributes=FALSE)
[1] TRUE
> vcomp(fit1, fit3, tol=1e-7)
[1] TRUE
> aeq(residuals(fit1), residuals(fit3))
[1] TRUE
> 
> # Now a model with random slopes and intercepts
> fit5 <- lme(y ~ trt, random= ~ trt|center, data=tdata,
+             method="REML")
> fit6 <- lmekin(y~ trt + (1+trt | center), data=tdata,
+                method="REML")
> aeq(fit5$logLik, fit6$loglik)
[1] TRUE
> aeq(fit5$sigma, fit6$sigma, tol=1e-4)
[1] TRUE
> aeq(fixef(fit5), fixef(fit6), tol=1e-4)
[1] TRUE
> 
> 
> #
> # Use an lme data set
> #
> mystool <- as.data.frame(ergoStool) #get rid of contrast attributes
> 
> efit1 <-  lme(effort ~ Type, data=mystool, random= ~1|Subject,
+             method="ML")
> efit2 <- lmekin(effort ~ Type + (1|Subject), mystool)
> aeq(efit1$logLik, efit2$loglik)
[1] TRUE
> aeq(fixef(efit1), fixef(efit2))
[1] TRUE
> vcomp(efit1, efit2, tol=1e-5)
[1] TRUE
> aeq(efit1$sigma, efit2$sigma, tol=1e-5)
[1] TRUE
> 
> efit3 <-lme(effort ~ Type, data=mystool, random= ~1|Subject,
+             method="REML")
> efit4 <- lmekin(effort ~ Type + (1|Subject), mystool, method="REML")
> 
> aeq(efit3$logLik, efit4$loglik)
[1] TRUE
> aeq(fixef(efit3), fixef(efit4))
[1] TRUE
> vcomp(efit3, efit4, tol=1e-3)
[1] TRUE
> aeq(efit3$sigma, efit4$sigma, tol=1e-4)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   2.92    0.32    3.25