Running 'RFsimulate.R' [9s/10s]
  Running 'krige.R' [7s/7s]
  Running 'lgcp.R' [7s/8s]
  Running 'lgm.R' [15s/16s]
  Running 'lgmRaster.R' [105s/70s]
  Running 'likfitLgm.R' [9s/10s]
  Running 'matern.R' [5s/6s]
  Running 'maternGmrfPrec.R' [8s/9s]
  Running 'profLlgm.R' [12s/11s]
  Running 'simLgcp.R' [6s/7s]
  Running 'stackRasterList.R' [5s/5s]
Running the tests in 'tests/lgmRaster.R' failed.
Complete output:
  > #+ setup
  > library('geostatsp')
  Loading required package: Matrix
  Loading required package: raster
  Loading required package: sp
  > #'
  > 
  > #' # simulated data
  > 
  > 
  > #+ simData
  > if(.Platform$OS.type == 'windows') {
  + 	Ncell =  30
  + } else {
  + 	Ncell = 40
  + }
  > myRaster = squareRaster(extent(0,6000,0,6000), Ncell)
  > 
  > myParam=c(oneminusar=0.1, conditionalVariance=2.5^2,shape=2)
  > myQ = maternGmrfPrec(myRaster, param=myParam)
  > attributes(myQ)$info$optimalShape
       shape   variance      range   cellSize 
    4.092496 110.524266 900.000000 150.000000 
  > set.seed(0)
  > mySim = RFsimulate(attributes(myQ)$info$optimalShape, myRaster)
  New output format of RFsimulate: S4 object of class 'RFsp';
  for a bare, but faster array format use 'RFoptions(spConform=FALSE)'.
  > 
  > otherPar = c(intercept=1, beta = 2, tau=10)
  > myCov = myRaster
  > values(myCov) = rep(seq(-1,1,len=ncol(myCov)), nrow(myCov))
  > 
  > myLambda = otherPar['intercept'] + otherPar['beta'] * myCov + mySim
  > myY = myLambda
  > values(myY)=rnorm(prod(dim(myLambda)),values(myLambda), sd=otherPar['tau']) 
  > 
  > names(myCov) = 'x'
  > names(myY) = gsub("^layer\\.","sim", names(mySim))
  > #'
  > 
  > 
  > 
  > #+ plotSim, fig.cap='simulated data', fig.subcap =c('random effect','observed')
  > plot(mySim)
  > 
  > plot(myY)
  > #'
  > 
  > #' grid search	
  > 
  > #+ simGrid
  > myResR = lgm(formula = sim ~ x, 
  +     data=raster::stack(myY, myCov), 
  +     oneminusar = exp(seq(log(0.05), log(0.2),len=25)),
  +     nugget = exp(seq(log(5), log(100),len=21)), shape=2, 
  +     adjustEdges=TRUE,
  +     mc.cores=1+(.Platform$OS.type=='unix') )		
  > #'
  > 
  > #+ simPlot
  > Sbreaks = c(-100000,-50,-20,-10, -5,  -2, -1,0)
  > 
  > myCol = mapmisc::colourScale(
  +     breaks = Sbreaks + max(myResR$array[,-1,'logLreml',]),
  +     style='fixed',
  +     col=terrain.colors
  + )
  > 
  > image(	
  +     as.numeric(dimnames(myResR$array)$propNugget)[-1], 
  +     as.numeric(dimnames(myResR$array)$oneminusar),
  +     myResR$array[1,-1,'logLreml',],
  +     xlab = 'propNugget', ylab='oneminusar',
  +     log='xy',
  +     col=myCol$col, breaks=myCol$breaks)
  > mapmisc::legendBreaks("topright", breaks = Sbreaks, col=myCol$col)
  > 
  > 
  > points(myResR$param['propNugget'], myResR$param['oneminusar'])
  > text(myResR$param['propNugget'], myResR$param['oneminusar'], 'mle', pos=3)
  > 
  > points(otherPar['tau']^2/myParam['conditionalVariance'], myParam['oneminusar'], pch=16, col='red')
  > text(otherPar['tau']^2/myParam['conditionalVariance'], myParam['oneminusar'], 'truth', pos=3)
  > #'
  > 
  > 
  > #+ results
  > myResR$param
                        mlDeviance                     remlDeviance 
                        0.00000000                       0.00000000 
                             tausq                             xisq 
                       93.90633978                       4.19961919 
                            boxcox                          sigmasq 
                        1.00000000                      92.44806887 
                       (Intercept)                                x 
                        1.97966863                       0.66607461 
              (Intercept)BetaHatSe                       xBetaHatSe 
                        2.85481766                       6.26365369 
                             range                            shape 
                      959.27881425                       2.00000000 
                       sigmasqXisq                     optimalShape 
                       22.01344093                       1.40743615 
     optimalVarianceAndShape.shape optimalVarianceAndShape.variance 
                        1.57574045                      28.22056200 
                        oneminusar                       propNugget 
                        0.08908987                      22.36067977 
  > attributes(myQ)$info$optimalShape
       shape   variance      range   cellSize 
    4.092496 110.524266 900.000000 150.000000 
  > otherPar
  intercept      beta       tau 
          1         2        10 
  > #'
  > 
  > #' checkResults
  > 
  > #+ checkResultsSetup
  > 
  > oneRes = c(
  +     oneminusar = as.numeric(dimnames(myResR$array)$oneminusar[5]), 
  +     shape=2, 
  +     propNugget = as.numeric(dimnames(myResR$array)$propNugget[4]))
  > 
  > Q = maternGmrfPrec(
  +     myRaster, 
  +     param=c(conditionalVariance=1, oneRes[c('oneminusar','shape')]),
  +     adjustEdges=TRUE)
  > 
  > Qinv = as.matrix(solve(Q))
  > 
  > obsCov = cbind(y=values(myY), intercept=1, x=values(myCov))
  > Xseq = 2:ncol(obsCov)
  > Nobs = nrow(obsCov)
  > #'
  > 
  > #' no nugget
  > 
  > #+ checkNoNugget
  > (xProdQ =  crossprod(obsCov, Q) %*% obsCov)
  3 x 3 Matrix of class "dgeMatrix"
                        y     intercept             x
  y         351652.940041  1.372903e+00 -1.706098e+00
  intercept      1.372903  6.182566e-01 -1.484923e-14
  x             -1.706098 -1.830480e-14  3.241081e-01
  > xProdQ - myResR$model$extras$ssq[,,1,'ssq',  as.character(oneRes['oneminusar'])]
  3 x 3 Matrix of class "dgeMatrix"
                        y     intercept             x
  y         -4.138565e-08 -4.449774e-13  2.026379e-12
  intercept  9.907630e-13 -2.087219e-14  1.465494e-14
  x          3.330669e-13  6.217249e-15 -1.104672e-14
  > 
  > myResR$array[,1,c('(Intercept)BetaHat','xBetaHat'),as.character(oneRes['oneminusar'])]
  (Intercept)BetaHat           xBetaHat 
            2.220604          -5.263978 
  > myResR$model$extras$ssq[Xseq,1,1,'beta',  as.character(oneRes['oneminusar'])]
  (Intercept)           x 
     2.220604   -5.263978 
  > (betahat=drop(solve(xProdQ[Xseq,Xseq]) %*% xProdQ[Xseq, -Xseq]))
  [1]  2.220604 -5.263978
  > 
  > myResR$array[,1,c('xisqHatMl','xisqHatReml'),as.character(oneRes['oneminusar'])]
    xisqHatMl xisqHatReml 
     219.7756    220.0506 
  > myResR$model$extras$ml[,1,c('profiledVarianceHatMl','profiledVarianceHatReml'),  as.character(oneRes['oneminusar'])]
    profiledVarianceHatMl profiledVarianceHatReml 
                 219.7756                220.0506 
  > resid = obsCov[,1] - betahat[1] - obsCov[,3]*betahat[2]
  > (varhat = diag(crossprod(resid, Q) %*% resid)/(nrow(obsCov) - c(ml=0,reml=length(Xseq))))
        ml     reml 
  219.7756 220.0506 
  > 
  > 
  > 
  > -0.5*myResR$model$extras$ml[,'0','m2logLml', as.character(oneRes['oneminusar'])]
  [1] -6856.989
  > myResR$array[,1,'logLml',as.character(oneRes['oneminusar'])]
  [1] -6856.989
  > 
  > if(requireNamespace("mvtnorm", quietly=TRUE))
  +   mvtnorm::dmvnorm(
  +       resid, 
  +       sigma = varhat['ml'] * Qinv, log=TRUE)
  [1] -6856.989
  > #'
  > 
  > #' with nugget
  > 
  > #+ checkWithNugget
  > V =  (1/oneRes['propNugget']) * Qinv + Diagonal(nrow(Q))
  > Vinv = solve(V)
  > 
  > (xProdQ =  crossprod(obsCov, Vinv) %*% obsCov)
  3 x 3 Matrix of class "dgeMatrix"
                        y     intercept             x
  y         122219.865304  8.199378e+00 -3.062111e+00
  intercept      8.199378  3.961904e+00 -3.047562e-13
  x             -3.062111 -3.102241e-13  1.999803e+00
  > myResR$model$extras$ssq[,,  as.character(oneRes['propNugget']),'ssq',  as.character(oneRes['oneminusar'])]
               datasetCol
  datasetRow              sim  (Intercept)             x
    sim         122219.865304 8.199378e+00 -3.062111e+00
    (Intercept)      8.199378 3.961904e+00  2.753353e-14
    x               -3.062111 2.753353e-14  1.999803e+00
  > 
  > 
  > (betahat= drop(solve(xProdQ[Xseq,Xseq]) %*% xProdQ[Xseq, -Xseq]))
  [1]  2.069555 -1.531206
  > myResR$array[,
  +     as.character(oneRes['propNugget']), 
  +     c('(Intercept)BetaHat','xBetaHat'),
  +     as.character(oneRes['oneminusar'])]
  (Intercept)BetaHat           xBetaHat 
            2.069555          -1.531206 
  > 
  > resid = obsCov[,1] - betahat[1] - obsCov[,3]*betahat[2]
  > (varhat = diag(crossprod(resid, Vinv) %*% resid)/(nrow(obsCov) - c(ml=0,reml=length(Xseq))))
        ml     reml 
  76.37388 76.46947 
  > myResR$array[,
  +     as.character(oneRes['propNugget']), 
  +     c('tausqHatMl','tausqHatReml'),
  +     as.character(oneRes['oneminusar'])]
    tausqHatMl tausqHatReml 
      76.37388     76.46947 
  > 
  > -0.5*myResR$model$extras$ml[,as.character(oneRes['propNugget']),'m2logLml', as.character(oneRes['oneminusar'])]
  [1] -6167.726
  > myResR$array[,as.character(oneRes['propNugget']),'logLml',as.character(oneRes['oneminusar'])]
  [1] -6167.726
  > 
  > if(requireNamespace("mvtnorm", quietly=TRUE))
  +   mvtnorm::dmvnorm(
  +       resid, 
  +       sigma = as.matrix(varhat['ml'] * V), log=TRUE)
  [1] -6167.726
  > 
  > #' 
  > 
  > #' # swiss rain
  > 
  > #+ swissRain
  > data('swissRainR')
  > 
  > anotherx = raster(swissRainR[['alt']])
  > values(anotherx) = seq(0,1,len=ncell(anotherx))
  > names(anotherx) = "myvar"
  > 
  > swissRainR2 = brick(swissRainR[['alt']], 
  +     sqrt(swissRainR[['prec1']]),
  +     anotherx)
  > #'
  > 
  > #+ aggregateIfWindows
  > 
  > if(.Platform$OS.type=='windows') {
  +   swissRainR2 = raster::aggregate(swissRainR2, fact=2)
  + }
  > #'
  > 
  > #+ swissRainFit
  > 
  > swissResR =  lgm(
  +     formula=layer ~ alt+ myvar, 
  +     data=swissRainR2, shape=2,
  +     oneminusar = exp(seq(log(0.05), log(0.1), len=11)),
  +     nugget = exp(seq(log(0.25), log(2.5), len=11)),
  +     adjustEdges=TRUE,
  +     mc.cores=1+(.Platform$OS.type=='unix') )		
  > #'
  > 
  > #+ swissRainPlot
  > 
  > myCol = mapmisc::colourScale(
  +     breaks = Sbreaks + max(swissResR$array[,-1,'logLreml',]),
  +     style='fixed',
  +     col=terrain.colors
  + )
  > 
  > image(	
  +     as.numeric(dimnames(swissResR$array)$propNugget)[-1], 
  +     as.numeric(dimnames(swissResR$array)$oneminusar), 
  +     swissResR$array[1,-1,'logLreml',],
  +     xlab = 'propNugget', ylab='oneminusar',
  +     log='xy',
  +     col=myCol$col, breaks=myCol$breaks)
  > mapmisc::legendBreaks("topright", breaks = Sbreaks, col=myCol$col)
  > #'
  > 
  > 
  > 
  > #' boxcox
  > #+ boxCox
  > 
  > 
  > yBC = sqrt(myY + 1 - minValue(myY))
  > names(yBC) = names(myY)
  > myResBC = lgm(
  +     formula = sim ~ x, 
  +     data=raster::stack(yBC, myCov), 
  +     oneminusar = exp(seq(log(0.05), log(0.15), len=11)),
  +     nugget = exp(seq(log(5), log(50), len=11)),
  +     shape=2, 
  +     mc.cores=1+(.Platform$OS.type=='unix'), 
  +     fixBoxcox=FALSE,
  +     adjustEdges=FALSE)
  Error in ml[bestBcList[1, 1], , , , drop = FALSE] : 
    invalid subscript type 'list'
  Calls: lgm ... callGeneric -> eval -> eval -> lgm -> lgm -> loglikGmrf
  Execution halted