#####################################################
############ Spatial regression modeling ############
#####################################################
require(spdep);require(Matrix)
data(boston)
y <- boston.c[, "CMEDV" ]
x <- boston.c[,c("CRIM","ZN","INDUS", "CHAS", "NOX","RM", "AGE",
"DIS" ,"RAD", "TAX", "PTRATIO", "B", "LSTAT")]
xgroup<- boston.c[,"TOWN"]
coords<- boston.c[,c("LON","LAT")]
meig <- meigen(coords=coords)
# meig<- meigen_f(coords=coords) ## for large samples
#####################################################
####### Gaussian regression #########################
res <- resf(y = y, x = x, meig = meig)
res
plot_s(res) ## spatially dependent component (intercept)
#### Group-wise random intercepts
#res2 <- resf(y = y, x = x, meig = meig, xgroup = xgroup)
#### Group-level spatial dependence (s_id) + random intercepts (xgroup)
#meig_g<- meigen(coords=coords, s_id = xgroup)
#res3 <- resf(y = y, x = x, meig = meig_g, xgroup = xgroup)
#### Coefficients varying depending on x
#res4 <- resf(y = y, x = x, meig = meig, nvc = TRUE)
#res4
#plot_s(res4) # spatially dependent component (intercept)
#plot_s(res4,5) # spatial plot of the 5-th NVC
#plot_s(res4,6) # spatial plot of the 6-th NVC
#plot_s(res4,13)# spatial plot of the 13-th NVC
#plot_n(res4,5) # 1D plot of the 5-th NVC
#plot_n(res4,6) # 1D plot of the 6-th NVC
#plot_n(res4,13)# 1D plot of the 13-th NVC
#####################################################
###### Non-Gaussian regression ######################
#### Model for non-Gaussian continuous data
# - Probability distribution is estimated from data
#ng5 <- nongauss_y( tr_num = 2 )# 2 SAL transformations to Gaussianize y
#res5 <- resf(y = y, x = x, meig = meig, nongauss = ng5)
#res5 ## tr_num may be selected by comparing BIC
#plot(res5$pdf,type="l") # Estimated probability density function
#res5$skew_kurt # Skew and kurtosis of the estimated PDF
#res5$pred_quantile[1:2,]# predicted value by quantile
#coef_marginal(res5) # Estimated marginal effects (dy/dx)
#### Model for non-Gaussian and non-negative continuous data
# - Probability distribution is estimated from data
#ng6 <- nongauss_y( tr_num = 2, y_nonneg = TRUE )
#res6 <- resf(y = y, x = x, meig = meig, nongauss = ng6 )
#coef_marginal(res6)
#### Overdispersed Poisson model for count data
# - y: count data
#ng7 <- nongauss_y( y_type = "count" )
#res7 <- resf(y = y, x = x, meig = meig, nongauss = ng7 )
#### General model for count data
# - y: count data
# - Probability distribution is estimated from data
#ng8 <- nongauss_y( y_type = "count", tr_num = 2 )
#res8 <- resf(y = y, x = x, meig = meig, nongauss = ng8 )
#####################################################
################ 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", "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)
#res9 <- resf(y=y,x=x,meig=meig )
#res9
Run the code above in your browser using DataLab