Read in data
cars = readOGR("usa48_usedcars.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/jessieshapiro/Documents/Documents/Education/University of Utah/GEOG6000/Labs/lab12/usa48_usedcars.shp", layer: "usa48_usedcars"
## with 48 features
## It has 3 fields
Visualize data
coords = coordinates(cars)
plot(cars, axes = T)
points(coords, pch=16, col=2)

Build a nieghborhood function

The Queen’s case is used to build a neighborhood function linking the 48 states. The Rooke’s case, Delaunay triangulation and K nearest neighbor were also considered to build a neighborhood matrix but either had too many links far away from eachother or had too few linkages.

carsqueen = poly2nb(cars)
plot(cars)
plot(carsqueen, coords, add=T, col=2, lwd=3)

Build a spatial weight matrix
dstscars = nbdists(carsqueen, coords)
invd = function(x) {1/(x/1000)}
idwcars = lapply(dstscars, invd)
swmcars = nb2listw(carsqueen, glist = idwcars, style = "B")
summary(swmcars)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 48 
## Number of nonzero links: 214 
## Percentage nonzero weights: 9.288194 
## Average number of links: 4.458333 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8 
##  1  4  9 11 10  9  2  2 
## 1 least connected region:
## 16 with 1 link
## 2 most connected regions:
## 22 39 with 8 links
## 
## Weights style: B 
## Weights constants summary:
##    n   nn       S0       S1        S2
## B 48 2304 59814.62 50413392 369595519
Test for spatial autocorrelation
cars.listw = nb2listw(carsqueen)
moran.test(cars$price_1960, cars.listw, alternative="two.sided", randomisation=TRUE)
## 
##  Moran I test under randomisation
## 
## data:  cars$price_1960  
## weights: cars.listw    
## 
## Moran I statistic standard deviate = 8.1752, p-value = 2.954e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.783561543      -0.021276596       0.009692214
moran.mc(cars$price_1960, cars.listw, nsim=999, alternative='greater')
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  cars$price_1960 
## weights: cars.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.78356, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Based on the p-value of 2.954-e16, Z-score of 8.1752 and Moran’s I statistic of 0.78356, we can reject the null hypothesis that there is no significant spatial autocorrelation.

Build simple linear model
lmcars = lm(cars$price_1960 ~ cars$tax_charge)
summary(lmcars)
## 
## Call:
## lm(formula = cars$price_1960 ~ cars$tax_charge)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -116.701  -45.053   -1.461   43.400  107.807 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1435.7506    27.5796  52.058  < 2e-16 ***
## cars$tax_charge    0.6872     0.1754   3.918 0.000294 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57.01 on 46 degrees of freedom
## Multiple R-squared:  0.2503, Adjusted R-squared:  0.234 
## F-statistic: 15.35 on 1 and 46 DF,  p-value: 0.0002939
moran.mc(residuals(lmcars), cars.listw, nsim=999)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  residuals(lmcars) 
## weights: cars.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.57482, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

Coefficients of the model results: 1435.7506 (Intercept), 0.6872(Tax Charge)

R2: 0.2503

The p-value of 0.001 and statistic of 0.57482 from the Monte-Carlo simulation indicates that spatial autocorrelation is present in the residuals.

Lagrange multiplier test
summary(lm.LMtests(lmcars, cars.listw, test=c("LMerr","LMlag")))
##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = cars$price_1960 ~ cars$tax_charge)
## weights: cars.listw
##  
##       statistic parameter   p.value    
## LMerr    31.793         1 1.715e-08 ***
## LMlag    40.664         1 1.808e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm.LMtests(lmcars, cars.listw, test=c("RLMerr","RLMlag")))
##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = cars$price_1960 ~ cars$tax_charge)
## weights: cars.listw
##  
##        statistic parameter  p.value   
## RLMerr  0.051751         1 0.820044   
## RLMlag  8.922943         1 0.002816 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The non-robust version had statistically significant results for both the error (p-value 1.715e-08) and lag model (p-value 1.808e-10). The robust version only showed significant results for the lag model (p-value 0.002816). The model of choice is the lag model.

Build spatial model
lmcarsLAG = lagsarlm(price_1960 ~ tax_charge, data=cars, cars.listw)
summary(lmcarsLAG)
## 
## Call:
## lagsarlm(formula = price_1960 ~ tax_charge, data = cars, listw = cars.listw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -77.6781 -16.9504   4.2497  19.5484  58.9811 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 309.42546  123.03167  2.5150   0.0119
## tax_charge    0.16711    0.10212  1.6364   0.1018
## 
## Rho: 0.78302, LR test value: 42.681, p-value: 6.4426e-11
## Asymptotic standard error: 0.081636
##     z-value: 9.5916, p-value: < 2.22e-16
## Wald statistic: 91.998, p-value: < 2.22e-16
## 
## Log likelihood: -239.8252 for lag model
## ML residual variance (sigma squared): 1036.7, (sigma: 32.197)
## Number of observations: 48 
## Number of parameters estimated: 4 
## AIC: 487.65, (AIC for lm: 528.33)
## LM test for residual autocorrelation
## test value: 2.1141, p-value: 0.14594

Coeffecients from lag model: 309.42546(Intercept), 0.16711(Tax Charge)

Rho: 0.78302

AIC: 478.65

AIC for lm: 528.33

Test residuals for remaining autocorrelation
moran.mc(residuals(lmcarsLAG), listw=cars.listw, nsim=999)
## Warning: Method residuals.sarlm moved to the spatialreg package
## Warning in residuals.sarlm(lmcarsLAG): install the spatialreg package
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  residuals(lmcarsLAG) 
## weights: cars.listw  
## number of simulations + 1: 1000 
## 
## statistic = -0.077924, observed rank = 299, p-value = 0.701
## alternative hypothesis: greater

The results of Moran’s I statistic for residuals of -0.0779 and p-value of 0.706, indicates the model accounts for spatial autocorrelation and there is no significant autocorrelation with the remaining in the residuals. The coefficient relating tax and delivery to car price is not that significant. This might be because supply and demand differ between states.