PenalizedSecondOrderCCL: Second order conditional composite likelihood for...

In kristianhessellund/Multilogreg: Semi-parametric models for multivariate point pattern data

Description

This is a conditional composite likelihood method for semi-parametric estimation of multivariate log Gaussian Cox processes (LGCP).

Usage

PenalizedSecondOrderCCL(
  X = NULL,
  psi0 = NULL,
  xi0 = NULL,
  sigma0 = NULL,
  phi0 = NULL,
  Beta,
  xibound = c(1e-05, Inf),
  sigmabound = c(1e-05, Inf),
  phibound = c(1e-05, Inf),
  lamb = 0,
  lat = 1,
  Rmax = NULL,
  covariate = NULL,
  ite = 100,
  tol = 1e-05,
  tol.alpha = 1e-10,
  tol.eta = 1e-10,
  mu = 1,
  prints = T
)

Arguments

X	A multivariate log Gaussian Cox process. Must be of class ppp.
psi0	Optional. Initial value for psi. If psi0 = NULL, then psi0 is simulated from Unif(-0.25,0.25).
xi0	Optional. Initial value xi. If xi0 = NULL, then xi0 is simulated using the size of the observation window.
sigma0	Optional. Initial value sigma. If sigma0 = NULL, then sigma0 is simulated from Unif(0.4,0.6).
phi0	Optional. Initial value phi. If phi0 = NULL, then phi0 is simulated using the size of the observation window.
Beta	Estimated first order parameters. Beta must be a matrix. The number of rows must correspond to the number of covariates. The number of columns must correspond to the number of point types.
xibound	Optional. Parameter space for xi. Default is xibound = (1e-5,Inf).
sigmabound	Optional. Boundary for sigma. Default is sigmabound = (1e-5,Inf).
phibound	Optional. Boundary for phi. Default is phi.bound = (1e-5,Inf).
lamb	Optional. lasso parameter. Default is lamb = 0. Can be a sequence of lambda-values.
lat	Number of common latent fields in the model. Default is lat = 1.
Rmax	Maximal distance between distinct pairs of points. If Rmax = NULL, then Rmax is determined by the observation window.
covariate	Covariates. Covariates must be a matrix. The rows corresponds to the points in the point pattern and the columns indicates the corresponding covariates that are observed at the location of the point pattern.
ite	Maximal number of iterations. Default is ite = 100.
tol	Convergence parameter for likelihood function. The default is tol = 1e-5.
tol.alpha	Convergence parameter for alpha. The default is tol.alpha = 1e-10.
tol.eta	Convergence parameter for eta. The default is tol.eta = 1e-10.
mu	Contraint parameter. The default is mu = 1.
prints	Optional. Default prints = TRUE. Printing value of likelihood function along with convergence rate for each iteration.

Details

To be updated.

Value

Return estimated parameters of multivariate LGCP X

Author(s)

Kristian Bjørn Hessellund, Ganggang Xu, Yongtao Guan and Rasmus Waagepetersen.

References

Hessellund, K. B., Xu, G., Guan, Y. and Waagepetersen, R. (2020) Second order semi-parametric inference for multivariate log Gaussian Cox processes.

Examples

### First order analysis of Washington D.C. street crimes ###

nspecies <- length(unique(markedprocess$marks))
ncovs <- dim(Z_process)[2]

# Initial parameters for street crime data
Beta0=matrix(runif(nspecies*ncovs,-0.5,0.5),nrow =nspecies)

# Choosing the last type of street crime as the baseline
Beta0 <- as.matrix(t(scale(Beta0,center=Beta0[nspecies,],scale=FALSE)))

# Parameter estimation
Betahat=FirstOrderCCL(X=markedprocess,Beta0 = Beta0,covariate = Z_process)

# Parameter estimates
colnames(Betahat$betahat)=c("Burglary", "Assault w. weapon", "Motor v. theft",
"Theft f. auto", "Robbery", "Other theft (baseline)")
rownames(Betahat$betahat)=c("Intercept", "%African", "%Hispanic", "%Males",
 "Median income", "%Household","%Bachelor", "Distance to police station")
Betahat

### Non-parametric estimates of cross PCF ratios ###

rseq = seq(0,3000,length=100)
cpcf=CrossPCF(X = markedprocess,covariate = Z_process,Beta = Betahat$betahat,r = rseq
,bwd = 200,cut = NULL)

plot(rseq,rep(1,100),type="l",col=1,ylim=c(0,2),xlab="r",ylab="",main="PCFs")
for(i in 1:(nspecies-1)){
  lines(rseq,cpcf$cPCF.m[i,i,],,col=i)
}

plot(rseq,rep(1,100),type="l",col=1,ylim=c(0,2),xlab="r",ylab="",main="PCFs")
for(i in 1:(nspecies-1)){
  for(j in 1:(nspecies-1)){
    if(i==j){next}
    lines(rseq,cpcf$cPCF.m[i,j,],,col=1)
  }
}

### Standard errors for Beta ###

var.est=Firstorder_var(X=markedprocess,covariate = Z_process,Beta=Betahat$betahat,
r = seq(0,3000,length=100),bwd = 200,cPCF = cpcf$cPCF.m,method = "pool")

std.error=matrix(sqrt(diag(var.est$COV)),nrow=ncovs)
colnames(std.error)=c("Burglary", "Assault w. weapon", "Motor v. theft","Theft f. auto", "Robbery")
rownames(std.error)=c("Intercept", "%African", "%Hispanic", "%Males",
 "Median income", "%Household","%Bachelor", "Distance to police station")
std.error

### Second order analysis of Washington D.C. street crimes ###

# Globals
len.x=markedprocess$window$xrange[2]-markedprocess$window$xrange[1]
len.y=markedprocess$window$yrange[2]-markedprocess$window$yrange[1]
lw=(len.x+len.y)/2
latent = 1
Rmax   = 1000
rl     = 100
r      = seq(0,Rmax,length=rl)

# Parameter spaces (upper,lower)
xib = c(lw/10,lw/500)
sib = c(10,0.01)
phb = c(lw/10,lw/1000)

notrun = TRUE
if(notrun ==FALSE){
# Parameter estimation
par.est = PenalizedSecondOrderCCL(X=markedprocess,covariate = Z_process,Beta = Betahat$betahat,
Rmax = Rmax,xibound=xib,sigmabound = sib,phibound=phb, lat = latent)

# Parameter estimates
alp.est = matrix(par.est[[1]]$alpha[,,1])
xi.est  = par.est[[1]]$xi
sig.est = par.est[[1]]$sigma
phi.est = par.est[[1]]$phi

## Pair correlation functions (PCF)
pcf.est=matrix(NA,nrow=nspecies,ncol=rl)
for(i in 1:nspecies){
  tmp=0
  for(k in 1:latent){tmp = tmp + alp.est[i,k]*alp.est[i,k]*exp(-r/xi.est[k])}
  pcf.est[i,]=exp(tmp+sig.est[i]^2*exp(-r/phi.est[i]))
}

## Plot of PCFs
lims=c(min(pcf.est,0.9),max(pcf.est,1.1))
plot(r,rep(1,rl),type="l",ylim=lims,main="PCFs",ylab="", lwd=2)
for(j in 1:nspecies){lines(r,pcf.est[j,],col=j+1,lwd=2)}
txt.legend=c(as.expression(bquote(g[11])),as.expression(bquote(g[22])),
as.expression(bquote(g[33])),as.expression(bquote(g[44])),
as.expression(bquote(g[55])),as.expression(bquote(g[66])))
legend("topright",legend =  txt.legend,lty=1,col=c(2,3,4,5,6,7), cex=1,lwd=2)

# cross-PCFs
cross.pcf.est=matrix(NA,nrow=nspecies*(nspecies-1)/2,ncol=rl)
l=0
for(i in 1:nspecies){
 for(j in i:nspecies){
   if(i==j){next}
   l=l+1
   tmp=0
   for(k in 1:latent){tmp = tmp + alp.est[i,k]*alp.est[j,k]*exp(-r/xi.est[k])}
   cross.pcf.est[l,]=exp(tmp)
  }
}

# Plot of cross-PCFs
ylims=c(min(0.9,cross.pcf.est),max(1.1,cross.pcf.est))
ltys=c(rep(1,8),rep(2,8))
plot(r,rep(1,rl),type="l",main="cross-PCFs",ylab="",ylim=ylims)
for(i in 1:(nspecies*(nspecies-1)/2)){lines(r,cross.pcf.est[i,],col=i+1,lty=ltys[i],lwd=2)}
txt.legend=c(as.expression(bquote(g[12])),as.expression(bquote(g[13])),
as.expression(bquote(g[14])),as.expression(bquote(g[15])),as.expression(bquote(g[16])),
as.expression(bquote(g[23])),as.expression(bquote(g[24])),as.expression(bquote(g[25])),
as.expression(bquote(g[26])),as.expression(bquote(g[34])),as.expression(bquote(g[35])),
as.expression(bquote(g[36])),as.expression(bquote(g[45])),as.expression(bquote(g[46])),
as.expression(bquote(g[56])))
legend("topright",legend =  txt.legend,lty=c(rep(1,8),rep(2,8)),col=2:17,cex=0.7,bg="white",lwd=2)
}

kristianhessellund/Multilogreg documentation built on Jan. 1, 2021, 7:23 a.m.