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@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.biv...@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
