Accessability Maps

df = bus.ca %>% group_by(caid) %>% summarise(bus = mean(bus)); df$bus = df$bus/sum(df$bus)
m = map = makemap(df,"bus"); m
df = tnc.ca %>% group_by(pickupca) %>% summarise(tnc = mean(tnc)); df$tnc = df$tnc/sum(df$tnc)
m = map = makemap(df,"tnc",by.y = "pickupca"); m
df = cta.l %>% group_by(pickupca) %>% summarise(train = mean(m)); df$train = df$train/sum(df$train)
m = map = makemap(df,"train",by.y = "pickupca"); m

We use similar to that proposed by Bernardinelli, Clayton, Pascutto, Montomoli, Ghislandi, and Songini (1995), which represents the spatio-temporal pattern in the mean response with spatially varying linear time trends. We assume the data is Gaussian. The model estimates autocorrelated linear time trends for each community area (areal unit). Thus it is appropriate if the goal of the analysis is to estimate which areas are exhibiting increasing or decreasing (linear) trends in the response over time. The full model specification is given below.

We assume we have a set of \(K\) non-overlaping areal units \(S_1,\ldots,S_K\) and data is recorded for each areal unit for \(N\) consecutive time steps \(t = 1,\ldots,N\), then the general hierarchical model is \[ \begin{aligned} Y_{kt} \sim & f(y_{kt} \mid \mu_{kt},\nu^2)\\ g(\mu_{kt}) = & x^T_{kt}\beta+O_{kt} + \psi_{kt}\\ \beta \sim & N(\mu_{\beta},\Sigma_{\beta}) \end{aligned} \] Given that we assume Normal (continious) data, we use \(f =\) Noraml density and link function \(g\) to be identity function.

The spatio-temporal correlations are modeled using \[ \psi_{kt} = \beta_1 + \phi_k + (\alpha + \delta_k)\left(\dfrac{t - \bar t}{N}\right) \] Thus, the temporal pattern is modeled by a local trend \(\delta_k\) and the spatial correlations are modeled via \(\phi_k\), both are random effects, which are conditionally normal \[ \begin{aligned} \phi_k \mid \phi_{-k}, W \sim & N\left(\dfrac{1}{c}\sum_{j=1}^Kw_{kj}\phi_j, \dfrac{\tau^2_{int}}{c}\right),~~~c = \sum_{j=1}^Kw_{kj}-1 +1/\rho_{int}\\ \delta_k \mid \delta_{-k},W \sim & N\left(\dfrac{1}{q}\sum_{j=1}^Kw_{kj}\delta_j, \dfrac{\tau^2_{slo}}{q}\right),~~~q = \sum_{j=1}^Kw_{kj}-1 +1/\rho^2_{slo}\\ \tau^2_{int},\tau^2_{slo} \sim & IG(a,b)\\ \tau^2_{int},\tau^2_{slo} \sim & Uniform(0,1)\\ \alpha \sim N(\mu_{\alpha},\sigma^2_{\alpha}) \end{aligned} \] Eech community area has its own linear trend with inctercept \(\beta_1 + \phi_k\) and slope \(\alpha + \delta_k\). Here \(\rho_{int},\rho_{slo}\) are spatial dependence parameters, with values of one corresponding to strong spatial smoothness that is equivalent to the intrinsic CAR prior proposed by Besag et al. (1991), while values of zero correspond to independence (for example if \(\rho_{slo}\) then \(\delta_k \sim N(0, \tau^2_{slo})\).

Now lets look at the results of the models

ca = sort(unique(bus.tnc.ca$caid))
m = modelmap(readRDS("model-cache/bus-tnc-1.Rda"),ca); m

Findings:

m = readRDS("model-cache/bus-tnc-1.Rda")
m$summary.results[,1:3]
##                 Median       2.5%      97.5%
## (Intercept)    15.6881     9.2965    21.9667
## bus             0.4761     0.3795     0.5744
## alpha          -0.7932    -2.5505     0.9689
## tau2.int    15250.8194 10497.7775 22154.7443
## tau2.slo        0.0084     0.0021     0.0925
## nu2            59.8091    54.6324    65.6019
## rho.int         0.7446     0.4245     0.9624
## rho.slo         0.3770     0.0167     0.9195
library(coda)
beta.samples <- mcmc.list(m$samples$beta)
plot(beta.samples)

m = readRDS("model-cache/bus-tnc-1.Rda")
apply(m$samples[["phi"]],2,quantile,probs=c(0.05,0.95))[,1:5]
##          var1      var2     var3       var4       var5
## 5%  -29.20275 -24.20779 39.98131 -17.584224 -16.614614
## 95% -21.49399 -16.38874 48.00248  -9.922877  -9.513528
phi.samples <- mcmc.list(m$samples[["phi"]][,1:10])
plot(phi.samples)

- There is no local effet on the trend

m = readRDS("model-cache/bus-tnc-1.Rda")
apply(m$samples[["delta"]],2,quantile,probs=c(0.05,0.95))[,1:5]
##           var1       var2       var3       var4       var5
## 5%  -0.1763081 -0.1569927 -0.1427589 -0.1231031 -0.1213945
## 95%  0.1691506  0.1498601  0.1530900  0.1201097  0.1182249
hist(m$residuals$pearson,breaks = 100, main="Standartized residuals")

# With L-System

model = readRDS("model-cache/tnc-train.Rda")
ca = sort(unique(bus.tnc.cta.l.ca$caid))
m = modelmap(model,ca); m
model$summary.results[,1:3]
##                  Median        2.5%       97.5%
## (Intercept)    127.5000     89.8817    164.2657
## train            0.4308      0.2895      0.5744
## alpha           11.5450     -7.8647     28.3170
## tau2.int    388636.9237 229171.2010 740640.3867
## tau2.slo         0.0080      0.0021      0.1035
## nu2           3288.6626   2913.6250   3711.6146
## rho.int          0.1601      0.0090      0.5657
## rho.slo          0.3755      0.0191      0.9173
hist(model$residuals$pearson,breaks = 100, main="Standartized residuals")

With L-System and bus

model = readRDS("model-cache/tnc-bus_train.Rda")
ca = sort(unique(bus.tnc.cta.l.ca$caid))
m = modelmap(model,ca); m
model$summary.results[,1:3]
##                  Median        2.5%       97.5%
## (Intercept)    -33.3342   -100.1132     32.4743
## bus              2.4717      1.5956      3.3383
## train            0.2772      0.1223      0.4227
## alpha            0.0080    -17.0117     17.4833
## tau2.int    328802.0112 200068.5778 628969.1932
## tau2.slo         0.0083      0.0021      0.0912
## nu2           3130.1253   2767.8649   3530.2913
## rho.int          0.1176      0.0053      0.4641
## rho.slo          0.3699      0.0162      0.9176
hist(model$residuals$pearson,breaks = 100, main="Standartized residuals")