Rで解析:データの分布をグラフで確認「ggridges」パッケージ

Rの解析に役に立つ記事
スポンサーリンク

データの分布を確認するのに便利なパッケージの紹介です。気象庁から22.05.03の各都道府県の最高気温を0:00から23:00まで1時間毎に取得し、joyplotを作成するコマンドも紹介します。

パッケージバージョンは0.5.3。実行コマンドはwindows 11のR version 4.1.2で確認しています。

スポンサーリンク

パッケージのインストール

下記コマンドを実行してください。

#パッケージのインストール
install.packages("ggridges")

コマンドの紹介

詳細はコマンド、パッケージのヘルプを確認してください。

#パッケージの読み込み:libraryコマンド
library("ggridges")

###データ例の作成#####
#tidyverseパッケージがなければインストール
if(!require("tidyverse", quietly = TRUE)){
  install.packages("tidyverse");require("tidyverse")
}
set.seed(1234)
n <- 500
TestData <- data.frame(Group = sample(paste0("Group", 1:5), n, replace = TRUE),
                          Time = sample(1:5, n, replace = TRUE),
                          height = rnorm(n),
                          Data2 = rnorm(n) + rnorm(n) + rnorm(n))
#geom_ridgeline&#29992;&#12398;&#12487;&#12540;&#12479;
TestRidgeLine <- TestData %>% group_by(Group, Time) %>%
  summarise_all(lst(mean)) %>% mutate(Yposi = recode(Group, Group1 = 0,
                                                     Group2 = .3, Group3 = .5,
                                                     Group4 = .7, Group5 = .9))
#######

#&#39640;&#12373;&#12434;&#25351;&#23450;&#12375;&#12383;&#12456;&#12522;&#12450;&#12503;&#12525;&#12483;&#12488;&#12434;&#20316;&#25104;:geom_ridgeline&#12467;&#12510;&#12531;&#12489;
#0&#20197;&#19979;&#12398;&#39640;&#12373;&#34920;&#31034;&#12377;&#12427;&#31684;&#22258;&#12434;&#25351;&#23450;:min_height&#12458;&#12503;&#12471;&#12519;&#12531;
ggplot(TestRidgeLine, aes(x = Time, y = Yposi,
                          height = height_mean, group = Yposi,
                          fill = Group)) +
  geom_ridgeline(show.legend = F, alpha = .5,
                 min_height = min(TestRidgeLine[, 3]))

#&#12487;&#12540;&#12479;&#20998;&#24067;&#12434;&#12503;&#12525;&#12483;&#12488;:geom_density_ridges&#12467;&#12510;&#12531;&#12489;
#&#12464;&#12521;&#12501;&#19979;&#37096;&#32218;&#20184;&#12365;&#12503;&#12525;&#12483;&#12488;:geom_density_ridges2&#12467;&#12510;&#12531;&#12489;
ggplot(TestData, aes(x = height, y = Group, fill = Group)) +
  geom_density_ridges2(scale = 1) + facet_wrap(~Group)

#&#22615;&#33394;&#12434;&#25351;&#23450;:scale_fill_cyclical&#12467;&#12510;&#12531;&#12489;
#&#26528;&#32218;&#12434;&#25351;&#23450;:scale_color_cyclical&#12467;&#12510;&#12531;&#12489;
ggplot(TestData, aes(x = height, y = Group, fill = Group, color = Group)) +
  geom_density_ridges2(scale = 1, size = 1.5) + facet_wrap(~Group) +
  scale_fill_cyclical(values = c("blue", "green", "yellow")) +
  scale_color_cyclical(values = c("red", "black"))

出力例

・geom_ridgelineコマンド

・geom_density_ridgesコマンド

・scale_fill_cyclicalコマンド

気象庁から最高気温を取得

最高気温はNewMaxTempに格納しています。

#&#12300;tidyverse&#12301;&#12497;&#12483;&#12465;&#12540;&#12472;&#12434;&#35501;&#12415;&#36796;&#12415;
if(!require("tidyverse", quietly = TRUE)){
  install.packages("tidyverse");require("tidyverse")
}
#&#37117;&#36947;&#24220;&#30476;&#12434;&#28310;&#20633;
JpanPref <- c("&#21271;&#28023;&#36947;", "&#38738;&#26862;&#30476;", "&#23721;&#25163;&#30476;", "&#23470;&#22478;&#30476;", "&#31119;&#23798;&#30476;", "&#33576;&#22478;&#30476;", "&#21315;&#33865;&#30476;",
              "&#31179;&#30000;&#30476;", "&#23665;&#24418;&#30476;", "&#26032;&#28511;&#30476;", "&#26627;&#26408;&#30476;", "&#22524;&#29577;&#30476;", "&#26481;&#20140;&#37117;", "&#32676;&#39340;&#30476;",
              "&#23665;&#26792;&#30476;", "&#31070;&#22856;&#24029;&#30476;", "&#23500;&#23665;&#30476;", "&#38263;&#37326;&#30476;", "&#38745;&#23713;&#30476;", "&#30707;&#24029;&#30476;", "&#31119;&#20117;&#30476;",
              "&#23696;&#38428;&#30476;", "&#24859;&#30693;&#30476;", "&#28363;&#36032;&#30476;", "&#19977;&#37325;&#30476;", "&#20140;&#37117;&#24220;", "&#22856;&#33391;&#30476;", "&#21644;&#27468;&#23665;&#30476;",
              "&#20853;&#24235;&#30476;", "&#22823;&#38442;&#24220;", "&#40165;&#21462;&#30476;", "&#23713;&#23665;&#30476;", "&#23798;&#26681;&#30476;", "&#24195;&#23798;&#30476;", "&#39321;&#24029;&#30476;",
              "&#24499;&#23798;&#30476;", "&#24859;&#23195;&#30476;", "&#39640;&#30693;&#30476;", "&#23665;&#21475;&#30476;", "&#31119;&#23713;&#30476;", "&#22823;&#20998;&#30476;", "&#23470;&#23822;&#30476;",
              "&#20304;&#36032;&#30476;", "&#29066;&#26412;&#30476;", "&#40575;&#20816;&#23798;&#30476;", "&#38263;&#23822;&#30476;", "&#27798;&#32260;&#30476;")

#&#26178;&#38291;&#25991;&#23383;&#21015;&#12434;&#20316;&#25104;
Hour <- paste0(formatC(0:23, width = 2, flag = "0"), "00")

#&#12487;&#12540;&#12479;&#20445;&#31649;&#29992;&#22793;&#25968;
NewMaxTemp <- data.frame()

for(i in seq(Hour)){
  ###&#27671;&#35937;&#24193;&#12424;&#12426;20220503&#12398;&#27598;&#26178;&#12398;&#26368;&#39640;&#27671;&#28201;&#12434;&#21462;&#24471;#####
  #&#21442;&#32771;:https://www.data.jma.go.jp/obd/stats/data/mdrr/docs/csv_dl_readme.html
  MaxTemp <- read.csv(paste0("https://www.data.jma.go.jp/obd/stats/data/mdrr/tem_rct/alltable/mxtemsadext00_20220503", Hour[i], ".csv"),
                         header = T, fileEncoding = "cp932")
  #&#26368;&#39640;&#27671;&#28201;&#20966;&#29702;
  GetMaxTemp <- NULL
  for(n in 1:47){
    #&#37117;&#36947;&#24220;&#30476;&#12434;&#25277;&#20986;
    GetPrefData <- MaxTemp[which(MaxTemp[, 2] %in% grep(JpanPref[n], MaxTemp[, 2], value = TRUE)),]
    #&#26368;&#39640;&#27671;&#28201;&#12434;&#38477;&#38918;&#12391;&#20006;&#12403;&#26367;&#12360;
    GetPrefData <- GetPrefData[order(GetPrefData[, 10], decreasing = TRUE),]
    #&#26368;&#39640;&#27671;&#28201;&#12434;&#21462;&#24471;
    GetMaxTemp <- c(GetMaxTemp, GetPrefData[1, 10])
  }
  
  HourTemp <- cbind(Hour[i], JpanPref, GetMaxTemp)
  NewMaxTemp <- rbind(NewMaxTemp, HourTemp)
  
}

#&#21015;&#21517;&#12434;&#20184;&#19982;
colnames(NewMaxTemp) <- c("Hour", "Pref", "MaxTemp")

#&#26368;&#39640;&#27671;&#28201;&#12434;&#25968;&#20516;&#21270;
NewMaxTemp[, 3] <- type.convert(NewMaxTemp[, 3], as.is = TRUE)

#&#26368;&#39640;&#27671;&#28201;&#12391;&#37117;&#36947;&#24220;&#30476;&#12434;&#20006;&#12403;&#26367;&#12360;
#&#28310;&#20633;
NewMaxTemp %>%
  group_by(Pref) %>%
  summarise(Max = max(MaxTemp)) %>%
  arrange(Max) %>%
  mutate(Pref = factor(Pref)) %>%
  select(Pref) -> OrderVecPref
#&#20006;&#12403;&#26367;&#12360;&#12392;Hour&#12434;factor&#21270;
NewMaxTemp %>%
  mutate(Hour = factor(Hour),
         Pref = factor(Pref, levels = OrderVecPref$Pref)) -> NewMaxTemp

joyplotの作成

geom_density_ridgesコマンドを使う

#joyplotの作成:geom_density_ridgesコマンドを使う
ggplot(NewMaxTemp, aes(x = MaxTemp, y = Pref, fill = Pref)) +
  geom_density_ridges(show.legend = F) +
  geom_vline(xintercept = 26.6, col = "#ffc1c1") +
  theme(axis.text.y = element_text(size = 7.5),
        axis.text = element_text(colour = "#ffffe0"),
        panel.grid = element_blank(),
        panel.background = element_rect(fill = "#0a0a0a"),
        plot.background = element_rect(fill = "#0a0a0a"),
        plot.margin = unit(rep(0.2, 4), "cm"))

出力例


少しでも、あなたの解析が楽になりますように!!

タイトルとURLをコピーしました