Thank you, Achim.
So, then a follow up question to the community: is there a package out
there that allows one to calculate a correct DW for a regression with
weights?
Dimitri
On Fri, Aug 12, 2011 at 2:42 PM, Achim Zeileis <achim.zeil...@uibk.ac.at> wrote:
> 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.
>>
>
--
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.
