FAQ
Dear Roger (and any other interested parties),

Thanks very much for responding. I tried adding the scale() argument you suggested, but that didn't seem to make any difference. I've set a seed, and to make it even more easily reproducible, I've loosely followed code from: http://www.mail-archive.com/r-sig-geo at stat.math.ethz.ch/msg00799.html. The contiguities (in addition to the data) are now generated directly using the code below.

I'm still getting consistent upward bias for the Intercept, and otherwise perfect recoveries of the data-generating parameters. I tried reversing the sign of rho, and that didn't make any difference.

Any ideas? Sorry to bother you with this, but I'd like to know why the simulation is generating this result.

Thanks again,
Malcolm

library(spdep)
sims <- 1000 # set number of simulations
rho <- 0.2 # set autocorrelation coefficient
Bs <- c(2, 5, 3, -2) # set a vector of betas
nb7rt <- cell2nb(7, 7, torus=TRUE) # generate contiguities
listw <- nb2listw(nb7rt)
mat <- nb2mat(nb7rt) # create contiguity matrix
set.seed(20100817)
e <- matrix(rnorm(sims*length(nb7rt)), nrow=length(nb7rt)) # create random errors
e <- scale(e, center=TRUE, scaleúLSE) # constrain mean of errors to zero
X0 <- matrix(1, ncol=sims, nrow=length(nb7rt)) # create Intercept
X1 <- matrix(rnorm(sims*length(nb7rt), 4, 2), nrow=length(nb7rt)) # generate some covariates
X2 <- matrix(rnorm(sims*length(nb7rt), 2, 2), nrow=length(nb7rt))
X3 <- matrix(rnorm(sims*length(nb7rt), -3, 1), nrow=length(nb7rt))
Xbe <- X0*Bs[1]+X1*Bs[2]+X2*Bs[3]+X3*Bs[4]+e
y <- solve(diag(length(nb7rt)) - rho*mat) %*% Xbe # generate lagged y data
lag_res1 <- lapply(1:sims, function(i) lagsarlm(y[,i] ~ X1[,i] + X2[,i] + X3[,i], listw=listw)) # fit the model
apply(do.call("rbind", lapply(lag_res1, coefficients)), 2, mean) # mean estimates are excellent, except Int
apply(do.call("rbind", lapply(lag_res1, coefficients)), 2, var) # variance is also a lot higher for Int

On 16 Aug 2010, at 19:53, Roger Bivand wrote:
On Mon, 16 Aug 2010, Malcolm Fairbrother wrote:

Dear list,

I am running some simulations, trying to use lagsarlm (from the spdep package) to recover the parameters used to generate the data. In a basic simulation I am running, I am finding that I am able to recover rho almost perfectly, and all but one of the betas perfectly. However, the beta attached to the constant is substantially biased, for reasons I cannot understand.

The code I am using is below. The spatial weights matrix is from the 48 contiguous U.S. states (rows sum to 1). Can anyone see where I am going wrong? Or is the biased B0 coefficient somehow a consequence of the particular neighbourhood structure I'm using? Any help or tips would be much appreciated.
Malcolm,

You didn't set a seed, so I can't reproduce this exactly, but I think that the mean of e will be added to your constant, won't it?

n <- 48
e <- matrix(rnorm(n*1000, mean=0, sd=4), 1000, n)
summary(apply(e, 1, mean))

Could you do a scale(e, center=TRUE, scaleúLSE) to force it to mean zero?

Hope this helps,

Roger
- Malcolm

n <- dim(W)[1] # sample size
Bs <- c(2, 5, 3, -2) # vector of Beta coefficients
rho <- 0.2 # set autocorrelation coefficient
bres <- matrix(NA, nrow00, ncol=5)
for (i in 1:1000) {
+ e <- rnorm(n, mean=0, sd=4)
+ X1 <- rnorm(n, 4, 2) # create some independent variables
+ X2 <- rnorm(n, 2, 2)
+ X3 <- rnorm(n, -3, 1)
+ X <- cbind(rep(1, n), X1, X2, X3)
+ y <- (solve(diag(n)-rho*W)) %*% ((X%*%Bs)+e) # generate lagged Ys
+ data <- as.data.frame(cbind(y, X))
+ lagmod <- lagsarlm(y ~ X1 + X2 + X3, dataÚta, oxw.listw1)
+ bres[i,] <- coefficients(lagmod)
+ }
apply(bres, 2, mean)
[1] 0.1876292 2.6275783 4.9897827 3.0060831 -1.9883803
apply(bres, 2, median)
[1] 0.1907057 2.4000895 4.9887496 3.0043258 -2.0076960

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

## Related Discussions

 view thread | post posts ‹ prev | 5 of 6 | next ›
Discussion Overview
 group r-sig-geo categories r posted Aug 15, '10 at 11:19p active Aug 17, '10 at 5:54p posts 6 users 2 website r-project.org irc #r

### 2 users in discussion

Content

People

Support

Translate

site design / logo © 2022 Grokbase