On Fri, 12 Aug 2011, Dimitri Liakhovitski wrote:

Hello!

I have a data frame mysample (sorry for a long way of creating it
below - but I need it in this form, and it works). I regress Y onto X1
through X11 - first without weights, then with weights:

regtest1<-lm(Y~., data=mysample[-13]))
regtest2<-lm(Y~., data=mysample[-13]),weights=mysample$weight)
summary(regtest1)
summary(regtest2)

Then I calculate Durbin-Watson for both regressions using 2 different packages:

library(car)
library(lmtest)

durbinWatsonTest(regtest1)[2]
dwtest(regtest1)$stat

durbinWatsonTest(regtest2)[2]
dwtest(regtest2)$stat

When there are no weights, the Durbin-Watson statistic is the same.
But when there are weights, 2 packages give Durbin-Watson different
statistics. Anyone knows why?

The result of dwtest() is wrong. Internally, dwtest() extracts the model matrix and response (but no weights) and does all processing based on these. Thus, it computes the DW statistic for regtest1 not regtest2.

I've just added a patch to my source code which catches this problem and throws a meaningful error message. It will be part of the next release (0.9-29) in due course.

Of course, this doesn't help you with computing the DW statistic for the weighted regression but hopefully it reduces the confusion about the different behaviors...
Z

Also, it's interesting that both of them are also different from what
SPSS spits out...

Thank you!
Dimitri


############################################
### Run the whole code below to create mysample:

intercor<-0.3        # intercorrelation among all predictors
k<-10                # number of predictors
sigma<-matrix(intercor,nrow=k,ncol=k) # matrix of intercorrelations
among predictors
diag(sigma)<-1

require(mvtnorm)
set.seed(123)
mypop<-as.data.frame(rmvnorm(n=100000, mean=rep(0,k), sigma=sigma,
method="chol"))
names(mypop)<-paste("x",1:k,sep="")
set.seed(123)
mypop$x11<-sample(c(0,1),100000,replace=T)

set.seed(123)
betas<-round(abs(rnorm(k+1)),2) # desired betas
Y<-as.matrix(mypop) %*% betas
mypop<-cbind(mypop, Y)
rSQR<-.5
VARofY<- mean(apply(as.data.frame(mypop$Y),2,function(x){x^2})) -
mean(mypop$Y)^2
mypop$Y<-mypop$Y + rnorm(100000, mean=0, sd=sqrt(VARofY/rSQR-VARofY))

n<-200
set.seed(123)
cases.for.sample<-sample(100000,n,replace=F)
mysample<-mypop[cases.for.sample,]
mysample<-cbind(mysample[k+2],mysample[1:(k+1)])  #dim(sample)
weight<-rep(1:10,20);weight<-weight[order(weight)]
mysample$weight<-weight



--
Dimitri Liakhovitski
marketfusionanalytics.com

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to