# NOT RUN {
data(boston)
hr0 <- lm(log(MEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) +
AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT), data=boston.c)
summary(hr0)
logLik(hr0)
gp0 <- lm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) +
AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT), data=boston.c)
summary(gp0)
logLik(gp0)
lm.morantest(hr0, nb2listw(boston.soi))
# }
# NOT RUN {
require(maptools)
boston.tr <- readShapePoly(system.file("etc/shapes/boston_tracts.shp",
package="spdep")[1], ID="poltract",
proj4string=CRS(paste("+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66",
"+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")))
boston_nb <- poly2nb(boston.tr)
# }
# NOT RUN {
gp1 <- errorsarlm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2)
+ I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT),
data=boston.c, nb2listw(boston.soi), method="Matrix",
control=list(tol.opt = .Machine$double.eps^(1/4)))
summary(gp1)
gp2 <- lagsarlm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2)
+ AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT),
data=boston.c, nb2listw(boston.soi), method="Matrix")
summary(gp2)
# }
# NOT RUN {
## Conversion table 1980/1970
# ICPSR_07913.zip
# 07913-0001-Data.txt
# http://dx.doi.org/10.3886/ICPSR07913.v1
# Provider: ICPSR
# Content: text/plain; charset="us-ascii"
#
# TY - DATA
# T1 - Census of Population and Housing 1980 [United States]:
# 1970-Pre 1980 Tract Relationships
# AU - United States Department of Commerce. Bureau of the Census
# DO - 10.3886/ICPSR07913.v1
# PY - 1984-06-28
# UR - http://dx.doi.org/10.3886/ICPSR07913.v1
# PB - Inter-university Consortium for Political and Social Research
# (ICPSR) [distributor]
# ER -
# widths <- c(ID=5L, FIPS70State=2L, FIPS70cty=3L, Tract70=6L, FIPS80State=2L,
# FIPS80cty=3L, f1=7L, CTC=6L, f2=2L, intersect1=3L, intersect2=3L, name=30L)
# dta0 <- read.fwf("07913-0001-Data.txt", unname(widths),
# col.names=names(widths), colClasses=rep("character", 12), as.is=TRUE)
# sub <- grep("25", dta0$FIPS80State)
# MA <- dta0[sub,]
## match against boston data set
# library(spdep)
# data(boston)
# bTR <- boston.c$TRACT
# x1 <- match(as.integer(MA$Tract70), bTR)
# BOSTON <- MA[!is.na(x1),]
## MA 1990 tracts
# library(rgdal)
# MAtr90 <- readOGR(".", "tr25_d90")
## counties in the BOSTON SMSA
## https://www.census.gov/population/metro/files/lists/historical/90nfips.txt
## 1123 Boston-Lawrence-Salem-Lowell-Brockton, MA NECMA
## 1123 25 009 Essex County
## 1123 25 017 Middlesex County
## 1123 25 021 Norfolk County
## 1123 25 023 Plymouth County
## 1123 25 025 Suffolk County
# BOSTON_SMSA <- MAtr90[MAtr90$CO <!-- %in% c("009", "017", "021", "023", "025"),] -->
# proj4string(BOSTON_SMSA) <- CRS(paste("+proj=longlat +datum=NAD27 +no_defs",
# "+ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat"))
# CTC4 <- substring(BOSTON$CTC, 1, 4)
# CTC4u <- unique(CTC4)
# TB_CTC4u <- match(BOSTON_SMSA$TRACTBASE, CTC4u)
## match 1980 tracts with 1990
# BOSTON_SMSA1 <- BOSTON_SMSA[!is.na(TB_CTC4u),]
## union Polygons objects with same 1970 tract code
#library(rgeos)
# BOSTON_SMSA2 <- gUnaryUnion(BOSTON_SMSA1,
# id=as.character(BOSTON_SMSA1$TRACTBASE))
## reorder data set
# mm <- match(as.integer(as.character(row.names(BOSTON_SMSA2))), boston.c$TRACT)
# df <- boston.c[mm,]
# row.names(df) <- df$TRACT
# row.names(BOSTON_SMSA2) <- as.character(as.integer(row.names(BOSTON_SMSA2)))
## create SpatialPolygonsDataFrame
# BOSTON_SMSA3 <- SpatialPolygonsDataFrame(BOSTON_SMSA2,
# data=data.frame(poltract=row.names(BOSTON_SMSA2),
# row.names=row.names(BOSTON_SMSA2)))
# BOSTON_SMSA4 <- spCbind(BOSTON_SMSA3, df)
# mm1 <- match(boston.c$TRACT, row.names(BOSTON_SMSA4))
# BOSTON_SMSA5 <- BOSTON_SMSA4[mm1,]
#writeOGR(BOSTON_SMSA5, ".", "boston_tracts", driver="ESRI Shapefile",
# overwrite_layer=TRUE)
# moran.test(boston.c$CMEDV, nb2listw(boston.soi))
# moran.test(BOSTON_SMSA5$CMEDV, nb2listw(boston.soi))
# }
Run the code above in your browser using DataLab