تحلیل داده‌های زلزله با R-بخش اول
تحلیل داده‌های زلزله با R - بخش اول
در این بخش از سری آموزش‌های تحلیل داده‌های زلزله، قصد داریم تا به معرفی تجزیه و تحلیل اکتشافی یک مجموعه داده بپردازیم. همچنین به گزارش وقوع زلزله و حوادث مشابه در یک پنجره زمانی خاص را نشان می‌دهیم همچنین یک تحلیل اکتشافی اساسی و برخی از اصول و ابزار‌های اصلی را در این مورد بررسی میکنیم. سپس نحوه نشان دادن محتوای مجموعه داده‌های زمین لرزه‌ها را بروی نقشه جغرافیایی به همراه وضعیت و نمودارهای تعاملی و متحرک نشان خواهیم داد. سرانجام، تجزیه و تحلیل خوشه‌ای از زمین لرزه واقع در مناطق جغرافیایی خاص ارائه خواهد شد.
- بخش اول : تحلیل اکتشافی متغیرهای کمی
- بخش دوم : تحلیل اکتشافی متغیرهای طبقه‌بندی شده
- بخش سوم : نقشه‌های استاتیک، تعاملی و متحرک
- بخش چهارم : تجزیه و تحلیل خوشه‌ای
همه انالیزهای گفته شده با پکیج‌های زیر قابل انجام هست.
library(ggplot2) library(ggmap) library(dplyr) library(factoextra) library(Hmisc) library(gridExtra) library(lubridate) library(corrplot) library(vcd) library(gmodels)
install.packages(c("ggplot2",
"ggmap",
"dplyr",
"factoextra",
"Hmisc",
"gridExtra",
"lubridate",
"corrplot",
"vcd",
"gmodels"))
دریافت داده:
داده‌های زمین لرزه را از سایت لرزه نگاری usgs دریافت میکنیم مخصوصا داده‌های ۳۰ روز اخیر. توجه داشته باشید که بعضی مجموعه داده‌ها روز به روز ذخیره شده و در مجموع بازه ۳۰ روزه را پوشش می‌دهد. که در حال حاضر با این نوع داده‌ها کار نمی‌کنیم. داده‌ها را دانلود میکنیم و برای استفاده بارگذاری می‌نماییم. توجه کنید که داده‌های زمین لرزه متغیرهای محلی هستند.
با استفاده از دستور زیر داده‌ها را مستقیما از طریق نرم افزار دانلود و در نرم افزار لود میکند توجه کنید که قبل از انجام این کار ابتدا پکیج‌های بالا نصب و فراخوانی کنید.
if("all_week.csv" %in% dir(".") == FALSE) {
url="https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.csv"
download.file(url = url, destfile = "all_week.csv")
}
quakes = read.csv("all_week.csv", header=TRUE, sep=',', stringsAsFactors = FALSE)
با استفاده از دستور زیر میتوان ساختار dataframe را برحسب متغیرهای factor و time تغییر داد تا در ادامه دچار مشکل نشویم.
quakes$time = ymd_hms(quakes$time) quakes$updated = ymd_hms(quakes$updated) quakes$magType = as.factor(quakes$magType) quakes$net = as.factor(quakes$net) quakes$type = as.factor(quakes$type) quakes$status = as.factor(quakes$status) quakes$locationSource = as.factor(quakes$locationSource) quakes$magSource = as.factor(quakes$magSource)
با استفاده از دستور زیر می‌توانید تعداد سطر و ستون مجموعه داده خود را مشاهده کنید (ابعاد داده را نشان می‌دهد)
dim(quakes) ## [1] 8407 22
با دستور زیر می‌توانیم داده ها را طوری بازچینی بکنیم که ترتیب سطر ها به شکل نزولی باشد.
quakes = arrange(quakes, -row_number())
با دستور زیر می‌توانیم تعداد سطرهایی که داده کامل دارند(داده گمشده در بین آنها نیست) را مشاهده کنیم. توجه کنید که داده‌های شما ممکن هست تعداد متفاوتی را نمایش بدهد بدلیل اینکه تعداد زلزله ها هر لحظه بروزرسانی می‌شود.
sum(complete.cases(quakes)) ## [1] 3957
کسر زیر درصد تعداد سطرهای کامل و بدون داده گمشده را نشان می‌دهد.
sum(complete.cases(quakes))/nrow(quakes) ## [1] 0.4706792
برای مشاهده نگاه کلی به مجموعه داده‌ها از دستورات زیر می‌توانید استفاده کنید
head(quakes) ## time latitude longitude depth mag magType nst gap ## 1 2019-02-14 16:35:36 53.1050 -164.7006 20.20 2.30 ml NA NA ## 2 2019-02-14 16:42:10 33.5345 -116.7100 2.46 1.52 ml 37 35 ## 3 2019-02-14 16:52:34 19.2010 -155.5120 34.29 2.34 ml 57 81 ## 4 2019-02-14 17:00:44 60.5937 -147.8132 13.20 2.10 ml NA NA ## 5 2019-02-14 17:02:05 64.3332 -148.4033 11.50 0.90 ml NA NA ## 6 2019-02-14 17:06:09 65.4253 -153.3378 14.10 1.00 ml NA NA ## dmin rms net id updated ## 1 NA 0.30 ak ak01922ox4si 2019-03-01 21:51:59 ## 2 0.03833 0.20 ci ci37531042 2019-02-14 17:00:40 ## 3 0.03755 0.09 hv hv70811146 2019-02-14 19:21:24 ## 4 NA 0.42 ak ak01922pb2w7 2019-03-01 21:51:36 ## 5 NA 0.65 ak ak01922pbfay 2019-03-01 21:51:37 ## 6 NA 0.83 ak ak01922pcak4 2019-03-01 21:52:00 ## place type horizontalError depthError ## 1 134km SSE of Akutan, Alaska earthquake NA 27.00 ## 2 4km SW of Anza, CA earthquake 0.22 0.26 ## 3 3km W of Pahala, Hawaii earthquake 0.42 0.53 ## 4 51km ESE of Whittier, Alaska earthquake NA 0.20 ## 5 43km SE of North Nenana, Alaska earthquake NA 0.40 ## 6 65km WNW of Tanana, Alaska earthquake NA 0.10 ## magError magNst status locationSource magSource ## 1 NA NA reviewed ak ak ## 2 0.159 27 reviewed ci ci ## 3 0.177 21 reviewed hv hv ## 4 NA NA reviewed ak ak ## 5 NA NA reviewed ak ak ## 6 NA NA reviewed ak ak tail(quakes) # run it to show results
برای مشاهده ساختار هر متغیر تشکیل دهنده مجموعه داده و اینکه این متغیرها از چه نوعی هستند وسایر اطلاعات کلی از دستور زیر استفاده کنید.
str(quakes) ## 'data.frame': 8407 obs. of 22 variables: ## $ time : POSIXct, format: "2019-02-14 16:35:36" "2019-02-14 16:42:10" ... ## $ latitude : num 53.1 33.5 19.2 60.6 64.3 ... ## $ longitude : num -165 -117 -156 -148 -148 ... ## $ depth : num 20.2 2.46 34.29 13.2 11.5 ... ## $ mag : num 2.3 1.52 2.34 2.1 0.9 1 0.3 0.4 2.4 1.82 ... ## $ magType : Factor w/ 10 levels "","mb","mb_lg",..: 6 6 6 6 6 6 6 6 6 4 ... ## $ nst : int NA 37 57 NA NA NA NA NA NA 5 ... ## $ gap : num NA 35 81 NA NA NA NA NA NA 141 ... ## $ dmin : num NA 0.0383 0.0376 NA NA ... ## $ rms : num 0.3 0.2 0.09 0.42 0.65 0.83 0.48 0.7 0.55 0.09 ... ## $ net : Factor w/ 14 levels "ak","ci","hv",..: 1 2 3 1 1 1 1 1 1 13 ... ## $ id : chr "ak01922ox4si" "ci37531042" "hv70811146" "ak01922pb2w7" ... ## $ updated : POSIXct, format: "2019-03-01 21:51:59" "2019-02-14 17:00:40" ... ## $ place : chr "134km SSE of Akutan, Alaska" "4km SW of Anza, CA" "3km W of Pahala, Hawaii" "51km ESE of Whittier, Alaska" ... ## $ type : Factor w/ 7 levels "chemical explosion",..: 2 2 2 2 2 2 2 2 2 2 ... ## $ horizontalError: num NA 0.22 0.42 NA NA NA NA NA NA 1.06 ... ## $ depthError : num 27 0.26 0.53 0.2 0.4 0.1 0.4 0.5 0.6 3.08 ... ## $ magError : num NA 0.159 0.177 NA NA NA NA NA NA 0.244 ... ## $ magNst : int NA 27 21 NA NA NA NA NA NA 4 ... ## $ status : Factor w/ 2 levels "automatic","reviewed": 2 2 2 2 2 2 2 2 2 2 ... ## $ locationSource : Factor w/ 15 levels "ak","ci","hv",..: 1 2 3 1 1 1 1 1 1 14 ... ## $ magSource : Factor w/ 15 levels "ak","ci","hv",..: 1 2 3 1 1 1 1 1 1 14 ...
توضیحات مختصری درباره متغیرهای تشکیل دهنده این مجموعه داده در زیر بیان میکنیم که هر متغیر بیانگر چه ویژگی می‌باشد.
time
latitude
longitude
depth
mag
magType
nst
gap
dmin
rms
net
id
update
place
type
horizontalError
depthError
magError
magNst
status
locationSource
magSource
– زمانی که حادثه (زلزله) رخداده است.
– بیانگر عرض جغرافیایی می‌باشد که دامنه تغییر ان از -90 درجه تا +90 درجه می‌باشد.
– بیانگر طول جغرافیایی می‌باشد که دامنه تغییر ان از -180 درجه تا +180 درجه می‌باشد.
– عمق رخداد زمین لرزه برحسب کیلومتر.
– بزرگای زلزله از دامنه 1.0- تا 10.0 می‌باشد.
– روش و یا الگوریتم مورد استفاده برای محاسبه بزرگای هر رخداد.
– تعداد ایستگاه‌های لرزه شناسی برای تعیین مکان زلزله
– بزرگترین زاویه در بین زوایای ایستگاه‌های گیرنده زلزله
– فاصله افقی از مرکز زلزله را به نزدیکترین ایستگاه
– مربع میانگین ریشه (RMS) باقی‌مانده زمان پیمودن، برحسب ثانیه، با استفاده از تمام وزن‌ها
– شناسه کمک کننده داده، مشخصه شبکه‌ای که با استفاده از داده‌های آن این زمین لرزه و اطلاعات آن به ثبت رسیده است.
– یک مشخصه یکتا برای رخداد، این شناسه مورد نظر در حال حاضر برای این ویژگی هست که ممکن است در طول زمان تغییر کند.
– زمانی که رویداد اخیرا به روز شده است.
– توضیحات متنی از نام منطقه جغرافیایی نزدیک به رخداد. که ممکن هست اسم شهر باشد.
– نوع رخداد لرزه‌ای.
– عدم اطمینان از محل گزارش شده از این رویداد برحسب کیلومتر.
– عدم اطمینان از عمق گزارش شده از این رویداد برحسب کیلومتر.
– عدم اطمینان از بزرگای گزارش شده از این رویداد .
– مجموع ایستگاه‌های لرزه‌نگاری مورد استفاده شده برای محاسبه بزرگای این رخداد زلزله.
– مشخص کننده اینکه آیا این رویداد به‌وسیله نیروی انسانی بازبینی شده است.
– شبکه ای که در اصل نویسنده محل گزارش شده از این رویداد است.
– شبکه ای که در اصل نویسنده بزرگای گزارش شده از این رویداد است.
با دستورات زیر می‌توانید فاکتورهای ذخیره شده در متغیر type را مشاهده کنید و سطح‌بندی نوع زلزله را تعیین کنید.
sort(unique(quakes$type)) ## [1] chemical explosion earthquake explosion ## [4] ice quake other event quarry blast ## [7] rock burst ## 7 Levels: chemical explosion earthquake explosion ... rock burst
برای مشاهده مدت زمان مجموعه داده خود می‌توانید از دستورات زیر استفاده کنید. بدین منظور ابتدا تعداد سطرهای مجموعه داده خود را در یک متغیر ذخیره کرده و سپس از طریق فراخوانی داده متغیر زمان را ایندکس اخرین سطر که قبلا در متغیر nr ذخیره کرده بودیم را چاپ میکنیم.
nr = nrow(quakes) c(quakes$time[1], quakes$time[nr]) ## [1] "2019-02-14 16:35:36 UTC" "2019-03-16 16:06:56 UTC"
تجزیه و تحلیل اکتشافی
در این بخش به آنالیز آماری تک تک متغیرها می‌پردازیم و استباط های آماری را انجام می‌دهیم. در زیر نگاه مختصری به نوع متغیرهای می‌اندازیم.
(numeric_vars = names(which(sapply(quakes, class) == "numeric"))) ## [1] "latitude" "longitude" "depth" ## [4] "mag" "gap" "dmin" ## [7] "rms" "horizontalError" "depthError" ## [10] "magError" (integer_vars = names(which(sapply(quakes, class) == "integer"))) ## [1] "nst" "magNst" (factor_vars = names(which(sapply(quakes, class) == "factor"))) ## [1] "magType" "net" "type" "status" ## [5] "locationSource" "magSource" (character_vars = names(which(sapply(quakes, class) == "character"))) ## [1] "id" "place" setdiff(names(quakes), c(numeric_vars, integer_vars, factor_vars, character_vars)) ## [1] "time" "updated"
متغیرهای کمی
متغیرهای کمی شامل عددی (پیوسته) و عدد صحیح ( گسسته) می‌باشند که قبلا معرفی کرده‌ایم. خلاصه کوتاهی را برای هر یک در زیر نشان می‌دهیم
quantitative_vars = c(numeric_vars, integer_vars) sapply(quakes[, quantitative_vars], summary) ## $latitude ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -59.96 33.84 38.82 40.94 60.03 72.32 ## ## $longitude ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -180.0 -149.3 -120.6 -113.1 -115.7 180.0 ## ## $depth ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -3.45 3.00 8.58 23.21 20.10 619.21 ## ## $mag ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## -0.710 0.900 1.400 1.651 2.000 7.500 2 ## ## $gap ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 12.0 69.0 103.0 122.3 158.6 353.0 2472 ## ## $dmin ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.0003 0.0219 0.0638 0.5698 0.2102 32.4620 2659 ## ## $rms ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.0000 0.1000 0.2000 0.3386 0.5500 1.8400 2 ## ## $horizontalError ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.080 0.250 0.440 1.865 1.290 27.650 3073 ## ## $depthError ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 0.380 0.700 2.934 1.900 745.500 ## ## $magError ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.0000 0.0960 0.1500 0.1835 0.2280 5.3600 2793 ## ## $nst ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 3.00 8.00 14.00 18.85 24.00 199.00 3368 ## ## $magNst ## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ## 0.00 4.00 8.00 15.35 17.00 415.00 2488
تابع describe() در پکیج Hmisc توضیحات خوبی را در این زمینه ارائه می‌دهد.
describe(quakes[, quantitative_vars]) ## quakes[, quantitative_vars] ## ## 12 Variables 8407 Observations ## --------------------------------------------------------------------------- ## latitude ## n missing distinct Info Mean Gmd .05 .10 ## 8407 0 7416 1 40.94 19.67 13.49 19.16 ## .25 .50 .75 .90 .95 ## 33.84 38.82 60.03 63.08 64.98 ## ## lowest : -59.9641 -59.9052 -59.8885 -59.8790 -59.8610 ## highest: 71.7325 71.7340 71.7684 71.8899 72.3184 ## --------------------------------------------------------------------------- ## longitude ## n missing distinct Info Mean Gmd .05 .10 ## 8407 0 7617 1 -113.1 47.68 -156.15 -152.83 ## .25 .50 .75 .90 .95 ## -149.33 -120.56 -115.67 -67.39 22.34 ## ## lowest : -179.9963 -179.8762 -179.8751 -179.8321 -179.7528 ## highest: 179.0796 179.1853 179.3115 179.8629 179.9946 ## --------------------------------------------------------------------------- ## depth ## n missing distinct Info Mean Gmd .05 .10 ## 8407 0 2755 1 23.21 32.26 0.00 0.97 ## .25 .50 .75 .90 .95 ## 3.00 8.58 20.10 59.00 98.84 ## ## lowest : -3.45 -3.44 -3.41 -3.40 -3.39, highest: 610.04 612.77 614.79 616.79 619.21 ## --------------------------------------------------------------------------- ## mag ## n missing distinct Info Mean Gmd .05 .10 ## 8405 2 412 1 1.651 1.225 0.270 0.470 ## .25 .50 .75 .90 .95 ## 0.900 1.400 2.000 3.216 4.500 ## ## lowest : -0.71 -0.64 -0.50 -0.40 -0.37, highest: 6.20 6.30 6.40 7.00 7.50 ## --------------------------------------------------------------------------- ## gap ## n missing distinct Info Mean Gmd .05 .10 ## 5935 2472 919 1 122.3 77.88 38.0 48.0 ## .25 .50 .75 .90 .95 ## 69.0 103.0 158.6 224.9 285.9 ## ## lowest : 12.00 13.00 14.00 15.00 16.00, highest: 349.33 349.73 350.00 352.00 353.00 ## --------------------------------------------------------------------------- ## dmin ## n missing distinct Info Mean Gmd .05 .10 ## 5748 2659 4519 1 0.5698 0.9607 0.005299 0.008580 ## .25 .50 .75 .90 .95 ## 0.021943 0.063785 0.210200 1.449100 3.035300 ## ## lowest : 0.0002619 0.0003986 0.0004850 0.0007119 0.0007695 ## highest: 24.2990000 26.6220000 30.2310000 31.8220000 32.4620000 ## --------------------------------------------------------------------------- ## rms ## n missing distinct Info Mean Gmd .05 .10 ## 8405 2 676 1 0.3386 0.3293 0.03 0.04 ## .25 .50 .75 .90 .95 ## 0.10 0.20 0.55 0.79 0.95 ## ## lowest : 0.0000 0.0032 0.0046 0.0090 0.0098, highest: 1.6000 1.7700 1.7800 1.8200 1.8400 ## --------------------------------------------------------------------------- ## horizontalError ## n missing distinct Info Mean Gmd .05 .10 ## 5334 3073 529 1 1.865 2.616 0.15 0.18 ## .25 .50 .75 .90 .95 ## 0.25 0.44 1.29 7.10 9.40 ## ## lowest : 0.08 0.09 0.10 0.11 0.12, highest: 23.60 24.76 26.58 27.60 27.65 ## --------------------------------------------------------------------------- ## depthError ## n missing distinct Info Mean Gmd .05 .10 ## 8407 0 756 0.999 2.934 4.387 0.180 0.200 ## .25 .50 .75 .90 .95 ## 0.380 0.700 1.900 6.894 12.470 ## ## Value 0 10 20 30 40 80 90 120 180 260 ## Frequency 7304 751 105 234 6 1 1 1 1 1 ## Proportion 0.869 0.089 0.012 0.028 0.001 0.000 0.000 0.000 0.000 0.000 ## ## Value 660 750 ## Frequency 1 1 ## Proportion 0.000 0.000 ## --------------------------------------------------------------------------- ## magError ## n missing distinct Info Mean Gmd .05 .10 ## 5614 2793 664 1 0.1835 0.1391 0.037 0.056 ## .25 .50 .75 .90 .95 ## 0.096 0.150 0.228 0.350 0.430 ## ## lowest : 0.000000000 0.001000000 0.002000000 0.003000000 0.004180669 ## highest: 3.670000000 3.950000000 4.140000000 4.990000000 5.360000000 ## --------------------------------------------------------------------------- ## nst ## n missing distinct Info Mean Gmd .05 .10 ## 5039 3368 104 0.999 18.85 15.39 4 6 ## .25 .50 .75 .90 .95 ## 8 14 24 39 49 ## ## lowest : 3 4 5 6 7, highest: 126 129 144 145 199 ## --------------------------------------------------------------------------- ## magNst ## n missing distinct Info Mean Gmd .05 .10 ## 5919 2488 171 0.997 15.35 18.02 1 2 ## .25 .50 .75 .90 .95 ## 4 8 17 30 49 ## ## lowest : 0 1 2 3 4, highest: 355 367 370 371 415 ## ---------------------------------------------------------------------------
یک تابع کمکی برای قطعه‌بندی با استفاده از پکیج Kmisc که در CRAN نگهداری می‌شود را در زیر آورده‌ایم
## from Kmisc
chunk = function(min, max, size, by=1) {
if( missing(max) ) {
max = min
min = 1
}
n = ceiling( (max-min) / (size*by) )
out = vector("list", n)
lower = min
upper = by * (min + size - 1)
for( i in 1:length(out) ) {
out[[i]] = seq(lower, min(upper, max), by=by)
lower = lower + size * by
upper = upper + size * by
}
names(out) = paste(sep="", "chunk", 1:length(out))
return(out)
}
نمودار جعبه‌ای برای هر متغیر کمی در زیر نشان داده شده است.
gp = invisible(lapply(quantitative_vars, function(x) {
ggplot(data = quakes, aes(y=eval(parse(text=x)))) + geom_boxplot(fill = "palegreen4") + ylab(x) + ggtitle(x) + theme(legend.position="none") + coord_flip()}))
grob_plots = invisible(lapply(chunk(1, length(gp), 4), function(x) {
marrangeGrob(grobs=lapply(gp[x], ggplotGrob), nrow=2, ncol=2)}))
grob_plots$chunk1
grob_plots$chunk2
grob_plots$chunk3
نمودار چگالی برای تصویر سازی اینکه تراکم داده در کدام نقطه بیشتر هست کمک میکند. با استفاده کد زیر می‌توانید نمودارهای مشابه را رسم کنید که ما فقط یک نمونه را می‌آوریم.
gp = invisible(lapply(quantitative_vars, function(x) {
ggplot(data = quakes, aes(x=eval(parse(text=x)))) + geom_density(fill = "palegreen4") + ylab(x) + xlab("") + ggtitle(x) + theme(legend.position="none")}))
grob_plots = invisible(lapply(chunk(1, length(gp), 4), function(x) {
marrangeGrob(grobs=lapply(gp[x], ggplotGrob), nrow=2, ncol=2)}))
grob_plots$chunk1
grob_plots$chunk2
grob_plots$chunk3
حال نمودار جعبه‌ای که در بالا رسم شد را برحسب نوع رخداد هرکدام رسم میکنیم
فقط یک نمونه از نمودار را به عنوان نمونه می‌آوریم
gp = invisible(lapply(quantitative_vars, function(x) {
ggplot(data = quakes, aes(x = type, y = eval(parse(text=x)))) + geom_boxplot(fill = "palegreen4") + ylab(x) + ggtitle(x) + theme(legend.position="none") + coord_flip()}))
grob_plots = invisible(lapply(chunk(1, length(gp), 4), function(x) {
marrangeGrob(grobs=lapply(gp[x], ggplotGrob), nrow=2, ncol=2)}))
grob_plots$chunk1
grob_plots$chunk1
grob_plots$chunk2
grob_plots$chunk3
به‌همین ترتیب می‌توانید تابع چگالی را نیز برحسب نوع رخداد رسم کنید
فقط یک نمونه از نمودار را به عنوان نمونه می‌آوریم
gp = invisible(lapply(quantitative_vars, function(x) {
ggplot(data = quakes, aes(x=eval(parse(text=x)))) + geom_density(fill = "palegreen4") + ylab(x) + xlab("") + ggtitle(x) + theme(legend.position="none") + facet_wrap( ~ type)}))
grob_plots = invisible(lapply(chunk(1, length(gp), 4), function(x) {
marrangeGrob(grobs=lapply(gp[x], ggplotGrob), nrow=2, ncol=2)}))
grob_plots$chunk1
grob_plots$chunk1
grob_plots$chunk2
grob_plots$chunk3
برای تحقیقات بهتر و نمایش کامل اطلاعات میتوان از نمودار پراکنش استفاده کرد.
با استفاده از کد زیر این نمودار را می‌کشیم.
earthquakes = quakes %>% filter(type == "earthquake") pairs(earthquakes %>% select(quantitative_vars)) hist_data = earthquakes %>% select(mag, depth) hist_data = hist_data[complete.cases(hist_data),] p = ggplot(data = hist_data, aes(x = depth, y = mag)) + geom_point() + guides(fill=FALSE) p
بالای ۳۰۰ کیلومتر از این داده‌ها در بزرگای بیشتر از ۴ اتفاق افتاده‌اند، برای مشاهده روابط بین کیلومتر و بزرگای زلزله از مدل خطی زیر استفاده میکنیم.
با استفاده از کد زیر این نمودار را می‌کشیم.
mag_depth_lr = lm(mag ~ poly(depth, 3), data = hist_data) summary(mag_depth_lr) ## ## Call: ## lm(formula = mag ~ poly(depth, 3), data = hist_data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2.1845 -0.7333 -0.2310 0.4361 4.7112 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.65543 0.01177 140.625 < 2e-16 *** ## poly(depth, 3)1 40.93311 1.06794 38.329 < 2e-16 *** ## poly(depth, 3)2 -18.69337 1.06794 -17.504 < 2e-16 *** ## poly(depth, 3)3 6.64030 1.06794 6.218 5.29e-10 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.068 on 8226 degrees of freedom ## Multiple R-squared: 0.1807, Adjusted R-squared: 0.1804 ## F-statistic: 604.7 on 3 and 8226 DF, p-value: less than 2.2e-16
تمامی ضرایب رگرسیونی معنی‌دار هستند در زیر نمودار مقادیر اصلی در مقابل مدل برازش داده شده را رسم میکنیم.
تحلیل همبستگی
با استفاده از کد زیر همبستگی محاسبه شده و در جدول زیر آورده شده است.
(cor_mat = cor(quakes[, quantitative_vars], use = "pairwise.complete.obs")) ## latitude longitude depth mag ## latitude 1.00000000 -0.52875417 -0.14867697 -0.38540388 ## longitude -0.52875417 1.00000000 0.14551789 0.55605921 ## depth -0.14867697 0.14551789 1.00000000 0.38277138 ## mag -0.38540388 0.55605921 0.38277138 1.00000000 ## gap -0.03794518 0.02975063 0.05196967 0.02808442 ## dmin -0.47773533 0.44108043 0.30279438 0.57339025 ## rms 0.13784391 0.19709935 0.33147665 0.60872325 ## horizontalError -0.53481018 0.61381314 0.54880413 0.76651963 ## depthError -0.06707420 0.06460787 0.03415940 0.08011669 ## magError 0.02470576 -0.09122383 -0.06033724 -0.10270236 ## nst -0.11490376 -0.27741899 -0.04172062 0.28344750 ## magNst -0.19571866 0.26728165 0.26750233 0.46199100 ## gap dmin rms horizontalError ## latitude -0.03794518 -0.47773533 0.13784391 -0.53481018 ## longitude 0.02975063 0.44108043 0.19709935 0.61381314 ## depth 0.05196967 0.30279438 0.33147665 0.54880413 ## mag 0.02808442 0.57339025 0.60872325 0.76651963 ## gap 1.00000000 0.02719460 0.01882822 0.24321328 ## dmin 0.02719460 1.00000000 0.53628845 0.65243264 ## rms 0.01882822 0.53628845 1.00000000 0.75539036 ## horizontalError 0.24321328 0.65243264 0.75539036 1.00000000 ## depthError 0.32860091 0.08272005 0.03538935 0.26089694 ## magError 0.10447472 -0.10511159 -0.12456313 -0.09591121 ## nst -0.49361171 -0.19298946 0.06461759 -0.24813050 ## magNst -0.18564633 0.28321115 0.42430650 0.29769313 ## depthError magError nst magNst ## latitude -0.06707420 0.02470576 -0.11490376 -0.19571866 ## longitude 0.06460787 -0.09122383 -0.27741899 0.26728165 ## depth 0.03415940 -0.06033724 -0.04172062 0.26750233 ## mag 0.08011669 -0.10270236 0.28344750 0.46199100 ## gap 0.32860091 0.10447472 -0.49361171 -0.18564633 ## dmin 0.08272005 -0.10511159 -0.19298946 0.28321115 ## rms 0.03538935 -0.12456313 0.06461759 0.42430650 ## horizontalError 0.26089694 -0.09591121 -0.24813050 0.29769313 ## depthError 1.00000000 -0.01048942 -0.20439173 -0.01479938 ## magError -0.01048942 1.00000000 -0.07789094 -0.14474547 ## nst -0.20439173 -0.07789094 1.00000000 0.72021297 ## magNst -0.01479938 -0.14474547 0.72021297 1.00000000
با استفاده از نمودارهای زیر همبستگی بین داده ها به‌طور قابل فهم تری نمایش داده می‌شود.
نقشه
با استفاده از دستور زیر نقشه زلزله های رخداده در فاصله زمانی داده های برداشت شده را می‌توانیم رسم کنیم.
world = map_data('world')
title = paste("Earthquake map from ", paste(earthquakes$time[1], earthquakes$time[nrow(earthquakes)], sep = " to "))
p = ggplot() + geom_map(data = world, map = world, aes(x = long, y=lat, group=group, map_id=region), fill="white", colour="#7f7f7f", size=0.5)
p = p + geom_point(data = earthquakes, aes(x=longitude, y = latitude, colour = mag)) + scale_colour_gradient(low = "#00AA00",high = "#FF00AA") + ggtitle(title)
p
