Earthquake depth and the number of stations reporting
Introduction
When looking into research on earthquake depth, Earthquakes can occur anywhere between the Earth’s surface and about 700 kilometers below the surface. Earthquake depth range of 0 - 700 km is divided into three zones: shallow, intermediate, and deep. The strength of shaking from an earthquake diminishes with increasing distance from the earthquake’s source, so the strength of shaking at the surface from an earthquake that occurs at 500km deep is considerably less than if the same earthquake had occurred at 20 km depth.
Packages Used
library(Zelig)
library(dplyr)
library(tidyr)
library(survival)
Data
data(quakes)
head(quakes)
Explaining Data
The data above shows the number of stations reporting after different level earthquakes.A data frame with 1000 observations on 5 variables. My hypothesis for this analysis is that how many times stations will be reporting is directly effected by the depth of the earth quake. I believe that earthquakes with a larger depth cause more damage, therefore causing more stations to report on it.
Poisson Model
Does depth of the earthquake effect how many stations are reporting?
z.out <- zelig(stations ~ mag + depth + long + lat, model = "poisson", data = quakes, cite = F)
summary(z.out)
Model:
Call:
z5$zelig(formula = stations ~ mag + depth + long + lat, data = quakes)
Deviance Residuals:
Min 1Q Median 3Q Max
-7.3543 -1.1201 -0.1238 0.9457 5.9071
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.9057762 0.1848110 -21.134 < 2e-16
mag 1.2088383 0.0118700 101.840 < 2e-16
depth 0.0002722 0.0000257 10.591 < 2e-16
long 0.0098097 0.0009717 10.095 < 2e-16
lat 0.0068245 0.0011614 5.876 4.2e-09
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 12198.5 on 999 degrees of freedom
Residual deviance: 2764.3 on 995 degrees of freedom
AIC: 7950.4
Number of Fisher Scoring iterations: 4
Next step: Use 'setx' method
A Poisson model is useful for when you have a count variable (how many times does something occur). For this data I am using the poisson model because of how many times stations are reporting. I am Using the “Count Poisson Model” because the main response variable is an integer. I will estimate the depth of an earthquake given its physical measurements.
Magnitude Range
a.range = min(quakes$mag):max(quakes$mag)
x <- setx(z, mag = a.range)
s <- sim(z, x = x)
ci.plot(s)

Depth Range
a.range = min(quakes$depth):max(quakes$depth)
x <- setx(z, depth = a.range)
s <- sim(z, x = x)
ci.plot(s)

x.low <- setx(z.out, depth = 100)
x.high <- setx1(z.out, depth = 600)
x <- setx(z, depth = mean(quakes$depth))
x1 <- setx(z, depth = mean(quakes$depth) + sd(quakes$depth))
s <- sim(z, x = x, x1 = x1)
summary(s)
sim x :
-----
ev
mean sd 50% 2.5% 97.5%
[1,] 29.47034 0.1822166 29.47298 29.1255 29.81901
pv
mean sd 50% 2.5% 97.5%
[1,] 29.243 5.255623 29 19 40
sim x1 :
-----
ev
mean sd 50% 2.5% 97.5%
[1,] 31.25722 0.258944 31.262 30.75319 31.76446
pv
mean sd 50% 2.5% 97.5%
[1,] 31.467 5.664362 32 21 42
fd
mean sd 50% 2.5% 97.5%
[1,] 1.78688 0.1655012 1.785758 1.470837 2.109044
From these models, we can see that depth does effect the amount of stations that are reporting. The higher (or deeper) the depth of the earthquakes, the more stations that are reporting.
Plots
plot(s.out)

Quartiles of Depth
xl <- setx(z, depth = quantile(quakes$depth, .25))
xl1 <- setx(z, depth = quantile(quakes$depth, .25) + sd(quakes$depth))
sl <- sim(z, x = xl, x1 = xl1)
fdl <- sl$get_qi(xvalue="x1", qi="fd")
xm <- setx(z, depth = quantile(quakes$depth, .50))
xm1 <- setx(z, depth = quantile(quakes$depth, .50) + sd(quakes$depth))
sm <- sim(z, x = xm, x1 = xm1)
fdm <- sm$get_qi(xvalue="x1", qi="fd")
xh <- setx(z, depth = quantile(quakes$depth, .75))
xh1 <- setx(z, depth = quantile(quakes$depth, .75) + sd(quakes$depth))
sh <- sim(z, x = xh, x1 = xh1)
fdh <- sh$get_qi(xvalue="x1", qi="fd")
Combining the results
d <- as.data.frame(cbind(fdl, fdm, fdh))
head(d)
tidd <- d %>%
gather(quantile, dif, 1:3)
head(tidd)
tidd %>%
group_by(quantile) %>%
summarise(mean = mean(dif), sd = sd(dif))
GG Plot
library(ggplot2)
ggplot(tidd, aes(dif)) + geom_histogram() + facet_grid(~ quantile)

Overall Thoughts
What we have learned overall: Based on this data set, we learned, that when it comes to earthquakes, stations are more interested in reporting on earthquakes that have a larger depth, as well as an overall larger magnitude.
LS0tDQp0aXRsZTogIkVhcnRocXVha2UgTWFnbml0dWRlIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCiMjIEVhcnRocXVha2UgZGVwdGggYW5kIHRoZSBudW1iZXIgb2Ygc3RhdGlvbnMgcmVwb3J0aW5nIA0KDQojIyBJbnRyb2R1Y3Rpb24gDQoqV2hlbiBsb29raW5nIGludG8gcmVzZWFyY2ggb24gZWFydGhxdWFrZSBkZXB0aCwgRWFydGhxdWFrZXMgY2FuIG9jY3VyIGFueXdoZXJlIGJldHdlZW4gdGhlIEVhcnRoJ3Mgc3VyZmFjZSBhbmQgYWJvdXQgNzAwIGtpbG9tZXRlcnMgYmVsb3cgdGhlIHN1cmZhY2UuIEVhcnRocXVha2UgZGVwdGggcmFuZ2Ugb2YgMCAtIDcwMCBrbSBpcyBkaXZpZGVkIGludG8gdGhyZWUgem9uZXM6IHNoYWxsb3csIGludGVybWVkaWF0ZSwgYW5kIGRlZXAuIFRoZSBzdHJlbmd0aCBvZiBzaGFraW5nIGZyb20gYW4gZWFydGhxdWFrZSBkaW1pbmlzaGVzIHdpdGggaW5jcmVhc2luZyBkaXN0YW5jZSBmcm9tIHRoZSBlYXJ0aHF1YWtlJ3Mgc291cmNlLCBzbyB0aGUgc3RyZW5ndGggb2Ygc2hha2luZyBhdCB0aGUgc3VyZmFjZSBmcm9tIGFuIGVhcnRocXVha2UgdGhhdCBvY2N1cnMgYXQgNTAwa20gZGVlcCBpcyBjb25zaWRlcmFibHkgbGVzcyB0aGFuIGlmIHRoZSBzYW1lIGVhcnRocXVha2UgaGFkIG9jY3VycmVkIGF0IDIwIGttIGRlcHRoLioNCg0KIyMgUGFja2FnZXMgVXNlZCANCmBgYHtyfQ0KbGlicmFyeShaZWxpZykNCmxpYnJhcnkoZHBseXIpDQpsaWJyYXJ5KHRpZHlyKQ0KbGlicmFyeShzdXJ2aXZhbCkNCmBgYA0KDQojIyBEYXRhIA0KYGBge3J9DQpkYXRhKHF1YWtlcykNCmhlYWQocXVha2VzKQ0KYGBgDQoNCiMjIEV4cGxhaW5pbmcgRGF0YSANCipUaGUgZGF0YSBhYm92ZSBzaG93cyB0aGUgbnVtYmVyIG9mIHN0YXRpb25zIHJlcG9ydGluZyBhZnRlciBkaWZmZXJlbnQgbGV2ZWwgZWFydGhxdWFrZXMuQSBkYXRhIGZyYW1lIHdpdGggMTAwMCBvYnNlcnZhdGlvbnMgb24gNSB2YXJpYWJsZXMuIE15IGh5cG90aGVzaXMgZm9yIHRoaXMgYW5hbHlzaXMgaXMgdGhhdCBob3cgbWFueSB0aW1lcyBzdGF0aW9ucyB3aWxsIGJlIHJlcG9ydGluZyBpcyBkaXJlY3RseSBlZmZlY3RlZCBieSB0aGUgZGVwdGggb2YgdGhlIGVhcnRoIHF1YWtlLiBJIGJlbGlldmUgdGhhdCBlYXJ0aHF1YWtlcyB3aXRoIGEgbGFyZ2VyIGRlcHRoIGNhdXNlIG1vcmUgZGFtYWdlLCB0aGVyZWZvcmUgY2F1c2luZyBtb3JlIHN0YXRpb25zIHRvIHJlcG9ydCBvbiBpdC4qDQoNCiMjIFBvaXNzb24gTW9kZWwgDQoqRG9lcyBkZXB0aCBvZiB0aGUgZWFydGhxdWFrZSBlZmZlY3QgaG93IG1hbnkgc3RhdGlvbnMgYXJlIHJlcG9ydGluZz8qDQpgYGB7cn0NCnoub3V0IDwtIHplbGlnKHN0YXRpb25zIH4gbWFnICsgZGVwdGggKyBsb25nICsgbGF0LCBtb2RlbCA9ICJwb2lzc29uIiwgZGF0YSA9IHF1YWtlcywgY2l0ZSA9IEYpDQpzdW1tYXJ5KHoub3V0KQ0KYGBgDQoqQSBQb2lzc29uIG1vZGVsIGlzIHVzZWZ1bCBmb3Igd2hlbiB5b3UgaGF2ZSBhIGNvdW50IHZhcmlhYmxlIChob3cgbWFueSB0aW1lcyBkb2VzIHNvbWV0aGluZyBvY2N1cikuIEZvciB0aGlzIGRhdGEgSSBhbSB1c2luZyB0aGUgcG9pc3NvbiBtb2RlbCBiZWNhdXNlIG9mIGhvdyBtYW55IHRpbWVzIHN0YXRpb25zIGFyZSByZXBvcnRpbmcuIEkgYW0gVXNpbmcgdGhlICJDb3VudCBQb2lzc29uIE1vZGVsIiBiZWNhdXNlIHRoZSBtYWluIHJlc3BvbnNlIHZhcmlhYmxlIGlzIGFuIGludGVnZXIuIEkgd2lsbCBlc3RpbWF0ZSB0aGUgZGVwdGggb2YgYW4gZWFydGhxdWFrZSBnaXZlbiBpdHMgcGh5c2ljYWwgbWVhc3VyZW1lbnRzLiogDQoNCiMjIE1hZ25pdHVkZSBSYW5nZSANCmBgYHtyfQ0KYS5yYW5nZSA9IG1pbihxdWFrZXMkbWFnKTptYXgocXVha2VzJG1hZykNCnggPC0gc2V0eCh6LCBtYWcgPSBhLnJhbmdlKQ0KcyA8LSBzaW0oeiwgeCA9IHgpDQpjaS5wbG90KHMpDQpgYGANCg0KIyMgRGVwdGggUmFuZ2UgDQpgYGB7cn0NCmEucmFuZ2UgPSBtaW4ocXVha2VzJGRlcHRoKTptYXgocXVha2VzJGRlcHRoKQ0KeCA8LSBzZXR4KHosIGRlcHRoID0gYS5yYW5nZSkNCnMgPC0gc2ltKHosIHggPSB4KQ0KY2kucGxvdChzKQ0KYGBgDQoNCmBgYHtyfQ0KeC5sb3cgPC0gc2V0eCh6Lm91dCwgZGVwdGggPSAxMDApDQp4LmhpZ2ggPC0gc2V0eDEoei5vdXQsIGRlcHRoID0gNjAwKQ0KYGBgDQoNCmBgYHtyfQ0KeCA8LSBzZXR4KHosIGRlcHRoID0gbWVhbihxdWFrZXMkZGVwdGgpKQ0KeDEgPC0gc2V0eCh6LCBkZXB0aCA9IG1lYW4ocXVha2VzJGRlcHRoKSArIHNkKHF1YWtlcyRkZXB0aCkpDQpzIDwtIHNpbSh6LCB4ID0geCwgeDEgPSB4MSkNCnN1bW1hcnkocykNCmBgYA0KKkZyb20gdGhlc2UgbW9kZWxzLCB3ZSBjYW4gc2VlIHRoYXQgZGVwdGggZG9lcyBlZmZlY3QgdGhlIGFtb3VudCBvZiBzdGF0aW9ucyB0aGF0IGFyZSByZXBvcnRpbmcuIFRoZSBoaWdoZXIgKG9yIGRlZXBlcikgdGhlIGRlcHRoIG9mIHRoZSBlYXJ0aHF1YWtlcywgdGhlIG1vcmUgc3RhdGlvbnMgdGhhdCBhcmUgcmVwb3J0aW5nLioNCg0KIyMgUGxvdHMgDQpgYGB7cn0NCnBsb3Qocy5vdXQpDQpgYGANCg0KIyMgUXVhcnRpbGVzIG9mIERlcHRoIA0KYGBge3J9DQp4bCA8LSBzZXR4KHosIGRlcHRoID0gcXVhbnRpbGUocXVha2VzJGRlcHRoLCAuMjUpKQ0KeGwxIDwtIHNldHgoeiwgZGVwdGggPSBxdWFudGlsZShxdWFrZXMkZGVwdGgsIC4yNSkgKyBzZChxdWFrZXMkZGVwdGgpKQ0Kc2wgPC0gc2ltKHosIHggPSB4bCwgeDEgPSB4bDEpDQpmZGwgPC0gc2wkZ2V0X3FpKHh2YWx1ZT0ieDEiLCBxaT0iZmQiKQ0KDQp4bSA8LSBzZXR4KHosIGRlcHRoID0gcXVhbnRpbGUocXVha2VzJGRlcHRoLCAuNTApKQ0KeG0xIDwtIHNldHgoeiwgZGVwdGggPSBxdWFudGlsZShxdWFrZXMkZGVwdGgsIC41MCkgKyBzZChxdWFrZXMkZGVwdGgpKQ0Kc20gPC0gc2ltKHosIHggPSB4bSwgeDEgPSB4bTEpDQpmZG0gPC0gc20kZ2V0X3FpKHh2YWx1ZT0ieDEiLCBxaT0iZmQiKQ0KDQp4aCA8LSBzZXR4KHosIGRlcHRoID0gcXVhbnRpbGUocXVha2VzJGRlcHRoLCAuNzUpKQ0KeGgxIDwtIHNldHgoeiwgZGVwdGggPSBxdWFudGlsZShxdWFrZXMkZGVwdGgsIC43NSkgKyBzZChxdWFrZXMkZGVwdGgpKQ0Kc2ggPC0gc2ltKHosIHggPSB4aCwgeDEgPSB4aDEpDQpmZGggPC0gc2gkZ2V0X3FpKHh2YWx1ZT0ieDEiLCBxaT0iZmQiKQ0KYGBgDQoNCiMjIENvbWJpbmluZyB0aGUgcmVzdWx0cyANCmBgYHtyfQ0KZCA8LSBhcy5kYXRhLmZyYW1lKGNiaW5kKGZkbCwgZmRtLCBmZGgpKQ0KaGVhZChkKQ0KYGBgDQoNCmBgYHtyfQ0KdGlkZCA8LSBkICU+JSANCiAgZ2F0aGVyKHF1YW50aWxlLCBkaWYsIDE6MykNCmhlYWQodGlkZCkNCmBgYA0KDQpgYGB7cn0NCnRpZGQgJT4lIA0KICBncm91cF9ieShxdWFudGlsZSkgJT4lIA0KICBzdW1tYXJpc2UobWVhbiA9IG1lYW4oZGlmKSwgc2QgPSBzZChkaWYpKQ0KYGBgDQoNCiMjIEdHIFBsb3QgDQpgYGB7cn0NCmxpYnJhcnkoZ2dwbG90MikNCg0KZ2dwbG90KHRpZGQsIGFlcyhkaWYpKSArIGdlb21faGlzdG9ncmFtKCkgKyBmYWNldF9ncmlkKH4gcXVhbnRpbGUpDQpgYGANCg0KIyMgT3ZlcmFsbCBUaG91Z2h0cyANCipXaGF0IHdlIGhhdmUgbGVhcm5lZCBvdmVyYWxsOiBCYXNlZCBvbiB0aGlzIGRhdGEgc2V0LCB3ZSBsZWFybmVkLCB0aGF0IHdoZW4gaXQgY29tZXMgdG8gZWFydGhxdWFrZXMsIHN0YXRpb25zIGFyZSBtb3JlIGludGVyZXN0ZWQgaW4gcmVwb3J0aW5nIG9uIGVhcnRocXVha2VzIHRoYXQgaGF2ZSBhIGxhcmdlciBkZXB0aCwgYXMgd2VsbCBhcyBhbiBvdmVyYWxsIGxhcmdlciBtYWduaXR1ZGUuKg0K