lagsarlm with simulated data
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=FALSE) # 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=FALSE) 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=1000, 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=data, 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