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 SVC models ##################
#####################################################
#### SVC or constant coefficients on x ##############
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 (n > 5,000), the following additional learning
### mitigates an degeneracy/over-smoothing problem in SVCs
# res_adj<- addlearn_local(res)
# res_adj
# plot_s(res_adj,0)
# plot_s(res_adj,1)
# plot_s(res_adj,2)
#### SVC on x #######################################
#res2 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_sel = FALSE )
#### Group-level SVC or constant coefficients on x ##
#### Group-wise random intercepts ###################
#meig_g <- meigen(coords, s_id=xgroup)
#res3 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig_g,xgroup=xgroup)
#####################################################
############## Gaussian SNVC models #################
#####################################################
#### SNVC, SVC, NVC, or constant coefficients on x ###
#res4 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_nvc =TRUE)
#### SNVC, SVC, NVC, or constant coefficients on x ###
#### NVC or Constant coefficients on xconst ##########
#res5 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_nvc =TRUE, xconst_nvc=TRUE)
#plot_s(res5,0) # Spatially varying intercept
#plot_s(res5,1) # Spatial plot of the SNVC (SVC + NVC) on x[,1]
#plot_s(res5,1,btype="svc")# Spatial plot of SVC in the SNVC
#plot_s(res5,1,btype="nvc")# Spatial plot of NVC in the SNVC
#plot_n(res5,1) # 1D plot of the NVC
#plot_s(res5,6,xtype="xconst")# Spatial plot of the NVC on xconst[,6]
#plot_n(res5,6,xtype="xconst")# 1D plot of the NVC on xconst[,6]
#####################################################
############## Non-Gaussian SVC models ##############
#####################################################
#### Generalized model for continuous data ##########
# - Probability distribution is estimated from data
#ng6 <- nongauss_y( tr_num = 2 )# 2 SAL transformations to Gaussianize y
#res6 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng6 )
#res6 # tr_num may be selected by comparing BIC (or AIC)
#coef_marginal_vc(res6) # marginal effects from x (dy/dx)
#plot(res6$pdf,type="l") # Estimated probability density function
#res6$skew_kurt # Skew and kurtosis of the estimated PDF
#res6$pred_quantile[1:2,]# predicted value by quantile
#### Generalized model for non-negative continuous data
# - Probability distribution is estimated from data
#ng7 <- nongauss_y( tr_num = 2, y_nonneg = TRUE )
#res7 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng7 )
#coef_marginal_vc(res7)
#### Overdispersed Poisson model for count data #####
# - y is assumed as a count data
#ng8 <- nongauss_y( y_type = "count" )
#res8 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng8 )
#### Generalized model for count data ###############
# - y is assumed as a count data
# - Probability distribution is estimated from data
#ng9 <- nongauss_y( y_type = "count", tr_num = 2 )
#res9 <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, nongauss = ng9 )
Run the code above in your browser using DataLab