[R-sig-Geo] Moran's I with Objects of Different Length

[Roger Bivand]

On Tue, 11 Aug 2020, Chanda Chiseni wrote:

> I am trying to perform moran's I test on residuals from by probit model
> using k-nearest neighbour weights, however i run into an error which i
> cant seem to find the solution anywhere online the error is
>
>> moran.test(residuals.glm(svyprobitest),knear2weight)
> Error in moran.test(residuals.glm(svyprobitest), knear2weight) :
>  objects of different length
>
>
> The codes I use to come up to this are below
>
>> library(foreign)
>> dhsanalysis2= read.dta("DHSSpatialreg11.dta")
> ##removing NAs in PSU
>> dhsanalysis3=subset(dhsanalysis2, !is.na(psu))
> ##generating a survey design
>> mydesighdhs= svydesign(ids =~psu, data =dhsanalysis3, weight = ~wgt,
> strata = ~v023,nest=TRUE)
> ##Defining coordinate colums
>> coordinates(dhsanalysis3)= c("longnum","latnum")
> #Defining Projection
>> proj4string(dhsanalysis3) <- CRS("+init=epsg:4326")
> ##binding longiude and latitude
>> lon2<- dhsanalysis3$longnum
>> lat2<- dhsanalysis3$latnum
>> coords<- cbind(lon2,lat2)
> ##Creating spatial weights based on the nearest neighbour
>> knear2= knearneigh(coords,k=2,longlat=T)
> Warning message:
> In knearneigh(coords, k = 2, longlat = T) :
>  knearneigh: identical points found ### I get a warning message on
> identical points found, is this a problem and how would i deal with this

It usually indicates muddled thinking. You have more than one observation 
at the same point, which could be an unintended duplicate (copy of the 
same observation), or it could indicate that any spatial process has 
multiple values at that point. In geostatistics, one would jitter the 
points to introduce distance. In your case, this probably isn't households 
living at one address point. All the observations at the same point will 
get the same sets of neighbours. If there are many co-located 
observations, k in k-nearest may be too small to include them all.

>> knear2.nb= knn2nb(knear2)
>> knear2weight= nb2listw(knear2.nb, style="W",zero.policy=T)
>
> ##declaring categorical variables as factor variables
>> dhsanalysis3$hivpositive.f <- factor(dhsanalysis3$hivpositive)
>> dhsanalysis3$protestant.f <- factor(dhsanalysis3$protestant2)
>> dhsanalysis3$Married.f <- factor(dhsanalysis3$Married)
>> dhsanalysis3$female.f <- factor(dhsanalysis3$female)
>> dhsanalysis3$urban1.f <- factor(dhsanalysis3$urban1)
>> dhsanalysis3$river10kmdum.f <- factor(dhsanalysis3$river10kmdum)
>> dhsanalysis3$Explorer50kmdum.f <- factor(dhsanalysis3$Explorer50kmdum)
>> dhsanalysis3$Rail50kmdum.f <- factor(dhsanalysis3$Rail50kmdum)
>> dhsanalysis3$Province1.f <- factor(dhsanalysis3$Province1)
>> dhsanalysis3$WealthIndex.f <- factor(dhsanalysis3$WealthIndex)
>> dhsanalysis3$occupation2.f <- factor(dhsanalysis3$occupation2)
>> dhsanalysis3$highested.f <- factor(dhsanalysis3$highested)
>> protestant.f= dhsanalysis3$protestant.f
>> hivpositive.f=dhsanalysi3$hivpositive.f
>> Married.f=dhsanalysis3$Married.f
>> female.f=dhsanalysis3$female.f
>> urban1.f=dhsanalysis3$urban1.f
>> river10kmdum.f=dhsanalysis3$river10kmdum.f
>> Explorer50kmdum.f=dhsanalysis2$Explorer50kmdum.f
>> Province1.f=dhsanalysis3$Province1.f
>> WealthIndex.f=dhsanalysis3$WealthIndex.f
>> occupation2.f=dhsanalysis3$occupation2.f
>> highested.f=dhsanalysis3$highested.f
>> Age=dhsanalysis3$Age
>> Age2=dhsanalysis3$Age2
>> HIVKnowledge=dhsanalysis3$HIVKnowledge
>> churchkm=dhsanalysis3$churchkm
>> lnhospitalkm=dhsanalysis3$lnhospitalkm
>> lnElevationMean=dhsanalysis3$lnElevationMean
>
>> ###VARIABLES IN SURVEY DESIGN
>> mydesighdhs$hivpositive.f <- factor(dhsanalysis3$hivpositive)
>> mydesighdhs$protestant.f <- factor(dhsanalysis3$protestant2)
>> mydesighdhs$Married.f <- factor(dhsanalysis3$Married)
>> mydesighdhs$female.f <- factor(dhsanalysis3$female)
>> mydesighdhs$urban1.f <- factor(dhsanalysis3$urban1)
>> mydesighdhs$river10kmdum.f <- factor(dhsanalysis3$river10kmdum)
>> mydesighdhs$Explorer50kmdum.f <- factor(dhsanalysis3$Explorer50kmdum)
>> mydesighdhs$Rail50kmdum.f <- factor(dhsanalysis3$Rail50kmdum)
>> mydesighdhs$Province1.f <- factor(dhsanalysis3$Province1)
>> mydesighdhs$WealthIndex.f <- factor(dhsanalysis3$WealthIndex)
>> mydesighdhs$occupation2.f <- factor(dhsanalysis3$occupation2)
>> mydesighdhs$highested.f <- factor(dhsanalysis3$highested)
>> mydesighdhs$Age=Age
>> mydesighdhs$Age2=Age2
>> mydesighdhs$HIVKnowledge=HIVKnowledge
>> mydesighdhs$churchkm=churchkm
>> mydesighdhs$lnhospitalkm=lnhospitalkm
>> mydesighdhs$lnElevationMean=lnElevationMean
>
>> svyprobitest= svyglm(hivpositive.f~ churchkm +lnhospitalkm +protestant.f+
> Age+
> Age2+Married.f+female.f+urban1.f+river10kmdum.f+Explorer50kmdum.f+Rail50kmdum.f+lnElevationMean+Province1.f
> + WealthIndex.f + HIVKnowledge+ occupation2.f+
> highested.f,design=mydesighdhs,family = quasibinomial(link =
> "probit"),data=dhsanalysis3)
> Error in .subset2(x, i, exact = exact) : subscript out of bounds

So you've an error anyway.

>
> ## Iran my probit model without implementing the survey design in the model
> just to see whether Moran test is working
>> svyprobitest2= glm(hivpositive.f~ churchkm +lnhospitalkm +protestant.f+
> Age+
> Age2+Married.f+female.f+urban1.f+river10kmdum.f+Explorer50kmdum.f+Rail50kmdum.f+lnElevationMean+Province1.f
> + WealthIndex.f + HIVKnowledge+ occupation2.f+ highested.f,family =
> quasibinomial(link = "probit"),data=dhsanalysis3)
> #Moran test
>> moran.test(residuals.glm(svyprobitest),knear2weight)
> Error in moran.test(residuals.glm(svyprobitest), knear2weight) :
>  objects of different length
>

Testing residuals using moran.test() is usually wrong, as the expectation 
and variance of the statistic are based on a null model (intercept only), 
not a model with covariates.

What are length(residuals.glm(svyprobitest)) and length(knear2.nb)? Did 
glm drop observations with missing values? If so, you should subset both 
the data submitted to glm() and the neighbour object so that they match.

Hope this helps,

Roger

>
> Kind Regards,
>
> Michael Chanda Chiseni
>
> Phd Candidate
>
> Department of Economic History
>
> Lund University
>
> Visiting address: Alfa 1, Scheelevägen 15 B, 22363 Lund
>
>
>
> *Africa is not poor, it is poorly managed (Ellen Johnson-Sirleaf ). *
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en