County Analysis
Doc Larry
as of 8/31/2020
Load Libraries
## Loading required package: carData## Loading required package: sp## Loading required package: raster##
## Attaching package: 'raster'## The following objects are masked from 'package:MASS':
##
## area, select## Warning: package 'leaflet' was built under R version 4.0.2## Checking rgeos availability: TRUE## Warning: package 'rgdal' was built under R version 4.0.2## rgdal: version: 1.5-10, (SVN revision 1006)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/lfult/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/lfult/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.## rgeos version: 0.5-3, (SVN revision 634)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-2
## Polygon checking: TRUE##
## Attaching package: 'shiny'## The following object is masked from 'package:ggExtra':
##
## runExample## Warning: package 'sf' was built under R version 4.0.2## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: Matrix## Registered S3 methods overwritten by 'spatialreg':
## method from
## residuals.stsls spdep
## deviance.stsls spdep
## coef.stsls spdep
## print.stsls spdep
## summary.stsls spdep
## print.summary.stsls spdep
## residuals.gmsar spdep
## deviance.gmsar spdep
## coef.gmsar spdep
## fitted.gmsar spdep
## print.gmsar spdep
## summary.gmsar spdep
## print.summary.gmsar spdep
## print.lagmess spdep
## summary.lagmess spdep
## print.summary.lagmess spdep
## residuals.lagmess spdep
## deviance.lagmess spdep
## coef.lagmess spdep
## fitted.lagmess spdep
## logLik.lagmess spdep
## fitted.SFResult spdep
## print.SFResult spdep
## fitted.ME_res spdep
## print.ME_res spdep
## print.lagImpact spdep
## plot.lagImpact spdep
## summary.lagImpact spdep
## HPDinterval.lagImpact spdep
## print.summary.lagImpact spdep
## print.sarlm spdep
## summary.sarlm spdep
## residuals.sarlm spdep
## deviance.sarlm spdep
## coef.sarlm spdep
## vcov.sarlm spdep
## fitted.sarlm spdep
## logLik.sarlm spdep
## anova.sarlm spdep
## predict.sarlm spdep
## print.summary.sarlm spdep
## print.sarlm.pred spdep
## as.data.frame.sarlm.pred spdep
## residuals.spautolm spdep
## deviance.spautolm spdep
## coef.spautolm spdep
## fitted.spautolm spdep
## print.spautolm spdep
## summary.spautolm spdep
## logLik.spautolm spdep
## print.summary.spautolm spdep
## print.WXImpact spdep
## summary.WXImpact spdep
## print.summary.WXImpact spdep
## predict.SLX spdep##
## Attaching package: 'spdep'## The following objects are masked from 'package:spatialreg':
##
## anova.sarlm, as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, bptest.sarlm, can.be.simmed,
## cheb_setup, coef.gmsar, coef.sarlm, coef.spautolm, coef.stsls,
## create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm,
## deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
## errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm,
## fitted.SFResult, fitted.spautolm, get.ClusterOption,
## get.coresOption, get.mcOption, get.VerboseOption,
## get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone,
## gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts,
## Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
## lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm,
## LR1.sarlm, LR1.spautolm, LU_prepermutate_setup, LU_setup,
## Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
## mom_calc_int2, moments_setup, powerWeights, predict.sarlm,
## predict.SLX, print.gmsar, print.ME_res, print.sarlm,
## print.sarlm.pred, print.SFResult, print.spautolm, print.stsls,
## print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
## print.summary.stsls, residuals.gmsar, residuals.sarlm,
## residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup,
## SE_interp_setup, SE_whichMin_setup, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup,
## SpatialFiltering, spautolm, spBreg_err, spBreg_lag, spBreg_sac,
## stsls, subgraph_eigenw, summary.gmsar, summary.sarlm,
## summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlm## The following objects are masked from 'package:elsa':
##
## geary, moran## -- Attaching packages ----------------------------------------------- tidyverse 1.3.0 --## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## v purrr 0.3.4## Warning: package 'tibble' was built under R version 4.0.2## -- Conflicts -------------------------------------------------- tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x tidyr::expand() masks Matrix::expand()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x dplyr::recode() masks car::recode()
## x dplyr::select() masks raster::select(), MASS::select()
## x purrr::some() masks car::some()
## x tidyr::unpack() masks Matrix::unpack()## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## lift## Warning: package 'glmnet' was built under R version 4.0.2## Loaded glmnet 4.0-2## Warning: package 'spatialEco' was built under R version 4.0.2##
## Attaching package: 'spatialEco'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:spData':
##
## elev## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:grid':
##
## explode## The following object is masked from 'package:elsa':
##
## correlogram## The following object is masked from 'package:raster':
##
## shift## Warning: package 'lmSupport' was built under R version 4.0.2## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car## Loaded lars 1.2## Warning: package 'glmpath' was built under R version 4.0.2## Loading required package: survival##
## Attaching package: 'survival'## The following object is masked from 'package:spatialEco':
##
## concordance## The following object is masked from 'package:caret':
##
## clusterLoad Geography and Flat Files
setwd("D:/Covid19")
#Load shape and merge with data
myshape=shapefile("cb_2018_us_county_500k")
countydata=read.csv("reducedcounty4.csv",fileEncoding="UTF-8-BOM")
countydata$M=as.numeric(countydata$countyFIPS)Merge Shape and Flat File
myshape$M=as.numeric(myshape$GEOID)
counties=merge(myshape, countydata, by="M",type="left")Load Online Data & Merge
We write the file for replication at a specific date.
#deaths=read.csv("https://usafactsstatic.blob.core.windows.net/public/data/covid-19/covid_deaths_usafacts.csv", fileEncoding="UTF-8-BOM")
#write.csv(deaths,"deaths.csv", row.names=FALSE)Select Date for Run
deaths=read.csv("deaths.csv")
deaths=deaths[, c(1:4, ncol(deaths))] #this selects the last observation
colnames(deaths)=c("GEOID","County.Name","State","stateFIPS","deaths")
deaths$M=as.numeric(deaths$GEOID)Merge Online Data with Flat Files
counties=merge(counties,deaths, by="M",type="left")
counties$Deaths.100K=counties$deaths/(counties$Pop2020/100000)
counties@data$Deaths=counties@data$GEOID.y=counties@data$stateFIPS=counties@data$deaths=
counties@data$M=counties@data$STATEFP=counties@data$COUNTYFP=counties@data$COUNTYNS=
counties@data$AFFGEOID=counties@data$GEOID.x=counties@data$LSAD=counties@data$ALAND=
counties@data$AWATER=counties@data$County.Name=
counties@data$State=NULL
colnames(counties@data)## [1] "NAME" "countyFIPS" "Pop2020" "PopDensity" "Native"
## [6] "Hispanic" "Black" "Asian" "Prop65" "UE2019"
## [11] "Poverty" "Smoke" "AdultObesity" "Alcohol" "Alzheimers"
## [16] "Asthma" "AtrialFib" "Cancer" "COPD" "Depression"
## [21] "Diabetes" "DrugAbuse" "HIV" "HeartFail" "HepB"
## [26] "Stroke" "AcuteBeds" "CMI" "WParty" "Deaths.100K"Omit Incomplete Cases and Geography (sp.na.omit)
library(Amelia)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.6, built: 2019-11-24)
## ## Copyright (C) 2005-2020 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(counties@data)
counties2=sp.na.omit((counties))## Deleting rows: 2728293032232332432532632732848348450250350450550650750850951051175284688288395710681069107612001201120212031210139614301431144214431475147614771487151515261527152815841592160316041615161616171618161916251649172617331777186018641882188519091910191319321933195420222041205420622063206820822112211721692182218521862258242725392541257725782579258825922623263526592660266727092716272527432748276327862792280828442890289128972899290129343054missmap(counties2@data)
Describe Rate & Density & Proportion Democratic
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:lars':
##
## error.bars## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:car':
##
## logitdescribe(counties2@data$Deaths.100K)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 3116 34.03 46.75 17.75 24.33 26.32 0 461.16 461.16 2.88 12.42
## se
## X1 0.84describe(counties2@data$PopDensity)## vars n mean sd median trimmed mad min max range skew
## X1 1 3116 106.34 696.94 17.5 29.6 20.62 0.09 27755.49 27755.4 26.31
## kurtosis se
## X1 897.45 12.49Establish Correlations
library(corrplot)## corrplot 0.84 loadedmycol=colorRampPalette(c("red","orange","yellow",
"white","green", "dark green"))(20)
mycor=cor(counties2@data[,c(13:ncol(counties2@data))])
corrplot(mycor, method="ellipse", type="upper",
addCoef.col=TRUE, tl.cex=.7, number.cex=.5, insig="blank",
order="hclust", hclust.method="centroid", number.digits=2,
col=mycol)
Scale
counties3=counties2
mytrain=counties2@data[,c(4:ncol(counties2@data))]
colnames(mytrain)## [1] "PopDensity" "Native" "Hispanic" "Black" "Asian"
## [6] "Prop65" "UE2019" "Poverty" "Smoke" "AdultObesity"
## [11] "Alcohol" "Alzheimers" "Asthma" "AtrialFib" "Cancer"
## [16] "COPD" "Depression" "Diabetes" "DrugAbuse" "HIV"
## [21] "HeartFail" "HepB" "Stroke" "AcuteBeds" "CMI"
## [26] "WParty" "Deaths.100K"for (i in 1:ncol(mytrain)){mytrain[,i]=as.numeric(scale(mytrain[,i]))}
for (i in 4: ncol(counties2@data)){
counties2@data[,i]=as.numeric(scale(counties2@data[,i]))}
x=model.matrix(Deaths.100K~., data=mytrain)[,-ncol(mytrain)]
y=mytrain$Deaths.100KInitial Regression Model With All Variables
library(car)
library(MASS)
fullOLS=lm(Deaths.100K~., data=mytrain)
temp1=summary(fullOLS)
VIF=c(0,vif(fullOLS))
final=cbind(temp1$coefficients,VIF)
round(final,3)## Estimate Std. Error t value Pr(>|t|) VIF
## (Intercept) 0.000 0.014 0.000 1.000 0.000
## PopDensity 0.163 0.017 9.716 0.000 1.397
## Native 0.090 0.018 4.951 0.000 1.635
## Hispanic 0.133 0.022 5.924 0.000 2.490
## Black 0.408 0.029 14.156 0.000 4.126
## Asian 0.008 0.019 0.409 0.682 1.799
## Prop65 0.022 0.019 1.164 0.245 1.734
## UE2019 0.079 0.018 4.337 0.000 1.639
## Poverty 0.018 0.027 0.673 0.501 3.706
## Smoke -0.061 0.026 -2.336 0.020 3.418
## AdultObesity -0.045 0.019 -2.406 0.016 1.761
## Alcohol 0.041 0.020 2.108 0.035 1.907
## Alzheimers 0.112 0.021 5.261 0.000 2.234
## Asthma -0.049 0.020 -2.410 0.016 2.029
## AtrialFib 0.017 0.021 0.842 0.400 2.094
## Cancer -0.010 0.020 -0.490 0.624 2.004
## COPD -0.074 0.027 -2.805 0.005 3.493
## Depression 0.036 0.023 1.594 0.111 2.525
## Diabetes 0.183 0.027 6.856 0.000 3.514
## DrugAbuse -0.027 0.022 -1.231 0.218 2.422
## HIV -0.074 0.021 -3.578 0.000 2.105
## HeartFail -0.027 0.021 -1.286 0.199 2.202
## HepB -0.048 0.021 -2.303 0.021 2.170
## Stroke 0.092 0.022 4.225 0.000 2.369
## AcuteBeds -0.006 0.018 -0.331 0.741 1.583
## CMI 0.038 0.017 2.192 0.028 1.481
## WParty 0.029 0.019 1.477 0.140 1.849temp1$r.squared## [1] 0.3768526myres=residuals(fullOLS)
#mystep=stepAIC(fullOLS, direction="both", trace=FALSE)
#anova(mystep)Plot Full OLS Residuals on Map
counties2@data$res0=myres
counties2@data$sd_breaks0=scale(counties2@data$res0)[,1]
mymin=min(counties2@data$sd0)## Warning in min(counties2@data$sd0): no non-missing arguments to min; returning
## Infmymax=max(counties2@data$sd0)## Warning in max(counties2@data$sd0): no non-missing arguments to max; returning -
## Infmy_breaks =seq(-5,5)
tm_shape(counties2, bbox=tmaptools::bb(xlim=c(-14, -5), ylim=c(2, 9),width=5, ext=6, relative = TRUE)) + tm_fill("sd_breaks0", title = "Residuals", style = "fixed", breaks = my_breaks, palette = "-RdBu") + tm_borders(alpha = 0.1) + tm_layout(main.title = "Residuals", main.title.size = 1,
legend.position = c("right", "bottom"), legend.title.size = 0.8)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Warning: Values have found that are higher than the highest break## Variable(s) "sd_breaks0" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# Determine Spatial Analysis to Run for the Full Model
mylm0=lm(Deaths.100K~., data=mytrain)
list.queen<-poly2nb(counties2, queen=TRUE)
W=rwm=nb2listw(list.queen, style="W", zero.policy=TRUE)
test1=lm.morantest(mylm0, rwm, alternative="two.sided", zero.policy=TRUE)
m=test1$estimate[1]
n=test1$p.value
temptest=lm.LMtests(mylm0, rwm, test=c("LMerr","LMlag","RLMerr","RLMlag","SARMA"),zero.policy=TRUE)
a=temptest$LMerr$p.value
b=temptest$LMlag$p.value
y=temptest$RLMerr$p.value
z=temptest$RLMlag$p.value
ans=c(m,n,a,b,y,z)
names(ans)=c("Moran's I","p-value","LMerror", "LMLag", "RLMerr", "RLMlag")
ans## Moran's I p-value LMerror LMLag RLMerr
## 2.524971e-01 9.393718e-127 0.000000e+00 0.000000e+00 6.238845e-02
## RLMlag
## 0.000000e+00Lasso Prep
x=as.matrix(mytrain[,1:ncol(mytrain)-1])
y=as.matrix(mytrain[, ncol(mytrain)])
mynames=colnames(mytrain)[-ncol(mytrain)]
colnames(x)=mynamesLasso with Tuned Lambda
set.seed(1234)
lambda <- 10^seq(-3, 3, length = 100)
model1a <- train(
Deaths.100K ~., data = mytrain,
method = "glmnet",
tuneGrid = expand.grid(alpha = 1, lambda = lambda))## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.head(model1a)## $method
## [1] "glmnet"
##
## $modelInfo
## $modelInfo$label
## [1] "glmnet"
##
## $modelInfo$library
## [1] "glmnet" "Matrix"
##
## $modelInfo$type
## [1] "Regression" "Classification"
##
## $modelInfo$parameters
## parameter class label
## 1 alpha numeric Mixing Percentage
## 2 lambda numeric Regularization Parameter
##
## $modelInfo$grid
## function(x, y, len = NULL, search = "grid") {
## if(search == "grid") {
## numLev <- if(is.character(y) | is.factor(y)) length(levels(y)) else NA
## if(!is.na(numLev)) {
## fam <- ifelse(numLev > 2, "multinomial", "binomial")
## } else fam <- "gaussian"
## init <- glmnet::glmnet(Matrix::as.matrix(x), y,
## family = fam,
## nlambda = len+2,
## alpha = .5)
## lambda <- unique(init$lambda)
## lambda <- lambda[-c(1, length(lambda))]
## lambda <- lambda[1:min(length(lambda), len)]
## out <- expand.grid(alpha = seq(0.1, 1, length = len),
## lambda = lambda)
## } else {
## out <- data.frame(alpha = runif(len, min = 0, 1),
## lambda = 2^runif(len, min = -10, 3))
## }
## out
## }
##
## $modelInfo$loop
## function(grid) {
## alph <- unique(grid$alpha)
## loop <- data.frame(alpha = alph)
## loop$lambda <- NA
## submodels <- vector(mode = "list", length = length(alph))
## for(i in seq(along = alph)) {
## np <- grid[grid$alpha == alph[i],"lambda"]
## loop$lambda[loop$alpha == alph[i]] <- np[which.max(np)]
## submodels[[i]] <- data.frame(lambda = np[-which.max(np)])
## }
## list(loop = loop, submodels = submodels)
## }
## <bytecode: 0x00000000067b2f70>
##
## $modelInfo$fit
## function(x, y, wts, param, lev, last, classProbs, ...) {
## numLev <- if(is.character(y) | is.factor(y)) length(levels(y)) else NA
##
## theDots <- list(...)
##
## if(all(names(theDots) != "family")) {
## if(!is.na(numLev)) {
## fam <- ifelse(numLev > 2, "multinomial", "binomial")
## } else fam <- "gaussian"
## theDots$family <- fam
## }
##
## ## pass in any model weights
## if(!is.null(wts)) theDots$weights <- wts
##
## if(!(class(x)[1] %in% c("matrix", "sparseMatrix")))
## x <- Matrix::as.matrix(x)
##
## modelArgs <- c(list(x = x,
## y = y,
## alpha = param$alpha),
## theDots)
##
## out <- do.call(glmnet::glmnet, modelArgs)
## if(!is.na(param$lambda[1])) out$lambdaOpt <- param$lambda[1]
## out
## }
## <bytecode: 0x00000000522eea50>
##
## $modelInfo$predict
## function(modelFit, newdata, submodels = NULL) {
## if(!is.matrix(newdata)) newdata <- Matrix::as.matrix(newdata)
## if(length(modelFit$obsLevels) < 2) {
## out <- predict(modelFit, newdata, s = modelFit$lambdaOpt)
## } else {
## out <- predict(modelFit, newdata, s = modelFit$lambdaOpt, type = "class")
## }
## if(is.matrix(out)) out <- out[,1]
##
## if(!is.null(submodels)) {
## if(length(modelFit$obsLevels) < 2) {
## tmp <- as.list(as.data.frame(predict(modelFit, newdata, s = submodels$lambda),
## stringsAsFactors = TRUE))
## } else {
## tmp <- predict(modelFit, newdata, s = submodels$lambda, type = "class")
## tmp <- if(is.matrix(tmp)) as.data.frame(tmp, stringsAsFactors = FALSE) else as.character(tmp)
## tmp <- as.list(tmp)
## }
## out <- c(list(out), tmp)
## }
## out
## }
## <bytecode: 0x00000000424e96a8>
##
## $modelInfo$prob
## function(modelFit, newdata, submodels = NULL) {
## obsLevels <- if("classnames" %in% names(modelFit)) modelFit$classnames else NULL
## probs <- predict(modelFit,
## Matrix::as.matrix(newdata),
## s = modelFit$lambdaOpt,
## type = "response")
## if(length(obsLevels) == 2) {
## probs <- as.vector(probs)
## probs <- as.data.frame(cbind(1-probs, probs), stringsAsFactors = FALSE)
## colnames(probs) <- modelFit$obsLevels
## } else {
## probs <- as.data.frame(probs[,,1,drop = FALSE], stringsAsFactors = FALSE)
## names(probs) <- modelFit$obsLevels
## }
## if(!is.null(submodels)) {
## tmp <- predict(modelFit,
## Matrix::as.matrix(newdata),
## s = submodels$lambda,
## type = "response")
## if(length(obsLevels) == 2) {
## tmp <- as.list(as.data.frame(tmp, stringsAsFactors = TRUE))
## tmp <- lapply(tmp,
## function(x, lev) {
## x <- as.vector(x)
## tmp <- data.frame(1-x, x)
## names(tmp) <- lev
## tmp
## },
## lev = modelFit$obsLevels)
## } else tmp <- apply(tmp, 3, function(x) data.frame(x))
## probs <- if(is.list(tmp)) c(list(probs), tmp) else list(probs, tmp)
## }
## probs
## }
##
## $modelInfo$predictors
## function(x, lambda = NULL, ...) {
## if(is.null(lambda))
## {
## if(length(lambda) > 1) stop("Only one value of lambda is allowed right now")
## if(!is.null(x$lambdaOpt)) {
## lambda <- x$lambdaOpt
## } else stop("must supply a value of lambda")
## }
## allVar <- if(is.list(x$beta)) rownames(x$beta[[1]]) else rownames(x$beta)
## out <- unlist(predict(x, s = lambda, type = "nonzero"))
## out <- unique(out)
## if(length(out) > 0) {
## out <- out[!is.na(out)]
## out <- allVar[out]
## }
## out
## }
##
## $modelInfo$varImp
## function(object, lambda = NULL, ...) {
## if(is.null(lambda)) {
## if(length(lambda) > 1) stop("Only one value of lambda is allowed right now")
## if(!is.null(object$lambdaOpt)) {
## lambda <- object$lambdaOpt
## } else stop("must supply a value of lambda")
## }
## beta <- predict(object, s = lambda, type = "coef")
## if(is.list(beta)) {
## out <- do.call("cbind", lapply(beta, function(x) x[,1]))
## out <- as.data.frame(out, stringsAsFactors = TRUE)
## } else out <- data.frame(Overall = beta[,1])
## out <- abs(out[rownames(out) != "(Intercept)",,drop = FALSE])
## out
## }
##
## $modelInfo$levels
## function(x) if(any(names(x) == "obsLevels")) x$obsLevels else NULL
##
## $modelInfo$tags
## [1] "Generalized Linear Model" "Implicit Feature Selection"
## [3] "L1 Regularization" "L2 Regularization"
## [5] "Linear Classifier" "Linear Regression"
##
## $modelInfo$sort
## function(x) x[order(-x$lambda, x$alpha),]
##
## $modelInfo$trim
## function(x) {
## x$call <- NULL
## x$df <- NULL
## x$dev.ratio <- NULL
## x
## }
##
##
## $modelType
## [1] "Regression"
##
## $results
## alpha lambda RMSE Rsquared MAE RMSESD RsquaredSD
## 1 1 1.000000e-03 0.8019046 0.3511695 0.5151952 0.04485129 0.02862406
## 2 1 1.149757e-03 0.8018485 0.3512109 0.5151332 0.04485889 0.02861459
## 3 1 1.321941e-03 0.8017824 0.3512581 0.5150611 0.04486558 0.02860211
## 4 1 1.519911e-03 0.8017078 0.3513146 0.5149809 0.04487146 0.02859093
## 5 1 1.747528e-03 0.8016239 0.3513766 0.5148911 0.04487895 0.02858305
## 6 1 2.009233e-03 0.8015313 0.3514444 0.5147923 0.04488979 0.02857464
## 7 1 2.310130e-03 0.8014311 0.3515152 0.5146820 0.04490056 0.02856069
## 8 1 2.656088e-03 0.8013253 0.3515848 0.5145605 0.04491194 0.02853915
## 9 1 3.053856e-03 0.8012124 0.3516568 0.5144311 0.04492172 0.02851177
## 10 1 3.511192e-03 0.8010884 0.3517363 0.5142866 0.04493361 0.02848474
## 11 1 4.037017e-03 0.8009552 0.3518177 0.5141341 0.04494714 0.02844961
## 12 1 4.641589e-03 0.8008162 0.3518983 0.5139760 0.04496742 0.02840986
## 13 1 5.336699e-03 0.8006780 0.3519660 0.5138172 0.04498949 0.02835059
## 14 1 6.135907e-03 0.8005431 0.3520211 0.5136675 0.04501209 0.02828369
## 15 1 7.054802e-03 0.8004219 0.3520501 0.5135447 0.04502644 0.02823126
## 16 1 8.111308e-03 0.8003353 0.3520271 0.5134615 0.04504203 0.02816539
## 17 1 9.326033e-03 0.8003014 0.3519257 0.5134218 0.04506796 0.02808044
## 18 1 1.072267e-02 0.8003301 0.3517327 0.5134290 0.04511107 0.02797146
## 19 1 1.232847e-02 0.8004031 0.3514953 0.5134666 0.04517434 0.02780192
## 20 1 1.417474e-02 0.8005577 0.3511644 0.5135497 0.04528387 0.02757213
## 21 1 1.629751e-02 0.8008274 0.3506956 0.5137222 0.04545648 0.02732448
## 22 1 1.873817e-02 0.8011674 0.3501615 0.5139737 0.04567067 0.02709056
## 23 1 2.154435e-02 0.8016521 0.3494715 0.5143201 0.04591416 0.02682238
## 24 1 2.477076e-02 0.8023298 0.3485795 0.5147631 0.04617863 0.02661283
## 25 1 2.848036e-02 0.8031877 0.3475507 0.5152764 0.04647696 0.02635794
## 26 1 3.274549e-02 0.8042583 0.3463346 0.5158998 0.04687216 0.02592635
## 27 1 3.764936e-02 0.8054844 0.3450274 0.5167119 0.04725175 0.02550025
## 28 1 4.328761e-02 0.8069823 0.3435205 0.5177737 0.04766852 0.02504306
## 29 1 4.977024e-02 0.8088146 0.3417410 0.5191693 0.04817410 0.02443798
## 30 1 5.722368e-02 0.8107774 0.3401697 0.5208266 0.04863015 0.02374716
## 31 1 6.579332e-02 0.8132584 0.3381946 0.5229878 0.04908885 0.02297484
## 32 1 7.564633e-02 0.8163410 0.3357839 0.5256815 0.04968359 0.02214844
## 33 1 8.697490e-02 0.8201456 0.3328751 0.5290537 0.05042253 0.02140393
## 34 1 1.000000e-01 0.8249981 0.3289432 0.5333495 0.05143959 0.02048753
## 35 1 1.149757e-01 0.8311550 0.3235622 0.5388883 0.05249435 0.02022823
## 36 1 1.321941e-01 0.8388837 0.3161844 0.5458607 0.05353525 0.02086977
## 37 1 1.519911e-01 0.8481626 0.3067794 0.5542395 0.05438206 0.02096215
## 38 1 1.747528e-01 0.8586348 0.2964443 0.5633996 0.05486304 0.02126914
## 39 1 2.009233e-01 0.8694514 0.2881858 0.5726452 0.05476729 0.02108798
## 40 1 2.310130e-01 0.8810366 0.2823942 0.5823216 0.05535954 0.02439549
## 41 1 2.656088e-01 0.8949625 0.2748842 0.5938898 0.05647102 0.02674784
## 42 1 3.053856e-01 0.9117752 0.2638591 0.6078513 0.05746037 0.02702881
## 43 1 3.511192e-01 0.9304662 0.2572342 0.6234794 0.05850477 0.02452861
## 44 1 4.037017e-01 0.9527710 0.2554854 0.6421094 0.06008698 0.02309465
## 45 1 4.641589e-01 0.9789552 0.2506544 0.6637686 0.05930496 0.02065851
## 46 1 5.336699e-01 0.9914412 NaN 0.6741097 0.05540203 NA
## 47 1 6.135907e-01 0.9914412 NaN 0.6741097 0.05540203 NA
## 48 1 7.054802e-01 0.9914412 NaN 0.6741097 0.05540203 NA
## 49 1 8.111308e-01 0.9914412 NaN 0.6741097 0.05540203 NA
## 50 1 9.326033e-01 0.9914412 NaN 0.6741097 0.05540203 NA
## 51 1 1.072267e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 52 1 1.232847e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 53 1 1.417474e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 54 1 1.629751e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 55 1 1.873817e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 56 1 2.154435e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 57 1 2.477076e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 58 1 2.848036e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 59 1 3.274549e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 60 1 3.764936e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 61 1 4.328761e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 62 1 4.977024e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 63 1 5.722368e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 64 1 6.579332e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 65 1 7.564633e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 66 1 8.697490e+00 0.9914412 NaN 0.6741097 0.05540203 NA
## 67 1 1.000000e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 68 1 1.149757e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 69 1 1.321941e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 70 1 1.519911e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 71 1 1.747528e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 72 1 2.009233e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 73 1 2.310130e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 74 1 2.656088e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 75 1 3.053856e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 76 1 3.511192e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 77 1 4.037017e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 78 1 4.641589e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 79 1 5.336699e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 80 1 6.135907e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 81 1 7.054802e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 82 1 8.111308e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 83 1 9.326033e+01 0.9914412 NaN 0.6741097 0.05540203 NA
## 84 1 1.072267e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 85 1 1.232847e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 86 1 1.417474e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 87 1 1.629751e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 88 1 1.873817e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 89 1 2.154435e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 90 1 2.477076e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 91 1 2.848036e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 92 1 3.274549e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 93 1 3.764936e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 94 1 4.328761e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 95 1 4.977024e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 96 1 5.722368e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 97 1 6.579332e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 98 1 7.564633e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 99 1 8.697490e+02 0.9914412 NaN 0.6741097 0.05540203 NA
## 100 1 1.000000e+03 0.9914412 NaN 0.6741097 0.05540203 NA
## MAESD
## 1 0.01749789
## 2 0.01750715
## 3 0.01751736
## 4 0.01752791
## 5 0.01754315
## 6 0.01755914
## 7 0.01757539
## 8 0.01759197
## 9 0.01760917
## 10 0.01762790
## 11 0.01765347
## 12 0.01768879
## 13 0.01772204
## 14 0.01775184
## 15 0.01775843
## 16 0.01776534
## 17 0.01775944
## 18 0.01776008
## 19 0.01775623
## 20 0.01775567
## 21 0.01776523
## 22 0.01776266
## 23 0.01772390
## 24 0.01770420
## 25 0.01767339
## 26 0.01768593
## 27 0.01766693
## 28 0.01763093
## 29 0.01761527
## 30 0.01762633
## 31 0.01763408
## 32 0.01767679
## 33 0.01777015
## 34 0.01798163
## 35 0.01831727
## 36 0.01868912
## 37 0.01883921
## 38 0.01866026
## 39 0.01825953
## 40 0.01858344
## 41 0.01934662
## 42 0.02003998
## 43 0.02078619
## 44 0.02200013
## 45 0.02079137
## 46 0.01753622
## 47 0.01753622
## 48 0.01753622
## 49 0.01753622
## 50 0.01753622
## 51 0.01753622
## 52 0.01753622
## 53 0.01753622
## 54 0.01753622
## 55 0.01753622
## 56 0.01753622
## 57 0.01753622
## 58 0.01753622
## 59 0.01753622
## 60 0.01753622
## 61 0.01753622
## 62 0.01753622
## 63 0.01753622
## 64 0.01753622
## 65 0.01753622
## 66 0.01753622
## 67 0.01753622
## 68 0.01753622
## 69 0.01753622
## 70 0.01753622
## 71 0.01753622
## 72 0.01753622
## 73 0.01753622
## 74 0.01753622
## 75 0.01753622
## 76 0.01753622
## 77 0.01753622
## 78 0.01753622
## 79 0.01753622
## 80 0.01753622
## 81 0.01753622
## 82 0.01753622
## 83 0.01753622
## 84 0.01753622
## 85 0.01753622
## 86 0.01753622
## 87 0.01753622
## 88 0.01753622
## 89 0.01753622
## 90 0.01753622
## 91 0.01753622
## 92 0.01753622
## 93 0.01753622
## 94 0.01753622
## 95 0.01753622
## 96 0.01753622
## 97 0.01753622
## 98 0.01753622
## 99 0.01753622
## 100 0.01753622
##
## $pred
## NULL
##
## $bestTune
## alpha lambda
## 17 1 0.009326033Use Tuning to Generate Estimate
set.seed(1234)
lambda=as.numeric(model1a$bestTune[2])
model1=glmnet(x,y, data=mytrain,alpha=1, lambda=lambda)
a=as.vector(round(coef(model1), 3)) #get the coefficients less interceptGet Coefficients from lars package
set.seed(1234)
mylars=lars(x,y,intercept=TRUE, max.steps=10000)
b=as.matrix(round(coefficients(mylars)[26,],3))
#temp=rep(0, length(x)*nrow(x))
#fakedata=matrix(data=temp, ncol=ncol(x))
#temp=predict(mylars, fakedata)
myint=-2.549066
mylars2=lars(x,y,intercept=TRUE, type="lasso")
b=c(-2.549,round(coef(mylars2)[mylars2$RSS==min(mylars2$RSS)],3))Get adjusted p-values from covTest
mycovTest=covTest(mylars2,x,y)
res=as.data.frame(mycovTest$results)
pvals=res[order(res$Predictor_Number), ]
#pvals=pvals[, -1]
pvals## Predictor_Number Drop_in_covariance P-value
## X.3 1 3.3247 0.0361
## X.8 2 4.3528 0.0130
## X.4 3 32.6801 0.0000
## X 4 457.8677 0.0000
## X.23 5 0.0355 0.9651
## X.19 6 0.1963 0.8218
## X.7 7 4.0297 0.0179
## X.18 8 0.1567 0.8550
## X.13 9 0.1107 0.8952
## X.11 10 0.1597 0.8524
## X.14 11 0.5182 0.5957
## X.2 12 9.6050 0.0001
## X.17 13 0.1410 0.8685
## X.22 14 0.1549 0.8565
## X.24 15 0.0035 0.9965
## X.9 16 12.4173 0.0000
## X.21 17 0.3514 0.7037
## X.1 18 56.2859 0.0000
## X.16 19 0.1222 0.8850
## X.12 20 2.7017 0.0673
## X.20 21 0.2473 0.7809
## X.10 22 1.6949 0.1838
## X.5 23 0.8393 0.4321
## X.25 24 0.1095 0.8963
## X.15 25 0.2208 0.8019
## X.6 26 2.4636 0.0853Compare Coefficients from glmnet and lars, p-values
newdata=data.frame(lasso_glmnet=a,
lasso_lars=b,
pvalue=c(0,pvals[,3]))
rownames(newdata)=c("Intercept",mynames)
#newdata
findata=newdata[newdata$pvalue<.10,]
findata[-1,]## lasso_glmnet lasso_lars pvalue
## PopDensity 0.146 0.163 0.0361
## Native 0.064 0.090 0.0130
## Hispanic 0.120 0.133 0.0000
## Black 0.376 0.408 0.0000
## UE2019 0.071 0.079 0.0179
## Alzheimers 0.112 0.112 0.0001
## COPD -0.073 -0.074 0.0000
## Diabetes 0.152 0.183 0.0000
## HIV -0.048 -0.074 0.0673
## WParty 0.026 0.029 0.0853Lasso Model
x2=as.matrix(mytrain[,c(1:4,6:7,12,16,18, 26)])
y2=as.vector(mytrain[, 27])
mynames=colnames(x2)
mylars3=lars(x2, y2, intercept=T, type="lasso")
larscoef=as.matrix(coef(mylars3)[max(which.min(mylars3$RSS)),])
temp=as.data.frame(covTest(mylars3, x2, y2)$results)
temp=temp[order(temp$Predictor_Number),]
temp$beta=larscoef
rownames(temp)=mynames
mypred=matrix(rep(0,ncol(x2)*nrow(x2)), ncol=ncol(x2))
myint=predict(mylars3, mypred)$fit[1]
temp## Predictor_Number Drop_in_covariance P-value beta
## PopDensity 1 3.2666 0.0383 0.13755784
## Native 2 3.2815 0.0377 0.05657627
## Hispanic 3 45.8823 0.0000 0.13183637
## Black 4 449.8628 0.0000 0.36899848
## Prop65 5 2.3870 0.0921 0.02507298
## UE2019 6 4.9396 0.0072 0.07543364
## Alzheimers 7 9.4370 0.0001 0.14897808
## COPD 8 18.5452 0.0000 -0.10427388
## Diabetes 9 55.3018 0.0000 0.16192376
## WParty 10 2.4184 0.0892 0.02432629myint## [1] 2.423441e-17Generate Predictions
predictions1 <- mylars3 %>% predict(x2)
res1=y2-predictions1$fit
hist(res1)
Lasso Residual Analysis
counties2@data$res1=res1
counties2@data$sd_breaks1=scale(counties2@data$res1)[,1]
mymin=min(counties2@data$sd)## Warning in min(counties2@data$sd): no non-missing arguments to min; returning
## Infmymax=max(counties2@data$sd)## Warning in max(counties2@data$sd): no non-missing arguments to max; returning -
## Infmyrange=max(abs(mymin), mymax)
my_breaks =seq(-5,5)
tm_shape(counties2, bbox=tmaptools::bb(xlim=c(-14, -5), ylim=c(2, 9),width=5, ext=5, relative = TRUE)) + tm_fill("sd_breaks1", title = "Residuals", style = "fixed", breaks = my_breaks, palette = "-RdBu") + tm_borders(alpha = 0.1) + tm_layout(main.title = "Residuals", main.title.size = 1,
legend.position = c("right", "bottom"), legend.title.size = 0.8)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Warning: Values have found that are higher than the highest break## Variable(s) "sd_breaks1" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Determine Spatial Analysis to Run, Reduced Model
mylm=lm(Deaths.100K~PopDensity+Native+Hispanic+Black+UE2019+Alzheimers+
COPD+Diabetes+WParty, data=mytrain)
test1=lm.morantest(mylm, rwm, alternative="two.sided", zero.policy=TRUE)
m=test1$estimate[1]
n=test1$p.value
temptest=lm.LMtests(mylm, rwm, test=c("LMerr","LMlag","RLMerr","RLMlag","SARMA"),zero.policy=TRUE)
a=temptest$LMerr$p.value
b=temptest$LMlag$p.value
y=temptest$RLMerr$p.value
z=temptest$RLMlag$p.value
ans=c(m,n,a,b,y,z)
names(ans)=c("Moran's I","p-value","LMerror", "LMLag", "RLMerr", "RLMlag")
ans## Moran's I p-value LMerror LMLag RLMerr
## 2.646846e-01 3.753957e-136 0.000000e+00 0.000000e+00 4.642088e-03
## RLMlag
## 0.000000e+00Run 1st Spatial Analysis, Scaled
mylag1=stsls(Deaths.100K~.,
zero.policy=TRUE, data=counties2@data[4:30], W)## Warning: Function stsls moved to the spatialreg packagecounties2@data$reslag1=mylag1$res
counties2@data$sd_breaks2=scale(counties2@data$reslag1)[,1]
moran1=moran.mc(counties2@data$reslag1, W, 50, zero.policy=TRUE) #LISA
mylisa1=lisa(counties2, d1=0, d2=2000, statistic="I", zcol="Deaths.100K")
summary(mylag1)##
## Call:spatialreg::stsls(formula = formula, data = data, listw = listw,
## zero.policy = zero.policy, na.action = na.action, robust = robust,
## HC = HC, legacy = legacy, W2X = W2X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.20336 -0.31444 -0.10689 0.16712 8.33438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Rho 0.6338717 0.0371946 17.0421 < 2.2e-16
## (Intercept) -0.0043200 0.0126343 -0.3419 0.7324087
## PopDensity 0.0663413 0.0159780 4.1521 3.295e-05
## Native 0.0695929 0.0161973 4.2966 1.735e-05
## Hispanic 0.0819591 0.0201559 4.0663 4.777e-05
## Black 0.1780831 0.0290048 6.1398 8.263e-10
## Asian -0.0093791 0.0169740 -0.5526 0.5805662
## Prop65 0.0222192 0.0166369 1.3355 0.1817017
## UE2019 0.0521403 0.0162481 3.2090 0.0013319
## Poverty 0.0120178 0.0243237 0.4941 0.6212517
## Smoke -0.0055221 0.0235848 -0.2341 0.8148768
## AdultObesity 0.0060735 0.0170369 0.3565 0.7214733
## Alcohol 0.0239738 0.0174748 1.3719 0.1700914
## Alzheimers 0.0734463 0.0190156 3.8624 0.0001123
## Asthma -0.0222880 0.0180645 -1.2338 0.2172761
## AtrialFib 0.0105867 0.0182875 0.5789 0.5626556
## Cancer -0.0157509 0.0178875 -0.8806 0.3785614
## COPD -0.0467567 0.0236687 -1.9755 0.0482151
## Depression 0.0426598 0.0200792 2.1246 0.0336221
## Diabetes 0.0779494 0.0244658 3.1861 0.0014422
## DrugAbuse -0.0327577 0.0196639 -1.6659 0.0957382
## HIV -0.0466848 0.0183977 -2.5375 0.0111636
## HeartFail -0.0088772 0.0187796 -0.4727 0.6364252
## HepB -0.0310786 0.0186398 -1.6673 0.0954491
## Stroke 0.0264753 0.0198258 1.3354 0.1817450
## AcuteBeds 0.0091541 0.0159198 0.5750 0.5652825
## CMI 0.0448932 0.0153793 2.9191 0.0035108
## WParty 0.0461883 0.0172089 2.6840 0.0072752
##
## Residual variance (sigma squared): 0.49719, (sigma: 0.70512)moran1##
## Monte-Carlo simulation of Moran I
##
## data: counties2@data$reslag1
## weights: W
## number of simulations + 1: 51
##
## statistic = -0.093994, observed rank = 1, p-value = 0.9804
## alternative hypothesis: greatermylisa1## class : SpatialPolygonsDataFrame
## features : 3116
## extent : -178.3347, -66.94989, 18.91036, 65.45352 (xmin, xmax, ymin, ymax)## Warning in proj4string(x): CRS object has comment, which is lost in output## crs : +proj=longlat +datum=NAD83 +no_defs
## variables : 2
## names : Ii, Z.Ii
## min values : -1.13779547311883e-15, 0.00031933549932932
## max values : 4.80093306596748e-15, 0.00031933549932932R2
totalss=var(counties2@data$Deaths.100K)*(3115)
1-mylag1$sse/totalss## [1] 0.5071169Evaluate Spatial Residuals
tm_shape(counties2, bbox=tmaptools::bb(xlim=c(-14, -5), ylim=c(2, 9),width=5, ext=5, relative = TRUE)) + tm_fill("sd_breaks2", title = "Residuals", style = "fixed", breaks = my_breaks, palette = "-RdBu") + tm_borders(alpha = 0.1) + tm_layout(main.title = "Residuals", main.title.size = 1,
legend.position = c("right", "bottom"), legend.title.size = 0.8)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Warning: Values have found that are higher than the highest break## Variable(s) "sd_breaks2" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Reduced GIS Model based on Lasso Coefficients, Scaled
mylag2=stsls(Deaths.100K~PopDensity+Native+Hispanic+Black+
UE2019+Alzheimers+COPD+Diabetes+WParty,
zero.policy=TRUE, data=counties2@data[4:31], W)## Warning: Function stsls moved to the spatialreg packagecounties2@data$reslag2=mylag2$res
counties2@data$sd_breaks3=scale(counties2@data$reslag2)[,1]
moran2=moran.mc(counties2@data$reslag2, W, 50, zero.policy=TRUE) #LISA
summary(mylag2)##
## Call:spatialreg::stsls(formula = formula, data = data, listw = listw,
## zero.policy = zero.policy, na.action = na.action, robust = robust,
## HC = HC, legacy = legacy, W2X = W2X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.16839 -0.31934 -0.11259 0.15880 8.47341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Rho 0.5885137 0.0381504 15.4262 < 2.2e-16
## (Intercept) -0.0040108 0.0126899 -0.3161 0.7519535
## PopDensity 0.0509149 0.0143005 3.5604 0.0003703
## Native 0.0592370 0.0136762 4.3314 1.482e-05
## Hispanic 0.0707999 0.0156422 4.5262 6.005e-06
## Black 0.1685428 0.0222440 7.5770 3.531e-14
## UE2019 0.0623166 0.0146281 4.2600 2.044e-05
## Alzheimers 0.0973908 0.0157923 6.1670 6.960e-10
## COPD -0.0525788 0.0191494 -2.7457 0.0060379
## Diabetes 0.0794890 0.0214076 3.7131 0.0002047
## WParty 0.0323228 0.0151588 2.1323 0.0329835
##
## Residual variance (sigma squared): 0.50157, (sigma: 0.70822)moran2##
## Monte-Carlo simulation of Moran I
##
## data: counties2@data$reslag2
## weights: W
## number of simulations + 1: 51
##
## statistic = -0.069989, observed rank = 1, p-value = 0.9804
## alternative hypothesis: greaterR2 Calculation
1-mylag2$sse/totalss## [1] 0.500038Reduced GIS Model based on Lasso Coefficients, Unscaled
mylag3=stsls(Deaths.100K~PopDensity+Native+Hispanic+Black+
UE2019+Alzheimers+COPD+Diabetes+WParty,
zero.policy=TRUE, data=counties2@data[4:31], W)## Warning: Function stsls moved to the spatialreg packagesummary(mylag3)##
## Call:spatialreg::stsls(formula = formula, data = data, listw = listw,
## zero.policy = zero.policy, na.action = na.action, robust = robust,
## HC = HC, legacy = legacy, W2X = W2X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.16839 -0.31934 -0.11259 0.15880 8.47341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Rho 0.5885137 0.0381504 15.4262 < 2.2e-16
## (Intercept) -0.0040108 0.0126899 -0.3161 0.7519535
## PopDensity 0.0509149 0.0143005 3.5604 0.0003703
## Native 0.0592370 0.0136762 4.3314 1.482e-05
## Hispanic 0.0707999 0.0156422 4.5262 6.005e-06
## Black 0.1685428 0.0222440 7.5770 3.531e-14
## UE2019 0.0623166 0.0146281 4.2600 2.044e-05
## Alzheimers 0.0973908 0.0157923 6.1670 6.960e-10
## COPD -0.0525788 0.0191494 -2.7457 0.0060379
## Diabetes 0.0794890 0.0214076 3.7131 0.0002047
## WParty 0.0323228 0.0151588 2.1323 0.0329835
##
## Residual variance (sigma squared): 0.50157, (sigma: 0.70822)Plot Death Rate by Percent Democrat
forplot=read.csv("partypolitics.csv")
boxplot(forplot$Death_Rate~forplot$Winning_Party, col=c("red", "blue"),
horizontal=TRUE, notch=TRUE, xlab="Death Rate", ylab="Party")