source("SemiparametricTest_revision_final.R") library("spatstat") load("dc.covars.Rda") load("DCcrime.Rda") DC.grids <- read.csv("DCgrids.csv",header=T) grid.mark_block <- ppp(x=DC.grids$Easting, y=DC.grids$Northing,marks=DC.grids$census_block,window=win) grid.mark_blockg <- ppp(x=DC.grids$Easting, y=DC.grids$Northing,marks=DC.grids$census_blockgroup,window=win) grid.mark_ctract <- ppp(x=DC.grids$Easting, y=DC.grids$Northing,marks=DC.grids$census_tract,window=win) rmax.pre <-600 r <- seq(0,rmax.pre,l=401) ###if the number of knots is too small, there maybe some numerical issues pd.K <- 100 degree <- 3 for(model in c("Full","Reduced","Improved")){ set.seed(1) if(model=="Full") ##Full model for Assult object <- kppm(Assult~black.per+bachelor.per+poverty.per+Onpremisedist+Liquordist+Grocerydist+Bdist1+pop.den_blockg+I(black.per_block*populus_block)+I(house.per_block*populus_block)+area_block+I(area_block^2)+I(pop.den_blockg^2),covariates=dc.covars, clusters="Thomas") if(model=="Reduced") ##Reduced model for Assult object <- kppm(Assult~black.per+bachelor.per+poverty.per+Onpremisedist+Liquordist+Grocerydist+Bdist1+pop.den_blockg+I(black.per_block*populus_block)+area_block,covariates=dc.covars, clusters="Thomas") if(model=="Improved") ##Improved model for Assult object <- kppm(Assult~black.per+bachelor.per+Licensedist+Bdist1+pop.den_blockg+I(black.per_block*populus_block)+area_block,covariates=dc.covars, clusters="Thomas") beta.pre <- object$po$coef wfun <- function(dist){dist} R_range <- c(10,pi/median(nndist(Assult))*rmax.pre) obj.cv <- NonPCF_pd_MCE_cv_parallel(object,beta.pre,r,degree,grid.num=NULL,rescale=TRUE,R_range,n.R=10,Iter=1,n.rcv=100,wfun,d.out=100,err=10^(-3),numCores=10) par(mfrow=c(2,2)) plot(obj.cv$R_seq,obj.cv$cv.score,type="l",lwd=2) plot(obj.cv$R_seq,obj.cv$lscv.score,type="l",lwd=2) pd.R <- obj.cv$R_opt.ls pd.K <- obj.cv$K_opt.ls obj.pcf.ls <- NonPCF_pd_MCE(object,beta.pre,r,R=pd.R,degree,K=pd.K,grid.num=NULL,rescale=TRUE,Iter=1,wfun,d.out=100,err=10^(-4)) plot(r,obj.pcf.ls$gfun.pd(r),lwd=2,col="red",type="l") object.ls <- improve.kppm_new(object,fast=TRUE,eps.rmax = 0.01,gfun=obj.pcf.ls$gfun.pd,cov=TRUE,rmax.pre=rmax.pre,covariate.center=TRUE) betathat.ls <- object.ls$po$coef se.ls <- sqrt(diag(object.ls$improve$vcov)) betathat.ls se.ls pvalue.ls <- 2*(1-pnorm(abs(betathat.ls)/se.ls)) pvalue.ls mu <- object.ls$improve$mu/object.ls$improve$wt U <- object.ls$improve$U if(model=="Reduced"){ mu_reduced <- as.im(data.frame("x"= U$x, "y" = U$y, "value" =log(mu))) gfun_ass_red <- obj.pcf.ls$gfun.pd } if(model=="Improved"){ mu_improved <- as.im(data.frame("x"= U$x, "y" = U$y, "value" =log(mu))) gfun_ass_imp <- obj.pcf.ls$gfun.pd } ##census block level pvlaues r1.ls <- c(LocalKppm_Firstorder_parallel(object.ls,numCores=50,grid.mark=grid.mark_block)) ##census block-group level pvlaues r2.ls <- c(LocalKppm_Firstorder_parallel(object.ls,numCores=50,grid.mark=grid.mark_blockg)) ##census tract level pvlaues r3.ls <- c(LocalKppm_Firstorder_parallel(object.ls,numCores=50,grid.mark=grid.mark_ctract)) names(r1.ls) <-names(r2.ls) <-names(r3.ls) <- c("Semi-L", "Semi-G", "Semi-C") if(model=="Improved") save.image("Assult_imp.RData") if(model=="Reduced") save.image("Assult_red.RData") if(model=="Full") save.image("Assult_full.RData") }