#####################################################
################ SVC modeling #######################
#####################################################
require(spdep)
data(boston)
y <- boston.c[, "CMEDV"]
x <- boston.c[,c("CRIM", "AGE")]
xconst <- boston.c[,c("ZN","DIS","RAD","NOX", "TAX","RM", "PTRATIO", "B")]
xgroup <- boston.c[,"TOWN"]
coords <- boston.c[,c("LON", "LAT")]
meig <- meigen(coords=coords)
# meig <- meigen_f(coords=coords) ## for large samples
#####################################################
####### Gaussian regression with SVC ################
res <- resf_vc(y=y,x=x,xconst=xconst,meig=meig )
res
plot_s(res,0) # Spatially varying intercept
plot_s(res,1) # 1st SVC (Not shown because the SVC is estimated constant)
plot_s(res,2) # 2nd SVC
#### For large samples (e.g., n > 5,000), the following
#### additional learning often improves the modeling accuracy
# res_adj<- addlearn_local(res)
# res_adj
# plot_s(res_adj,0)
# plot_s(res_adj,1)
# plot_s(res_adj,2)
#### Group-level SVC (s_id) + random intercepts (xgroup)
# meig_g <- meigen(coords, s_id=xgroup)
# res2 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig_g,xgroup=xgroup)
#####################################################
####### Gaussian regression with SVC + NVC ##########
# res3 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_nvc =TRUE)
# plot_s(res3,0) # Spatially varying intercept
# plot_s(res3,1) # Spatial plot of the varying coefficient (SVC + NVC) on x[,1]
# plot_s(res3,1,btype="svc")# Spatial plot of SVC in the coefficient
# plot_s(res3,1,btype="nvc")# Spatial plot of NVC in the coefficient
# plot_n(res3,1) # 1D plot of the NVC
#####################################################
######## Non-Gaussian regression with SVC ###########
#### Model for non-Gaussian continuous data
# - Probability distribution is estimated from data
# ng4 <- nongauss_y( tr_num = 2 )# 2 SAL transformations to Gaussianize y
# res4 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng4 )
# res4 # tr_num may be selected by comparing BIC
# coef_marginal_vc(res4) # marginal effects from x (dy/dx)
# plot(res4$pdf,type="l") # Estimated probability density function
# res4$skew_kurt # Skew and kurtosis of the estimated PDF
# res4$pred_quantile[1:2,]# predicted value by quantile
#### Model for non-Gaussian and non-negative continuous data
# - Probability distribution is estimated from data
## 2 SAL trans. + 1 Box-Cox trans. to Gaussianize y
# ng5 <- nongauss_y( tr_num = 2, y_nonneg = TRUE )
# res5 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng5 )
# coef_marginal_vc(res5)
#### Overdispersed Poisson model for count data
# - y: count data
#ng6 <- nongauss_y( y_type = "count" )
#res6 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng6 )
#### General model for count data
# - y: count data
# - Probability distribution is estimated from data
#ng7 <- nongauss_y( y_type = "count", tr_num = 2 )
#res7 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng7 )
#####################################################
################ STVC modeling ######################
#####################################################
# See \url{https://github.com/dmuraka/spmoran}
#require(spData)
#data(house)
#dat0 <- st_as_sf(house)
#dat <- data.frame(st_coordinates(dat0), dat0)
#y <- log(dat[,"price"])
#x <- dat[,c("lotsize","TLA")]
#xconst <- dat[,c("rooms","beds")]
#byear <- house$yrbuilt
#syear <- as.numeric(as.character(house$syear))#factor -> numeric
#coords_z<- cbind(byear,syear)
#meig <- meigen_f(coords=coords, coords_z=cbind(byear,syear),interact=TRUE)
#res8 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig )
#res8
## Varying intercept for byear <=1950 and syear==1998
#plot_s(res8,0, coords_z1_lim=c(-Inf, 1950),coords_z2_lim=1998)
## 1st STVCs which are significant at the 5 percent level, for byear <= 1950
#plot_s(res8,1, coords_z1_lim=c(-Inf, 1950), pmax=0.05)
## 2nd STVC for byear >= 1951
#plot_s(res8,2, coords_z1_lim=c(1951,Inf))
Run the code above in your browser using DataLab