For this homework, I used the 2016 Ethiopian Demographic and Health Survey(EDHS). The spatial data for Ethipia is downloaded from the Spatial Data Repository, DHS website (https://spatialdata.dhsprogram.com/home/).

The purpose of these homework to examine the association between (avg height for age standard deviation) which is a continous measure and other explanatory variables (proportion of poor mothers, proportion of non educated mothers, average children evern born and mean age of mothers), using various spatial Autoregression models to determine the best model.

library(maptools)
library(sf)

library(reticulate)
library(sp)
library(spdep)
library(RQGIS) #you need Qgis 2.18 installed, not 3!
library(mapview)
library(rgdal)
library(ggplot2)
library(dplyr)
library(gridExtra)
library(maps)
library(tibble)
library(knitr)
library(viridis)
library(rvest)
library(ggsn)
library(ggmap)
library(RColorBrewer)
library(haven)
library(survey)
library(dplyr)
library(foreign)
library(lme4)
library(car)
library(dplyr)
library(nortest)
library(lmtest)
library(classInt)
library(sandwich)
eth_16_points<-st_read("~/Google Drive/DHS Resources/ GPS data from DHS/ ET_2016_DHS_ETGE71FL",layer = "ETGE71FL")
## Reading layer `ETGE71FL' from data source `/Users/fikrewoldbitew/Google Drive/DHS Resources/ GPS data from DHS/ ET_2016_DHS_ETGE71FL' using driver `ESRI Shapefile'
## Simple feature collection with 643 features and 20 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5.684342e-14 ymin: 5.684342e-14 xmax: 47.00792 ymax: 14.50213
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
eths_16_points<-st_transform(eth_16_points, crs=20136 )
 
eths_border<-st_read("~/Google Drive/fikre/ETH_adm_DIVA_GIS", layer= "ETH_adm0")
## Reading layer `ETH_adm0' from data source `/Users/fikrewoldbitew/Google Drive/fikre/ETH_adm_DIVA_GIS' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 32.99994 ymin: 3.322099 xmax: 47.98618 ymax: 14.89996
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
eths_border<-st_transform(eths_border, crs=20136 )

eths_16_points<-st_crop(eths_16_points,eths_border)
 
 
#eths_von_16<-st_voronoi(st_union(eths_16_points))
#eth_vonmap_16<-st_intersection(st_cast(eths_von_16), st_union(eths_border))
 
 
#eth_vonmap_16<-st_intersection(eths_border,eths_von_16)
 
#eth_vonmap_16<-st_collection_extract(eth_vonmap_16,type="POLYGON")
 

#pv<-st_sf(eths_16_points, eth_vonmap_16)
library(SDraw)
library(rgeos)
out<-voronoi.polygons(as(eths_16_points, "Spatial"), bounding.polygon = as(eths_border, "Spatial"))
out@data<-as.data.frame(eths_16_points)
clip<- raster::intersect(out, eths_border)

eth_vonmap_16<-st_as_sf(clip)
plot(eth_vonmap_16["DHSCLUST"])

If we have the 2005 DHS data, we could generate estimates of some quality, say proportion of women underweight in Ethiopia

dhs_2016<- read.dta( "~/Google Drive/DHS Resources/ET_2016_DHS_08152017_2025_103714/etcr70dt/ETCR70FL.DTA", convert.factors = TRUE)

#dhs_2016<-read.dta("~/Google Drive/DHS Resources/ET_2016_DHS_08152017_2025_103714/etir70dt/ETIR70FL.dta", convert.factors = F)

edhs16<-subset(dhs_2016,select=c("v445","v001", "v440","v002", "v003", "v005", "v012", "v013", "v021","v022", "v024", "v025","v102", "v106", "v113", "v115","v116", "v119", "v120", "v121", "v122", "v123", "v124", "v125", "v127", "v129", "v153", "v190" , "v201","v213", "v217", "v302", "v440", "v445", "v501", "v511", "v525", "v623", "v701","v705", "v717"))
 
edhs16$SVYYEAR<-2016


#eth_16_points<-st_read("~/Google Drive/DHS Resources/ GPS data from DHS/ ET_2016_DHS_ETGE71FL",layer = "ETGE71FL")
mdat16<-merge(edhs16, eth_16_points, by.x="v021", by.y="DHSCLUST" ,all.x=T)
#tail(names(edhs11$mdat1), n=20)
mdat16<- edhs16[edhs16$v440!=9999&edhs16$v440!=9998,]

mdat16$stunt<-ifelse(mdat16$v440==9999 | mdat16$v440==9998 , NA, mdat16$v440/100)

mdat16$stunt_bin <- ifelse(mdat16$v106=='0',1,0)

#Women who are currently pregnant or pospartum/amenorrheic-Not Included  

#Education 
#mdat16$educ<-Recode(mdat16$v106, recodes="0='0Noeducation'; 1='1Primary'; 2='2Secodary'; 3='3Higher'; else=NA", as.factor=T)



mdat16$pno_educ<-Recode(mdat16$v106, recodes="'no education'=0;else=1", as.factor=T)


#Residence
mdat16$rural<-Recode(mdat16$v025, recodes="'rural'=1; else=0",  as.factor = T)


#mdat16$rur<-Recode(mdat16$v025, recodes='1=1;  else=0',  as.factor = F)
#mdat16$count<-1

#Age
mdat16$age<-Recode(mdat16$v013, recodes="1='15-19'; 2:3='20-29'; 4:5='30-39'; 6:7='40-49'; else='NA'", as.factor=T)

#Children Ever Born 
mdat16$parity<-Recode(mdat16$v201, recodes="0='0'; 1:2='1-2'; 3:4='3-4'; 5:18='5+'; else=NA", as.factor=T)

#wealth Index
mdat16$wealthindex<-Recode(mdat16$v190, recodes="1:2='1Poor'; 3='2Middle'; 4:5='3Rich'; else=NA", as.factor=T)

mdat16$ppoor<-Recode(mdat16$v190, recodes="'poorest'=1;'poorer'=1;else=0", as.factor=T)

mdat16$pwt <-mdat16$v005/1000000
mdat16$weight<-mdat16$v005/1000000
options(survey.lonely.psu = "adjust")
mdat16$yrstratum<-paste(mdat16$SVYYEAR, mdat16$v022, sep="-")

mdat16$yrpsu<-paste(mdat16$SVYYEAR, mdat16$v021, sep="-")

#table(mdat16$yrstratum)

des<-svydesign(ids=~yrpsu, strata=~yrstratum, weights=~weight, data = mdat16[is.na(mdat16$weight)==F,], nest=T)
no_educ<-svyby(~pno_educ, by=~v021, design = des, FUN = svymean, na.rm=T)

poor<-svyby(~ppoor, by=~v021, design = des, FUN = svymean, na.rm=T)
  1. Depedent Variable: Height for Age standard deviations (mean_stunt) Independent variables: mean_stunt= mean_stunt, wt = wt,mean_age, Avg childen ever born=avg_ceb, %noeduc= pno_educ, %poor = ppoor
mdat_1<-merge(mdat16, no_educ, by.x="v021", by.y="v021", sort=F)

mdat_new2<-merge(mdat_1, poor, by.x="v021", by.y="v021", sort=F)

datagroup<-mdat_new2%>%group_by(v021,pno_educ1, ppoor1) %>% 
summarize(mean_stunt = mean(stunt, na.rm = TRUE), 
          wt = mean(pwt, na.rm = TRUE),
          mean_age= mean(v012, na.rm=TRUE), 
          avg_ceb= mean(v201, na.rm=TRUE))

2.Ordinary OLS model for Mean Height for Age Standard deviations and a serious of weights using knearneigh weights(k=2 through k=6)

datagroup<-na.omit(datagroup)
fit.ols<-lm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1, data =datagroup)
datagroup$stunt_ests <- fitted(fit.ols)
datagroup$olsresid<-rstudent(fit.ols)

place_ests_16<-aggregate(stunt_ests~v021, data=datagroup, FUN=mean)

place_ests_residuals<-aggregate(olsresid~v021, data=datagroup, FUN=mean)
library(maptools)
eth_border=sf::st_read("~/Google Drive/fikre/ETH_adm_DIVA_GIS/ETH_adm1.shp")
## Reading layer `ETH_adm1' from data source `/Users/fikrewoldbitew/Google Drive/fikre/ETH_adm_DIVA_GIS/ETH_adm1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 32.99994 ymin: 3.322099 xmax: 47.98618 ymax: 14.89996
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##Map some estimates:



dhs_est_16<-left_join(eth_vonmap_16, place_ests_16, by=c("DHSCLUST"= "v021"))

dhs_est_resid<-left_join(eth_vonmap_16, place_ests_residuals, by=c("DHSCLUST"= "v021"))


#plot(dhs_est_16["stunt_ests"])
brks<-quantile(c(dhs_est_16$stunt_ests), probs = seq(0,1,length.out = 5), na.rm = T)

brks_res<-quantile(c(place_ests_residuals$olsresid), probs = seq(0,1,length.out = 5), na.rm = T)


dhs_est_16$stunt_quart_2016<-cut(dhs_est_16$stunt_ests, brks, include.lowest = T)
dhs_est_resid$res_quart<-cut(dhs_est_resid$olsresid, brks_res, include.lowest = T)

eth.geo.16<- st_as_sf(dhs_est_16)
eth.geo.resd<- st_as_sf(dhs_est_resid)

p1<-eth.geo.16%>%ggplot()+geom_sf(aes(fill=stunt_quart_2016))+geom_sf(data=eth_border,fill=NA, color="black",  size = 1)+
  scale_colour_brewer(palette = "YlOrRd",direction = -1 )+scale_fill_brewer(palette = "YlOrRd", na.value="grey", direction = -1)+guides(fill=guide_legend(title="Quartile"))+ylab(label = "")+xlab(label = "")+ggtitle(label = "Spatial Distribution of children's Avg. Height for Age SD, Ethiopia 2016", subtitle = "EDHS 2016")+theme(legend.text = element_text(size = 10))+theme(plot.title = element_text(color="Red", size=10, face="bold"), plot.subtitle = element_text(color="blue", size=10, face="bold"))+theme(legend.title = element_text(color="blue", size=10, face="bold"))


p1

#ggsave(plot=p1,filename = "~/Google Drive/PhD_Courses/Spring 2019/SpatialDemography_7263/HW/chilren_stunting.png", width = 30, height = 20,units = "in", dpi=100)




p2<-eth.geo.resd%>%ggplot()+geom_sf(aes(fill=res_quart))+geom_sf(data=eth_border,fill=NA, color="black",  size = 1)+
  scale_colour_brewer(palette = "PuBuGn", direction = -1)+scale_fill_brewer(palette = "PuBuGn", na.value="grey", direction = -1)+guides(fill=guide_legend(title="Quartile"))+ylab(label = "")+xlab(label = "")+ggtitle(label = "OLS Residuals for Avg. Children Height for Age SD, Ethiopia 2016", subtitle = "EDHS 2016")+theme(legend.text = element_text(size = 10))+theme(plot.title = element_text(color="Red", size=10, face="bold"), plot.subtitle = element_text(color="blue", size=10, face="bold"))+theme(legend.title = element_text(color="blue", size=10, face="bold"))



p2

#ggsave(plot=p1,filename = "~/Google Drive/PhD_Courses/Spring 2019/SpatialDemography_7263/HW/chilren_stunting.png", width = 30, height = 20,units = "in", dpi=100)
  1. Ordinary OLS model for Mean Height for Age Standard deviations and a serious of weights using knearneigh weights(k=2 through k=6)
#Fit  
eth_pts<-readOGR("~/Google Drive/DHS Resources/ GPS data from DHS/ ET_2016_DHS_ETGE71FL",layer = "ETGE71FL")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/fikrewoldbitew/Google Drive/DHS Resources/ GPS data from DHS/ ET_2016_DHS_ETGE71FL", layer: "ETGE71FL"
## with 643 features
## It has 20 fields
datagr_new<-datagroup[!is.na(datagroup$v021),]
eth_pts<-eth_pts[eth_pts$DHSCLUST!=149&eth_pts$DHSCLUST!=324&eth_pts$DHSCLUST!=500,]
m.geo<-merge(eth_pts@data, datagr_new, by.x="DHSCLUST", by.y="v021",  sort=F)
m.geo<- na.omit(m.geo)  
eth_pts@data<-m.geo


#mapview(eth_pts["DHSCLUST"])
fit.ols<-lm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1, data = eth_pts)
fit.ols.wt<-lm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1, data =eth_pts, weights = wt)
summary(fit.ols)
## 
## Call:
## lm(formula = mean_stunt ~ mean_age + avg_ceb + ppoor1 + pno_educ1, 
##     data = eth_pts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.03161 -0.42020 -0.06428  0.34888  2.50088 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.096678   0.256816  -4.270 2.25e-05 ***
## mean_age    -0.019366   0.008809  -2.199 0.028266 *  
## avg_ceb      0.104259   0.027561   3.783 0.000170 ***
## ppoor1       0.311025   0.095573   3.254 0.001197 ** 
## pno_educ1    0.436899   0.113805   3.839 0.000136 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6255 on 635 degrees of freedom
## Multiple R-squared:  0.06594,    Adjusted R-squared:  0.06006 
## F-statistic: 11.21 on 4 and 635 DF,  p-value: 8.608e-09
summary(fit.ols.wt)
## 
## Call:
## lm(formula = mean_stunt ~ mean_age + avg_ceb + ppoor1 + pno_educ1, 
##     data = eth_pts, weights = wt)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.55711 -0.19224  0.01128  0.22829  2.03638 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.897318   0.210111  -4.271 2.25e-05 ***
## mean_age    -0.017157   0.007176  -2.391 0.017102 *  
## avg_ceb      0.076210   0.020043   3.802 0.000157 ***
## ppoor1      -0.124867   0.073718  -1.694 0.090784 .  
## pno_educ1    0.200008   0.083614   2.392 0.017045 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4288 on 635 degrees of freedom
## Multiple R-squared:  0.02655,    Adjusted R-squared:  0.02042 
## F-statistic:  4.33 on 4 and 635 DF,  p-value: 0.001835
options(tigris_class = "sf")

#eths_pts<-as(eths_pts, "Spatial")

#nbs<-poly2nb(eths_pts, queen = T)
#wt<-nb2listw(nbs, style = "W")

et.nb6<-knearneigh(coordinates(eth_pts), k=6)
et.nb6<-knn2nb(et.nb6)
et.wt6<-nb2listw(et.nb6, style="W")

moran.test(eth_pts$mean_stunt, listw=et.wt6)
## 
##  Moran I test under randomisation
## 
## data:  eth_pts$mean_stunt  
## weights: et.wt6    
## 
## Moran I statistic standard deviate = 23.69, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.4946457718     -0.0015649452      0.0004387478
et.nb5<-knearneigh(coordinates(eth_pts), k=5)
et.nb5<-knn2nb(et.nb5)
et.wt5<-nb2listw(et.nb5, style="W")
moran.test(eth_pts$mean_stunt, listw=et.wt5)
## 
##  Moran I test under randomisation
## 
## data:  eth_pts$mean_stunt  
## weights: et.wt5    
## 
## Moran I statistic standard deviate = 21.178, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.483606739      -0.001564945       0.000524850
et.nb4<-knearneigh(coordinates(eth_pts), k=4)
et.nb4<-knn2nb(et.nb4)
et.wt4<-nb2listw(et.nb4, style="W")
moran.test(eth_pts$mean_stunt, listw=et.wt4)
## 
##  Moran I test under randomisation
## 
## data:  eth_pts$mean_stunt  
## weights: et.wt4    
## 
## Moran I statistic standard deviate = 18.874, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.481400822      -0.001564945       0.000654826
et.nb3<-knearneigh(coordinates(eth_pts), k=3)
et.nb3<-knn2nb(et.nb3)
et.wt3<-nb2listw(et.nb3,style="W")
moran.test(eth_pts$mean_stunt, listw=et.wt3)
## 
##  Moran I test under randomisation
## 
## data:  eth_pts$mean_stunt  
## weights: et.wt3    
## 
## Moran I statistic standard deviate = 15.986, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.4690059481     -0.0015649452      0.0008664981
et.nb2<-knearneigh(coordinates(eth_pts), k=2)
et.nb2<-knn2nb(et.nb2)
et.wt2<-nb2listw(et.nb2,style="W")
moran.test(eth_pts$mean_stunt, listw=et.wt2)
## 
##  Moran I test under randomisation
## 
## data:  eth_pts$mean_stunt  
## weights: et.wt2    
## 
## Moran I statistic standard deviate = 12.794, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.459402538      -0.001564945       0.001298163

The above results depicts that a positive autocorrelation in the residuals of mean height for age standard deviations is observed in every neighborhood specifications employed.

resi<-c(lm.morantest(fit.ols, listw=et.wt2)$estimate[1],
        lm.morantest(fit.ols, listw=et.wt3)$estimate[1],
        lm.morantest(fit.ols, listw=et.wt4)$estimate[1],
        lm.morantest(fit.ols, listw=et.wt5)$estimate[1],
        lm.morantest(fit.ols, listw=et.wt6)$estimate[1])
plot(resi, type="l")

b. After testing for autocorrelation in the residual using each neighborhood specification, the k=6 and k=5 has the largest redidual Moran’s I value, so we use k=6.

eth_pts$avg_stunt<-as.numeric(scale(datagroup$mean_stunt, center=T, scale=T))
library(sf)
qplot(eth_pts$avg_stunt,geom="histogram" )
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

According to WHO recommendation z-score less than -2 is most common criterion i.e 2 standard deviations (SD) below the median in reference population. , Height/Age z-score < -2 SD are considered as “stunted”

#locali<-localmoran(eth_pts$mean_stunt,  listw=et.wt6, p.adjust.method="fdr")
locali<-localmoran(datagroup$olsresid,  listw=et.wt6, p.adjust.method="fdr")

datagroup$locali<-locali[,1]
datagroup$localp<-locali[,5]

datagroup$locali <- fitted(fit.ols)

place_ests_mlocal<-aggregate(locali~v021, data=datagroup, FUN=mean)

dhs_est_mlocal<-left_join(eth_vonmap_16, place_ests_mlocal, by=c("DHSCLUST"= "v021"))

brk_lmor<-quantile(c(dhs_est_mlocal$locali), probs = seq(0,1,length.out = 5), na.rm = T)


dhs_est_mlocal$loc_quart<-cut(dhs_est_mlocal$locali, brk_lmor, include.lowest = T)

eth.geo.moran<- st_as_sf(dhs_est_mlocal)


p1<-eth.geo.moran%>%ggplot()+geom_sf(aes(fill=loc_quart))+geom_sf(data=eths_border,fill=NA, color="black",  size = 1)+
  scale_colour_brewer(palette = "RdBu" )+scale_fill_brewer(palette = "RdBu", na.value="grey")+guides(fill=guide_legend(title="Quartile"))+ylab(label = "")+xlab(label = "")+ggtitle(label = "Local Moran's of OLS Residuals", subtitle = "EDHS 2016")+theme(legend.text = element_text(size = 10))+theme(plot.title = element_text(color="Red", size=10, face="bold"), plot.subtitle = element_text(color="blue", size=10, face="bold"))+theme(legend.title = element_text(color="blue", size=10, face="bold"))

p1

  1. In the below section alternative Special Autoregressive Models (SAM) have been examined using k=6 neighborhood specification.
# Spatial Autoregressive Error Model

fit.err<-errorsarlm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1, data = eth_pts, listw=et.wt6, method = "MC")

#summary(fit.err,Nagelkerke=T)

#Spatial Autoregressive Lag Model
fit.lag<-lagsarlm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1, data = eth_pts, listw=et.wt6, type="lag", method = "MC")
#summary(fit.lag, Nagelkerke=T)


#Spatial Autoregression Model 

fit.spauto.w<-spautolm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1,data = eth_pts, listw=et.wt6, family="SAR", weights=wt, method = "MC")
#summary(fit.spauto.w, Nagelkerke=T)


#Spatial Autoregressive Lag Model 

fit.lag.d<-lagsarlm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1,data = eth_pts, listw=et.wt6, type="mixed", method = "MC")

#Spatial Autoregressive Error Model

fit.err.d<-errorsarlm(mean_stunt~mean_age+avg_ceb+ppoor1+pno_educ1,data = eth_pts, listw=et.wt6, etype = "emixed", method = "MC")

#summary(fit.err.d,Nagelkerke=T)
#As a pretty good indicator of which model is best, look at the AIC of each
AIC(fit.ols)
## [1] 1222.729
AIC(fit.ols.wt)
## [1] 1369.384
AIC(fit.lag)
## [1] 978.4153
AIC(fit.err)
## [1] 982.1506
AIC(fit.spauto.w)
## [1] 1265.569
AIC(fit.lag.d)
## [1] 968.9212
AIC(fit.err.d)
## [1] 973.9633

The above results clearly depicts that all spacial auregression models have lower AIC values when compared to the OLS models. However, the Spatial Autoregressive Lag Model was found to have a minimum AIC compared to other models. So, I chose Spatial Autoregressive Lag Model based on minimum AIC criteria.

summary(fit.lag.d, Nagelkerke=T)
## 
## Call:lagsarlm(formula = mean_stunt ~ mean_age + avg_ceb + ppoor1 + 
##     pno_educ1, data = eth_pts, listw = et.wt6, type = "mixed", 
##     method = "MC")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.635925 -0.322735 -0.020156  0.303995  1.971338 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   -0.21992287  0.43356867 -0.5072 0.611987
## mean_age       0.00037104  0.00744747  0.0498 0.960265
## avg_ceb        0.03861428  0.02385247  1.6189 0.105473
## ppoor1        -0.03331627  0.08317538 -0.4006 0.688748
## pno_educ1      0.26652908  0.10005754  2.6638 0.007727
## lag.mean_age  -0.01982427  0.01456613 -1.3610 0.173519
## lag.avg_ceb    0.01356235  0.04329551  0.3133 0.754090
## lag.ppoor1     0.33117909  0.14160577  2.3387 0.019349
## lag.pno_educ1  0.05924464  0.16832478  0.3520 0.724864
## 
## Rho: 0.63138, LR test value: 206.04, p-value: < 2.22e-16
## Asymptotic standard error: 0.037647
##     z-value: 16.771, p-value: < 2.22e-16
## Wald statistic: 281.27, p-value: < 2.22e-16
## 
## Log likelihood: -473.4606 for mixed model
## ML residual variance (sigma squared): 0.24071, (sigma: 0.49062)
## Nagelkerke pseudo-R-squared: 0.38147 
## Number of observations: 640 
## Number of parameters estimated: 11 
## AIC: 968.92, (AIC for lm: 1173)
## LM test for residual autocorrelation
## test value: 20.573, p-value: 5.7398e-06

The results above Spatial Autoregressive Lag Model therefore clearly depicts that: - The percentage of poor and percentage of non-educated mothers are positvely significantly associated with stunting (avg height for age standard deviations). The increament of proportion of poor and non educated mothers are significant contributing factors for the expansion of stunting in Ethiopia.