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

> options(na.action=na.omit)
> 
> # This data set comes from a family study, where there were doubts about
> #  the correctness of the lmekin fit. (Courtesy M deAndrade)
> # It's a good test of several things: there are missing values, and the
> #  data is not in family order
> load('brdat.rda')
> 
> library(kinship2)
Loading required package: Matrix
Loading required package: quadprog
> library(nlme)  
> bped <- with(brdat, pedigree(id, father, mother, sex, famid=family))
> kmat <- 2*kinship(bped)
> 
> fit1 <- lmekin(height ~ sex + (1|id), data= brdat, varlist=kmat)
> 
> # Compute the loglik of a fit
> blik <- function(fit, kin=kmat) {
+     keptrows <- as.numeric(names(fit$resid))  #not tossed away
+     indx <- match(dimnames(kin)[[1]], brdat$id[keptrows], nomatch=0)
+     
+     k2 <- kin[indx>0, indx>0]   # retained subjects
+     n <- nrow(k2)                # number who remain
+     vmat <- fit$sigma^2 * (diag(n) + unlist(VarCorr(fit))*k2)
+     
+     r2 <- fit$resid[indx]   #residuals, in kmat order
+     smat <- gchol(as(vmat, "bdsmatrix"))  #cholesky decomposition
+     
+     rsum <- sum(r2 * solve(smat, r2, full=TRUE))
+     csum <- n*log(2*pi) + sum(log(diag(smat)))
+ 
+     # Check using Matrix routines
+     rsum2  <- sum(r2 * solve(vmat, r2))
+     csum2  <- n*log(2*pi) + 2*sum(log(diag(chol(vmat))))
+     if (!all.equal(c(rsum, csum), c(rsum2, csum2))) cat("Matrix error")
+    -(rsum + csum)/2
+ }
> 
> # Test on a model I can verify with lme
> tfit1 <- lmekin(height ~ sex + (1|family), brdat)
> tfit2 <- lme(height ~ sex, random= ~1|family, brdat, method="ML",
+              na.action=na.omit)
> tkmat <- with(brdat, bdsBlock(id, family))
> blik(tfit1, as(tkmat, 'dsCMatrix'))
[1] -5733.688
>                                  
> 
> proc.time()
   user  system elapsed 
   3.20    0.39    3.53