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@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.biv...@nhh.no
>
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
