Dear all,
though this is an old issue by now (but still one of the first google results on the topic), here is some code answering one of the questions: should you demean the spatial lag or lag the demeaned variable. As far as I understand the following results, it does not matter. Both ways produce identical results, which both conform to the LSDV approach. Please find the example code below: library(spdep) library(splm) library(texreg) ### Demean function dm<-function(x, id){ res<-x-ave(x,id,FUN=function(x) mean(x, na.rm=T)) res } data(Produc, package="plm") data(usaww) usa.lw<-mat2listw(usaww) Produc.pd<-pdata.frame(Produc, index=c("state", "year")) Produc.pd$Wpcap<-slag(Produc.pd$pcap, usa.lw) Produc.pd$Wpc<-slag(Produc.pd$pc, usa.lw) Produc.pd$Wunemp<-slag(Produc.pd$unemp, usa.lw) #### LSDV #### lsdv.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp +state, data=Produc.pd) summary(lsdv.mod) #### Demean after lag #### dal.mod<-plm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd, effect="individual", model="within") summary(dal.mod) # Manual demeaning Produc.pd[,4:ncol(Produc.pd)]<-apply(Produc.pd[,4:ncol(Produc.pd)], 2, function(x) dm(x, Produc.pd$state)) dal2.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd) summary(dal2.mod) #### Lag after demean #### # Replace slags by demeaned slags Produc.pd$Wpcap<-slag(Produc.pd$pcap, usa.lw) Produc.pd$Wpc<-slag(Produc.pd$pc, usa.lw) Produc.pd$Wunemp<-slag(Produc.pd$unemp, usa.lw) lad.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd) summary(lad.mod) screenreg(l=list(lsdv.mod, dal.mod, dal2.mod, lad.mod), omit.coef = "state", digits=6) #### Test if both ways produce identical vectors #### Produc.pd<-pdata.frame(Produc, index=c("state", "year")) x<-Produc.pd$gsp # lag after demean dx<-dm(x, Produc.pd$state) wdx<-slag(dx, usa.lw) # deaman after lag wx<-slag(x, usa.lw) dwx<-dm(wx, Produc.pd$state) # Test for equality all.equal(wdx, dwx) Best, Tobi Tobias Rüttenauer TU Kaiserslautern Erwin-Schrödinger-Straße 57 67663 Kaiserslautern [hidden email] Tel.: 0631 205 5785 > -----Ursprüngliche Nachricht----- > Von: R-sig-Geo [mailto:[hidden email]] Im Auftrag von > Roger Bivand > Gesendet: 30 August 2016 08:44 > An: [hidden email] > Cc: [hidden email] > Betreff: Re: [R-sig-Geo] SLX model for splm package in R > > On Tue, 30 Aug 2016, [hidden email] wrote: > > > Dear all, > > > > Impacts in SLX model are calculated directly by parameters estimation: > > > > http://onlinelibrary.wiley.com/doi/10.1111/jors.12188/abstract > > > > Beta1 gives direct impact and Beta 2 the indirect: > > > > y = beta1*X + beta2*WX + e > > > > Yes, but you need the total impact beta1+beta2, but more work to infer > that sum, which is what you need to assess the impact of X. > > > I understand W must multiply T times the X matrix (T is the lenght of > > time series for a panel data). And variables are demeaned after WX is > > calculated. > > Kronecker product of sparse matrices. Please contribute a patch for > create_WX if desirable. > > Please motivate the demean-after-lag concusion. Will this not affect the > direct link between X and WX? > > Roger > > > > > Daniel. > > > > ----- Mensagem original ----- > > > > > > De: "Roger Bivand" <[hidden email]> > > Para: "Burcu Ozuduru" <[hidden email]> > > Cc: [hidden email], [hidden email] > > Enviadas: Segunda-feira, 29 de Agosto de 2016 4:13:26 > > Assunto: Re: [R-sig-Geo] SLX model for splm package in R > > > > This is far too ad-hoc a discussion. For Spatial panel SLX, please > > read > > carefully: > > > > http://dx.doi.org/10.1016/j.spasta.2014.02.002 > > DOI 10.1007/s10109-015-0217-3 > > DOI: 10.1111/jors.12188 > > > > for a view of the state of play. One question - should you demean the > > X variables before taking spatial lags? This does not apply to > > cross-section models for which spdep::create_WX() was written, for > > obvious reasons. A second question - how do you handle the calculation > > of total impacts in this setting? In cross-section models, this is > > reasonably easy, by linear combination. But in spatial panel SLX > > models, is it easy? What does the literature suggest? > > > > On Mon, 29 Aug 2016, Burcu Ozuduru wrote: > > > >> Hi, > >> > >> This is a great question because I faced the same problem a while back. > >> What I did is to obtain the spatial weight matrix as (have spdep > >> installed). > >> > >> S47matrix<-W1%*%S47 > >> > >> and then put it in my data frame and equation. Btw, above W1 is a > >> spatial weights matrix I obtained from a neighbor list > >> (W1=as.matrix(as_dgRMatrix_listw(ccListw))) and S47 is the dependent > >> variable. > >> > >> m2<- > data.frame(S41,S46,S47,S65,S68,S47matrix,MAD10000,NQPDA10000,TPBt > >> A10000,DivA10000,Lnk10000) > >> summary(m2<-zeroinfl(S47~MAD10000+NQPDA10000+S47matrix, > >> dist="negbin")) > >> > >> I also wanted to ask whether I was on the right track. Professor > >> Bivand, if this is a common problem could you direct us to a relevant > >> page -if available and if at all possible- please? Your advice would > >> be sincerely appreciated. > > > > This is an inappropriate question for this thread (SLX spatial panel > > models). You must refer to the literature, and understand what you are > > doing: > > > > zeroinfl(S47~MAD10000+NQPDA10000+S47matrix, dist="negbin") > > > > is completely wrong, as S47matrix is Wy, the spatial lag of the > > dependent variable, and this spatial lag model of response S47 cannot > > be estimated using pscl::zeroinfl (do say which package provides the > > functions you mention). Your only recourse is to write your own > > Bayesian model. It is even extremely likely that this model doesn't > > make any substantive sense, unless you have a good understanding about > > why S47 at i (for all i) should be simultaneously caused by the full > > set of values of S47. A model in the error might make sense (a zero > > inflated negbin model with an additive mrf random effect (maybe > R2BayesX 10.18637/jss.v014.i11). > > > > An SLX zero inflated negbin model could maybe apply, in which case > > pscl::zeroinfl might be applicable, but then you'd be including the > > lags of MAD10000 and NQPDA10000, not the response, and inference on > > total impacts of these variables could be challenging rather than > > obvious (MCMC summaries of the summed samples of the coefficients > > would work, but here no MCMC I think - R2BayesX might help). > > > > Most of these issues are far more complex than students or their > > supervisors understand, do ask an experienced applied statistician > > (I'm not an applied statistician, don't rely on my suggestions, you > > yourself are responsible for your results). > > > > Hope this clarifies, > > > > Roger > > > >> > >> Kind Regards- > >> > >> Thanks- > >> > >> Burcu > >> > >> --- > >> > >> Burcu H. Ozuduru, > >> Gazi University, Faculty of Architecture, Department of City and > >> Regional Planning Maltepe-Ankara Turkey > >> t: +90-312-5823701 > >> e: [hidden email] > >> http://websitem.gazi.edu.tr/bozuduru > >> > >> > >> On Mon, Aug 29, 2016 at 12:42 AM, Roger Bivand > <[hidden email]> wrote: > >> > >>> On Sun, 28 Aug 2016, [hidden email] wrote: > >>> > >>> Dear colleagues, > >>>> > >>>> Does anyone know how to implement Spatial Lag of X (SLX) model for > >>>> panel data using splm package? I tried the following to create a WX > >>>> matrix and then apply Fixed and Random Effects: > >>>> > >>>> tb <- read.csv2("panel.csv", header = TRUE, sep = ";", dec = ".") > >>>> time <- length(unique(tb$year)) X <- as.matrix(tb[ ,3:20]) > >>>> > >>> > >>> class(X) is not "pseries", is it? > >>> > >>> WX <- slag(X, listw, maxlag = 1) > >>>> > >>>> I see that slag is applicable for vectors of pseries, but this > >>>> led to the error message: "Error in UseMethod("slag"): no applicable > method > >>>> for 'slag' applied to an object of class "c('matrix', 'double', 'numeric')" > >>>> > >>> > >>> slag.default() takes a vector, slag.pseries() a "pseries" object. > >>> > >>> > >>>> Is there any way to implement this SLX model for panel series? > >>>> > >>> > >>> Did you look at spdep::create_WX? > >>> > >>> > >>>> Regards, > >>>> > >>>> Daniel. > >>>> > >>>> [[alternative HTML version deleted]] > >>>> > >>>> > >>> Please do not post HTML - like insects, HTML carries infectious > >>> and wastes bandwidth and server processing capacity. > >>> > >>> _______________________________________________ > >>>> R-sig-Geo mailing list > >>>> [hidden email] > >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >>>> > >>>> > >>> -- > >>> Roger Bivand > >>> Department of Economics, Norwegian School of Economics, > >>> Helleveien 30, N-5045 Bergen, Norway. > >>> voice: +47 55 95 93 55; fax +47 55 95 91 00 > >>> e-mail: [hidden email] > >>> http://orcid.org/0000-0003-2392-6140 > >>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en > >>> http://depsy.org/person/434412 > >>> > >>> _______________________________________________ > >>> R-sig-Geo mailing list > >>> [hidden email] > >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >>> > >> > > > > > > -- > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N-5045 Bergen, Norway. > voice: +47 55 95 93 55; fax +47 55 95 91 00 > e-mail: [hidden email] > http://orcid.org/0000-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en > http://depsy.org/person/434412 > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Dear Daniel,
yes, that is basically the same answer to the problem. With time-series SLX you need to create a pseries and then apply slag. A second solution (using a data.frame) would be to apply lag.listw to each cross-section. Using my example below, this would be: # Reorder data by cross-sections, then time oo<-order(Produc$year, Produc$state) Produc<-Produc[oo,] # Build indices indic<-seq(1,length(unique(Produc$year))) inde<-as.numeric(rep(indic, each=length(unique(Produc$state)))) # Construct WX Produc$Wpcap <- unlist(tapply(as.vector(Produc$pcap), inde, function(u) lag.listw(usa.lw, u, NAOK=T), simplify=TRUE)) Best, Tobi Von: Daniel Furlan Amaral [mailto:[hidden email]] Gesendet: 29 December 2017 17:27 An: Tobias Rüttenauer <[hidden email]> Betreff: Re: [R-sig-Geo] SLX model for splm package in R Dear Tobias, Thank you for the answer. Coincidently, I published my answer on Stack Overflow today and answered my own question. I think we came up with the same solution. Please see by the link what you think: https://stackoverflow.com/questions/46004670/slx-model-spatial-econometrics-with-panel-in-r-data-using-splm-package-and-sla Best regards, Daniel 2017-12-22 8:18 GMT-02:00 Tobias Rüttenauer <[hidden email] <mailto:[hidden email]> >: Dear all, though this is an old issue by now (but still one of the first google results on the topic), here is some code answering one of the questions: should you demean the spatial lag or lag the demeaned variable. As far as I understand the following results, it does not matter. Both ways produce identical results, which both conform to the LSDV approach. Please find the example code below: library(spdep) library(splm) library(texreg) ### Demean function dm<-function(x, id){ res<-x-ave(x,id,FUN=function(x) mean(x, na.rm=T)) res } data(Produc, package="plm") data(usaww) usa.lw<-mat2listw(usaww) Produc.pd<-pdata.frame(Produc, index=c("state", "year")) Produc.pd$Wpcap<-slag(Produc.pd$pcap, usa.lw) Produc.pd$Wpc<-slag(Produc.pd$pc, usa.lw) Produc.pd$Wunemp<-slag(Produc.pd$unemp, usa.lw) #### LSDV #### lsdv.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp +state, data=Produc.pd) summary(lsdv.mod) #### Demean after lag #### dal.mod<-plm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd, effect="individual", model="within") summary(dal.mod) # Manual demeaning Produc.pd[,4:ncol(Produc.pd)]<-apply(Produc.pd[,4:ncol(Produc.pd)], 2, function(x) dm(x, Produc.pd$state)) dal2.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd) summary(dal2.mod) #### Lag after demean #### # Replace slags by demeaned slags Produc.pd$Wpcap<-slag(Produc.pd$pcap, usa.lw) Produc.pd$Wpc<-slag(Produc.pd$pc, usa.lw) Produc.pd$Wunemp<-slag(Produc.pd$unemp, usa.lw) lad.mod<-lm(gsp~pcap+pc+unemp + Wpcap+Wpc+Wunemp, data=Produc.pd) summary(lad.mod) screenreg(l=list(lsdv.mod, dal.mod, dal2.mod, lad.mod), omit.coef = "state", digits=6) #### Test if both ways produce identical vectors #### Produc.pd<-pdata.frame(Produc, index=c("state", "year")) x<-Produc.pd$gsp # lag after demean dx<-dm(x, Produc.pd$state) wdx<-slag(dx, usa.lw) # deaman after lag wx<-slag(x, usa.lw) dwx<-dm(wx, Produc.pd$state) # Test for equality all.equal(wdx, dwx) Best, Tobi Tobias Rüttenauer TU Kaiserslautern Erwin-Schrödinger-Straße 57 67663 Kaiserslautern [hidden email] <mailto:[hidden email]> Tel.: 0631 205 5785 > -----Ursprüngliche Nachricht----- > Von: R-sig-Geo [mailto:[hidden email] <mailto:[hidden email]> ] Im Auftrag von > Roger Bivand > Gesendet: 30 August 2016 08:44 > An: [hidden email] <mailto:[hidden email]> > Cc: [hidden email] <mailto:[hidden email]> > Betreff: Re: [R-sig-Geo] SLX model for splm package in R > > On Tue, 30 Aug 2016, [hidden email] <mailto:[hidden email]> wrote: > > > Dear all, > > > > Impacts in SLX model are calculated directly by parameters estimation: > > > > http://onlinelibrary.wiley.com/doi/10.1111/jors.12188/abstract > > > > Beta1 gives direct impact and Beta 2 the indirect: > > > > y = beta1*X + beta2*WX + e > > > > Yes, but you need the total impact beta1+beta2, but more work to infer > that sum, which is what you need to assess the impact of X. > > > I understand W must multiply T times the X matrix (T is the lenght of > > time series for a panel data). And variables are demeaned after WX is > > calculated. > > Kronecker product of sparse matrices. Please contribute a patch for > create_WX if desirable. > > Please motivate the demean-after-lag concusion. Will this not affect the > direct link between X and WX? > > Roger > > > > > Daniel. > > > > ----- Mensagem original ----- > > > > > > De: "Roger Bivand" <[hidden email] <mailto:[hidden email]> > > > Para: "Burcu Ozuduru" <[hidden email] <mailto:[hidden email]> > > > Cc: [hidden email] <mailto:[hidden email]> , [hidden email] <mailto:[hidden email]> > > Enviadas: Segunda-feira, 29 de Agosto de 2016 4:13:26 > > Assunto: Re: [R-sig-Geo] SLX model for splm package in R > > > > This is far too ad-hoc a discussion. For Spatial panel SLX, please > > read > > carefully: > > > > http://dx.doi.org/10.1016/j.spasta.2014.02.002 > > DOI 10.1007/s10109-015-0217-3 > > DOI: 10.1111/jors.12188 > > > > for a view of the state of play. One question - should you demean the > > X variables before taking spatial lags? This does not apply to > > cross-section models for which spdep::create_WX() was written, for > > obvious reasons. A second question - how do you handle the calculation > > of total impacts in this setting? In cross-section models, this is > > reasonably easy, by linear combination. But in spatial panel SLX > > models, is it easy? What does the literature suggest? > > > > On Mon, 29 Aug 2016, Burcu Ozuduru wrote: > > > >> Hi, > >> > >> This is a great question because I faced the same problem a while back. > >> What I did is to obtain the spatial weight matrix as (have spdep > >> installed). > >> > >> S47matrix<-W1%*%S47 > >> > >> and then put it in my data frame and equation. Btw, above W1 is a > >> spatial weights matrix I obtained from a neighbor list > >> (W1=as.matrix(as_dgRMatrix_listw(ccListw))) and S47 is the dependent > >> variable. > >> > >> m2<- > data.frame(S41,S46,S47,S65,S68,S47matrix,MAD10000,NQPDA10000,TPBt > >> A10000,DivA10000,Lnk10000) > >> summary(m2<-zeroinfl(S47~MAD10000+NQPDA10000+S47matrix, > >> dist="negbin")) > >> > >> I also wanted to ask whether I was on the right track. Professor > >> Bivand, if this is a common problem could you direct us to a relevant > >> page -if available and if at all possible- please? Your advice would > >> be sincerely appreciated. > > > > This is an inappropriate question for this thread (SLX spatial panel > > models). You must refer to the literature, and understand what you are > > doing: > > > > zeroinfl(S47~MAD10000+NQPDA10000+S47matrix, dist="negbin") > > > > is completely wrong, as S47matrix is Wy, the spatial lag of the > > dependent variable, and this spatial lag model of response S47 cannot > > be estimated using pscl::zeroinfl (do say which package provides the > > functions you mention). Your only recourse is to write your own > > Bayesian model. It is even extremely likely that this model doesn't > > make any substantive sense, unless you have a good understanding about > > why S47 at i (for all i) should be simultaneously caused by the full > > set of values of S47. A model in the error might make sense (a zero > > inflated negbin model with an additive mrf random effect (maybe > R2BayesX 10.18637/jss.v014.i11). > > > > An SLX zero inflated negbin model could maybe apply, in which case > > pscl::zeroinfl might be applicable, but then you'd be including the > > lags of MAD10000 and NQPDA10000, not the response, and inference on > > total impacts of these variables could be challenging rather than > > obvious (MCMC summaries of the summed samples of the coefficients > > would work, but here no MCMC I think - R2BayesX might help). > > > > Most of these issues are far more complex than students or their > > supervisors understand, do ask an experienced applied statistician > > (I'm not an applied statistician, don't rely on my suggestions, you > > yourself are responsible for your results). > > > > Hope this clarifies, > > > > Roger > > > >> > >> Kind Regards- > >> > >> Thanks- > >> > >> Burcu > >> > >> --- > >> > >> Burcu H. Ozuduru, > >> Gazi University, Faculty of Architecture, Department of City and > >> Regional Planning Maltepe-Ankara Turkey > >> t: +90-312-5823701 > >> e: [hidden email] <mailto:[hidden email]> > >> http://websitem.gazi.edu.tr/bozuduru > >> > >> > >> On Mon, Aug 29, 2016 at 12:42 AM, Roger Bivand > <[hidden email] <mailto:[hidden email]> > wrote: > >> > >>> On Sun, 28 Aug 2016, [hidden email] <mailto:[hidden email]> wrote: > >>> > >>> Dear colleagues, > >>>> > >>>> Does anyone know how to implement Spatial Lag of X (SLX) model for > >>>> panel data using splm package? I tried the following to create a WX > >>>> matrix and then apply Fixed and Random Effects: > >>>> > >>>> tb <- read.csv2("panel.csv", header = TRUE, sep = ";", dec = ".") > >>>> time <- length(unique(tb$year)) X <- as.matrix(tb[ ,3:20]) > >>>> > >>> > >>> class(X) is not "pseries", is it? > >>> > >>> WX <- slag(X, listw, maxlag = 1) > >>>> > >>>> I see that slag is applicable for vectors of pseries, but this > >>>> led to the error message: "Error in UseMethod("slag"): no applicable > method > >>>> for 'slag' applied to an object of class "c('matrix', 'double', 'numeric')" > >>>> > >>> > >>> slag.default() takes a vector, slag.pseries() a "pseries" object. > >>> > >>> > >>>> Is there any way to implement this SLX model for panel series? > >>>> > >>> > >>> Did you look at spdep::create_WX? > >>> > >>> > >>>> Regards, > >>>> > >>>> Daniel. > >>>> > >>>> [[alternative HTML version deleted]] > >>>> > >>>> > >>> Please do not post HTML - like insects, HTML carries infectious > >>> and wastes bandwidth and server processing capacity. > >>> > >>> _______________________________________________ > >>>> R-sig-Geo mailing list > >>>> [hidden email] <mailto:[hidden email]> > >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >>>> > >>>> > >>> -- > >>> Roger Bivand > >>> Department of Economics, Norwegian School of Economics, > >>> Helleveien 30, N-5045 Bergen, Norway. > >>> voice: +47 55 95 93 55 <tel:%2B47%2055%2095%2093%2055> ; fax +47 55 95 91 00 <tel:%2B47%2055%2095%2091%2000> > >>> e-mail: [hidden email] <mailto:[hidden email]> > >>> http://orcid.org/0000-0003-2392-6140 > >>> https://scholar.google.no/citations?user=AWeghB0AAAAJ <https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en> &hl=en > >>> http://depsy.org/person/434412 > >>> > >>> _______________________________________________ > >>> R-sig-Geo mailing list > >>> [hidden email] <mailto:[hidden email]> > >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >>> > >> > > > > > > -- > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N-5045 Bergen, Norway. > voice: +47 55 95 93 55 <tel:%2B47%2055%2095%2093%2055> ; fax +47 55 95 91 00 <tel:%2B47%2055%2095%2091%2000> > e-mail: [hidden email] <mailto:[hidden email]> > http://orcid.org/0000-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0AAAAJ <https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en> &hl=en > http://depsy.org/person/434412 > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] <mailto:[hidden email]> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Free forum by Nabble | Edit this page |