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