の続きです。前回は温室効果ガスの排出量のヒストグラムを描きました。分布は対数正規分布に近い分布でした。対数正規分布は企業の売上高のような毎年毎年の成長率の積み重ねの性質のデータにみられる分布です。温室効果ガスの排出量も前年の排出量からの成長率で決まってくるのでしょうね。
今回は、肝心の温室効果ガス排出量は増えたのか？減ったのか？を統計学的に検証してみようと思います。
まずは、2023年と2021年の差分を計算するために、データフレームをいろいろ操作します。
まず、2023年のデータフレームを作成しました。
2021年のデータフレームも作ります。
この２つのデータフレームを結合します。codeとtypeをキーにして結合します。codeだけだと、日本製鉄のように、事業所としてと荷主としての両方に登場する企業があるからです。
inner_join()関数で結合します。
l_tCO2_2021とl_tCO2_2023の差分を計算します。
中央値、平均値がマイナスの値です。温室効果ガスの排出量は減ったようですね。
ヒストグラムを描いてみます。
0の位置の赤い垂線を引いてみました。わずかではありますが、ヒストグラムの中心は0よりも左側にありますね。
t検定をしてみます。
p値は4.584e-05と0.05よりも小さい値です。95%信頼区間は、-0.027から-0.0095と0を含んでいません。統計学的にも平均値は0よりも小さいと言えそうです。
こんどは、type別の値をみてみます。
事業所と荷主は中央値、平均値ともにマイナスの値ですが、輸送はともにプラスの値ですね。ヒストグラムを描いてみます。
このヒストグラムを見るかぎり、輸送セクターの企業が排出量が増えたと断言はできそうにないですね。
t検定をしてみます。
p値は0.2837と0.05よりも大きく、95%信頼区間は、-0.014から0.046と0を含んでいます。排出量が増加したとは断言できませんね。
輸送セクターの観測数は24個と少ないので、t.test()関数だけでなく、inferパッケージを使ってシミュレーションベースでの検定もしてみましょう。
まず、inferパッケージを読み込みます。
平均値を計算して保存しておきます。
平均値は0.0159ですね。これはsummary()関数やt.test()関数での結果と同じです。(あたりまえですね)
続いて、シミュレーションします。
この1000回のシミュレーションの結果を視覚化します。
赤い垂線が0.159の位置です。
p値を計算します。
p値は0.228でした。t.test()関数で算出された値よりは小さいですが、それでも0.05よりは十分に大きいです。
輸送セクターの排出量が増えたとは言い切れないですね。
今回の結論としては、2023年のほうが2021年よりも排出量が減った、ということですね。
今回は以上です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# 2023年のデータフレーム
df2023 <- df |>
filter(year ==
2023) |>
rename(tCO2_2023 = tCO2,
l_tCO2_2023 = l_tCO2) |>
select(-year, -sector)
df2023
#
# 2021年のデータフレーム
df2021 <- df |>
filter(year == 2021) |>
rename(tCO2_2021 = tCO2,
l_tCO2_2021 = l_tCO2) |>
select(code, type, tCO2_2021, l_tCO2_2021)
df2021
#
# 2023年と2021年を横に結合
df_wide <- inner_join(df2021, df2023, by = c("type", "code"))
df_wide
#
# l_tCO2の差分
df_wide <- df_wide |>
mutate(dl_tCO2 = l_tCO2_2023 - l_tCO2_2021)
summary(df_wide$dl_tCO2)
#
# dl_tCO2のヒストグラム
df_wide |>
ggplot(aes(x = dl_tCO2)) +
geom_histogram(color = "white", bins = 100) +
geom_vline(xintercept = 0, color = "red")
#
# t検定
t.test(df_wide$dl_tCO2)
#
# type別のdl_tCO2
df_wide |>
group_by(type) |>
summarize(
n = n(),
min = min(dl_tCO2),
median = median(dl_tCO2),
mean = mean(dl_tCO2),
max = max(dl_tCO2),
sd = sd(dl_tCO2),
cv = sd(dl_tCO2) / mean(dl_tCO2)
)
#
# 輸送のdl_tCO2のヒストグラム
df_wide |>
filter(type == "輸送") |>
ggplot(aes(x = dl_tCO2)) +
geom_histogram(color = "white", bins = 7, boundary = 0)
#
# 輸送のdl_tCO2は増加したのか？
t.test(df_wide$dl_tCO2[df_wide$type == "輸送"])
#
# inferパッケージによる検定
library(infer)
#
# 平均値を計算
observed_mean <- df_wide |>
filter(type == "輸送") |>
specify(response = dl_tCO2) |>
calculate(stat = "mean")
observed_mean
#
# null distributionの生成
set.seed(999)
null_dist <- df_wide |>
filter(type == "輸送") |>
specify(response = dl_tCO2) |>
hypothesize(null = "point", mu = 0) |>
generate(reps = 1000, type = "bootstrap") |>
calculate(stat = "mean")
#
# null distributionを視覚化
null_dist |>
visualize() +
shade_p_value(observed_mean,
direction = "two-sided")
#
# p値の計算
null_dist |>
get_p_value(obs_stat = observed_mean,
direction = "two-sided")
#
(冒頭の画像は Bing Image Crator で生成しました。プロンプトは、A desert in a remote land, with an oasis, and a close-up view of a cactus flower Photograph です。)
