[R-sig-Geo] lagsarlm with simulated data

[Malcolm Fairbrother]

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@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
>