R: Estimate spatial or spatio-temporal semiparametric regression...

pspatfit {pspatreg}

R Documentation

Estimate spatial or spatio-temporal semiparametric regression models from a spatial econometric perspective.

Description

Estimate geoadditive spatial or spatio-temporal semiparametric regression models of type ps-sar, ps-sem, ps-sarar, ps-sdm, ps-sdem or ps-slx. These type of specifications are very general and they can include parametric and non-parametric covariates, spatial or spatio-temporal non-parametric trends and spatial lags of the dependent and independent variables and/or the noise of the model. The non-parametric terms (either trends or covariates) are modeled using P-Splines. The non-parametric trend can be decomposed in an ANOVA way including main and interactions effects of 2nd and 3rd order. The estimation method can be restricted maximum likelihood (REML) or maximum likelihood (ML).

Usage

pspatfit(
  formula,
  data,
  na.action,
  listw = NULL,
  type = "sim",
  method = "eigen",
  Durbin = NULL,
  zero.policy = NULL,
  interval = NULL,
  trs = NULL,
  cor = "none",
  dynamic = FALSE,
  demean = FALSE,
  eff_demean = "individual",
  index = NULL,
  control = list()
)

Arguments

`formula`	A formula similar to GAM specification including parametric and non-parametric terms. Parametric covariates are included in the usual way and non-parametric P-spline smooth terms are specified using `pspl(.)` and `pspt(.)` for the non-parametric covariates and spatial or spatio-temporal trend, respectively. More details in Details and Examples.
`data`	A data frame containing the parametric and non-parametric covariates included in the model. Also, if a `pspt(.)` term is included in formula, the data frame must include the spatial and temporal coordinates specified in `pspt(.)`. In this case coordinates must be ordered choosing time as fast index and spatial coordinates as low indexes. See `head(unemp_it)` as an example.
`na.action`	A function (default `options("na.action")`), can also be 'na.omit' or 'na.exclude' with consequences for residuals and fitted values. It may be necessary to set 'zero.policy' to 'TRUE' because this subsetting may create no-neighbour observations.
`listw`	Default = 'NULL' This will create a model with no spatial dependency. To include spatial dependency, `listw` should be a spatial neighbours list object created for example by `nb2listw` from spdep package; if `nb2listw` not given, set to the same spatial weights as the `listw` argument. It can also be a spatial weighting matrix of order (NxN) instead of a `listw` neighbours list object.
`type`	Type of spatial model specification following the usual spatial econometric terminology. Default = `"sim"` this creates a model with no type of spatial dependency. Types of spatial models available (similar to spsur package): `"sar"`, `"sem"`, `"sdm"`, `"sdem"`, `"sarar"`, or `"slx"`. When creating a `"slx"`, `"sdem"` or `"sdm"` model, it is necessary to include the formula of the Durbin part in the `Durbin` argument in the same way than spsur or spatialreg packages. There are examples on how to create these models in Examples section.
`method`	Similar to the corresponding parameter of `lagsarlm` function in spatialreg package. "eigen" (default) - the Jacobian is computed as the product of (1 - rho*eigenvalue) using `eigenw` from package spatialreg. For big samples (> 500) method = "eigen" is not recommended. Use "spam" or "Matrix_J" for strictly symmetric weights lists of styles "B" and "C", or made symmetric by similarity (Ord, 1975, Appendix C) if possible for styles "W" and "S", using code from the spam or Matrix packages to calculate the determinant; "Matrix" and "spam_update" provide updating Cholesky decomposition methods; "LU" provides an alternative sparse matrix decomposition approach. In addition, there are "Chebyshev" and Monte Carlo "MC" approximate log-determinant methods; the Smirnov/Anselin (2009) trace approximation is available as "moments". Three methods: "SE_classic", "SE_whichMin", and "SE_interp" are provided experimentally, the first to attempt to emulate the behaviour of Spatial Econometrics toolbox ML fitting functions. All use grids of log determinant values, and the latter two attempt to ameliorate some features of "SE_classic".
`Durbin`	Default = 'NULL'. If model is of `type = "sdm"`, `"sdem"` or `"slx"` then this argument should be a formula of the subset of explanatory variables to be spatially lagged in the right hand side part of the model. See `spsurml` for a similar argument.
`zero.policy`	Similar to the corresponding parameter of `lagsarlm` function in spatialreg package. If 'TRUE' assign zero to the lagged value of zones without neighbours, if 'FALSE' assign 'NA' - causing `pspatfit()` to terminate with an error. Default = 'NULL'.
`interval`	Search interval for autoregressive parameter. Default = 'NULL'.
`trs`	Similar to the corresponding parameter of `lagsarlm` function in spatialreg package. Default 'NULL', if given, a vector of powered spatial weights matrix traces output by `trW`.
`cor`	Type of temporal correlation for temporal data. Possible values are `"none"` (default) or `"ar1"`.
`dynamic`	Logical value to set a dynamic model. Dynamic models include a temporal lag of the dependent variable in the right-hand side of the equation. Default = 'FALSE'.
`demean`	Logical value to include a demeaning for panel data. Default = 'FALSE'. The demeaning is done previously to the estimation for both parametric and nonparametric terms. It is not possible to set `demean = TRUE` when spatio-temporal trends are included.
`eff_demean`	Type of demeaning for panel data. Possible values are `"individual"` (default), `"time"` or `"twoways"`.
`index`	Vector of variables indexing panel data. First variable corresponds to individuals and second variable corresponds to temporal coordinate (fast index). It follows the same rules than `plm` function in package plm.
`control`	List of extra control arguments. See Control Arguments section below.

Details

Function to estimate the model:

y = (\rho*W_{N} \otimes I_T) y + f(s_1,s_2,\tau_{t}) + X \beta + (W_{N} \otimes I_T) X \theta + \sum_{i = 1}^k g(z_i) + \sum_{i = 1}^k g((\gamma_i*W_{N} \otimes I_T) z_i) + \epsilon

where:

f(s_1,s_2,\tau_t) is a smooth spatio-temporal trend of the spatial coordinates s1,s_2 and of the temporal coordinates \tau_t.
X is a matrix including values of parametric covariates.
g(z_i) are non-parametric smooth functions of the covariates z_i.
W_N is the spatial weights matrix.
\rho is the spatial spillover parameter.
I_T is an identity matrix of order T (T=1 for pure spatial data).
\epsilon ~ N(0,R) where R = \sigma^2 I_T if errors are uncorrelated or it follows an AR(1) temporal autoregressive structure for serially correlated errors.

Including non-parametric terms

The non-parametric terms are included in formula using pspt(.) for spatial or spatio-temporal trends and pspl(.) for other non-parametric smooth additive terms. For example, if a model includes:

An spatio-temporal trend with variables long and lat as spatial coordinates,and year as temporal coordinate.
Two non-parametric covariates named empgrowth and serv.
Three parametric covariates named partrate, agri and cons.

Then, the formula should be written as (choosing default values for each term):

unrate ~ partrate + agri + cons + pspl(serv) + pspl(empgrowth) + pspt(long,lat,year)

For a spatial trend case, the term pspt(.) does not include a temporal coordinate, that is, in the previous example would be specified as pspt(long,lat).

How to use pspl() and pspt()

Note that both in pspl(.) and pspt(.), we have to include the number of knots, named nknots, which is the dimension of the basis used to represent the smooth term. The value of nknots should not be less than the dimension of the null space of the penalty for the term, see null.space.dimension and choose.k from mgcv package to know how to choose nknots.

In pspl(.) the default is nknots = 10, see the help of pspl function. In this term we can only include single variables, so if we want more than one non-parametric variable we will use a pspl(.) term for each nonparametric variable.

On the other hand, pspt(.) is used for spatial smoothing (when temporal coordinate is 'NULL') or spatio-temporal smoothing (when a variable is provided for the temporal coordinate). The default for the temporal coordinate is time = NULL, see the help of pspt, and the default number of knots are nknots = c(10, 10, 5). If only include spatial smoothing, nknots will be a length 2 vector indicating the basis for each spatial coordinate. For spatio-temporal smoothing, it will be a length 3 vector.

ANOVA descomposition

In many situations the spatio-temporal trend, given by f(s_1,s_2,\tau_t), can be very complex and the use of a multidimensional smooth function may not be flexible enough to capture the structure in the data. Furthermore, the estimation of this trend can become computationally intensive especially for large databases.
To solve this problem, Lee and Durban (2011) proposed an ANOVA-type decomposition of this spatio-temporal trend where spatial and temporal main effects, and second- and third-order interaction effects can be identified as:

f(s_1, s_2, \tau_t) = f_1(s_1) + f_2(s_2) + f_t(\tau_t) + f_{1,2}(s_1, s_2) + f_{1,t}(s_1, \tau_t) + f_{2,t}(s_2, \tau_t) + f_{1,2,t}(s_1, s_2, \tau_t)

In this equation the decomposition of the spatio-temporal trend is as follows:

Main effects given by the functions f_1(s_1), f_2(s_2) and f_t(\tau_t).
Second-order interaction effects given by the functions f_{1,2}(s_1,s_2), f_{1,t}(s_1,\tau_t) and f_{2,t}(s_2,\tau_t).
Third-order interaction effect given by the function f_{1,2,t}(s_1,s_2,\tau_t).

In this case, each effect can have its own degree of smoothing allowing a greater flexibility for the spatio-temporal trend. The ANOVA decomposition of the trend can be set as an argument in pspt(.) terms choosing psanova = TRUE.

For example to choose an ANOVA decomposition in the previous case we can set:

pspt(long, lat, year, nknots = c(18,18,8), psanova = TRUE)

In most empirical cases main effects functions are more flexible than interaction effects functions and therefore, the number of knots in B-Spline bases for interaction effects do not need to be as big as the number of knots for main effects. Lee et al., (2013) proposed a nested basis procedure in which the number of knots for the interaction effects functions are reduced using divisors such that the space spanned by B-spline bases used for interaction effects are a subset of the space spanned by B-spline bases used for main effects. The divisors can be specified as an argument in pspt(.) terms.
To do this, there are three arguments available inside pspt() to define the divisors. These arguments are named nest_sp1, nest_sp2 and nest_time, respectively. The value for these arguments are vector parameters including divisors of the nknots values.

For example, if we set nest_sp1 = c(1,2,2) between the arguments of pspl(.), we will have all knots for main effect of s_1, 18/2=9 knots for each second-order effect including s_1, and 8/2=4 knots for the third order effect including s_1. It is important that the vector of numbers will be integer divisors of the values in nknots. See section Examples for more details.

Eventually, any effect function can be excluded of the ps-anova spatio-temporal trend. To exclude main effects, the arguments f1_main, f2_main or ft_main have to be set to 'FALSE' (default='TRUE'). We can also exclude the second- and third-order effects functions setting to 'FALSE' the arguments f12_int, f1t_int, f2t_int or f12t_int in pspl(.).

All the terms included in the model are jointly fitted using Separation of Anisotropic Penalties (SAP) algorithm (see Rodriguez-Alvarez et al., (2015)) which allows to the mixed model reparameterization of the model. For type of models "sar", "sem", "sdm", "sdem", "sarar" or cor = "ar1", the parameters \rho, \lambda and \phi are numerically estimated using bobyqa function implemented in package minqa. In these cases, an iterative process between SAP and numerical optimization of \rho, \lambda and \phi is applied until convergence. See details in Minguez et al., (2018).

Plotting non-parametric terms

To plot the non-linear functions corresponding to non-parametric terms we need to compute the fitted values, and standard erros, using fit_terms() function and, afterwards, use plot_terms() function to plot the non-linear functions.
An example of how plot the functions of non-parametric terms given by "var1" and "var2" variables is given by the next lines of code (it is assumed that a previous model has been fitted using pspatfit(.) and saved as an object named model):

list_varnopar <- c("var1", "var2")
terms_nopar <- fit_terms(model, list_varnopar)
plot_terms(terms_nopar, data)

The data argument of plot_terms() usually corresponds to the dataframe used to fitted the model although a different database can be used to plot the non-parametric terms.

Spatial impacts

For the spatial models given by type = "sar", "sdm", "sdem", "sarar" or "slx" it is possible to compute spatial spillovers as usual in spatial econometric specifications. Nevertheless, in this case we need to distinguish between parametric and non-parametric covariates when computing spatial impacts.

spatial impacts for parametric covariates
In this case, the spatial impacts are computed in the usual way using simulation. See LeSage and Page (2009) for computational details. The function impactspar() computes the direct, indirect and total impacts for parametric covariates and return and object similar to the case of spatialreg and spsur packages. The inference for "sar", "sdm", and "sarar" types is based on simulations and for "slx" and "sdem" types the standard errors or total impacts are computed using the variance-covariance matrix of the fitted model. The summary() method can be used to present the the complete table of spatial impacts in this parametric case. See the help of impactspar to know the additional arguments of the function. A little example is given in the next lines of code:

imp_parvar <- impactspar(MODEL, listw = W)
summary(imp_parvar)
spatial impacts for non-parametric covariates
In this case direct, indirect and total spatial impacts functions are obtained using impactsnopar. The details of computation and inference can be obtained from the help of impactsnopar. The argument viewplot of impactsnopar have to be set as 'TRUE' to plot the spatial impacts functions. Another way to get the same plots is using plot_impactsnopar function with the output of impactsnopar. Next lines give an example of both cases:

imp_nparvar <- impactsnopar(MODEL, listw = W, viewplot = TRUE)
imp_nparvar <- impactsnopar(MODEL, listw = W, viewplot = FALSE)
plot_impactsnopar(imp_nparvar, data = DATA)

Value

A list object of class pspatreg

`call`	Matched call.
`terms`	The terms object used.
`contrasts`	(only where relevant) the contrasts used for parametric covariates.
`xlevels`	(only where relevant) a record of the levels of the parametric factors used in fitting.
`data`	dataframe used as database.
`nsp`	number of spatial observations.
`nt`	number of temporal observations. It is set to `nt=1` for spatial data.
`nfull`	total number of observations.
`edftot`	Equivalent degrees of freedom for the whole model.
`edfspt`	Equivalent degrees of freedom for smooth spatio-temporal or spatial trend.
`edfnopar`	Equivalent degrees of freedom for non-parametric covariates.
`psanova`	TRUE if spatio-temporal or spatial trend is PS-ANOVA.
`type`	Value of `type` argument in the call to `pspatfit`.
`listw`	Value of `listw` argument in the call to `pspatfit`.
`Durbin`	Value of `Durbin` argument in the call to `pspatfit`.
`cor`	Value of `cor` argument in the call to `pspatfit`.
`dynamic`	Value of `dynamic` argument in the call to `pspatfit`.
`demean`	Value of `demean` argument in the call to `pspatfit`.
`eff_demean`	Value of `eff_demean` argument in the call to `pspatfit`.
`index`	Value of `index` argument in the call to `pspatfit`.
`bfixed`	Estimated betas corresponding to fixed effects in mixed model.
`se_bfixed`	Standard errors of fixed betas.
`brandom`	Estimated betas corresponding to random effects in mixed model.
`se_brandom`	Standard errors of random betas.
`vcov_fr`	Covariance matrix of fixed and random effects using frequentist or sandwich method.
`vcov_by`	Covariance matrix of fixed and random effects using bayesian method.
`rho`	Estimated rho for spatial lag of the dependent variable.
`se_rho`	Standard error of `rho`.
`delta`	Estimated delta for spatial error models.
`se_delta`	Standard error of `delta`.
`phi`	Estimated phi. If `cor="none"` always `phi=0`.
`se_phi`	Standard error of `phi`.
`fitted.values`	Vector of fitted values of the dependent variable.
`se_fitted.values`	Vector of standard errors of `fitted.values`.
`fitted.values_Ay`	Vector of fitted values of the spatial lag of dependent variable: `(\rho*W_N \otimes I_T) y`.
`se_fitted.values_Ay`	Vector of standard errors of `fitted.values_Ay`.
`residuals`	Vector of residuals.
`df.residual`	Equivalent degrees of freedom for `residuals`.
`sig2`	Residual variance computed as SSR/df.residual.
`llik`	Log-likelihood value.
`llik_reml`	Restricted log-likelihood value.
`aic`	Akaike information criterion.
`bic`	Bayesian information criterion.
`sp1`	First spatial coordinate.
`sp2`	Second spatial coordinate.
`time`	Time coordinate.
`y`	Dependent variable.
`X`	Model matrix for fixed effects.
`Z`	Model matrix for random effects.

Control Arguments

`optim`	method of estimation: `"llik_reml"` (default) or `"llik"`.
`typese`	method to compute standard errors. `"sandwich"` or `"bayesian"` (default). See Fahrmeir et al, pp. 375 for details of computations.
`vary_init`	Initial value of the noise variance in the model. Default = `NULL`.
`trace`	A logical value set to TRUE to show intermediate results during the estimation process. Default = FALSE.
`tol1`	Numerical value for the tolerance of convergence of penalization parameters during the estimation process. Default 1e-3. This tolerance is only used for small samples (<= 500 observations).
`tol2`	Numerical value for the tolerance of convergence of total estimated degrees of freedom ("edftot") during the estimation process. Default 1e-1. This tolerance is used for medium or big samples (> 500 observations).
`tol3`	Numerical value for the tolerance of convergence of spatial and correlation parameters during the estimation process. Default 1e-2.
`maxit`	An integer value for the maximum number of iterations until convergence. Default = 200.
`rho_init`	An initial value for `rho` parameter. Default 0.
`delta_init`	An initial value for `delta` parameter. Default 0.
`phi_init`	An initial value for `phi` parameter. Default 0.
`Imult`	default 2; used for preparing the Cholesky decompositions for updating in the Jacobian function
`super`	if `NULL` (default), set to `FALSE` to use a simplicial decomposition for the sparse Cholesky decomposition and method "Matrix_J", set to `as.logical(NA)` for method "Matrix", if `TRUE`, use a supernodal decomposition
`cheb_q`	default 5; highest power of the approximating polynomial for the Chebyshev approximation
`MC_p`	default 16; number of random variates
`MC_m`	default 30; number of products of random variates matrix and spatial weights matrix
`spamPivot`	default "MMD", alternative "RCM"
`in_coef`	default 0.1, coefficient value for initial Cholesky decomposition in "spam_update"
`type`	default "MC", used with method "moments"; alternatives "mult" and "moments", for use if trs is missing
`correct`	default `TRUE`, used with method "moments" to compute the Smirnov/Anselin correction term
`trunc`	default `TRUE`, used with method "moments" to truncate the Smirnov/Anselin correction term
`SE_method`	default "LU", may be "MC"
`nrho`	default 200, as in SE toolbox; the size of the first stage lndet grid; it may be reduced to for example 40
`interpn`	default 2000, as in SE toolbox; the size of the second stage lndet grid
`SElndet`	default `NULL`, may be used to pass a pre-computed SE toolbox style matrix of coefficients and their lndet values to the "SE_classic" and "SE_whichMin" methods
`LU_order`	default `FALSE`; used in "LU_prepermutate", note warnings given for lu method
`pre_eig`	default `NULL`; may be used to pass a pre-computed vector of eigenvalues

Author(s)

Roman Minguez	roman.minguez@uclm.es
Roberto Basile	roberto.basile@univaq.it
Maria Durban	mdurban@est-econ.uc3m.es
Gonzalo Espana-Heredia	gehllanza@gmail.com

References

Basile, R.; Durban, M.; Minguez, R.; Montero, J. M.; and Mur, J. (2014). Modeling regional economic dynamics: Spatial dependence, spatial heterogeneity and nonlinearities. Journal of Economic Dynamics and Control, (48), 229-245. <doi:10.1016/j.jedc.2014.06.011>
Eilers, P. and Marx, B. (1996). Flexible Smoothing with B-Splines and Penalties. Statistical Science, (11), 89-121.
Eilers, P. and Marx, B. (2021). Practical Smoothing. The Joys of P-Splines. Cambridge University Press.
Fahrmeir, L.; Kneib, T.; Lang, S.; and Marx, B. (2021). Regression. Models, Methods and Applications (2nd Ed.). Springer.
Lee, D. and Durban, M. (2011). P-Spline ANOVA Type Interaction Models for Spatio-Temporal Smoothing. Statistical Modelling, (11), 49-69. <doi:10.1177/1471082X1001100104>
Lee, D. J., Durban, M., and Eilers, P. (2013). Efficient two-dimensional smoothing with P-spline ANOVA mixed models and nested bases. Computational Statistics & Data Analysis, (61), 22-37. <doi:10.1016/j.csda.2012.11.013>
LeSage, J. and Pace, K. (2009). Introduction to Spatial Econometrics. CRC Press, Boca Raton.
Minguez, R.; Basile, R. and Durban, M. (2020). An Alternative Semiparametric Model for Spatial Panel Data. Statistical Methods and Applications, (29), 669-708. <doi:10.1007/s10260-019-00492-8>
Montero, J., Minguez, R., and Durban, M. (2012). SAR models with nonparametric spatial trends: A P-Spline approach. Estadistica Espanola, (54:177), 89-111.
Rodriguez-Alvarez, M. X.; Kneib, T.; Durban, M.; Lee, D.J. and Eilers, P. (2015). Fast smoothing parameter separation in multidimensional generalized P-splines: the SAP algorithm. Statistics and Computing 25 (5), 941-957. <doi:10.1007/s11222-014-9464-2>
Wood, S.N. (2017). Generalized Additive Models. An Introduction with R (second edition). CRC Press, Boca Raton.

Examples

################################################

##########################
library(pspatreg)
###############################################
# Examples using spatial data of Ames Houses.
###############################################
# Getting and preparing the data
library(spdep)
library(sf)
ames <- AmesHousing::make_ames() # Raw Ames Housing Data
ames_sf <- st_as_sf(ames, coords = c("Longitude", "Latitude"))
ames_sf$Longitude <- ames$Longitude
ames_sf$Latitude <- ames$Latitude
ames_sf$lnSale_Price <- log(ames_sf$Sale_Price)
ames_sf$lnLot_Area <- log(ames_sf$Lot_Area)
ames_sf$lnTotal_Bsmt_SF <- log(ames_sf$Total_Bsmt_SF+1)
ames_sf$lnGr_Liv_Area <- log(ames_sf$Gr_Liv_Area)
ames_sf1 <- ames_sf[(duplicated(ames_sf$Longitude) == FALSE), ]

####  GAM pure with pspatreg
form1 <- lnSale_Price ~ Fireplaces + Garage_Cars +
          pspl(lnLot_Area, nknots = 20) + 
          pspl(lnTotal_Bsmt_SF, nknots = 20) +
          pspl(lnGr_Liv_Area, nknots = 20)    
gampure <- pspatfit(form1, data = ames_sf1)
summary(gampure)

######################  Get Non-parametric terms of GAM with pspatreg
list_varnopar <- c("lnLot_Area", "lnTotal_Bsmt_SF", 
"lnGr_Liv_Area")
terms_nopar <- fit_terms(gampure, list_varnopar, intercept = TRUE)
######################  Plot non-parametric terms
plot_terms(terms_nopar, ames_sf1)

 
########### Constructing the spatial weights matrix
coord_sf1 <- cbind(ames_sf1$Longitude, ames_sf1$Latitude)
k5nb <- knn2nb(knearneigh(coord_sf1, k = 5, 
                          longlat = TRUE, use_kd_tree = FALSE), sym = TRUE)
lw_ames <- nb2listw(k5nb, style = "W", 
                  zero.policy = FALSE)
                  
#####################  GAM + SAR Model
gamsar <- pspatfit(form1, data = ames_sf1, 
                   type = "sar", listw = lw_ames,
                   method = "Chebyshev")
summary(gamsar)
                  
######### Non-Parametric Total, Direct and Indirect impacts
### with impactsnopar(viewplot = TRUE)
nparimpacts <- impactsnopar(gamsar, 
                            listw = lw_ames, 
                            viewplot = TRUE)
############ Non-Parametric Total, Direct and Indirect impacts
### with impactsnopar(viewplot = FALSE) and using plot_impactsnopar()
nparimpacts <- impactsnopar(gamsar, listw = lw_ames, viewplot = FALSE)
plot_impactsnopar(nparimpacts, data = ames_sf1, smooth = TRUE)

###################### Parametric Total, Direct and Indirect impacts
parimpacts <- impactspar(gamsar, listw = lw_ames)
summary(parimpacts)

###############################################
### Models with 2d spatial trend
form2 <- lnSale_Price ~ Fireplaces + Garage_Cars +
          pspl(lnLot_Area, nknots = 20) + 
          pspl(lnTotal_Bsmt_SF, nknots = 20) +
          pspl(lnGr_Liv_Area, nknots = 20) +
          pspt(Longitude, Latitude, 
               nknots = c(10, 10), 
               psanova = FALSE)
#####################  GAM + GEO Model
gamgeo2d <- pspatfit(form2, data = ames_sf1)
summary(gamgeo2d)

gamgeo2dsar <- pspatfit(form2, data = ames_sf1,
                        type = "sar", 
                        listw = lw_ames, 
                        method = "Chebyshev")
summary(gamgeo2dsar)
####### plot spatial trend for spatial point coordinate
plot_sp2d(gamgeo2dsar, data = ames_sf1)
### Models with psanova 2d spatial trend
form3 <- lnSale_Price ~ Fireplaces + Garage_Cars +
          pspl(lnLot_Area, nknots = 20) + 
          pspl(lnTotal_Bsmt_SF, nknots = 20) +
          pspl(lnGr_Liv_Area, nknots = 20) +
          pspt(Longitude, Latitude, 
               nknots = c(10, 10), 
               psanova = TRUE)
gamgeo2danovasar <- pspatfit(form3, data = ames_sf1,
                        type = "sar", 
                        listw = lw_ames, method = "Chebyshev")
summary(gamgeo2danovasar)
####### plot spatial trend for spatial point coordinate
plot_sp2d(gamgeo2danovasar, data = ames_sf1, 
addmain = TRUE, addint = TRUE)

## Comparison between models
anova(gampure, gamsar, gamgeo2d, gamgeo2dsar,
gamgeo2danovasar, lrtest = FALSE)

###############################################
###################### Examples using a panel data of rate of
###################### unemployment for 103 Italian provinces in 1996-2019.
###############################################
## load spatial panel and Wsp_it
## 103 Italian provinces. Period 1996-2019
data(unemp_it, package = "pspatreg")
## Wsp_it is a matrix. Create a neighboord list 
lwsp_it <- spdep::mat2listw(Wsp_it)
 ### Models with spatio-temporal trend
 ### Spatio-temporal semiparametric ANOVA model without spatial lag
 ### Interaction terms f12,f1t,f2t and f12t with nested basis
 ### Remark: nest_sp1, nest_sp2 and nest_time must be divisors of nknots
 form4 <- unrate ~ partrate + agri + cons +
                   pspl(serv, nknots = 15) + 
                   pspl(empgrowth, nknots = 20) +
                   pspt(long, lat, year, 
                        nknots = c(18, 18, 12),
                        psanova = TRUE, 
                        nest_sp1 = c(1, 2, 3), 
                        nest_sp2 = c(1, 2, 3),
                        nest_time = c(1, 2, 2))
 sptanova <- pspatfit(form4, data = unemp_it)
 summary(sptanova)
 
### Create sf object to make the plot 
### of spatio-temporal trends
library(sf)
unemp_it_sf <- st_as_sf(dplyr::left_join(
                              unemp_it, 
                              map_it,  
                        by = c("prov" = "COD_PRO")))
###### Plot spatio-temporal trends for different years
plot_sp3d(sptanova, data = unemp_it_sf, 
          time_var = "year", 
          time_index = c(1996, 2005, 2019),
          addmain = FALSE, addint = FALSE)
###### Plot of spatio-temporal trend, main effects 
######      and interaction effect for a year
plot_sp3d(sptanova, data = unemp_it_sf, 
          time_var = "year", 
          time_index = c(2019),
          addmain = TRUE, addint = TRUE)
          
###### Plot of temporal trends for each province
plot_sptime(sptanova, 
            data = unemp_it, 
            time_var = "year", 
            reg_var = "prov")

 
 ###############################################
 ### Spatio-temporal semiparametric ANOVA model without spatial lag
 ### Now we repeat previous spatio-temporal model but 
 ### restricting some interactions
 ### Interaction terms f12,f1t and f12t with nested basis
 ### Interaction term f2t restricted to 0
 
  form5 <- unrate ~ partrate + agri + cons + empgrowth +
                  pspl(serv, nknots = 15) + 
                  pspt(long, lat, year, 
                       nknots = c(18, 18, 6), 
                       psanova = TRUE,
                       nest_sp1 = c(1, 2, 3), 
                       nest_sp2 = c(1, 2, 3),
                       nest_time = c(1, 2, 2), 
                       f2t_int = FALSE)
 ## Add sar specification and ar1 temporal correlation 
 sptanova2_sar_ar1 <- pspatfit(form5, data = unemp_it, 
                              listw = lwsp_it, 
                              type = "sar",
                              cor = "ar1")
summary(sptanova2_sar_ar1)                     
################ Comparison with parametric panels            
######################  Demeaning (Within Estimators)
formpar <- unrate ~ partrate + agri + cons
# Not demeaning model
param <- pspatfit(formpar, data = unemp_it, listw = lwsp_it)
summary(param)
# Demeaning model
param_dem <- pspatfit(formpar, data = unemp_it,
                      demean = TRUE,
                      index = c("prov", "year"),
                      eff_demean = "individual" )
summary(param_dem)
# Compare results with plm package
param_plm <- plm::plm(formula = formpar,
                      data = unemp_it,
                      index = c("prov", "year"),
                      effect = "individual",
                      model = "within")
summary(param_plm)                                              
param_dem_time <- pspatfit(formpar, 
                      data = unemp_it, 
                      listw = lwsp_it,
                      demean = TRUE,
                      eff_demean = "time",
                      index = c("prov", "year"))
summary(param_dem_time)
param_plm_time <- plm::plm(formula = formpar,
                      data = unemp_it,
                      index = c("prov", "year"),
                      effect = "time",
                      model = "within")
summary(param_plm_time)
param_dem_twoways <- pspatfit(formpar, data = unemp_it,
                      demean = TRUE,
                      eff_demean = "twoways",
                      index = c("prov", "year") )
summary(param_dem_twoways)
param_plm_twoways <- plm::plm(formula = formpar,
                      data = unemp_it,
                      index = c("prov", "year"),
                      effect = "twoways",
                      model = "within")
summary(param_plm_twoways)
##### Demeaning with nonparametric covariates
formgam <- unrate ~ partrate + agri + cons +  
                    pspl(serv, nknots = 15) +
                    pspl(empgrowth, nknots = 20)
                    
gam_dem <- pspatfit(formula = formgam,
                      data = unemp_it,
                      demean = TRUE,
                      index = c("prov", "year"))
summary(gam_dem)   
# Compare with GAM pure without demeaning                    
gam <- pspatfit(formula = formgam,
                 data = unemp_it)
summary(gam)

## Demeaning with type = "sar" model
gamsar_dem <- pspatfit(formula = formgam,
                      data = unemp_it,
                      type = "sar", 
                      listw = lwsp_it,
                      demean = TRUE,
                      index = c("prov", "year"))
summary(gamsar_dem)

[Package pspatreg version 1.1.2 Index]