R Under development (unstable) (2024-02-19 r85946 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.

> #+ setup
> library('geostatsp')
Loading required package: Matrix
Loading required package: terra
terra 1.7.71
> #'
> 
> #' # simulated data
> 
> # exclude this line to use the RandomFields package
> options(useRandomFields = FALSE)
> 
> Ncell =  40
> 
> myRaster = squareRaster(ext(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)
install the RandomFields package for faster simulations
> 
> 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=c(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)		
> #'
> 
> #+ simPlot
> Sbreaks = c(-100000,-50,-20,-10, -5,  -2, -1,0)
> 
> if(requireNamespace("mapmisc", quietly=TRUE)) {
+ myCol = mapmisc::colourScale(
+     breaks = Sbreaks + max(myResR$array[,-1,'logLreml',], na.rm=TRUE),
+     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 
                     96.36887135                       3.71022520 
                          boxcox                          sigmasq 
                      1.00000000                     101.81739187 
                     (Intercept)                                x 
                     -0.20960843                       0.65228343 
            (Intercept)BetaHatSe                       xBetaHatSe 
                      3.47935734                       7.47260548 
                           range                            shape 
                   1021.72841054                       2.00000000 
                     sigmasqXisq                     optimalShape 
                     27.44237517                       1.43855082 
   optimalVarianceAndShape.shape optimalVarianceAndShape.variance 
                      1.59322794                      34.94381807 
                      oneminusar                       propNugget 
                      0.07937005                      25.97386040 
> 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"
                    sim     intercept            x
sim       353839.892739 -2.021419e+00 2.074457e+00
intercept     -2.021419  6.182566e-01 4.476280e-14
x              2.074457  3.394507e-14 3.241081e-01
> xProdQ - myResR$model$extras$ssq[,,1,'ssq',  as.character(oneRes['oneminusar'])]
3 x 3 Matrix of class "dgeMatrix"
                    sim     intercept             x
sim       -5.809125e-08 -3.050893e-13  1.159961e-12
intercept -9.503509e-14  4.440892e-15  8.722190e-15
x         -3.463896e-14 -1.364880e-14 -2.181588e-14
> 
> myResR$array[,1,c('(Intercept)BetaHat','xBetaHat'),as.character(oneRes['oneminusar'])]
(Intercept)BetaHat           xBetaHat 
         -3.269546           6.400511 
> myResR$model$extras$ssq[Xseq,1,1,'beta',  as.character(oneRes['oneminusar'])]
(Intercept)           x 
  -3.269546    6.400511 
> (betahat=drop(solve(xProdQ[Xseq,Xseq]) %*% xProdQ[Xseq, -Xseq]))
intercept         x 
-3.269546  6.400511 
> 
> myResR$array[,1,c('xisqHatMl','xisqHatReml'),as.character(oneRes['oneminusar'])]
  xisqHatMl xisqHatReml 
   221.1375    221.4143 
> myResR$model$extras$ml[,1,c('profiledVarianceHatMl','profiledVarianceHatReml'),  as.character(oneRes['oneminusar'])]
  profiledVarianceHatMl profiledVarianceHatReml 
               221.1375                221.4143 
> 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 
221.1375 221.4143 
> 
> 
> 
> -0.5*myResR$model$extras$ml[,'0','m2logLml', as.character(oneRes['oneminusar'])]
[1] -6861.931
> myResR$array[,1,'logLml',as.character(oneRes['oneminusar'])]
[1] -6861.931
> 
> if(requireNamespace("mvtnorm", quietly=TRUE))
+   mvtnorm::dmvnorm(
+       resid, 
+       sigma = varhat['ml'] * Qinv, log=TRUE)
> #'
> 
> #' 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"
                    sim     intercept            x
sim       124321.060517 -1.342479e+00 2.559467e+00
intercept     -1.342479  3.961904e+00 3.621256e-12
x              2.559467  3.623296e-12 1.999803e+00
> myResR$model$extras$ssq[,,  as.character(oneRes['propNugget']),'ssq',  as.character(oneRes['oneminusar'])]
             datasetCol
datasetRow              sim   (Intercept)             x
  sim         124321.060517 -1.342479e+00  2.559467e+00
  (Intercept)     -1.342479  3.961904e+00 -9.295934e-13
  x                2.559467 -9.295934e-13  1.999803e+00
> 
> 
> (betahat= drop(solve(xProdQ[Xseq,Xseq]) %*% xProdQ[Xseq, -Xseq]))
 intercept          x 
-0.3388471  1.2798594 
> myResR$array[,
+     as.character(oneRes['propNugget']), 
+     c('(Intercept)BetaHat','xBetaHat'),
+     as.character(oneRes['oneminusar'])]
(Intercept)BetaHat           xBetaHat 
        -0.3388471          1.2798594 
> 
> 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 
77.69833 77.79558 
> myResR$array[,
+     as.character(oneRes['propNugget']), 
+     c('tausqHatMl','tausqHatReml'),
+     as.character(oneRes['oneminusar'])]
  tausqHatMl tausqHatReml 
    77.69833     77.79558 
> 
> -0.5*myResR$model$extras$ml[,as.character(oneRes['propNugget']),'m2logLml', as.character(oneRes['oneminusar'])]
[1] -6181.48
> myResR$array[,as.character(oneRes['propNugget']),'logLml',as.character(oneRes['oneminusar'])]
[1] -6181.48
> 
> if(requireNamespace("mvtnorm", quietly=TRUE))
+   mvtnorm::dmvnorm(
+       resid, 
+       sigma = as.matrix(varhat['ml'] * V), log=TRUE)
> 
> #' 
> 
> #' # swiss rain
> 
> #+ swissRain
> data('swissRainR')
> swissRainR= unwrap(swissRainR)
> 
> anotherx = rast(swissRainR[['alt']])
> values(anotherx) = seq(0,1,len=ncell(anotherx))
> names(anotherx) = "myvar"
> 
> swissRainR2 = c(swissRainR[['alt']], 
+     sqrt(swissRainR[['prec1']]),
+     anotherx)
> #'
> 
> #+ aggregateIfWindows
> 
> if(.Platform$OS.type=='windows') {
+   swissRainR2 = aggregate(swissRainR2, fact=2)
+ }
> #'
> 
> #+ swissRainFit
> 
> swissResR =  lgm(
+     formula=prec1 ~ alt+ myvar, 
+     data=swissRainR2, shape=2,
+     oneminusar = exp(seq(log(0.01), log(0.1), len=11)),
+     nugget = exp(seq(log(0.25), log(2.5), len=11)),
+     adjustEdges=TRUE )		
Warning messages:
1: In maternGmrfPrec.dgCMatrix(NN, param = c(shape = as.vector(shape),  :
  grid is 45 by 16which may be too small for range 19.8997487421323 cells, 1- Ar param0.01
2: In maternGmrfPrec.dgCMatrix(NN, param = c(shape = as.vector(shape),  :
  grid is 45 by 16which may be too small for range 17.7124615423636 cells, 1- Ar param0.0125892541179417
3: In maternGmrfPrec.dgCMatrix(NN, param = c(shape = as.vector(shape),  :
  grid is 45 by 16which may be too small for range 15.7601693452856 cells, 1- Ar param0.0158489319246111
> #'
> 
> #+ swissRainPlot
> if(requireNamespace("mapmisc", quietly=TRUE)) {
+ myCol = mapmisc::colourScale(
+     breaks = Sbreaks + max(swissResR$array[,-1,'logLreml',], na.rm=TRUE),
+     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 - min(minmax(myY)))
> names(yBC) = names(myY)
> myResBC = lgm(
+     formula = sim ~ x, 
+     data=c(yBC, myCov), 
+     oneminusar = exp(seq(log(0.05), log(0.15), len=11)),
+     nugget = exp(seq(log(5), log(50), len=11)),
+     shape=2, reml=FALSE, 
+     mc.cores=1, 
+     fixBoxcox=FALSE,
+     adjustEdges=FALSE)
> #'
> 
> #+ boxCoxPlot
> 
> plot(myResBC$profL$boxcox,type='o', 
+   ylim=max(c(-10000,myResBC$profL$boxcox[,2]), na.rm=TRUE)-c(3,0))
> 
> myResBC$param
                      mlDeviance                     remlDeviance 
                       0.0000000                        0.0000000 
                           tausq                             xisq 
                      27.9378913                        2.8004229 
                          boxcox                          sigmasq 
                       2.0500000                       34.9486223 
                     (Intercept)                                x 
                      26.2773154                       -0.3980429 
            (Intercept)BetaHatSe                       xBetaHatSe 
                       0.3032953                        0.5147615 
                           range                            shape 
                     810.8258529                        2.0000000 
                     sigmasqXisq                     optimalShape 
                      12.4797660                        1.3122568 
   optimalVarianceAndShape.shape optimalVarianceAndShape.variance 
                       1.5269332                       16.3399839 
                      oneminusar                       propNugget 
                       0.1204112                        9.9763116 
> if(requireNamespace("mapmisc", quietly=TRUE)){
+ myCol = mapmisc::colourScale(
+     breaks = Sbreaks,
+     style='fixed',
+     col=terrain.colors
+ )
+ 
+ image(	
+     myResBC$profL$twoDim$x[-1], 
+     myResBC$profL$twoDim$y,
+     myResBC$profL$twoDim$z[-1,],
+     log='xy', 
+     ylab = 'range', xlab='propNugget',
+     col=myCol$col, breaks=myCol$breaks+max(myResBC$array[,,'logLreml',], na.rm=TRUE))
+ mapmisc::legendBreaks("topright",  myCol)
+ points(myResBC$param['propNugget'], myResBC$param['oneminusar'])
+ }
> #'
> 
> #' optimizing, doesn't work
> 
> #+ optimizepropNugget
> 
> if(Sys.info()['user'] =='patrick' & FALSE) {
+   myResRopt = lgm(
+       formula = sim ~ x, 
+       data=c(myY, myCov), 
+       oneminusar = seq(0.05, 0.2,len=4),
+       shape=2)		
+   
+   if(!interactive()) pdf("doesntwork.pdf")
+   plot(myResRopt$array[,,'oneminusar',], myResRopt$array[,,'propNugget',])
+   
+   if(!interactive()) dev.off()	
+   
+   swissResRoptAr =  lgm(
+       formula=layer ~ alt+ myvar, 
+       data=swissRainR2, shape=2,
+       oneminusar = seq(0.1, 0.5, len=4),
+       adjustEdges=FALSE
+   )
+   
+   swissResRopt =  lgm(
+       formula=layer ~ alt+ myvar, 
+       data=swissRainR2, shape=2,
+       adjustEdges=FALSE
+   )
+   
+   
+   swissResRopt$summary
+   
+ # with edge correction.  
+ # time consuming, only run this if Patrick is checking
+   
+   
+   
+ # optimize only nugget
+   swissResROptNug =  lgm(
+       formula=prec1 ~ alt+ myvar, 
+       data=swissRainR2, shape=2,
+       oneminusar=seq(0.05, 0.1, len=4),
+       adjustEdges=FALSE,fixNugget=TRUE,
+       mc.cores=1
+   )
+   
+   plot(swissResROptNug$profL$range, type='l')
+ }
> 
> 
> #'
> 
> proc.time()
   user  system elapsed 
  52.92    1.32   54.26