tests/testthat/test-predVar-Matern-corrMatrix.R

In spaMM: Mixed-Effect Models, with or without Spatial Random Effects

cat(crayon::yellow("test-predVar-Matern-corrMatrix"))
if (spaMM.getOption("example_maxtime")>0.7) { ##  not based on real timing
  cat(crayon::yellow(': THREE messages "spaMM is not able... " expected\n    (when corrnames are "Gibraltar" ... and ZAnames are "-5.3469:36.1291" ...)\n'))
  
  # checks predVar w/o permutation, + Matern vs corrMatrix, + spprec T/F  ... + perm_Q
  
  data("blackcap")
  MLcorMat <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
                         nu=4,rho=0.4)
  blackcap$name <- as.factor(rownames(blackcap))
  bidon <- list(T=TRUE)
  old_perm_Q <- spaMM.getOption("perm_Q")
  for (nam in c("T","N")) { # so it ends on default, NULL
    spaMM.options(perm_Q=bidon[[nam]])
    ## (1) Single variable present in the data 
    (f1 <- HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
                 corrMatrix=MLcorMat,method="ML", 
                 control.HLfit=list(sparse_precision=TRUE)))
    (f2 <- corrHLfit(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
                    ranFix=list(corrPars=list("1"=list(nu=4,rho=0.4))),method="ML"))
    (f3 <- HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
                corrMatrix=MLcorMat,HLmethod="ML", 
                control.HLfit=list(sparse_precision=FALSE)))
    (f4 <- corrHLfit(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
                    ranFix=list(corrPars=list("1"=list(nu=4,rho=0.4))),method="ML"))
    (f5 <- HLCor(migStatus ~ means+ corrMatrix(1|longitude+latitude),data=blackcap,
                corrMatrix=MLcorMat,method="ML")) # Check that order of data is respected in the Zmatrix for this "unsafe" input.
    # imput as precision matrix
    (f6 <- fitme(migStatus ~ means + corrMatrix(1|name), data=blackcap,
                covStruct=list(precision=as_precision(MLcorMat))))
    # Manual version of the same:
    as_mat <- proxy::as.matrix(MLcorMat, diag=1) 
    prec_mat <- solve(as_mat) ## precision factor matrix
    f7 <- fitme(migStatus ~ means + corrMatrix(1|name), data=blackcap,
                covStruct=list(precision=prec_mat))
    # the last two fitme spprec fits ssslightly differ from previous ones, even from f1,f2 spprec fits => its a fitme vs others difference. 
    testthat::expect_true(diff(range(c(logLik(f1),logLik(f2),logLik(f3),logLik(f4),logLik(f5),logLik(f6),logLik(f7))))<1e-8) # 
    testthat::expect_true(diff(range(c(predict(f1, newdata=f1$data[1,]), predict(f2, newdata=f1$data[1,]), predict(f3, newdata=f1$data[1,]),
                                       predict(f4, newdata=f1$data[1,]), predict(f5, newdata=f1$data[1,]), predict(f6, newdata=f1$data[1,]), 
                                       predict(f7, newdata=f1$data[1,]))))<1e-5)  
    #  I once had a problem with perm_Q=TRUE in test-predVar-Matern-corrMatrix -> predict(f2,.) 
    # so make sure to check this when changing relevant code. Seems OK now.
    # New ZA was Z_1x1 * A_nxn[identified row] -> 1xn and cov_newLv_oldv_list was 1xn
    # so the product failed. The new ZA should instead be 1x1.
    (c4a <- get_predVar(f4,variances=list(cov=TRUE))[3:2,3:2])
    (c4b <- get_predVar(f4,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c1 <- get_predVar(f1,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c2 <- get_predVar(f2,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c3 <- get_predVar(f3,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c5 <- get_predVar(f5,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c6 <- get_predVar(f6,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    (c7 <- get_predVar(f7,newdata=f4$data[c(3,2),],variances=list(cov=TRUE)))
    testthat::expect_true(diff(range(c(c1[1],c2[1],c3[1],c4a[1],c4b[1],c5[1],c6[1],c7[1])))<1e-5)  
    
  }
  spaMM.options(perm_Q=old_perm_Q)
} else {
  cat(crayon::yellow(': ONE message "spaMM is not able... " expected\n    (when corrnames are "Gibraltar" ... and ZAnames are "-5.3469:36.1291" ...)\n'))
}



data("blackcap") 
## Here we manually reconstruct the correlation matrix 
##  of the ML fit produced by corrHLfit:
MLcorMat <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
                       nu=0.6285603,rho=0.0544659)
blackcap$name <- as.factor(rownames(blackcap))                
## (1) Single variable present in the data 
HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
      corrMatrix=MLcorMat,method="ML")
## (2) Same, permuted: still gives correct result
perm <- sample(14)
# Permuted matrix (with permuted names)
pmat <- as.matrix(MLcorMat)[perm,perm]
HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap,
      corrMatrix=as.dist(pmat),method="ML")
## (3) Other grouping terms (note the messages):
HLCor(migStatus ~ means+ corrMatrix(1|longitude+latitude),data=blackcap,
      corrMatrix=MLcorMat,method="ML")

# incidentally, check that prediction runs for levels new relative to the data, not reltive to the corrMatrix
subfit <- HLCor(migStatus ~ means+ corrMatrix(1|name),data=blackcap[1:12,],
                corrMatrix=MLcorMat,method="ML")
predict(subfit,newdata=blackcap[13:14,])