R/parameterEstimation.R
In bentaylor1/lgcp: Log-Gaussian Cox Process
###
# Functions for parameter estimation
###
##' ginhomAverage function
##'
##' A function to estimate the inhomogeneous pair correlation function for a spatiotemporal point process. See equation (8) of Diggle P, Rowlingson B, Su T (2005).
##'
##' @param xyt an object of class stppp
##' @param spatial.intensity A spatialAtRisk object
##' @param temporal.intensity A temporalAtRisk object
##' @param time.window time interval contained in the interval xyt$tlim over which to compute average. Useful if there is a lot of data over a lot of time points.
##' @param rvals Vector of values for the argument r at which g(r) should be evaluated (see ?pcfinhom). There is a sensible default.
##' @param correction choice of edge correction to use, see ?pcfinhom, default is Ripley isotropic correction
##' @param suppresswarnings Whether or not to suppress warnings generated by pcfinhom
##' @param ... other parameters to be passed to pcfinhom, see ?pcfinhom
##' @return time average of inhomogenous pcf, equation (13) of Brix and Diggle 2001.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
ginhomAverage <- function(xyt,spatial.intensity,temporal.intensity,time.window=xyt$tlim,rvals=NULL,correction="iso",suppresswarnings=FALSE,...){
time.window <- sort(as.integer(time.window))
verifyclass(spatial.intensity,"spatialAtRisk")
verifyclass(temporal.intensity,"temporalAtRisk")
density <- as.im(spatial.intensity)
approxscale <- mean(sapply(xyt$t[xyt$t>=time.window[1] & xyt$t<=time.window[2]],temporal.intensity))
den <- density
den$v <- den$v * approxscale
xyt$t <- as.integer(xyt$t)
numsamp <- min(c(200,xyt$n))
if(!suppresswarnings){
if (is.null(rvals)){
gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,...)
}
else{
gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals,...)
}
}
else{
if (is.null(rvals)){
suppressWarnings(gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,...))
}
else{
suppressWarnings(gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals,...))
}
}
r <- gin$r
nam <- names(gin)
nam <- nam[nam!="r"]
tls <- sapply(time.window[1]:time.window[2],function(x){sum(xyt$t==x)})
tls <- (time.window[1]:time.window[2])[tls!=0]
ntls <- length(tls)
pb <- txtProgressBar(min=tls[1],max=rev(tls)[1],style=3)
if(!suppresswarnings){
pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(pcfinhom(xyt[xyt$t==t],lambda=den,r=r,...))})
}
else{
suppressWarnings(pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(pcfinhom(xyt[xyt$t==t],lambda=den,r=r,...))}))
}
close(pb)
li <- as.list(pcf[[1]])
ct <- 1
if(ntls>1){
for(i in 2:ntls){
if (!class(pcf[[i]])[1]=="try-error"){
li <- add.list(li,as.list(pcf[[i]]))
ct <- ct + 1
}
}
}
li <- smultiply.list(li,1/ct)
ginhom <- gin
for (n in nam){
idx <- which(names(ginhom)==n)
ginhom[[idx]] <- li[[idx]]
}
attr(ginhom,"correction") <- correction
cat("Returning an average of",ct,"curves\n")
return(ginhom)
}
##' KinhomAverage function
##'
##' A function to estimate the inhomogeneous K function for a spatiotemporal point process. The method of computation is similar to
##' \link{ginhomAverage}, see eq (8) Diggle P, Rowlingson B, Su T (2005) to see how this is computed.
##'
##' @param xyt an object of class stppp
##' @param spatial.intensity A spatialAtRisk object
##' @param temporal.intensity A temporalAtRisk object
##' @param time.window time interval contained in the interval xyt$tlim over which to compute average. Useful if there is a lot of data over a lot of time points.
##' @param rvals Vector of values for the argument r at which the inhmogeneous K function should be evaluated (see ?Kinhom). There is a sensible default.
##' @param correction choice of edge correction to use, see ?Kinhom, default is Ripley isotropic correction
##' @param suppresswarnings Whether or not to suppress warnings generated by Kinhom
##' @return time average of inhomogenous K function.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
KinhomAverage <- function(xyt,spatial.intensity,temporal.intensity,time.window=xyt$tlim,rvals=NULL,correction="iso",suppresswarnings=FALSE){
time.window <- sort(as.integer(time.window))
verifyclass(spatial.intensity,"spatialAtRisk")
verifyclass(temporal.intensity,"temporalAtRisk")
density <- as.im(spatial.intensity)
approxscale <- mean(sapply(xyt$t[xyt$t>=time.window[1] & xyt$t<=time.window[2]],temporal.intensity))
den <- density
den$v <- den$v * approxscale
xyt$t <- as.integer(xyt$t)
numsamp <- min(c(200,xyt$n))
if(!suppresswarnings){
if (is.null(rvals)){
Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den)
}
else{
Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals)
}
}
else{
if (is.null(rvals)){
suppressWarnings(Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den))
}
else{
suppressWarnings(Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals))
}
}
r <- Kin$r
nam <- names(Kin)
nam <- nam[nam!="r"]
tls <- sapply(time.window[1]:time.window[2],function(x){sum(xyt$t==x)})
tls <- (time.window[1]:time.window[2])[tls!=0]
ntls <- length(tls)
pb <- txtProgressBar(min=tls[1],max=rev(tls)[1],style=3)
if(!suppresswarnings){
pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(Kinhom(xyt[xyt$t==t],lambda=den,r=r))})
}
else{
suppressWarnings(pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(Kinhom(xyt[xyt$t==t],lambda=den,r=r))}))
}
close(pb)
li <- as.list(pcf[[1]])
ct <- 1
if(ntls>1){
for(i in 2:ntls){
if (!class(pcf[[i]])[1]=="try-error"){
li <- add.list(li,as.list(pcf[[i]]))
ct <- ct + 1
}
}
}
li <- smultiply.list(li,1/ct)
Kinhom <- Kin
for (n in nam){
idx <- which(names(Kinhom)==n)
Kinhom[[idx]] <- li[[idx]]
}
attr(Kinhom,"correction") <- correction
cat("Returning an average of",ct,"curves\n")
return(Kinhom)
}
##' spatialparsEst function
##'
##' Having estimated either the pair correlation or K functions using respectively \link{ginhomAverage} or \link{KinhomAverage}, the spatial
##' parameters sigma and phi can be estimated. This function provides a visual tool for this estimation procedure.
##'
##' To get a good choice of parameters, it is likely that the routine will have to be called several times in order to refine
##' the choice of sigma.range and phi.range.
##'
##' @param gk an R object; output from the function KinhomAverage or ginhomAverage
##' @param sigma.range range of sigma values to consider
##' @param phi.range range of phi values to consider
##' @param spatial.covmodel correlation type see ?CovarianceFct
##' @param covpars vector of additional parameters for certain classes of covariance function (eg Matern), these must be supplied in the order given in ?CovarianceFct
##' @param guess logical. Perform an initial guess at paramters? Alternative (the default) sets initial values in the middle of sigma.range and phi.range. NOTE: automatic parameter estimation can be can be unreliable.
##' @return rpanel function to help choose sigma nad phi by eye
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{KinhomAverage}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
spatialparsEst <- function(gk,sigma.range,phi.range,spatial.covmodel,covpars=c(),guess=FALSE){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
corchoice <- c() # to get rid of 'no visible binding' messages on checking
sigma <- c()
phi <- c()
sigmamin <- sigma.range[1]
sigmamax <- sigma.range[2]
phimin <- phi.range[1]
phimax <- phi.range[2]
idx <- which(names(gk)==attr(gk,"correction"))
env <- NULL
if (attr(gk,"fname") == "g[inhom]"){
panfun <- function(p){
env <<- environment()
r <- gk$r
egr <- suppressWarnings(exp(sapply(r,gu,sigma=p$sigma,phi=p$phi,model=spatial.covmodel,additionalparameters=covpars))-1)
gvals <- gk[[idx]]
if (p$transform!="none"){
if (p$transform=="log"){
gvals <- log(gk[[idx]])
egr <- log(egr)
}
}
plot(NULL,main="g[inhom]",xlim=c(0,max(r,na.rm=TRUE)),ylim=c(0,max(c(gvals[!is.infinite(gvals)&!is.nan(gk[[idx]])],egr[!is.infinite(egr)&!is.nan(egr)]))),xlab="r",ylab="g_inhom(r)",sub=paste(spatial.covmodel,"covariance function"))
lines(r,gvals)
lines(r,egr,col="orange")
legend("topright",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
return(p)
}
}
else if (attr(gk,"fname") == "K[inhom]"){
panfun <- function(p){
env <<- environment()
r <- gk$r
egr <- suppressWarnings(exp(sapply(r,gu,sigma=p$sigma,phi=p$phi,model=spatial.covmodel,additionalparameters=covpars))-1)
rdiff <- diff(r[1:2]) # do the integral on the discrete partition of r given by Kinhom
kest <- rdiff * egr[1] * 2 * pi * r[1]
for(i in 2:length(r)){
kest[i] <- kest[i-1] + rdiff * egr[i] * 2 * pi * r[i]
}
gvals <- gk[[idx]]
if (p$transform!="none"){
if (p$transform=="^1/4"){
gvals <- (gk[[idx]])^(1/4)
kest <- kest^(1/4)
}
}
plot(NULL,main="K[inhom]",xlim=c(0,max(r,na.rm=TRUE)),ylim=c(0,max(c(gvals[!is.infinite(gvals)],kest[!is.infinite(kest)]))),xlab="r",ylab="K_inhom(r)",sub=paste(spatial.covmodel,"covariance function"))
lines(r,gvals)
lines(r,kest,col="orange")
legend("topleft",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
return(p)
}
}
else{
stop("Argument gk is not of correct class, see ?ginhomAverage, ?KinhomAverage")
}
pancontrol <- rp.control("Parameter Estimation",aschar=FALSE)
##rp.radiogroup(pancontrol,var=corchoice,values=c("exponential","matern","whittle"),action=panfun,initval="exponential")
if (attr(gk,"fname") == "g[inhom]"){
rp.radiogroup(pancontrol,variable=transform,vals=c("none","log"),action=panfun,initval="none")
}
else if (attr(gk,"fname") == "K[inhom]"){
rp.radiogroup(pancontrol,variable=transform,vals=c("none","^1/4"),action=panfun,initval="^1/4")
}
# "quick" optimiser to get poor initial values for sigma and phi
r <- gk$r
if(guess){
if (attr(gk,"fname") == "g[inhom]"){
spfun <- function(sigmaphi){
egr <- suppressWarnings(exp(sapply(r,gu,sigma=sigmaphi[1],phi=sigmaphi[2],model=spatial.covmodel,additionalparameters=covpars))-1)
S <- suppressWarnings((egr-gk[[idx]])^2)
S[is.infinite(S)] <- NA
W <- 1/r^2
W[is.infinite(W)] <- NA
sidx <- !(is.na(W) | is.na(S))
ans <- sum((W*S)[sidx],na.rm=TRUE)/sum(W[sidx])
return(ans)
}
ops <- optim(c(sigmamin + (sigmamax-sigmamin)/2,phimin + (phimax-phimin)/2),spfun)
sigmainit <- ops$par[1]
phiinit <- ops$par[2]
}
else if (attr(gk,"fname") == "K[inhom]"){
sigphifun <- function(sigphi){
egr <- suppressWarnings(exp(sapply(r,gu,sigma=sigphi[1],phi=sigphi[2],model=spatial.covmodel,additionalparameters=covpars))-1)
rdiff <- diff(r[1:2]) # do the integral on the discrete partition of r given by Kinhom
kest <- rdiff * egr[1] * 2 * pi * r[1]
for(i in 2:length(r)){
kest[i] <- kest[i-1] + rdiff * egr[i] * 2 * pi * r[i]
}
gvals <- gk[[idx]]
gvals <- (gk[[idx]])^(1/4)
kest <- kest^(1/4)
gvals[is.infinite(gvals)] <- NA
kest[is.infinite(kest)] <- NA
S <- suppressWarnings((kest-gvals)^2)
S[is.infinite(S)] <- NA
W <- 1/r^2
W[is.infinite(W)] <- NA
sidx <- !(is.na(W) | is.na(S))
ans <- sum((W*S)[sidx],na.rm=TRUE)/sum(W[sidx])
return(ans)
}
ops <- optim(c(sigmamin + (sigmamax-sigmamin)/2,phimin + (phimax-phimin)/2),sigphifun)
sigmainit <- ops$par[1]
phiinit <- ops$par[2]
}
}
else{
sigmainit <- sigmamin + (sigmamax-sigmamin)/2
phiinit <- phimin + (phimax-phimin)/2
}
rp.slider(pancontrol,variable=sigma,from=sigmamin,to=sigmamax,initval=sigmainit,action=panfun,showvalue=TRUE)
rp.slider(pancontrol,variable=phi,from=phimin,to=phimax,initval=phiinit,action=panfun,showvalue=TRUE)
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
# construct a geometry string with the current height and new width:
hhh = as.numeric(tkwinfo("height",pancontrol$.handle))+50
www = 300
ggg = sprintf("%dx%d",www,hhh)
# resize
tkwm.geometry(pancontrol$.handle,ggg)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
## get the variables. Two points:
sigma <- get("p",envir=env)$sigma #as.numeric(tclvalue(pancontrol$sigma.tcl))
phi <- get("p",envir=env)$phi #as.numeric(tclvalue(pancontrol$phi.tcl))
return(list(sigma=sigma,phi=phi))
}
##' Cvb function
##'
##' This function is used in \code{thetaEst} to estimate the temporal correlation parameter, theta.
##'
##' @param xyt object of class stppp
##' @param spatial.intensity bivariate density estimate of lambda, an object of class im (produced from density.ppp for example)
##' @param N number of integration points
##' @param spatial.covmodel spatial covariance model
##' @param covpars additional covariance parameters
##' @return a function, see below.
##' Computes Monte carlo estimate of function C(v;beta) in Brix and Diggle 2001 pp 829 (... note later corrigendum to paper (2003) corrects the expression given in this paper)
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' }
##' @seealso \link{thetaEst}
##' @export
Cvb <- function(xyt,spatial.intensity,N=100,spatial.covmodel,covpars){
verifyclass(spatial.intensity,"im")
sar <- spatialAtRisk(list(X=spatial.intensity$xcol,Y=spatial.intensity$yrow,Zm=t(spatial.intensity$v)))
gsx <- length(xvals(sar))
gsy <- length(yvals(sar))
xy <- cbind(rep(xvals(sar),gsy),rep(yvals(sar),each=gsx))
wt <- as.vector(zvals(sar))
wt[is.na(wt)] <- 0
sidx <- sample(1:(gsx*gsy),N,prob=wt)
xy <- xy[sidx,]
pd <- as.vector(pairdist(xy))
cvb <- function(nu,sigma,phi,theta){
return(mean(exp(exp(-nu*theta)*gu(pd,sigma=sigma,phi=phi,model=spatial.covmodel,additionalparameters=covpars))-1))
}
return(cvb)
}
##' thetaEst function
##'
##' A tool to visually estimate the temporal correlation parameter theta; note that sigma and phi must have first been estiamted.
##'
##' @param xyt object of class stppp
##' @param spatial.intensity A spatial at risk object OR a bivariate density estimate of lambda, an object of class im (produced from density.ppp for example),
##' @param temporal.intensity either an object of class temporalAtRisk, or one that can be coerced into that form. If NULL (default), this is estimated from the data, seee ?muEst
##' @param sigma estimate of parameter sigma
##' @param phi estimate of parameter phi
##' @param theta.range range of theta values to consider
##' @param N number of integration points in computation of C(v,beta) (see Brix and Diggle 2003, corrigendum to Brix and Diggle 2001)
##' @param spatial.covmodel spatial covariance model
##' @param covpars additional covariance parameters
##' @return An r panel tool for visual estimation of temporal parameter theta
##' NOTE if lambdaEst has been invoked to estimate lambda, then the returned density should be passed to thetaEst as the argument spatial.intensity
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{lambdaEst}, \link{muEst}
##' @export
thetaEst <- function(xyt,spatial.intensity=NULL,temporal.intensity=NULL,sigma,phi,theta.range=c(0,10),N=100,spatial.covmodel="exponential",covpars=c()){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
if (inherits(spatial.intensity,"spatialAtRisk")){
spatial.intensity <- as.im(spatial.intensity)
}
if (is.null(spatial.intensity)){
spatial.intensity <- density(xyt)
}
if (is.null(temporal.intensity)){
temporal.intensity <- muEst(xyt)
}
else{
if (!inherits(temporal.intensity,"temporalAtRisk")){
temporal.intensity <- temporalAtRisk(temporal.intensity,tlim=xyt$tlim,xyt=xyt)
}
if(!all(xyt$tlim==attr(temporal.intensity,"tlim"))){
stop("Incompatible temporal.intensity, integer time limits (tlim and temporal.intensity$tlim) do not match")
}
}
thetamin <- theta.range[1]
thetamax <- theta.range[2]
tlim <- as.integer(xyt$tlim)
cvb <- Cvb(xyt=xyt,spatial.intensity=spatial.intensity,N=N,spatial.covmodel=spatial.covmodel,covpars=covpars)
xytlim <- sort(as.integer(xyt$tlim))
tvec <- xytlim[1]:xytlim[2]
if(length(tvec)==1){
stop("Insufficiently many time intervals: type as.integer(xyt$tlim) into console")
}
theta <- c() # to get rid of 'no visible binding' messages on checking
panfun <- function(p){
sigma <- as.numeric(p$sigma)
phi <- as.numeric(p$phi)
theta <- p$theta
uqt <- as.numeric(names(table(as.integer(xyt$t))))
tfit <- sapply(uqt,temporal.intensity)
autocov <- acf(table(as.integer(xyt$t)) - tfit,type="covariance",plot=FALSE)
r <- 0:(length(autocov$acf)-1)
plot(r,autocov$acf,type="l",main="Autocovariance",xlab="lag",ylab="acov(lag)")
legend("topright",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
abline(h=0)
rseq <- seq(r[1],r[length(r)],length.out=100)
theo <- sapply(rseq,cvb,sigma=sigma,phi=phi,theta=theta)
lines(rseq,theo*(autocov$acf[1]/theo[1]),type="l",col="orange")
return(p)
}
pancontrol <- rp.control("Parameter Estimation",aschar=FALSE)
rp.textentry(pancontrol,variable=sigma,initval=sigma,action=panfun)
rp.textentry(pancontrol,variable=phi,initval=phi,action=panfun)
rp.slider(pancontrol,variable=theta,from=thetamin,to=thetamax,initval=(thetamax-thetamin)/2,action=panfun,showvalue=TRUE)
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
## get the variables. Two points:
theta <- as.numeric(tclvalue(pancontrol$theta.tcl))
return(theta)
}
##' lambdaEst function
##'
##' Generic function for estimating bivariate densities by eye. Specific methods exist for stppp objects and ppp objects.
##'
##' @param xyt an object
##' @param ... additional arguments
##' @return method lambdaEst
##' @seealso \link{lambdaEst.stppp}, \link{lambdaEst.ppp}
##' @export
lambdaEst <- function(xyt,...){
UseMethod("lambdaEst")
}
##' lambdaEst.stppp function
##'
##' A tool for the visual estimation of lambda(s) via a 2 dimensional smoothing of the case locations. For parameter estimation, the alternative is
##' to estimate lambda(s) by some other means, convert it into a spatialAtRisk object and then into a pixel image object using the build in coercion
##' methods, this \code{im} object can then be fed to \link{ginhomAverage}, \link{KinhomAverage} or \link{thetaEst} for instance.
##'
##' The function lambdaEst is built directly on the density.ppp function and as such, implements a bivariate
##' Gaussian smoothing kernel. The bandwidth is initially that which is automatically chosen by the default method
##' of density.ppp. Since image plots of these kernel density estimates may not have appropriate
##' colour scales, the ability to adjust this is given with the slider 'colour adjustment'. With colour adjustment set
##' to 1, the default image.plot for the equivalent pixel image object is shown and for values less than 1, the colour
##' scheme is more spread out, allowing the user to get a better feel for the density that is being fitted. NOTE: colour
##' adjustment does not affect the returned density and the user should be aware that the returned density will 'look like'
##' that displayed when colour adjustment is set equal to 1.
##'
##'
##' @method lambdaEst stppp
##' @param xyt object of class stppp
##' @param weights Optional vector of weights to be attached to the points. May include negative values. See ?density.ppp.
##' @param edge Logical flag: if TRUE, apply edge correction. See ?density.ppp.
##' @param bw optional bandwidth. Set to NULL by default, which calls teh resolve.2D.kernel function for computing an initial value of this
##' @param ... arguments to be passed to plot
##' @return This is an rpanel function for visual choice of lambda(s), the output is a variable, varname, with the density *per unit time*
##' the variable varname can be fed to the function ginhomAverage or KinhomAverage as the argument density (see for example ?ginhomAverage), or into the
##' function thetaEst as the argument spatial.intensity.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{spatialAtRisk}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{muEst}
##' @export
lambdaEst.stppp <- function(xyt,weights=c(),edge=TRUE,bw=NULL,...){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
adjust <- c() # to get rid of 'no visible binding' messages on checking
plotdata <- c()
pow <- c()
if(is.null(bw)){
bw <- resolve.2D.kernel(NULL,x = xyt,adjust = 1)$sigma
}
varlist <- list(bw=bw,adjust=1,plotdata=TRUE,pow=1)
env <- NULL
panfun <- function(p){
env <<- environment()
bw <- as.numeric(p$bw)
adjust <- as.numeric(p$adjust)
d <- density(xyt,bandwidth=bw,weights=weights,edge=edge)
plotden <- d
plotden$v <- (d$v / diff(xyt$tlim))^p$pow
plot(plotden,main="Density",...)
if(p$plotdata){
points(xyt,pch="+",cex=0.5)
}
return(p)
}
panfun(varlist)
pancontrol <- rp.control("Density Estimation")
rp.textentry(pancontrol,variable=bw,initval=bw,action=panfun,title="bandwidth")
rp.checkbox(pancontrol,variable=plotdata,action=panfun,labels = "Plot Data",initval = TRUE)
rp.slider(pancontrol,variable=pow,from=0,to=1,initval=1,action=panfun,showvalue=TRUE,title="colour adjustment")
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
return(get("d",envir=env))
}
##' lambdaEst.ppp function
##'
##' A tool for the visual estimation of lambda(s) via a 2 dimensional smoothing of the case locations. For parameter estimation, the alternative is
##' to estimate lambda(s) by some other means, convert it into a spatialAtRisk object and then into a pixel image object using the build in coercion
##' methods, this \code{im} object can then be fed to \link{ginhomAverage}, \link{KinhomAverage} or \link{thetaEst} for instance.
##'
##' The function lambdaEst is built directly on the density.ppp function and as such, implements a bivariate
##' Gaussian smoothing kernel. The bandwidth is initially that which is automatically chosen by the default method
##' of density.ppp. Since image plots of these kernel density estimates may not have appropriate
##' colour scales, the ability to adjust this is given with the slider 'colour adjustment'. With colour adjustment set
##' to 1, the default image.plot for the equivalent pixel image object is shown and for values less than 1, the colour
##' scheme is more spread out, allowing the user to get a better feel for the density that is being fitted. NOTE: colour
##' adjustment does not affect the returned density and the user should be aware that the returned density will 'look like'
##' that displayed when colour adjustment is set equal to 1.
##'
##'
##' @method lambdaEst ppp
##' @param xyt object of class stppp
##' @param weights Optional vector of weights to be attached to the points. May include negative values. See ?density.ppp.
##' @param edge Logical flag: if TRUE, apply edge correction. See ?density.ppp.
##' @param bw optional bandwidth. Set to NULL by default, which calls teh resolve.2D.kernel function for computing an initial value of this
##' @param ... arguments to be passed to plot
##' @return This is an rpanel function for visual choice of lambda(s), the output is a variable, varname, with the density *per unit time*
##' the variable varname can be fed to the function ginhomAverage or KinhomAverage as the argument density (see for example ?ginhomAverage), or into the
##' function thetaEst as the argument spatial.intensity.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{spatialAtRisk}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{muEst}
##' @export
lambdaEst.ppp <- function(xyt,weights=c(),edge=TRUE,bw=NULL,...){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
adjust <- c() # to get rid of 'no visible binding' messages on checking
plotdata <- c()
pow <- c()
if(is.null(bw)){
bw <- resolve.2D.kernel(NULL,x = xyt,adjust = 1)$sigma
}
varlist <- list(bw=bw,adjust=1,plotdata=TRUE,pow=1)
env <- NULL
panfun <- function(p){
env <<- environment()
bw <- as.numeric(p$bw)
adjust <- as.numeric(p$adjust)
d <- density(xyt,sigma=bw,weights=weights,edge=edge)
plotden <- d
plotden$v <- d$v^p$pow
plot(plotden,main="Density",...)
if(p$plotdata){
points(xyt,pch="+",cex=0.5)
}
return(p)
}
panfun(varlist)
pancontrol <- rp.control("Density Estimation")
rp.textentry(pancontrol,variable=bw,initval=bw,action=panfun,title="bandwidth")
rp.checkbox(pancontrol,variable=plotdata,action=panfun,labels = "Plot Data",initval = TRUE)
rp.slider(pancontrol,variable=pow,from=0,to=1,initval=1,action=panfun,showvalue=TRUE,title="colour adjustment")
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
return(get("d",envir=env))
}
##' muEst function
##'
##' Computes a non-parametric estimate of mu(t). For the purposes of performing prediction, the alternatives are: (1) use a parameteric model as in Diggle P, Rowlingson B, Su T (2005),
##' or (2) a \link{constantInTime} model.
##'
##' @param xyt an stppp object
##' @param ... additional arguments to be passed to lowess
##' @return object of class temporalAtRisk giving the smoothed mut using the lowess function
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{temporalAtRisk}, \link{constantInTime}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}
##' @export
muEst <- function(xyt,...){
verifyclass(xyt,"stppp")
t <- as.integer(xyt$t)
tlim <- as.integer(xyt$tlim)
tseq <- tlim[1]:tlim[2]
counts <- sapply(tseq,function(x){sum(t==x)})
sm <- lowess(tseq,sqrt(counts),...)
return(temporalAtRisk(sm$y^2,tlim=tlim,xyt=xyt))
}
##' density.stppp function
##'
##' A wrapper function for \link{density.ppp}.
##'
##' @method density stppp
##' @param x an stppp object
##' @param bandwidth 'bandwidth' parameter, equivanent to parameter sigma in ?density.ppp ie standard deviation of isotropic Gaussian smoothing kernel.
##' @param ... additional arguments to be passed to density.ppp
##' @return bivariate density estimate of xyt; not this is a wrapper function for density.ppp
##' @seealso \link{density.ppp}
##' @export
density.stppp <- function(x,bandwidth=NULL,...){
density.ppp(x,...,sigma=bandwidth)
}
bentaylor1/lgcp documentation built on May 12, 2019, 2:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
###
# Functions for parameter estimation
###
##' ginhomAverage function
##'
##' A function to estimate the inhomogeneous pair correlation function for a spatiotemporal point process. See equation (8) of Diggle P, Rowlingson B, Su T (2005).
##'
##' @param xyt an object of class stppp
##' @param spatial.intensity A spatialAtRisk object
##' @param temporal.intensity A temporalAtRisk object
##' @param time.window time interval contained in the interval xyt$tlim over which to compute average. Useful if there is a lot of data over a lot of time points.
##' @param rvals Vector of values for the argument r at which g(r) should be evaluated (see ?pcfinhom). There is a sensible default.
##' @param correction choice of edge correction to use, see ?pcfinhom, default is Ripley isotropic correction
##' @param suppresswarnings Whether or not to suppress warnings generated by pcfinhom
##' @param ... other parameters to be passed to pcfinhom, see ?pcfinhom
##' @return time average of inhomogenous pcf, equation (13) of Brix and Diggle 2001.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
ginhomAverage <- function(xyt,spatial.intensity,temporal.intensity,time.window=xyt$tlim,rvals=NULL,correction="iso",suppresswarnings=FALSE,...){
time.window <- sort(as.integer(time.window))
verifyclass(spatial.intensity,"spatialAtRisk")
verifyclass(temporal.intensity,"temporalAtRisk")
density <- as.im(spatial.intensity)
approxscale <- mean(sapply(xyt$t[xyt$t>=time.window[1] & xyt$t<=time.window[2]],temporal.intensity))
den <- density
den$v <- den$v * approxscale
xyt$t <- as.integer(xyt$t)
numsamp <- min(c(200,xyt$n))
if(!suppresswarnings){
if (is.null(rvals)){
gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,...)
}
else{
gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals,...)
}
}
else{
if (is.null(rvals)){
suppressWarnings(gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,...))
}
else{
suppressWarnings(gin <- pcfinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals,...))
}
}
r <- gin$r
nam <- names(gin)
nam <- nam[nam!="r"]
tls <- sapply(time.window[1]:time.window[2],function(x){sum(xyt$t==x)})
tls <- (time.window[1]:time.window[2])[tls!=0]
ntls <- length(tls)
pb <- txtProgressBar(min=tls[1],max=rev(tls)[1],style=3)
if(!suppresswarnings){
pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(pcfinhom(xyt[xyt$t==t],lambda=den,r=r,...))})
}
else{
suppressWarnings(pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(pcfinhom(xyt[xyt$t==t],lambda=den,r=r,...))}))
}
close(pb)
li <- as.list(pcf[[1]])
ct <- 1
if(ntls>1){
for(i in 2:ntls){
if (!class(pcf[[i]])[1]=="try-error"){
li <- add.list(li,as.list(pcf[[i]]))
ct <- ct + 1
}
}
}
li <- smultiply.list(li,1/ct)
ginhom <- gin
for (n in nam){
idx <- which(names(ginhom)==n)
ginhom[[idx]] <- li[[idx]]
}
attr(ginhom,"correction") <- correction
cat("Returning an average of",ct,"curves\n")
return(ginhom)
}
##' KinhomAverage function
##'
##' A function to estimate the inhomogeneous K function for a spatiotemporal point process. The method of computation is similar to
##' \link{ginhomAverage}, see eq (8) Diggle P, Rowlingson B, Su T (2005) to see how this is computed.
##'
##' @param xyt an object of class stppp
##' @param spatial.intensity A spatialAtRisk object
##' @param temporal.intensity A temporalAtRisk object
##' @param time.window time interval contained in the interval xyt$tlim over which to compute average. Useful if there is a lot of data over a lot of time points.
##' @param rvals Vector of values for the argument r at which the inhmogeneous K function should be evaluated (see ?Kinhom). There is a sensible default.
##' @param correction choice of edge correction to use, see ?Kinhom, default is Ripley isotropic correction
##' @param suppresswarnings Whether or not to suppress warnings generated by Kinhom
##' @return time average of inhomogenous K function.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
KinhomAverage <- function(xyt,spatial.intensity,temporal.intensity,time.window=xyt$tlim,rvals=NULL,correction="iso",suppresswarnings=FALSE){
time.window <- sort(as.integer(time.window))
verifyclass(spatial.intensity,"spatialAtRisk")
verifyclass(temporal.intensity,"temporalAtRisk")
density <- as.im(spatial.intensity)
approxscale <- mean(sapply(xyt$t[xyt$t>=time.window[1] & xyt$t<=time.window[2]],temporal.intensity))
den <- density
den$v <- den$v * approxscale
xyt$t <- as.integer(xyt$t)
numsamp <- min(c(200,xyt$n))
if(!suppresswarnings){
if (is.null(rvals)){
Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den)
}
else{
Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals)
}
}
else{
if (is.null(rvals)){
suppressWarnings(Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den))
}
else{
suppressWarnings(Kin <- Kinhom(xyt[sample(1:xyt$n,numsamp,replace=FALSE)],lambda=den,r=rvals))
}
}
r <- Kin$r
nam <- names(Kin)
nam <- nam[nam!="r"]
tls <- sapply(time.window[1]:time.window[2],function(x){sum(xyt$t==x)})
tls <- (time.window[1]:time.window[2])[tls!=0]
ntls <- length(tls)
pb <- txtProgressBar(min=tls[1],max=rev(tls)[1],style=3)
if(!suppresswarnings){
pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(Kinhom(xyt[xyt$t==t],lambda=den,r=r))})
}
else{
suppressWarnings(pcf <- lapply(tls,function(t){setTxtProgressBar(pb,t);den<-density;den$v<-den$v*temporal.intensity(t);try(Kinhom(xyt[xyt$t==t],lambda=den,r=r))}))
}
close(pb)
li <- as.list(pcf[[1]])
ct <- 1
if(ntls>1){
for(i in 2:ntls){
if (!class(pcf[[i]])[1]=="try-error"){
li <- add.list(li,as.list(pcf[[i]]))
ct <- ct + 1
}
}
}
li <- smultiply.list(li,1/ct)
Kinhom <- Kin
for (n in nam){
idx <- which(names(Kinhom)==n)
Kinhom[[idx]] <- li[[idx]]
}
attr(Kinhom,"correction") <- correction
cat("Returning an average of",ct,"curves\n")
return(Kinhom)
}
##' spatialparsEst function
##'
##' Having estimated either the pair correlation or K functions using respectively \link{ginhomAverage} or \link{KinhomAverage}, the spatial
##' parameters sigma and phi can be estimated. This function provides a visual tool for this estimation procedure.
##'
##' To get a good choice of parameters, it is likely that the routine will have to be called several times in order to refine
##' the choice of sigma.range and phi.range.
##'
##' @param gk an R object; output from the function KinhomAverage or ginhomAverage
##' @param sigma.range range of sigma values to consider
##' @param phi.range range of phi values to consider
##' @param spatial.covmodel correlation type see ?CovarianceFct
##' @param covpars vector of additional parameters for certain classes of covariance function (eg Matern), these must be supplied in the order given in ?CovarianceFct
##' @param guess logical. Perform an initial guess at paramters? Alternative (the default) sets initial values in the middle of sigma.range and phi.range. NOTE: automatic parameter estimation can be can be unreliable.
##' @return rpanel function to help choose sigma nad phi by eye
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Baddeley AJ, Moller J, Waagepetersen R (2000). Non-and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329-350.
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{KinhomAverage}, \link{thetaEst}, \link{lambdaEst}, \link{muEst}
##' @export
spatialparsEst <- function(gk,sigma.range,phi.range,spatial.covmodel,covpars=c(),guess=FALSE){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
corchoice <- c() # to get rid of 'no visible binding' messages on checking
sigma <- c()
phi <- c()
sigmamin <- sigma.range[1]
sigmamax <- sigma.range[2]
phimin <- phi.range[1]
phimax <- phi.range[2]
idx <- which(names(gk)==attr(gk,"correction"))
env <- NULL
if (attr(gk,"fname") == "g[inhom]"){
panfun <- function(p){
env <<- environment()
r <- gk$r
egr <- suppressWarnings(exp(sapply(r,gu,sigma=p$sigma,phi=p$phi,model=spatial.covmodel,additionalparameters=covpars))-1)
gvals <- gk[[idx]]
if (p$transform!="none"){
if (p$transform=="log"){
gvals <- log(gk[[idx]])
egr <- log(egr)
}
}
plot(NULL,main="g[inhom]",xlim=c(0,max(r,na.rm=TRUE)),ylim=c(0,max(c(gvals[!is.infinite(gvals)&!is.nan(gk[[idx]])],egr[!is.infinite(egr)&!is.nan(egr)]))),xlab="r",ylab="g_inhom(r)",sub=paste(spatial.covmodel,"covariance function"))
lines(r,gvals)
lines(r,egr,col="orange")
legend("topright",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
return(p)
}
}
else if (attr(gk,"fname") == "K[inhom]"){
panfun <- function(p){
env <<- environment()
r <- gk$r
egr <- suppressWarnings(exp(sapply(r,gu,sigma=p$sigma,phi=p$phi,model=spatial.covmodel,additionalparameters=covpars))-1)
rdiff <- diff(r[1:2]) # do the integral on the discrete partition of r given by Kinhom
kest <- rdiff * egr[1] * 2 * pi * r[1]
for(i in 2:length(r)){
kest[i] <- kest[i-1] + rdiff * egr[i] * 2 * pi * r[i]
}
gvals <- gk[[idx]]
if (p$transform!="none"){
if (p$transform=="^1/4"){
gvals <- (gk[[idx]])^(1/4)
kest <- kest^(1/4)
}
}
plot(NULL,main="K[inhom]",xlim=c(0,max(r,na.rm=TRUE)),ylim=c(0,max(c(gvals[!is.infinite(gvals)],kest[!is.infinite(kest)]))),xlab="r",ylab="K_inhom(r)",sub=paste(spatial.covmodel,"covariance function"))
lines(r,gvals)
lines(r,kest,col="orange")
legend("topleft",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
return(p)
}
}
else{
stop("Argument gk is not of correct class, see ?ginhomAverage, ?KinhomAverage")
}
pancontrol <- rp.control("Parameter Estimation",aschar=FALSE)
##rp.radiogroup(pancontrol,var=corchoice,values=c("exponential","matern","whittle"),action=panfun,initval="exponential")
if (attr(gk,"fname") == "g[inhom]"){
rp.radiogroup(pancontrol,variable=transform,vals=c("none","log"),action=panfun,initval="none")
}
else if (attr(gk,"fname") == "K[inhom]"){
rp.radiogroup(pancontrol,variable=transform,vals=c("none","^1/4"),action=panfun,initval="^1/4")
}
# "quick" optimiser to get poor initial values for sigma and phi
r <- gk$r
if(guess){
if (attr(gk,"fname") == "g[inhom]"){
spfun <- function(sigmaphi){
egr <- suppressWarnings(exp(sapply(r,gu,sigma=sigmaphi[1],phi=sigmaphi[2],model=spatial.covmodel,additionalparameters=covpars))-1)
S <- suppressWarnings((egr-gk[[idx]])^2)
S[is.infinite(S)] <- NA
W <- 1/r^2
W[is.infinite(W)] <- NA
sidx <- !(is.na(W) | is.na(S))
ans <- sum((W*S)[sidx],na.rm=TRUE)/sum(W[sidx])
return(ans)
}
ops <- optim(c(sigmamin + (sigmamax-sigmamin)/2,phimin + (phimax-phimin)/2),spfun)
sigmainit <- ops$par[1]
phiinit <- ops$par[2]
}
else if (attr(gk,"fname") == "K[inhom]"){
sigphifun <- function(sigphi){
egr <- suppressWarnings(exp(sapply(r,gu,sigma=sigphi[1],phi=sigphi[2],model=spatial.covmodel,additionalparameters=covpars))-1)
rdiff <- diff(r[1:2]) # do the integral on the discrete partition of r given by Kinhom
kest <- rdiff * egr[1] * 2 * pi * r[1]
for(i in 2:length(r)){
kest[i] <- kest[i-1] + rdiff * egr[i] * 2 * pi * r[i]
}
gvals <- gk[[idx]]
gvals <- (gk[[idx]])^(1/4)
kest <- kest^(1/4)
gvals[is.infinite(gvals)] <- NA
kest[is.infinite(kest)] <- NA
S <- suppressWarnings((kest-gvals)^2)
S[is.infinite(S)] <- NA
W <- 1/r^2
W[is.infinite(W)] <- NA
sidx <- !(is.na(W) | is.na(S))
ans <- sum((W*S)[sidx],na.rm=TRUE)/sum(W[sidx])
return(ans)
}
ops <- optim(c(sigmamin + (sigmamax-sigmamin)/2,phimin + (phimax-phimin)/2),sigphifun)
sigmainit <- ops$par[1]
phiinit <- ops$par[2]
}
}
else{
sigmainit <- sigmamin + (sigmamax-sigmamin)/2
phiinit <- phimin + (phimax-phimin)/2
}
rp.slider(pancontrol,variable=sigma,from=sigmamin,to=sigmamax,initval=sigmainit,action=panfun,showvalue=TRUE)
rp.slider(pancontrol,variable=phi,from=phimin,to=phimax,initval=phiinit,action=panfun,showvalue=TRUE)
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
# construct a geometry string with the current height and new width:
hhh = as.numeric(tkwinfo("height",pancontrol$.handle))+50
www = 300
ggg = sprintf("%dx%d",www,hhh)
# resize
tkwm.geometry(pancontrol$.handle,ggg)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
## get the variables. Two points:
sigma <- get("p",envir=env)$sigma #as.numeric(tclvalue(pancontrol$sigma.tcl))
phi <- get("p",envir=env)$phi #as.numeric(tclvalue(pancontrol$phi.tcl))
return(list(sigma=sigma,phi=phi))
}
##' Cvb function
##'
##' This function is used in \code{thetaEst} to estimate the temporal correlation parameter, theta.
##'
##' @param xyt object of class stppp
##' @param spatial.intensity bivariate density estimate of lambda, an object of class im (produced from density.ppp for example)
##' @param N number of integration points
##' @param spatial.covmodel spatial covariance model
##' @param covpars additional covariance parameters
##' @return a function, see below.
##' Computes Monte carlo estimate of function C(v;beta) in Brix and Diggle 2001 pp 829 (... note later corrigendum to paper (2003) corrects the expression given in this paper)
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' }
##' @seealso \link{thetaEst}
##' @export
Cvb <- function(xyt,spatial.intensity,N=100,spatial.covmodel,covpars){
verifyclass(spatial.intensity,"im")
sar <- spatialAtRisk(list(X=spatial.intensity$xcol,Y=spatial.intensity$yrow,Zm=t(spatial.intensity$v)))
gsx <- length(xvals(sar))
gsy <- length(yvals(sar))
xy <- cbind(rep(xvals(sar),gsy),rep(yvals(sar),each=gsx))
wt <- as.vector(zvals(sar))
wt[is.na(wt)] <- 0
sidx <- sample(1:(gsx*gsy),N,prob=wt)
xy <- xy[sidx,]
pd <- as.vector(pairdist(xy))
cvb <- function(nu,sigma,phi,theta){
return(mean(exp(exp(-nu*theta)*gu(pd,sigma=sigma,phi=phi,model=spatial.covmodel,additionalparameters=covpars))-1))
}
return(cvb)
}
##' thetaEst function
##'
##' A tool to visually estimate the temporal correlation parameter theta; note that sigma and phi must have first been estiamted.
##'
##' @param xyt object of class stppp
##' @param spatial.intensity A spatial at risk object OR a bivariate density estimate of lambda, an object of class im (produced from density.ppp for example),
##' @param temporal.intensity either an object of class temporalAtRisk, or one that can be coerced into that form. If NULL (default), this is estimated from the data, seee ?muEst
##' @param sigma estimate of parameter sigma
##' @param phi estimate of parameter phi
##' @param theta.range range of theta values to consider
##' @param N number of integration points in computation of C(v,beta) (see Brix and Diggle 2003, corrigendum to Brix and Diggle 2001)
##' @param spatial.covmodel spatial covariance model
##' @param covpars additional covariance parameters
##' @return An r panel tool for visual estimation of temporal parameter theta
##' NOTE if lambdaEst has been invoked to estimate lambda, then the returned density should be passed to thetaEst as the argument spatial.intensity
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{lambdaEst}, \link{muEst}
##' @export
thetaEst <- function(xyt,spatial.intensity=NULL,temporal.intensity=NULL,sigma,phi,theta.range=c(0,10),N=100,spatial.covmodel="exponential",covpars=c()){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
if (inherits(spatial.intensity,"spatialAtRisk")){
spatial.intensity <- as.im(spatial.intensity)
}
if (is.null(spatial.intensity)){
spatial.intensity <- density(xyt)
}
if (is.null(temporal.intensity)){
temporal.intensity <- muEst(xyt)
}
else{
if (!inherits(temporal.intensity,"temporalAtRisk")){
temporal.intensity <- temporalAtRisk(temporal.intensity,tlim=xyt$tlim,xyt=xyt)
}
if(!all(xyt$tlim==attr(temporal.intensity,"tlim"))){
stop("Incompatible temporal.intensity, integer time limits (tlim and temporal.intensity$tlim) do not match")
}
}
thetamin <- theta.range[1]
thetamax <- theta.range[2]
tlim <- as.integer(xyt$tlim)
cvb <- Cvb(xyt=xyt,spatial.intensity=spatial.intensity,N=N,spatial.covmodel=spatial.covmodel,covpars=covpars)
xytlim <- sort(as.integer(xyt$tlim))
tvec <- xytlim[1]:xytlim[2]
if(length(tvec)==1){
stop("Insufficiently many time intervals: type as.integer(xyt$tlim) into console")
}
theta <- c() # to get rid of 'no visible binding' messages on checking
panfun <- function(p){
sigma <- as.numeric(p$sigma)
phi <- as.numeric(p$phi)
theta <- p$theta
uqt <- as.numeric(names(table(as.integer(xyt$t))))
tfit <- sapply(uqt,temporal.intensity)
autocov <- acf(table(as.integer(xyt$t)) - tfit,type="covariance",plot=FALSE)
r <- 0:(length(autocov$acf)-1)
plot(r,autocov$acf,type="l",main="Autocovariance",xlab="lag",ylab="acov(lag)")
legend("topright",lty=c("solid","solid"),col=c("black","orange"),legend=c("Empirical","Theoretical"))
abline(h=0)
rseq <- seq(r[1],r[length(r)],length.out=100)
theo <- sapply(rseq,cvb,sigma=sigma,phi=phi,theta=theta)
lines(rseq,theo*(autocov$acf[1]/theo[1]),type="l",col="orange")
return(p)
}
pancontrol <- rp.control("Parameter Estimation",aschar=FALSE)
rp.textentry(pancontrol,variable=sigma,initval=sigma,action=panfun)
rp.textentry(pancontrol,variable=phi,initval=phi,action=panfun)
rp.slider(pancontrol,variable=theta,from=thetamin,to=thetamax,initval=(thetamax-thetamin)/2,action=panfun,showvalue=TRUE)
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
## get the variables. Two points:
theta <- as.numeric(tclvalue(pancontrol$theta.tcl))
return(theta)
}
##' lambdaEst function
##'
##' Generic function for estimating bivariate densities by eye. Specific methods exist for stppp objects and ppp objects.
##'
##' @param xyt an object
##' @param ... additional arguments
##' @return method lambdaEst
##' @seealso \link{lambdaEst.stppp}, \link{lambdaEst.ppp}
##' @export
lambdaEst <- function(xyt,...){
UseMethod("lambdaEst")
}
##' lambdaEst.stppp function
##'
##' A tool for the visual estimation of lambda(s) via a 2 dimensional smoothing of the case locations. For parameter estimation, the alternative is
##' to estimate lambda(s) by some other means, convert it into a spatialAtRisk object and then into a pixel image object using the build in coercion
##' methods, this \code{im} object can then be fed to \link{ginhomAverage}, \link{KinhomAverage} or \link{thetaEst} for instance.
##'
##' The function lambdaEst is built directly on the density.ppp function and as such, implements a bivariate
##' Gaussian smoothing kernel. The bandwidth is initially that which is automatically chosen by the default method
##' of density.ppp. Since image plots of these kernel density estimates may not have appropriate
##' colour scales, the ability to adjust this is given with the slider 'colour adjustment'. With colour adjustment set
##' to 1, the default image.plot for the equivalent pixel image object is shown and for values less than 1, the colour
##' scheme is more spread out, allowing the user to get a better feel for the density that is being fitted. NOTE: colour
##' adjustment does not affect the returned density and the user should be aware that the returned density will 'look like'
##' that displayed when colour adjustment is set equal to 1.
##'
##'
##' @method lambdaEst stppp
##' @param xyt object of class stppp
##' @param weights Optional vector of weights to be attached to the points. May include negative values. See ?density.ppp.
##' @param edge Logical flag: if TRUE, apply edge correction. See ?density.ppp.
##' @param bw optional bandwidth. Set to NULL by default, which calls teh resolve.2D.kernel function for computing an initial value of this
##' @param ... arguments to be passed to plot
##' @return This is an rpanel function for visual choice of lambda(s), the output is a variable, varname, with the density *per unit time*
##' the variable varname can be fed to the function ginhomAverage or KinhomAverage as the argument density (see for example ?ginhomAverage), or into the
##' function thetaEst as the argument spatial.intensity.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{spatialAtRisk}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{muEst}
##' @export
lambdaEst.stppp <- function(xyt,weights=c(),edge=TRUE,bw=NULL,...){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
adjust <- c() # to get rid of 'no visible binding' messages on checking
plotdata <- c()
pow <- c()
if(is.null(bw)){
bw <- resolve.2D.kernel(NULL,x = xyt,adjust = 1)$sigma
}
varlist <- list(bw=bw,adjust=1,plotdata=TRUE,pow=1)
env <- NULL
panfun <- function(p){
env <<- environment()
bw <- as.numeric(p$bw)
adjust <- as.numeric(p$adjust)
d <- density(xyt,bandwidth=bw,weights=weights,edge=edge)
plotden <- d
plotden$v <- (d$v / diff(xyt$tlim))^p$pow
plot(plotden,main="Density",...)
if(p$plotdata){
points(xyt,pch="+",cex=0.5)
}
return(p)
}
panfun(varlist)
pancontrol <- rp.control("Density Estimation")
rp.textentry(pancontrol,variable=bw,initval=bw,action=panfun,title="bandwidth")
rp.checkbox(pancontrol,variable=plotdata,action=panfun,labels = "Plot Data",initval = TRUE)
rp.slider(pancontrol,variable=pow,from=0,to=1,initval=1,action=panfun,showvalue=TRUE,title="colour adjustment")
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
return(get("d",envir=env))
}
##' lambdaEst.ppp function
##'
##' A tool for the visual estimation of lambda(s) via a 2 dimensional smoothing of the case locations. For parameter estimation, the alternative is
##' to estimate lambda(s) by some other means, convert it into a spatialAtRisk object and then into a pixel image object using the build in coercion
##' methods, this \code{im} object can then be fed to \link{ginhomAverage}, \link{KinhomAverage} or \link{thetaEst} for instance.
##'
##' The function lambdaEst is built directly on the density.ppp function and as such, implements a bivariate
##' Gaussian smoothing kernel. The bandwidth is initially that which is automatically chosen by the default method
##' of density.ppp. Since image plots of these kernel density estimates may not have appropriate
##' colour scales, the ability to adjust this is given with the slider 'colour adjustment'. With colour adjustment set
##' to 1, the default image.plot for the equivalent pixel image object is shown and for values less than 1, the colour
##' scheme is more spread out, allowing the user to get a better feel for the density that is being fitted. NOTE: colour
##' adjustment does not affect the returned density and the user should be aware that the returned density will 'look like'
##' that displayed when colour adjustment is set equal to 1.
##'
##'
##' @method lambdaEst ppp
##' @param xyt object of class stppp
##' @param weights Optional vector of weights to be attached to the points. May include negative values. See ?density.ppp.
##' @param edge Logical flag: if TRUE, apply edge correction. See ?density.ppp.
##' @param bw optional bandwidth. Set to NULL by default, which calls teh resolve.2D.kernel function for computing an initial value of this
##' @param ... arguments to be passed to plot
##' @return This is an rpanel function for visual choice of lambda(s), the output is a variable, varname, with the density *per unit time*
##' the variable varname can be fed to the function ginhomAverage or KinhomAverage as the argument density (see for example ?ginhomAverage), or into the
##' function thetaEst as the argument spatial.intensity.
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{spatialAtRisk}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{muEst}
##' @export
lambdaEst.ppp <- function(xyt,weights=c(),edge=TRUE,bw=NULL,...){
## dummy function callback needed for OK button:
ok <- function(panel){
return(panel)
}
adjust <- c() # to get rid of 'no visible binding' messages on checking
plotdata <- c()
pow <- c()
if(is.null(bw)){
bw <- resolve.2D.kernel(NULL,x = xyt,adjust = 1)$sigma
}
varlist <- list(bw=bw,adjust=1,plotdata=TRUE,pow=1)
env <- NULL
panfun <- function(p){
env <<- environment()
bw <- as.numeric(p$bw)
adjust <- as.numeric(p$adjust)
d <- density(xyt,sigma=bw,weights=weights,edge=edge)
plotden <- d
plotden$v <- d$v^p$pow
plot(plotden,main="Density",...)
if(p$plotdata){
points(xyt,pch="+",cex=0.5)
}
return(p)
}
panfun(varlist)
pancontrol <- rp.control("Density Estimation")
rp.textentry(pancontrol,variable=bw,initval=bw,action=panfun,title="bandwidth")
rp.checkbox(pancontrol,variable=plotdata,action=panfun,labels = "Plot Data",initval = TRUE)
rp.slider(pancontrol,variable=pow,from=0,to=1,initval=1,action=panfun,showvalue=TRUE,title="colour adjustment")
## add our OK button. Make it a quitbutton:
rp.button(pancontrol,action=ok,title="OK",quitbutton=TRUE)
## now wait until our panel quits.
rp.block(pancontrol)
dev.off()
return(get("d",envir=env))
}
##' muEst function
##'
##' Computes a non-parametric estimate of mu(t). For the purposes of performing prediction, the alternatives are: (1) use a parameteric model as in Diggle P, Rowlingson B, Su T (2005),
##' or (2) a \link{constantInTime} model.
##'
##' @param xyt an stppp object
##' @param ... additional arguments to be passed to lowess
##' @return object of class temporalAtRisk giving the smoothed mut using the lowess function
##' @references
##' \enumerate{
##' \item Benjamin M. Taylor, Tilman M. Davies, Barry S. Rowlingson, Peter J. Diggle (2013). Journal of Statistical Software, 52(4), 1-40. URL http://www.jstatsoft.org/v52/i04/
##' \item Brix A, Diggle PJ (2001). Spatiotemporal Prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society, Series B, 63(4), 823-841.
##' \item Diggle P, Rowlingson B, Su T (2005). Point Process Methodology for On-line Spatio-temporal Disease Surveillance. Environmetrics, 16(5), 423-434.
##' }
##' @seealso \link{temporalAtRisk}, \link{constantInTime}, \link{ginhomAverage}, \link{KinhomAverage}, \link{spatialparsEst}, \link{thetaEst}, \link{lambdaEst}
##' @export
muEst <- function(xyt,...){
verifyclass(xyt,"stppp")
t <- as.integer(xyt$t)
tlim <- as.integer(xyt$tlim)
tseq <- tlim[1]:tlim[2]
counts <- sapply(tseq,function(x){sum(t==x)})
sm <- lowess(tseq,sqrt(counts),...)
return(temporalAtRisk(sm$y^2,tlim=tlim,xyt=xyt))
}
##' density.stppp function
##'
##' A wrapper function for \link{density.ppp}.
##'
##' @method density stppp
##' @param x an stppp object
##' @param bandwidth 'bandwidth' parameter, equivanent to parameter sigma in ?density.ppp ie standard deviation of isotropic Gaussian smoothing kernel.
##' @param ... additional arguments to be passed to density.ppp
##' @return bivariate density estimate of xyt; not this is a wrapper function for density.ppp
##' @seealso \link{density.ppp}
##' @export
density.stppp <- function(x,bandwidth=NULL,...){
density.ppp(x,...,sigma=bandwidth)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.