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:

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

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

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,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.

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

+ 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)

+ }

_______________________________________________

R-sig-Geo mailing list

R-sig-Geo at stat.math.ethz.ch

https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--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, nrow00, ncol=5)

for (i in 1:1000) {

+ e <- rnorm(n, mean=0, sd=4)Bs <- c(2, 5, 3, -2) # vector of Beta coefficients

rho <- 0.2 # set autocorrelation coefficient

bres <- matrix(NA, nrow00, ncol=5)

for (i in 1:1000) {

+ 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[1] 0.1876292 2.6275783 4.9897827 3.0060831 -1.9883803

apply(bres, 2, median)

_______________________________________________

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