イントロダクション

Rのデータセットであるairqualityを用いて、欠測がある変数のバウンドを求める。
データセットは下記のデータを使用する。

df = airquality

データの概要

airquality datasetは1973年5月から9月にかけてニューヨークで行われた毎日の空気質の測定値(計153日分)である。 本分析では変数Ozoneの欠測値のバウンドを求める。

データの前処理

変数Ozoneにおいて、欠測指標を作成する。 なお、本分析においてはパイプ演算子を使つためtidyverseパッケージをインストールする。

# install.packages("tidyverse")
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#欠測なら1,観測なら0の列を追加
df2 = df %>% 
    transmute(R = if_else(is.na(Ozone), 1,0), Ozone)

バウンドの設定

本分析では何も仮定しない場合のバウンドを求める。

\[ LOWER = E[Y|D=1]P(D=1) + Y_L P(D=0) \]

上記はバウンドの下限の数式である。

同様にバウンドの上限については下記の数式であらわされる。

\[ UPPER = E[Y|D=1]P(D=1) + Y_H P(D=0) \]

記号の意味として、

とする。

このとき、推定する変数のバウンドは

\[ LOWER <= E[Y] <= UPPER \]

となる。

バウンドの推定

上記の内容を踏まえ、実際にairqualityデータのOzoneのバウンドを推定する。

#観測率
obs_rate = sum(df2$R == 0) / length(df2$R)

#欠測率
mis_rate = 1 - obs_rate

#下限値と上限値を設定する
low_y = 0
high_y = 200



#下限
low = mean(df$Ozone, na.rm = T)*obs_rate + low_y*mis_rate

#上限
high = mean(df$Ozone, na.rm = T)*obs_rate + high_y*mis_rate

バウンドの下限

low
## [1] 31.94118

バウンドの上限

high
## [1] 80.30719

ここで、下限値Y_Lと上限値Y_Hに関しては、変数の性質や特徴を踏まえて、適切に設定する必要がある。

変数Ozoneはニューヨークで観測されたオゾン濃度の平均値であるから、値が0を下回ることは考えにくいとし、下限値を0とした。

しかし、上限値に関しては観測データの最大値である168を基準に設定したため、適切な値ではない可能性がある。

分析の結果

本分析の結論として何も仮定しない場合のOzoneのバウンドは

\[ 31.9 <= E[Y] <=80.3 \] となった。 したがって、Ozoneの欠測値はこの間に含まれると考えられる。

結論

分析の結果、Ozoneの欠測値は31.9 ~ 80.3の間に含まれると推定できた。
本分析の問題点としては、上限値の設定において、分析者の主観的要素が含まれており、適切な値でない可能性がある。
したがって、今後の目標として上限値が設定できない場合の部分識別の方法を用いて、推定を行う。

参考文献