spmoran
### https://cran.r-project.org/web/packages/spmoran/vignettes/boston_sample.pdf
library(spdep)
library(spmoran)
data(boston)
y <- boston.c[, "CMEDV" ]
x <- boston.c[,c("CRIM","ZN","INDUS", "CHAS", "NOX","RM", "AGE")]
coords<- boston.c[,c("LON","LAT")]
#########Distance-based ESF
meig <- meigen(coords=coords)
res <- esf(y=y,x=x,meig=meig, vif=10)
res
## Call:
## esf(y = y, x = x, vif = 10, meig = meig)
##
## ----Coefficients------------------------------
## Estimate SE t_value p_value
## (Intercept) 11.34040959 3.91692274 2.8952344 3.968277e-03
## CRIM -0.20942091 0.03048530 -6.8695702 2.089395e-11
## ZN 0.02322000 0.01384823 1.6767492 9.426799e-02
## INDUS -0.15063613 0.06823776 -2.2075188 2.776856e-02
## CHAS 0.15172838 0.93842988 0.1616832 8.716260e-01
## NOX -38.02167637 4.79403898 -7.9310320 1.651338e-14
## RM 6.33316024 0.36887955 17.1686403 1.842211e-51
## AGE -0.07820247 0.01564970 -4.9970593 8.274067e-07
##
## ----Spatial effects (residuals)---------------
## Estimate
## SE 6.8540461
## Moran.I/max(Moran.I) 0.6701035
##
## ----Error statistics--------------------------
## stat
## resid_SE 4.476459
## adjR2 0.762328
## logLik -1453.376154
## AIC 2996.752308
## BIC 3186.946458
meig_f<- meigen_f(coords)
res <- esf(y=y, x=x, meig=meig_f,vif=10, fn="all")
res <- resf(y = y, x = x, meig = meig)
res
## Call:
## resf(y = y, x = x, meig = meig)
##
## ----Coefficients------------------------------
## Estimate SE t_value p_value
## (Intercept) 6.63220230 3.94484187 1.6812340 9.340113e-02
## CRIM -0.19815202 0.03126666 -6.3374866 5.608680e-10
## ZN 0.01453736 0.01591772 0.9132813 3.615764e-01
## INDUS -0.15560251 0.06842940 -2.2739131 2.343446e-02
## CHAS 0.51046260 0.92329946 0.5528679 5.806244e-01
## NOX -31.26689872 5.02069114 -6.2276085 1.075127e-09
## RM 6.33993149 0.36671337 17.2885203 0.000000e+00
## AGE -0.06351411 0.01526957 -4.1595217 3.810684e-05
##
## ----Variance parameter------------------------
##
## Spatial effects (residuals):
## (Intercept)
## random_SE 6.7424423
## Moran.I/max(Moran.I) 0.6648678
##
## ----Error statistics--------------------------
## stat
## resid_SE 4.3515212
## adjR2(cond) 0.7735912
## rlogLik -1540.3812428
## AIC 3102.7624855
## BIC 3149.2543889
meig_f<- meigen_f(coords)
res <- resf(y = y, x = x, meig = meig_f)
res <- resf(y = y, x = x, meig = meig, nvc=TRUE)
res
## Call:
## resf(y = y, x = x, nvc = TRUE, meig = meig)
##
## ----Non-spatially varying coefficients (summary)----
##
## Coefficients:
## Intercept CRIM ZN INDUS
## Min. :25.41 Min. :-0.1822 Min. :0.02042 Min. :-0.2119
## 1st Qu.:25.41 1st Qu.:-0.1822 1st Qu.:0.02042 1st Qu.:-0.2119
## Median :25.41 Median :-0.1822 Median :0.02042 Median :-0.2119
## Mean :25.41 Mean :-0.1822 Mean :0.02042 Mean :-0.2119
## 3rd Qu.:25.41 3rd Qu.:-0.1822 3rd Qu.:0.02042 3rd Qu.:-0.2119
## Max. :25.41 Max. :-0.1822 Max. :0.02042 Max. :-0.2119
## CHAS NOX RM AGE
## Min. :1.375 Min. :-0.463 Min. :-0.78043 Min. :-0.06742
## 1st Qu.:1.375 1st Qu.: 6.083 1st Qu.:-0.40834 1st Qu.:-0.06742
## Median :1.375 Median : 7.792 Median :-0.16098 Median :-0.06742
## Mean :1.375 Mean : 7.074 Mean : 0.03975 Mean :-0.06742
## 3rd Qu.:1.375 3rd Qu.: 8.654 3rd Qu.: 0.19417 3rd Qu.:-0.06742
## Max. :1.375 Max. :11.517 Max. : 2.49406 Max. :-0.06742
##
## Statistical significance:
## Intercept CRIM ZN INDUS CHAS NOX RM AGE
## Not significant 0 0 506 0 0 506 472 0
## Significant (10% level) 0 0 0 0 506 0 7 0
## Significant ( 5% level) 0 0 0 0 0 0 10 0
## Significant ( 1% level) 506 506 0 506 0 0 17 506
##
## ----Variance parameter------------------------
##
## Spatial effects (residuals):
## (Intercept)
## random_SE 3.6981527
## Moran.I/max(Moran.I) 0.4490228
##
## Non-spatially varying coefficients:
## CRIM ZN INDUS CHAS NOX RM AGE
## random_SE 0 0 0 0 1.850518 0.2459548 0
##
## ----Error statistics--------------------------
## stat
## resid_SE 3.7949128
## adjR2(cond) 0.8271073
## rlogLik -1478.6128728
## AIC 2983.2257457
## BIC 3038.1707224
res <- resf(y = y, x = x, meig = meig, nvc=TRUE, nvc_sel=1:3)
res
## Call:
## resf(y = y, x = x, nvc = TRUE, nvc_sel = 1:3, meig = meig)
##
## ----Non-spatially varying coefficients (summary)----
##
## Coefficients:
## Intercept CRIM ZN INDUS
## Min. :8.04 Min. :-0.1978 Min. :-0.02646 Min. :-0.1618
## 1st Qu.:8.04 1st Qu.:-0.1978 1st Qu.: 0.02423 1st Qu.:-0.1618
## Median :8.04 Median :-0.1978 Median : 0.02423 Median :-0.1618
## Mean :8.04 Mean :-0.1978 Mean : 0.02047 Mean :-0.1618
## 3rd Qu.:8.04 3rd Qu.:-0.1978 3rd Qu.: 0.02423 3rd Qu.:-0.1618
## Max. :8.04 Max. :-0.1978 Max. : 0.07651 Max. :-0.1618
## CHAS NOX RM AGE
## Min. :0.5596 Min. :-32.04 Min. :6.218 Min. :-0.06464
## 1st Qu.:0.5596 1st Qu.:-32.04 1st Qu.:6.218 1st Qu.:-0.06464
## Median :0.5596 Median :-32.04 Median :6.218 Median :-0.06464
## Mean :0.5596 Mean :-32.04 Mean :6.218 Mean :-0.06464
## 3rd Qu.:0.5596 3rd Qu.:-32.04 3rd Qu.:6.218 3rd Qu.:-0.06464
## Max. :0.5596 Max. :-32.04 Max. :6.218 Max. :-0.06464
##
## Statistical significance:
## Intercept CRIM ZN INDUS CHAS NOX RM AGE
## Not significant 0 0 496 0 506 0 0 0
## Significant (10% level) 0 0 0 0 0 0 0 0
## Significant ( 5% level) 506 0 5 506 0 0 0 0
## Significant ( 1% level) 0 506 5 0 0 506 506 506
##
## ----Variance parameter------------------------
##
## Spatial effects (residuals):
## (Intercept)
## random_SE 6.6961726
## Moran.I/max(Moran.I) 0.6708207
##
## Non-spatially varying coefficients:
## CRIM ZN INDUS CHAS NOX RM AGE
## random_SE 2.923838e-08 0.008130434 2.056973e-07 0 0 0 0
##
## ----Error statistics--------------------------
## stat
## resid_SE 4.2790184
## adjR2(cond) 0.7797353
## rlogLik -1537.6449527
## AIC 3103.2899053
## BIC 3162.4614187
y <- boston.c[, "CMEDV"]
x <- boston.c[,c("CRIM", "AGE")]
xconst <- boston.c[,c("ZN","DIS","RAD","NOX", "TAX","RM", "PTRATIO", "B")]
coords <- boston.c[,c("LON","LAT")]
meig <- meigen(coords=coords)
res <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_sel = FALSE )
## [1] "------- Iteration 1 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3120.605"
## [1] "------- Iteration 2 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3114.252"
## [1] "------- Iteration 3 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3114.139"
## [1] "------- Iteration 4 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3114.138"
## Call:
## resf_vc(y = y, x = x, xconst = xconst, x_sel = FALSE, meig = meig)
##
## ----Spatially varying coefficients on x (summary)----
##
## Coefficient estimates:
## (Intercept) CRIM AGE
## Min. :12.03 Min. :-3.29294 Min. :-0.14986
## 1st Qu.:13.99 1st Qu.:-0.19941 1st Qu.:-0.08377
## Median :15.06 Median : 0.04993 Median :-0.06780
## Mean :15.70 Mean : 0.05902 Mean :-0.06582
## 3rd Qu.:17.31 3rd Qu.: 0.36587 3rd Qu.:-0.04710
## Max. :20.46 Max. : 1.83866 Max. : 0.04298
##
## Statistical significance:
## Intercept CRIM AGE
## Not significant 0 416 147
## Significant (10% level) 0 27 40
## Significant ( 5% level) 190 17 99
## Significant ( 1% level) 316 46 220
##
## ----Constant coefficients on xconst----------------------------
## Estimate SE t_value p_value
## ZN 0.03202068 0.013219003 2.422322 1.582817e-02
## DIS -1.47514930 0.334360238 -4.411856 1.292875e-05
## RAD 0.36064288 0.090818317 3.971037 8.368693e-05
## NOX -36.21088315 5.134427150 -7.052565 6.925571e-12
## TAX -0.01242296 0.003502523 -3.546862 4.320840e-04
## RM 6.49212566 0.326197980 19.902409 0.000000e+00
## PTRATIO -0.52573980 0.151594626 -3.468064 5.762765e-04
## B 0.02091202 0.003094117 6.758638 4.477529e-11
##
## ----Variance parameters----------------------------------
##
## Spatial variation (coefficients on x):
## (Intercept) CRIM AGE
## random_SE 3.9039831 1.59443322 0.05746111
## Moran.I/max(Moran.I) 0.6627375 0.04502003 0.06267778
##
## ----Error statistics-------------------------------------
## stat
## resid_SE 3.6706778
## adjR2(cond) 0.8375658
## rlogLik -1501.0302460
## AIC 3038.0604921
## BIC 3114.1381521
plot_s(res,0) # Spatially varying intercept
res <- resf_vc(y=y,x=x,xconst=xconst,meig=meig )
## [1] "------- Iteration 1 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3120.605"
## [1] "------- Iteration 2 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3107.452"
## [1] "------- Iteration 3 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3106.939"
## [1] "------- Iteration 4 -------"
## [1] "1/3"
## [1] "2/3"
## [1] "3/3"
## [1] "BIC: 3106.939"
## Call:
## resf_vc(y = y, x = x, xconst = xconst, meig = meig)
##
## ----Spatially varying coefficients on x (summary)----
##
## Coefficient estimates:
## (Intercept) CRIM AGE
## Min. :11.17 Min. :-0.1814 Min. :-0.14114
## 1st Qu.:12.96 1st Qu.:-0.1814 1st Qu.:-0.07938
## Median :14.24 Median :-0.1814 Median :-0.06451
## Mean :14.75 Mean :-0.1814 Mean :-0.06198
## 3rd Qu.:16.64 3rd Qu.:-0.1814 3rd Qu.:-0.04730
## Max. :19.40 Max. :-0.1814 Max. : 0.05117
##
## Statistical significance:
## Intercept CRIM AGE
## Not significant 0 0 123
## Significant (10% level) 0 0 48
## Significant ( 5% level) 255 0 100
## Significant ( 1% level) 251 506 235
##
## ----Constant coefficients on xconst----------------------------
## Estimate SE t_value p_value
## ZN 0.03473113 0.013895397 2.499470 1.278897e-02
## DIS -1.34121749 0.325351805 -4.122361 4.457033e-05
## RAD 0.29200513 0.082384859 3.544403 4.341876e-04
## NOX -29.36631205 4.942673715 -5.941382 5.622099e-09
## TAX -0.01371011 0.003512961 -3.902723 1.094893e-04
## RM 6.26622242 0.340562626 18.399619 0.000000e+00
## PTRATIO -0.53923932 0.151877936 -3.550478 4.245538e-04
## B 0.01971973 0.003090429 6.380904 4.338061e-10
##
## ----Variance parameters----------------------------------
##
## Spatial variation (coefficients on x):
## (Intercept) CRIM AGE
## random_SE 3.715171 0 0.05169206
## Moran.I/max(Moran.I) 0.789340 NA 0.04975523
##
## ----Error statistics-------------------------------------
## stat
## resid_SE 4.0311657
## adjR2(cond) 0.8048959
## rlogLik -1503.6570917
## AIC 3039.3141834
## BIC 3106.9387701
res <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_nvc =TRUE)
## [1] "------- Iteration 1 -------"
## [1] "1/5"
## [1] "2/5"
## [1] "3/5"
## [1] "4/5"
## [1] "5/5"
## [1] "BIC: 3118.893"
## [1] "------- Iteration 2 -------"
## [1] "1/5"
## [1] "2/5"
## [1] "3/5"
## [1] "4/5"
## [1] "5/5"
## [1] "BIC: 3110.52"
## [1] "------- Iteration 3 -------"
## [1] "1/5"
## [1] "2/5"
## [1] "3/5"
## [1] "4/5"
## [1] "5/5"
## [1] "BIC: 3110.519"
## [1] "------- Iteration 4 -------"
## [1] "1/5"
## [1] "2/5"
## [1] "3/5"
## [1] "4/5"
## [1] "5/5"
## [1] "BIC: 3110.519"
## Call:
## resf_vc(y = y, x = x, xconst = xconst, x_nvc = TRUE, meig = meig)
##
## ----Spatially and non-spatially varying coefficients on x (summary)----
##
## Coefficient estimates:
## (Intercept) CRIM AGE
## Min. :13.7 Min. :-0.1837 Min. :-0.16218
## 1st Qu.:13.7 1st Qu.:-0.1837 1st Qu.:-0.07425
## Median :13.7 Median :-0.1837 Median :-0.05491
## Mean :13.7 Mean :-0.1837 Mean :-0.04870
## 3rd Qu.:13.7 3rd Qu.:-0.1837 3rd Qu.:-0.02589
## Max. :13.7 Max. :-0.1837 Max. : 0.08386
##
## Statistical significance:
## Intercept CRIM AGE
## Not significant 0 0 169
## Significant (10% level) 0 0 45
## Significant ( 5% level) 506 0 85
## Significant ( 1% level) 0 506 207
##
## ----Constant coefficients on xconst----------------------------
## Estimate SE t_value p_value
## ZN 0.03621115 0.013711134 2.641003 8.549301e-03
## DIS -1.65624971 0.259776835 -6.375664 4.462573e-10
## RAD 0.30482416 0.081633873 3.734040 2.122319e-04
## NOX -27.93545504 4.891161245 -5.711416 2.021062e-08
## TAX -0.01337477 0.003493264 -3.828731 1.467697e-04
## RM 6.37243847 0.343764351 18.537229 0.000000e+00
## PTRATIO -0.56324946 0.150692561 -3.737739 2.092264e-04
## B 0.01926817 0.003112574 6.190429 1.336719e-09
##
## ----Variance parameters----------------------------------
##
## Spatial variation (coefficients on x):
## (Intercept) CRIM AGE
## random_SE 0.002753139 0 0.06316541
## Moran.I/max(Moran.I) 0.342830186 NA 0.23319003
##
## Non-spatial variation (coefficients on x):
## CRIM AGE
## random_SE 0 0
##
## ----Error statistics-------------------------------------
## stat
## resid_SE 4.0639128
## adjR2(cond) 0.8017132
## rlogLik -1505.4474458
## AIC 3042.8948916
## BIC 3110.5194783
res <- resf_vc(y=y,x=x,xconst=xconst,meig=meig, x_nvc =TRUE,
xconst_nvc=TRUE)
## [1] "------- Iteration 1 -------"
## [1] "1/13"
## [1] "2/13"
## [1] "3/13"
## [1] "4/13"
## [1] "5/13"
## [1] "7/13"
## [1] "8/13"
## [1] "9/13"
## [1] "10/13"
## [1] "11/13"
## [1] "12/13"
## [1] "13/13"
## [1] "BIC: 3023.44"
## [1] "------- Iteration 2 -------"
## [1] "1/13"
## [1] "2/13"
## [1] "3/13"
## [1] "4/13"
## [1] "5/13"
## [1] "7/13"
## [1] "8/13"
## [1] "9/13"
## [1] "10/13"
## [1] "11/13"
## [1] "12/13"
## [1] "13/13"
## [1] "BIC: 3013.009"
## [1] "------- Iteration 3 -------"
## [1] "1/13"
## [1] "2/13"
## [1] "3/13"
## [1] "4/13"
## [1] "5/13"
## [1] "7/13"
## [1] "8/13"
## [1] "9/13"
## [1] "10/13"
## [1] "11/13"
## [1] "12/13"
## [1] "13/13"
## [1] "BIC: 3012.86"
## [1] "------- Iteration 4 -------"
## [1] "1/13"
## [1] "2/13"
## [1] "3/13"
## [1] "4/13"
## [1] "5/13"
## [1] "7/13"
## [1] "8/13"
## [1] "9/13"
## [1] "10/13"
## [1] "11/13"
## [1] "12/13"
## [1] "13/13"
## [1] "BIC: 3012.858"
## [1] "------- Iteration 5 -------"
## [1] "1/13"
## [1] "2/13"
## [1] "3/13"
## [1] "4/13"
## [1] "5/13"
## [1] "7/13"
## [1] "8/13"
## [1] "9/13"
## [1] "10/13"
## [1] "11/13"
## [1] "12/13"
## [1] "13/13"
## [1] "BIC: 3012.857"
## Call:
## resf_vc(y = y, x = x, xconst = xconst, x_nvc = TRUE, xconst_nvc = TRUE,
## meig = meig)
##
## ----Spatially and non-spatially varying coefficients on x (summary)----
##
## Coefficient estimates:
## (Intercept) CRIM AGE
## Min. :34.99 Min. :-2.1670 Min. :-0.07495
## 1st Qu.:40.95 1st Qu.:-0.6135 1st Qu.:-0.07495
## Median :42.29 Median :-0.4158 Median :-0.07495
## Mean :42.44 Mean :-0.4289 Mean :-0.07495
## 3rd Qu.:43.78 3rd Qu.:-0.2156 3rd Qu.:-0.07495
## Max. :49.95 Max. : 0.5207 Max. :-0.07495
##
## Statistical significance:
## Intercept CRIM AGE
## Not significant 0 394 0
## Significant (10% level) 0 15 0
## Significant ( 5% level) 0 29 0
## Significant ( 1% level) 506 68 506
##
## ----Non-spatially varying coefficients on xconst (summary)----
##
## Coefficient estimates:
## ZN DIS RAD NOX
## Min. :0.02512 Min. :-1.107 Min. :0.6287 Min. :-23.30
## 1st Qu.:0.02512 1st Qu.:-1.107 1st Qu.:0.6287 1st Qu.:-19.37
## Median :0.02512 Median :-1.107 Median :0.6287 Median :-18.48
## Mean :0.02512 Mean :-1.107 Mean :0.6287 Mean :-18.55
## 3rd Qu.:0.02512 3rd Qu.:-1.107 3rd Qu.:0.6287 3rd Qu.:-17.57
## Max. :0.02512 Max. :-1.107 Max. :0.6287 Max. :-14.47
## TAX RM PTRATIO B
## Min. :-0.01512 Min. :0.5988 Min. :-0.6371 Min. :0.01371
## 1st Qu.:-0.01512 1st Qu.:0.8372 1st Qu.:-0.6371 1st Qu.:0.01371
## Median :-0.01512 Median :1.0394 Median :-0.6371 Median :0.01371
## Mean :-0.01512 Mean :1.2054 Mean :-0.6371 Mean :0.01371
## 3rd Qu.:-0.01512 3rd Qu.:1.3012 3rd Qu.:-0.6371 3rd Qu.:0.01371
## Max. :-0.01512 Max. :3.2979 Max. :-0.6371 Max. :0.01371
##
## Statistical significance:
## ZN DIS RAD NOX TAX RM PTRATIO B
## Not significant 0 0 0 185 0 414 0 0
## Significant (10% level) 506 0 0 217 0 27 0 0
## Significant ( 5% level) 0 0 0 40 0 23 0 0
## Significant ( 1% level) 0 506 506 64 506 42 506 506
##
## ----Variance parameters----------------------------------
##
## Spatial variation (coefficients on x):
## (Intercept) CRIM AGE
## random_SE 4.0639969 0.99802716 0
## Moran.I/max(Moran.I) 0.3274852 0.07446611 NA
##
## Non-spatial variation (coefficients on x):
## CRIM AGE
## random_SE 0.03403638 0
##
## Non-spatial variation (coefficients on xconst):
## ZN DIS RAD NOX TAX RM PTRATIO B
## random_SE 0 0 0 1.496749 0 0.2001897 0 0
##
## ----Error statistics-------------------------------------
## stat
## resid_SE 3.1950502
## adjR2(cond) 0.8766801
## rlogLik -1447.2765888
## AIC 2932.5531776
## BIC 3012.8573743
plot_n(res,4,xtype="xconst")#NVC on xconst[,4]
plot_n(res,6,xtype="xconst")#NVC on xconst[,6]
plot_s(res,4,xtype="xconst")#NVC on xconst[,4]
plot_s(res,6,xtype="xconst")#NVC on xconst[,6]
xgroup<- boston.c[,"TOWNNO"]
coords<- boston.c[,c("LON","LAT")]
meig_g<- meigen(coords=coords, s_id=xgroup)
x <- boston.c[,c("CRIM","ZN","INDUS", "CHAS", "NOX","RM", "AGE")]
res <- resf(y = y, x = x, meig = meig_g, xgroup = xgroup)
res
## Call:
## resf(y = y, x = x, xgroup = xgroup, meig = meig_g)
##
## ----Coefficients------------------------------
## Estimate SE t_value p_value
## (Intercept) -0.81545943 3.23135854 -0.2523581 8.008871e-01
## CRIM -0.04596392 0.02505503 -1.8345188 6.728064e-02
## ZN 0.03285021 0.02313784 1.4197611 1.564153e-01
## INDUS 0.03549188 0.11980486 0.2962474 7.671869e-01
## CHAS -0.62561231 0.72381491 -0.8643264 3.878995e-01
## NOX -26.38632672 3.88238119 -6.7964286 3.668488e-11
## RM 6.30273567 0.29409796 21.4307357 0.000000e+00
## AGE -0.06730232 0.01048068 -6.4215611 3.637544e-10
##
## ----Variance parameter------------------------
##
## Spatial effects (residuals):
## (Intercept)
## random_SE 5.074794
## Moran.I/max(Moran.I) 0.812936
##
## Group effects:
## xgroup
## ramdom_SE 4.4404
##
## ----Error statistics--------------------------
## stat
## resid_SE 3.2429178
## adjR2(cond) 0.8740022
## rlogLik -1465.8450362
## AIC 2955.6900724
## BIC 3006.4085124
res$b_g[[1]][1:5,]# Estimates in the first 5 groups
## Estimate SE t_value
## xgroup_0 2.165726 2.061093 1.0507657
## xgroup_1 3.747633 1.783543 2.1012294
## xgroup_2 6.544205 1.659184 3.9442318
## xgroup_3 2.431558 1.431325 1.6988163
## xgroup_4 1.036033 1.181672 0.8767521
r_res <-resf(y=y, x=x, meig=meig_g, xgroup=xgroup)
pred <-predict0(r_res, x0=x, meig0=meig_g, xgroup0=xgroup)
pred$pred[1:5,]
## pred xb sf_residual xgroup
## 1 23.70932 22.71407 -1.170482 2.165726
## 2 24.57615 22.21874 -1.390220 3.747633
## 3 30.58942 28.23201 -1.390220 3.747633
## 4 33.24998 28.19959 -1.493814 6.544205
## 5 33.62206 28.57167 -1.493814 6.544205
adat <- aggregate(data.frame(y, pred$pred),by=list(xgroup),mean)
adat[1:5,]
## Group.1 y pred xb sf_residual xgroup
## 1 0 24.00000 23.70932 22.71407 -1.170482 2.165726
## 2 1 28.15000 27.58279 25.22537 -1.390220 3.747633
## 3 2 32.76667 31.89132 26.84093 -1.493814 6.544205
## 4 3 19.42857 19.36679 18.51187 -1.576641 2.431558
## 5 4 16.71364 16.72781 17.10793 -1.416151 1.036033
require(rgdal)
require(rgeos)
require(dplyr)
boston.tr <- readOGR(system.file("shapes/boston_tracts.shp",package="spData")[1])
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\subas\Documents\R\win-library\4.0\spData\shapes\boston_tracts.shp", layer: "boston_tracts"
## with 506 features
## It has 36 fields
b1 <- st_as_sf(boston.tr)
b1_dissolve <- b1 %>% group_by(TOWNNO) %>% summarize() #dissolve
boston.tr2 <- as_Spatial(b1_dissolve)
boston.tr2@data$id<- 1:(dim(boston.tr2)[1])
b2_dat <- boston.tr2@data
b2_dat2 <- merge(b2_dat, adat,by.x="TOWNNO",by.y="Group.1",all.x=TRUE)
boston.tr2@data<- b2_dat2[order(b2_dat2$id),]
spplot(boston.tr2,c("y","pred"), lwd=0.3)
spplot(boston.tr2,c("xgroup","sf_residual"), lwd=0.3)
spplot(boston.tr2,"xb", lwd=0.3)
rv_res <- resf_vc(y=y, x=x, meig=meig_g, xgroup=xgroup,
x_sel=FALSE)
## [1] "------- Iteration 1 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3101.038"
## [1] "------- Iteration 2 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3098.223"
## [1] "------- Iteration 3 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3040.184"
## [1] "------- Iteration 4 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3039.573"
## [1] "------- Iteration 5 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3039.571"
## [1] "------- Iteration 6 -------"
## [1] "1/9"
## [1] "2/9"
## [1] "3/9"
## [1] "4/9"
## [1] "5/9"
## [1] "6/9"
## [1] "7/9"
## [1] "8/9"
## [1] "9/9"
## [1] "BIC: 3039.571"
pred_vc <- predict0_vc(rv_res, x0=x, meig0=meig_g, xgroup0=xgroup)
adat_vc <- aggregate(data.frame(y, pred_vc$pred), by=list(xgroup), mean)
adat_vc[1:5,]
## Group.1 y pred xb sf_residual xgroup
## 1 0 24.00000 23.67942 23.12339 -1.123210 1.679234
## 2 1 28.15000 27.81223 27.43916 -1.962508 2.335580
## 3 2 32.76667 32.28764 31.09015 -2.546255 3.743741
## 4 3 19.42857 19.25746 18.45659 -2.500409 3.301278
## 5 4 16.71364 16.68370 15.40354 -1.024050 2.304206
require(plm)
require(spData)
data(Produc, package = "plm")
data(us_states)
us_states2 <- data.frame(us_states$GEOID,us_states$NAME,st_coordinates(st_centroid(us_states)))
names(us_states2)[3:4]<- c("px","py")
us_states3 <- us_states2[order(us_states2[,1]),][-8,]
us_states3$state<- unique(Produc[,1])
pdat0 <- na.omit(merge(Produc,us_states3[,c(3:5)],by="state",all.x=TRUE,sort=FALSE))
pdat <- pdat0[order(pdat0$state,pdat0$year),]
pdat[1:5,]
## state year region pcap hwy water util pc gsp emp
## 1 ALABAMA 1970 6 15032.67 7325.80 1655.68 6051.20 35793.80 28418 1010.5
## 2 ALABAMA 1971 6 15501.94 7525.94 1721.02 6254.98 37299.91 29375 1021.9
## 3 ALABAMA 1972 6 15972.41 7765.42 1764.75 6442.23 38670.30 31303 1072.3
## 4 ALABAMA 1973 6 16406.26 7907.66 1742.41 6756.19 40084.01 33430 1135.5
## 5 ALABAMA 1974 6 16762.67 8025.52 1734.85 7002.29 42057.31 33749 1169.8
## unemp px py
## 1 4.7 -86.82645 32.7926
## 2 5.2 -86.82645 32.7926
## 3 4.7 -86.82645 32.7926
## 4 3.9 -86.82645 32.7926
## 5 5.5 -86.82645 32.7926
y <- log(pdat$gsp)
x <- data.frame(log_pcap=log(pdat$pcap), log_pc=log(pdat$pc),
log_emp=log(pdat$emp), unemp=pdat$unemp)
coords<- pdat[,c("px", "py")]
s_id <- pdat$state
meig_p<- meigen(coords,s_id=s_id)# Moran eigenvectors by states
pmod0 <- resf(y=y,x=x,meig=meig_p) # pooling model
xgroup<- pdat[,c("state")]
pmod1 <- resf(y=y,x=x,meig=meig_p,xgroup=xgroup)# individual model
xgroup<- pdat[,c("year")]
pmod2 <- resf(y=y,x=x,meig=meig_p,xgroup=xgroup)# time model
xgroup<- pdat[,c("state","year")]
pmod3 <- resf(y=y,x=x,meig=meig_p,xgroup=xgroup)# two-way model
s_g <- pmod3$b_g[[1]]
s_g[1:5,]# State-level group effects
## Estimate SE t_value
## state_ALABAMA -0.07162777 0.01390150 -5.152520
## state_ARIZONA -0.04406765 0.01668096 -2.641794
## state_ARKANSAS -0.07255311 0.01471153 -4.931718
## state_CALIFORNIA 0.24008052 0.01967542 12.202053
## state_COLORADO -0.11495815 0.01232160 -9.329805
t_g <- pmod3$b_g[[2]]
t_g[1:5,]# Year-level group effects
## Estimate SE t_value
## year_1970 -0.006015617 0.011091180 -0.5423784
## year_1971 0.002902572 0.010569185 0.2746259
## year_1972 0.013282446 0.010416809 1.2750974
## year_1973 0.021949814 0.010280020 2.1351917
## year_1974 -0.009852342 0.009679287 -1.0178789
pm0 <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
data = pdat, effect="twoways",model="random")
pm0
##
## Model Formula: log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
##
## Coefficients:
## (Intercept) log(pcap) log(pc) log(emp) unemp
## 2.3634993 0.0178529 0.2655895 0.7448989 -0.0045755
s_g_plm<- ranef(pm0,"individual")
t_g_plm<- ranef(pm0,"time")
plot(s_g_plm,s_g[,1],xlab="plm",ylab="resf")
abline(0,1,col="red")
plot(t_g_plm,t_g[,1],xlab="plm",ylab="resf")
abline(0,1,col="red")
y <- boston.c[, "CMEDV" ]
x <- boston.c[,c("CRIM","ZN","INDUS", "CHAS", "NOX","RM", "AGE")]
coords<- boston.c[,c("LON","LAT")]
meig <- meigen(coords=coords)
res <- resf_qr(y=y,x=x,meig=meig, boot=TRUE)
## [1] "------- Complete: tau=0.1 -------"
## [1] "------- Complete: tau=0.2 -------"
## [1] "------- Complete: tau=0.3 -------"
## [1] "------- Complete: tau=0.4 -------"
## [1] "------- Complete: tau=0.5 -------"
## [1] "------- Complete: tau=0.6 -------"
## [1] "------- Complete: tau=0.7 -------"
## [1] "------- Complete: tau=0.8 -------"
## [1] "------- Complete: tau=0.9 -------"
## Call:
## resf_qr(y = y, x = x, meig = meig, boot = TRUE)
##
## ----Coefficients------------------------------
## tau=0.1 tau=0.2 tau=0.3 tau=0.4 tau=0.5
## (Intercept) 23.86841970 29.16185736 26.550125353 21.16263694 17.151053980
## CRIM -0.36845124 -0.21172051 -0.106949379 -0.08357496 -0.070290258
## ZN -0.01169653 -0.01627637 -0.009652286 -0.01947512 -0.008198579
## INDUS 0.25009373 0.03992002 -0.111010420 -0.01521113 -0.096468769
## CHAS 0.98647836 1.28770409 0.438428954 0.26777796 -0.048278485
## NOX -32.89857783 -23.60303480 -15.109338348 -12.70090129 -11.263158727
## RM 0.71728433 0.49201634 1.169115918 2.21382993 3.004059676
## AGE 0.01977978 -0.05087471 -0.082548477 -0.11192561 -0.105681036
## tau=0.6 tau=0.7 tau=0.8 tau=0.9
## (Intercept) 13.999671526 11.28433168 -23.3939330 -57.24239068
## CRIM -0.064412593 -0.07823561 -0.1876252 -0.18934294
## ZN 0.007962903 0.01009742 0.1635369 0.03890142
## INDUS -0.167039581 -0.30344029 -0.9074079 -0.49797629
## CHAS -1.665298913 -1.51518801 -3.8773852 -0.04635798
## NOX -11.405913169 -20.36309658 -39.1980207 -41.26421537
## RM 3.730680883 5.25253569 13.7698457 19.62200618
## AGE -0.092068861 -0.07567382 -0.0587608 -0.03904752
##
## ----Spatial effects (residuals)---------------
## tau=0.1 tau=0.2 tau=0.3 tau=0.4 tau=0.5
## spcomp_SE 7.1522586 8.1254770 5.7952363 4.4135132 4.7198329
## spcomp_Moran.I/max(Moran.I) 0.2375865 0.3228553 0.3239407 0.3650454 0.5096847
## tau=0.6 tau=0.7 tau=0.8 tau=0.9
## spcomp_SE 4.8818059 6.3633073 16.9989855 16.3826940
## spcomp_Moran.I/max(Moran.I) 0.5690447 0.6935049 0.6757823 0.7203891
##
## ----Error statistics--------------------------
## tau=0.1 tau=0.2 tau=0.3 tau=0.4 tau=0.5 tau=0.6
## resid_SE 6.4395412 6.2086846 5.169030 4.7999618 4.5977255 4.8160068
## quasi_adjR2(cond) 0.6007294 0.6828421 0.666506 0.6183801 0.6229795 0.6121279
## tau=0.7 tau=0.8 tau=0.9
## resid_SE 5.6288391 12.2961444 18.6716254
## quasi_adjR2(cond) 0.6153019 0.6741455 0.4582676
