R Under development (unstable) (2023-10-18 r85349 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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(spatsurv) Welcome to 'spatsurv': Spatial Survival Analysis B. M. Taylor & B. S. Rowlingson. Type 'spatsurvVignette()' to view the package vignette. Type 'citation("spatsurv")' to view the citation for this package. Please see the spatsurv package NEWS file for latest additions, changes and bug fixes. > library(sp) > library(spatstat.explore) Loading required package: spatstat.data Loading required package: spatstat.geom spatstat.geom 3.2-6 Loading required package: spatstat.random spatstat.random 3.1-6 Loading required package: nlme spatstat.explore 3.2-3 > library(survival) > > par(mfrow=c(1,1)) > > set.seed(10) > > n <- 100 > DIST <- exponentialHaz() > > > OMEGA <- 1 > > > > # Generate spatially correlated survival data ... > dat <- simsurv(X=cbind( age=runif(n,5,50),sex=rbinom(n,1,0.5),cancer=rbinom(n,1,0.2)), + dist=DIST, + omega=OMEGA, + cov.model=ExponentialCovFct(), + mcmc.control=mcmcpars(nits=100,burn=10,thin=10)) | |===== | 11% : Burn-in | |========== | 22% : Burn-in | |=============== | 33% : Burn-in | |==================== | 44% : Burn-in | |========================= | 56% : Burn-in | |============================== | 67% : Burn-in | |=================================== | 78% : Burn-in | |======================================== | 89% : Burn-in | |=============================================| 100% : Burn-in | |= | 1% : Run [12/100] | |= | 2% : Run [13/100] | |== | 3% : Run [14/100] | |== | 4% : Run [15/100] | |=== | 6% : Run [16/100] | |=== | 7% : Run [17/100] | |==== | 8% : Run [18/100] | |==== | 9% : Run [19/100] | |===== | 10% : Run [20/100] | |===== | 11% : Run [21/100] | |====== | 12% : Run [22/100] | |====== | 13% : Run [23/100] | |======= | 15% : Run [24/100] | |======= | 16% : Run [25/100] | |======== | 17% : Run [26/100] | |======== | 18% : Run [27/100] | |========= | 19% : Run [28/100] | |========= | 20% : Run [29/100] | |========== | 21% : Run [30/100] | |========== | 22% : Run [31/100] | |=========== | 24% : Run [32/100] | |=========== | 25% : Run [33/100] | |============ | 26% : Run [34/100] | |============ | 27% : Run [35/100] | |============= | 28% : Run [36/100] | |============= | 29% : Run [37/100] | |============== | 30% : Run [38/100] | |============== | 31% : Run [39/100] | |=============== | 33% : Run [40/100] | |=============== | 34% : Run [41/100] | |================ | 35% : Run [42/100] | |================ | 36% : Run [43/100] | |================= | 37% : Run [44/100] | |================= | 38% : Run [45/100] | |================== | 39% : Run [46/100] | |================== | 40% : Run [47/100] | |=================== | 42% : Run [48/100] | |=================== | 43% : Run [49/100] | |==================== | 44% : Run [50/100] | |==================== | 45% : Run [51/100] | |===================== | 46% : Run [52/100] | |===================== | 47% : Run [53/100] | |====================== | 48% : Run [54/100] | |====================== | 49% : Run [55/100] | |======================= | 51% : Run [56/100] | |======================= | 52% : Run [57/100] | |======================== | 53% : Run [58/100] | |======================== | 54% : Run [59/100] | |========================= | 55% : Run [60/100] | |========================= | 56% : Run [61/100] | |========================== | 57% : Run [62/100] | |========================== | 58% : Run [63/100] | |=========================== | 60% : Run [64/100] | |=========================== | 61% : Run [65/100] | |============================ | 62% : Run [66/100] | |============================ | 63% : Run [67/100] | |============================= | 64% : Run [68/100] | |============================= | 65% : Run [69/100] | |============================== | 66% : Run [70/100] | |============================== | 67% : Run [71/100] | |=============================== | 69% : Run [72/100] | |=============================== | 70% : Run [73/100] | |================================ | 71% : Run [74/100] | |================================ | 72% : Run [75/100] | |================================= | 73% : Run [76/100] | |================================= | 74% : Run [77/100] | |================================== | 75% : Run [78/100] | |================================== | 76% : Run [79/100] | |=================================== | 78% : Run [80/100] | |=================================== | 79% : Run [81/100] | |==================================== | 80% : Run [82/100] | |==================================== | 81% : Run [83/100] | |===================================== | 82% : Run [84/100] | |===================================== | 83% : Run [85/100] | |====================================== | 84% : Run [86/100] | |====================================== | 85% : Run [87/100] | |======================================= | 87% : Run [88/100] | |======================================= | 88% : Run [89/100] | |======================================== | 89% : Run [90/100] | |======================================== | 90% : Run [91/100] | |========================================= | 91% : Run [92/100] | |========================================= | 92% : Run [93/100] | |========================================== | 93% : Run [94/100] | |========================================== | 94% : Run [95/100] | |=========================================== | 96% : Run [96/100] | |=========================================== | 97% : Run [97/100] | |============================================ | 98% : Run [98/100] | |============================================ | 99% : Run [99/100] | |=============================================| 100% : Run [100/100] Returning last set of simulated survival times. Mean acceptance: 0.6087682 > > coords <- dat$coords > SIGMA <- dat$cov.parameters[1] > PHI <- dat$cov.parameters[2] > > par(mfrow=c(2,2)) > plot(coords,col=grey(1-dat$survtimes/max(dat$survtimes)),pch=19) > > X <- as.data.frame(dat$X) # covariates > > survtimes <- dat$survtimes > censtimes <- runif(n,min(survtimes),max(survtimes)) > survdat <- gencens(survtimes,censtimes) > > > # priors > betaprior <- betapriorGauss(mean=0,sd=10) > omegaprior <- omegapriorGauss(mean=0,sd=10) > etaprior <- etapriorGauss(mean=log(c(SIGMA,PHI)),sd=c(0.3,0.3)) > priors <- mcmcPriors( betaprior=betaprior, + omegaprior=omegaprior, + etaprior=etaprior, + call=indepGaussianprior, + derivative=derivindepGaussianprior) > > # create SpatialPointsDataFrame containing the covariate data and coordinates of the survival data > spatdat <- SpatialPointsDataFrame(coords,data=as.data.frame(X)) > spatdat$ss <- survdat > > #browser() > > if(TRUE){ + ss <- survspat( formula=ss~age+sex+cancer, + data=spatdat, + dist=DIST, + #cov.model=covmodel(model="exponential",pars=NULL), + cov.model=ExponentialCovFct(), + mcmc.control=mcmcpars(nits=100,burn=10,thin=9), + priors=priors) + } Getting initial estimates of model parameters using maximum likelihood on non-spatial version of the model Initial optimisation via BFGS ... Refining optimum via Nelder Mead ... Refining optimum via BFGS ... $par [1] 0.04948869 -0.05338080 0.25718756 -1.35087123 $value [1] 88.71317 $counts function gradient 10 1 $convergence [1] 0 $message NULL Done. Calibrating MCMC algorithm and finding initial values ... Constructing quadratic approximation to posterior (this can take some time) ... Done. Running MCMC ... [,1] [,2] [,3] [,4] [,5] [1,] 0.0001323030 0.0002318211 -0.0001067002 -0.003441417 0.000000000 [2,] 0.0002318211 0.0735997075 0.0066363724 -0.045407068 0.000000000 [3,] -0.0001067002 0.0066363724 0.0896745628 -0.026319323 0.000000000 [4,] -0.0034414167 -0.0454070682 -0.0263193232 0.136134422 0.000000000 [5,] 0.0000000000 0.0000000000 0.0000000000 0.000000000 0.001545814 [6,] 0.0000000000 0.0000000000 0.0000000000 0.000000000 0.002018676 [7,] 0.0000000000 0.0000000000 0.0000000000 0.000000000 0.000000000 [8,] 0.0000000000 0.0000000000 0.0000000000 0.000000000 0.000000000 [,6] [,7] [,8] [1,] 0.000000000 0.00000000 0.00000000 [2,] 0.000000000 0.00000000 0.00000000 [3,] 0.000000000 0.00000000 0.00000000 [4,] 0.000000000 0.00000000 0.00000000 [5,] 0.002018676 0.00000000 0.00000000 [6,] 0.005276104 0.00000000 0.00000000 [7,] 0.000000000 0.03576937 0.00000000 [8,] 0.000000000 0.00000000 0.04722017 | |===== | 11% : Burn-in | |========== | 22% : Burn-in | |=============== | 33% : Burn-in | |==================== | 44% : Burn-in | |========================= | 56% : Burn-in | |============================== | 67% : Burn-in | |=================================== | 78% : Burn-in | |======================================== | 89% : Burn-in | |=============================================| 100% : Burn-in | |= | 1% : Run [12/100] | |= | 2% : Run [13/100] | |== | 3% : Run [14/100] | |== | 4% : Run [15/100] | |=== | 6% : Run [16/100] | |=== | 7% : Run [17/100] | |==== | 8% : Run [18/100] | |==== | 9% : Run [19/100] | |===== | 10% : Run [20/100] | |===== | 11% : Run [21/100] | |====== | 12% : Run [22/100] | |====== | 13% : Run [23/100] | |======= | 15% : Run [24/100] | |======= | 16% : Run [25/100] | |======== | 17% : Run [26/100] | |======== | 18% : Run [27/100] | |========= | 19% : Run [28/100] | |========= | 20% : Run [29/100] | |========== | 21% : Run [30/100] | |========== | 22% : Run [31/100] | |=========== | 24% : Run [32/100] | |=========== | 25% : Run [33/100] | |============ | 26% : Run [34/100] | |============ | 27% : Run [35/100] | |============= | 28% : Run [36/100] | |============= | 29% : Run [37/100] | |============== | 30% : Run [38/100] | |============== | 31% : Run [39/100] | |=============== | 33% : Run [40/100] | |=============== | 34% : Run [41/100] | |================ | 35% : Run [42/100] | |================ | 36% : Run [43/100] | |================= | 37% : Run [44/100] | |================= | 38% : Run [45/100] | |================== | 39% : Run [46/100] | |================== | 40% : Run [47/100] | |=================== | 42% : Run [48/100] | |=================== | 43% : Run [49/100] | |==================== | 44% : Run [50/100] | |==================== | 45% : Run [51/100] | |===================== | 46% : Run [52/100] | |===================== | 47% : Run [53/100] | |====================== | 48% : Run [54/100] | |====================== | 49% : Run [55/100] | |======================= | 51% : Run [56/100] | |======================= | 52% : Run [57/100] | |======================== | 53% : Run [58/100] | |======================== | 54% : Run [59/100] | |========================= | 55% : Run [60/100] | |========================= | 56% : Run [61/100] | |========================== | 57% : Run [62/100] | |========================== | 58% : Run [63/100] | |=========================== | 60% : Run [64/100] | |=========================== | 61% : Run [65/100] | |============================ | 62% : Run [66/100] | |============================ | 63% : Run [67/100] | |============================= | 64% : Run [68/100] | |============================= | 65% : Run [69/100] | |============================== | 66% : Run [70/100] | |============================== | 67% : Run [71/100] | |=============================== | 69% : Run [72/100] | |=============================== | 70% : Run [73/100] | |================================ | 71% : Run [74/100] | |================================ | 72% : Run [75/100] | |================================= | 73% : Run [76/100] | |================================= | 74% : Run [77/100] | |================================== | 75% : Run [78/100] | |================================== | 76% : Run [79/100] | |=================================== | 78% : Run [80/100] | |=================================== | 79% : Run [81/100] | |==================================== | 80% : Run [82/100] | |==================================== | 81% : Run [83/100] | |===================================== | 82% : Run [84/100] | |===================================== | 83% : Run [85/100] | |====================================== | 84% : Run [86/100] | |====================================== | 85% : Run [87/100] | |======================================= | 87% : Run [88/100] | |======================================= | 88% : Run [89/100] | |======================================== | 89% : Run [90/100] | |======================================== | 90% : Run [91/100] | |========================================= | 91% : Run [92/100] | |========================================= | 92% : Run [93/100] | |========================================== | 93% : Run [94/100] | |========================================== | 94% : Run [95/100] | |=========================================== | 96% : Run [96/100] | |=========================================== | 97% : Run [97/100] | |============================================ | 98% : Run [98/100] | |============================================ | 99% : Run [99/100] h = 0.8344136 | |=============================================| 100% : Run [100/100] Time taken: 1.52251 > > par(mfrow=c(3,3)) > plot(ss$etasamp[,1],type="s", main="eta 1") > plot(ss$etasamp[,2],type="s", main="eta 2") > plot(ss$omegasamp[,1],type="s", main="omega") > plot(ss$betasamp[,1],type="s", main="beta 1") > plot(ss$betasamp[,2],type="s", main="beta 2") > plot(ss$Ysamp[,1],type="s", main="Y 1") > plot(ss$Ysamp[,floor(n/2)],type="s", main="Y floor(n/2)") > plot(ss$tarrec,type="s", main="log posterior") > plot(dat$Y,colMeans(ss$Ysamp),xlab="TRUE Y",ylab="Estimated Y") > abline(0,1) > > # > #save(list=ls(),file="simout1.RData") > # > > proc.time() user system elapsed 6.79 0.59 7.37