Julia で統計解析　第 4 章　一変量統計¶

In [1]:
Version()

最新バージョン 2024-10-29 13:50

一変量統計¶

サンプルサイズ，（算術）平均値，不偏分散，標準偏差は基本統計量と呼ばれる。その他に五数要約値である最小値，第一四分位数，中央値（第二四分位数），第三四分位数，最大値も必要である。

R に用意されている airquality データセットを用いて例示する。

In [2]:
using RCall
airquality = rcopy(R"airquality");
println("size: ", size(airquality))
first(airquality, 5)

size: (153, 6)
Out[2]:
5×6 DataFrame
RowOzoneSolar_RWindTempMonthDay
Int64?Int64?Float64Int64Int64Int64
1411907.46751
2361188.07252
31214912.67453
41831311.56254
5missingmissing14.35655

Ozone, Solar_R は欠損値を含むので取り扱いに注意が必要である。dropmissing() により欠損値をリストワイズ除去したデータフレームも作っておく。同時に，分析に必要のない変数を「選択しない」（それ以外を選択する）ために select() を使う。

In [3]:
using DataFrames
airquality2 = dropmissing(select(airquality, Not([:Month, :Day])))
# airquality2 = dropmissing(select(airquality, [1, 2, 3, 4])) などと同じ
println("size: ", size(airquality2))
first(airquality2, 5)

size: (111, 4)
Out[3]:
5×4 DataFrame
RowOzoneSolar_RWindTemp
Int64Int64Float64Int64
1411907.467
2361188.072
31214912.674
41831311.562
5232998.665

一変数の場合¶

基礎統計量¶

単変量の基礎統計量を求める関数として，sum(), mean(), var(), std(), median(), quantile() がある。

個別の変数に対して適用する場合は以下のようにする。

In [4]:
using Statistics
mean(airquality.Ozone)
Out[4]:
missing

Ozone には欠損値 missing が含まれているので，平均値を求めるためには前もって skipmissing() を適用しなければならない。

In [5]:
mean(skipmissing(airquality.Ozone))
Out[5]:
42.12931034482759

何回も Ozone について統計処理をする場合，使用するたびに skipmissing() するのは無駄なので，結果を変数に代入しておけばよい。

なお，以下の使用例で最後に |> plintln をつけているのはこの資料作成ツールのくせに対応するためなので，REPL モードで実行する場合は不要である。

In [6]:
x = skipmissing(airquality.Ozone);
mean(x) |> println
var(x)  |> println
std(x)  |> println

42.12931034482759
1088.2005247376317
32.98788451443396

なお，airquality2 は一つでも欠損値を含む行は削除しているので，airquality2 においては skipmissing() を使う必要はないが，サンプルサイズが異なるので，結果が違うことに注意する必要がある。

In [7]:
mean(airquality2.Ozone)
Out[7]:
42.0990990990991

データの個数は length() で確認できるが skipmissing() で得られるデータは使用できない。以下で定義するような filter() を使う必要がある。filter() は，欠損値を除くデータを作る。skipmissing() は欠損値を含むデータというデータ型を作るだけで，実際に missing が除去されるわけではない。しかし，他の統計関数では skipmissing() による場合も filter(!ismissing, ) による場合も結果は同じである。

In [8]:
y = filter(!ismissing, airquality.Ozone)
mean(y)   |> println
var(y)    |> println
std(y)    |> println
length(y) |> println

42.12931034482759
1088.200524737631
32.987884514433944
116

パーセンタイル値¶

パーセンタイル値は quantile() で求める。第 2 引数で 任意の q 値をベクトルで指定できる。

In [9]:
using StatsBase
quantile(skipmissing(airquality.Ozone), [0, 0.25, 0.5, 0.75, 1])
Out[9]:
5-element Vector{Float64}:
   1.0
  18.0
  31.5
  63.25
 168.0
In [10]:
quantile(airquality2.Ozone, [0, 0.25, 0.5, 0.75, 1])
Out[10]:
5-element Vector{Float64}:
   1.0
  18.0
  31.0
  62.0
 168.0

度数分布¶

カテゴリー変数の度数分布は freqtable() で求める。freqtable() は欠損値 missing の個数もカウントする。

chickwts は，餌の添加物（feed）の種類による雛鳥の体重（weight）のデータセットである。

In [11]:
using RCall
chickwts = rcopy(R"chickwts")

using FreqTables
freqtable(chickwts.feed)
Out[11]:
6-element Named Vector{Int64}
Dim1      │ 
──────────┼───
casein    │ 12
horsebean │ 10
linseed   │ 12
meatmeal  │ 11
soybean   │ 14
sunflower │ 12

複数の変数の場合¶

分析対象が複数の変数の場合に，変数を 1 つずつ分析するのは効率が悪い。

まとめて処理するために，eachcol(), describe(), combine() 関数がある。

eachcol() を使う¶

eachcol() はデータフレームの各列に対して関数を適用し，結果をベクトルとして返す。

mean の後に . がついているのは，mean が複数の要素（データフレームの各列）に対して処理を行うようにするためである。

In [12]:
airquality = rcopy(R"airquality")
first(airquality, 10)
Out[12]:
10×6 DataFrame
RowOzoneSolar_RWindTempMonthDay
Int64?Int64?Float64Int64Int64Int64
1411907.46751
2361188.07252
31214912.67453
41831311.56254
5missingmissing14.35655
628missing14.96656
7232998.66557
8199913.85958
981920.16159
10missing1948.669510
In [13]:
M = mean.(eachcol(airquality));
M
Out[13]:
6-element Vector{Union{Missing, Float64}}:
   missing
   missing
  9.957516339869281
 77.88235294117646
  6.993464052287582
 15.803921568627452

REPL での出力と，println() での出力形式が異なるので，DataFrame() にしてから println() すればよい。

なお，そのままでは各行の変数名が表示されないので，DataFrame() で，variable=names(データフレーム名) とすれば 1 列目に変数名が付加される。

DataFrame() の中での，a=b という指定は，「変数 b を列名 a としてデータフレームに含める」という意味である（a と b が同じでも構わない）。

In [14]:
DataFrame(variable=names(airquality), Mean=M) |> println

6×2 DataFrame
 Row │ variable  Mean          
     │ String    Float64?      
─────┼─────────────────────────
   1 │ Ozone     missing       
   2 │ Solar_R   missing       
   3 │ Wind            9.95752
   4 │ Temp           77.8824
   5 │ Month           6.99346
   6 │ Day            15.8039

各変数に欠損値が含まれる場合には skipmissing() した後に統計関数を適用しなければならないので，以下のようにする。 map() は，第 2 引数のそれぞれの要素に対して 第 1 引数で定義される関数を適用する。関数は普通の関数と無名関数の場合がある。無名関数は「仮引数 -> 関数定義」の形をしている。

In [15]:
M = map(x -> mean(skipmissing(x)), eachcol(airquality))
M
Out[15]:
6-element Vector{Float64}:
  42.12931034482759
 185.93150684931507
   9.95751633986928
  77.88235294117646
   6.993464052287582
  15.803921568627452

上の例で使用した無名関数と同じことをする普通の関数として定義し，それを map() で使用する場合は以下のようになる。

In [16]:
mean2(x) = mean(skipmissing(x))
M = map(mean2, eachcol(airquality))
Out[16]:
6-element Vector{Float64}:
  42.12931034482759
 185.93150684931507
   9.95751633986928
  77.88235294117646
   6.993464052287582
  15.803921568627452

複数の統計量を計算してそれらを DataFtame() でまとめて一覧表にすることができる。

In [17]:
M = map(x -> mean(skipmissing(x)), eachcol(airquality))
V = map(x -> var(skipmissing(x)), eachcol(airquality))
S = map(x -> std(skipmissing(x)), eachcol(airquality))
n = map(x -> length(filter(!ismissing, x)), eachcol(airquality))
DataFrame(variable=names(airquality), n=n, Mean=M, Var=V, SD=S)
Out[17]:
6×5 DataFrame
RowvariablenMeanVarSD
StringInt64Float64Float64Float64
1Ozone11642.12931088.232.9879
2Solar_R146185.9328110.5290.0584
3Wind1539.9575212.41153.523
4Temp15377.882489.59139.46527
5Month1536.993462.006541.41652
6Day15315.803978.57978.86452

このように map() を使う方法は，同じような指定を何回も繰り返さなければならないが，必要な処理をまとめて関数化しておけば 1 行の関数呼び出しで実行できるので簡単に使用できる。

In [18]:
using Statistics
function describe2(df)
    M = map(x -> mean(skipmissing(x)), eachcol(df))
    V = map(x -> var(skipmissing(x)), eachcol(df))
    S = map(x -> std(skipmissing(x)), eachcol(df))
    n = map(x -> length(filter(!ismissing, x)), eachcol(df))
    DataFrame(variable=names(df), n=n, Mean=M, Var=V, SD=S)
end
Out[18]:
describe2 (generic function with 1 method)
In [19]:
describe2(airquality) |> println

6×5 DataFrame
 Row │ variable  n      Mean       Var         SD       
     │ String    Int64  Float64    Float64     Float64  
─────┼──────────────────────────────────────────────────
   1 │ Ozone       116   42.1293   1088.2      32.9879
   2 │ Solar_R     146  185.932    8110.52     90.0584
   3 │ Wind        153    9.95752    12.4115    3.523
   4 │ Temp        153   77.8824     89.5913    9.46527
   5 │ Month       153    6.99346     2.00654   1.41652
   6 │ Day         153   15.8039     78.5797    8.86452
In [20]:
describe2(airquality2) |> println

4×5 DataFrame
 Row │ variable  n      Mean       Var        SD       
     │ String    Int64  Float64    Float64    Float64  
─────┼─────────────────────────────────────────────────
   1 │ Ozone       111   42.0991   1107.29    33.276
   2 │ Solar_R     111  184.802    8308.74    91.1523
   3 │ Wind        111    9.93964    12.6573   3.55771
   4 │ Temp        111   77.7928     90.8203   9.52997

describe() を使う¶

describe() では，より簡単に複数の変数に対して複数の統計量をまとめて求めることができる。

必要な統計量は第 2 引数以降に任意の順序でシンボルによって列挙する。指定できるのは :mean, :std, :min, :q25, :median, :q75, :max, :nunique, :nmissing, :first, :last, :eltype である。

すべての統計量を計算するときは :all を指定すればよい。

describe() には総和 sum()，不偏分散 var() は実装されていない。

describe() ではパーセンタイル値は q0(min), q25, q50(median), q75, q100(max) を求めることができる。

なお，REPL では結果の出力桁数（幅）により表示が打ち切られる。 以下のように printlins() で明示的に出力すれば打ち切られることはない。 ただし，物理的な出力桁数の制限により折り返し表示されるので見にくくなることはある。

In [21]:
println(describe(airquality, :all))

6×16 DataFrame
 Row │ variable  mean       std       min   q25      median   q75      max    sum      nunique  nuniqueall  nmissing  nnonmissing  first  last   eltype                
     │ Symbol    Float64    Float64   Real  Float64  Float64  Float64  Real   Real     Nothing  Int64       Int64     Int64        Real   Real   Type                  
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Ozone      42.1293   32.9879    1      18.0      31.5    63.25  168     4887                     67        37          116   41     20    Union{Missing, Int64}
   2 │ Solar_R   185.932    90.0584    7     115.75    205.0   258.75  334    27146                    117         7          146  190    223    Union{Missing, Int64}
   3 │ Wind        9.95752   3.523     1.7     7.4       9.7    11.5    20.7   1523.5                   31         0          153    7.4   11.5  Float64
   4 │ Temp       77.8824    9.46527  56      72.0      79.0    85.0    97    11916                     40         0          153   67     68    Int64
   5 │ Month       6.99346   1.41652   5       6.0       7.0     8.0     9     1070                      5         0          153    5      9    Int64
   6 │ Day        15.8039    8.86452   1       8.0      16.0    23.0    31     2418                     31         0          153    1     30    Int64

sum(), var() に対応する関数 describe2() を定義しておく。指定できるのは :sum, :mean, :var, :std, :min, :q25, :median, :q75, :max, :nmissing である。

すべての統計量を計算するときは，統計量の指定をしないようにするか :all を指定すればよい。

In [22]:
function describe2(df, arg...=:all)
    q25(x) = quantile(x, 0.25)
    q75(x) = quantile(x, 0.75)
    nmissing(x) = nrow(df) - length(ismissing.(x))
    arg0 = [:sum, :mean, :var, :std, :min, :q25, :median, :q75, :max, :nmissing]
    func = [sum, mean, var, std, minimum, q25, median, q75, maximum, nmissing]
    funcname = ["sum", "mean", "var", "std", "minimum", "q25", "median", "q75", "maximum", "nmissing"]
    arg[1] == :all && (arg = arg0)
    res = DataFrame(variable = Symbol.(names(df)))
    n = ncol(df)
    m = length(arg)
    body = Array{Float64, 2}(undef, n, m)
    loc = indexin(arg, arg0)
    for i in 1:n
        x = skipmissing(df[!, i])
        body[i, :] = [func[j](x) for j in loc]
    end
    bodydf = DataFrame(body, funcname[loc])
    pos = indexin(["nmissing"], funcname[loc])[1]
    isnothing(pos) || (bodydf[!, pos] = Int.(bodydf[:, pos]))
    hcat(res, bodydf)
end

describe2(airquality) |> println; println()
describe2(airquality, :all) |> println; println()
describe2(airquality, :q25, :q75, :nmissing) |> println; println()
describe2(airquality2, :q25, :q75, :nmissing) |> println; println()

6×11 DataFrame
 Row │ variable  sum      mean       var         std       minimum  q25      median   q75      maximum  nmissing 
     │ Symbol    Float64  Float64    Float64     Float64   Float64  Float64  Float64  Float64  Float64  Int64    
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Ozone      4887.0   42.1293   1088.2      32.9879       1.0    18.0      31.5    63.25    168.0        37
   2 │ Solar_R   27146.0  185.932    8110.52     90.0584       7.0   115.75    205.0   258.75    334.0         7
   3 │ Wind       1523.5    9.95752    12.4115    3.523        1.7     7.4       9.7    11.5      20.7         0
   4 │ Temp      11916.0   77.8824     89.5913    9.46527     56.0    72.0      79.0    85.0      97.0         0
   5 │ Month      1070.0    6.99346     2.00654   1.41652      5.0     6.0       7.0     8.0       9.0         0
   6 │ Day        2418.0   15.8039     78.5797    8.86452      1.0     8.0      16.0    23.0      31.0         0

6×11 DataFrame
 Row │ variable  sum      mean       var         std       minimum  q25      median   q75      maximum  nmissing 
     │ Symbol    Float64  Float64    Float64     Float64   Float64  Float64  Float64  Float64  Float64  Int64    
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Ozone      4887.0   42.1293   1088.2      32.9879       1.0    18.0      31.5    63.25    168.0        37
   2 │ Solar_R   27146.0  185.932    8110.52     90.0584       7.0   115.75    205.0   258.75    334.0         7
   3 │ Wind       1523.5    9.95752    12.4115    3.523        1.7     7.4       9.7    11.5      20.7         0
   4 │ Temp      11916.0   77.8824     89.5913    9.46527     56.0    72.0      79.0    85.0      97.0         0
   5 │ Month      1070.0    6.99346     2.00654   1.41652      5.0     6.0       7.0     8.0       9.0         0
   6 │ Day        2418.0   15.8039     78.5797    8.86452      1.0     8.0      16.0    23.0      31.0         0

6×4 DataFrame
 Row │ variable  q25      q75      nmissing 
     │ Symbol    Float64  Float64  Int64    
─────┼──────────────────────────────────────
   1 │ Ozone       18.0     63.25        37
   2 │ Solar_R    115.75   258.75         7
   3 │ Wind         7.4     11.5          0
   4 │ Temp        72.0     85.0          0
   5 │ Month        6.0      8.0          0
   6 │ Day          8.0     23.0          0

4×4 DataFrame
 Row │ variable  q25      q75      nmissing 
     │ Symbol    Float64  Float64  Int64    
─────┼──────────────────────────────────────
   1 │ Ozone        18.0     62.0         0
   2 │ Solar_R     113.5    255.5         0
   3 │ Wind          7.4     11.5         0
   4 │ Temp         71.0     84.5         0

元の describe() に :sum と :var　にも対応する別バージョンの関数を呈示しておく。

In [23]:
function getcol(df, sym, newsym)
    df2 = describe(df, sym, :nmissing)
    if sym == :mean
        n = nrow(df)
        df2[!, sym] .*= (n .- df2.nmissing)
    else
        df2[!, sym] .*= df2[:, sym]
    end
    rename!(df2, sym => newsym)
end

function describe3(df, arg...)
    length(arg) == 1 && arg[1] == :all && (arg = (:sum, :mean, :var, :std, :min, :q25, :median, :q75, :max, :nunique, :nmissing, :first, :last, :eltype))
    df2 = DataFrame(variable = Symbol.(names(df)))
    for a in arg
        if a == :sum
            df3 = getcol(df, :mean, :sum)
        elseif a == :var
            df3 = getcol(df, :std, :var)
        else
            df3 = describe(df, a)
        end
        df2 = hcat(df2, select(df3, 2))
    end
    df2
end

describe3(airquality, :all)
Out[23]:
6×15 DataFrame
Rowvariablesummeanvarstdminq25medianq75maxnuniquenmissingfirstlasteltype
SymbolFloat64Float64Float64Float64RealFloat64Float64Float64RealNothingInt64RealRealType
1Ozone4887.042.12931088.232.9879118.031.563.25168374120Union{Missing, Int64}
2Solar_R27146.0185.9328110.5290.05847115.75205.0258.753347190223Union{Missing, Int64}
3Wind1523.59.9575212.41153.5231.77.49.711.520.707.411.5Float64
4Temp11916.077.882489.59139.465275672.079.085.09706768Int64
5Month1070.06.993462.006541.4165256.07.08.09059Int64
6Day2418.015.803978.57978.8645218.016.023.0310130Int64

複数の変数に対して，任意の q 値に対するパーセンタイル値を求めるには，以下のような関数を定義すればよい。missing を除去して計算される。

In [24]:
function quantile2(df, qs=[0, 0.25, 0.5, 0.75, 1])
    Names =["q$(Int(100q))" for q in qs]
    Names = replace(Names, "q0" => "min", "q50" => "median", "q100" => "max")
    qtile = map(x -> quantile(skipmissing(x), qs), eachcol(df))
    df2 = DataFrame(Matrix(reshape(vcat(qtile...), length(qs), :)'), Names)
    n = map(x -> length(filter(!ismissing, x)), eachcol(df))
    insertcols!(df2, 1, :variable => names(df), :n => n)
    df2
end

quantile2(airquality) |> println; println()
quantile2(airquality, [0.05, 0.25, 0.75, 0.95]) |> println; println()

6×7 DataFrame
 Row │ variable  n      min      q25      median   q75      max     
     │ String    Int64  Float64  Float64  Float64  Float64  Float64 
─────┼──────────────────────────────────────────────────────────────
   1 │ Ozone       116      1.0    18.0      31.5    63.25    168.0
   2 │ Solar_R     146      7.0   115.75    205.0   258.75    334.0
   3 │ Wind        153      1.7     7.4       9.7    11.5      20.7
   4 │ Temp        153     56.0    72.0      79.0    85.0      97.0
   5 │ Month       153      5.0     6.0       7.0     8.0       9.0
   6 │ Day         153      1.0     8.0      16.0    23.0      31.0

6×6 DataFrame
 Row │ variable  n      q5       q25      q75      q95     
     │ String    Int64  Float64  Float64  Float64  Float64 
─────┼─────────────────────────────────────────────────────
   1 │ Ozone       116     7.75    18.0     63.25    108.5
   2 │ Solar_R     146    24.25   115.75   258.75    311.5
   3 │ Wind        153     4.6      7.4     11.5      15.5
   4 │ Temp        153    60.2     72.0     85.0      92.0
   5 │ Month       153     5.0      6.0      8.0       9.0
   6 │ Day         153     2.0      8.0     23.0      29.4
In [25]:
function quantile3(df, q=[0, 0.25, 0.5, 0.75, 1])
    res = DataFrame(variable = Symbol.(names(df)))
    n = ncol(df)
    m = length(q)
    body = Array{Float64, 2}(undef, n, m)
    for i in 1:n
        body[i, :] = quantile(skipmissing(df[!, i]), q)
    end
    colnames = ["q$(Int(q0*100))" for q0 in q]
    replace!(colnames, "q0" =>  "min", "q50" => "median", "q100" => "max")
    hcat(res, DataFrame(body, colnames))
end

quantile3(airquality) |> println; println()
quantile3(airquality2) |> println; println()

6×6 DataFrame
 Row │ variable  min      q25      median   q75      max     
     │ Symbol    Float64  Float64  Float64  Float64  Float64 
─────┼───────────────────────────────────────────────────────
   1 │ Ozone         1.0    18.0      31.5    63.25    168.0
   2 │ Solar_R       7.0   115.75    205.0   258.75    334.0
   3 │ Wind          1.7     7.4       9.7    11.5      20.7
   4 │ Temp         56.0    72.0      79.0    85.0      97.0
   5 │ Month         5.0     6.0       7.0     8.0       9.0
   6 │ Day           1.0     8.0      16.0    23.0      31.0

4×6 DataFrame
 Row │ variable  min      q25      median   q75      max     
     │ Symbol    Float64  Float64  Float64  Float64  Float64 
─────┼───────────────────────────────────────────────────────
   1 │ Ozone         1.0     18.0     31.0     62.0    168.0
   2 │ Solar_R       7.0    113.5    207.0    255.5    334.0
   3 │ Wind          2.3      7.4      9.7     11.5     20.7
   4 │ Temp         57.0     71.0     79.0     84.5     97.0

パーセンタイル値を求める q 値を指定できる。

In [26]:
quantile2(airquality2, (1:2:9)./10)
Out[26]:
4×7 DataFrame
Rowvariablenq10q30medianq70q90
StringInt64Float64Float64Float64Float64Float64
1Ozone11111.020.031.049.089.0
2Solar_R11137.0137.0207.0248.0285.0
3Wind1115.78.09.711.514.9
4Temp11164.073.079.083.090.0

分析に使用する変数は以下のように指定できるので，新しくデータフレームを作る必要はない。

In [27]:
quantile2(airquality[:, [1, 3, 4]])
Out[27]:
3×7 DataFrame
Rowvariablenminq25medianq75max
StringInt64Float64Float64Float64Float64Float64
1Ozone1161.018.031.563.25168.0
2Wind1531.77.49.711.520.7
3Temp15356.072.079.085.097.0

combine(), select()/select!(), transform()/transform!() を使う¶

Hadley Wickham が書いた論文 "The Split-Apply-Combine Strategy for Data Analysis" にあるように，統計解析は「データをグループ化し(split)，それぞれのグループに関数を適用し(apply)，結果をまとめ上げる(combine)」という三段階からなる。

Julia ではこれを，combin2, select/select!, transform/transform! で行う。

combine() は，第 1 引数で渡されたデータフレームについて，第 2 引数以降で指示された関数を適用し，結果をデータフレームにまとめる。

本来は複数のサブグループに対して用いるものであるが，データフレーム全体も 1 個のサブグループに過ぎないので，適用してもなんの問題もない。

第 2 引数は何通りかの書き方がある。

  • nrow のようにデータフレームを引数とする関数名だけを書く場合。nrow() はデータフレームに作用する関数で，データフレームの行数を返す。結果のデータフレームで，nrow 列に 153 が書かれている。結果は 1 個の数値に限らない。たとえば size() はデータフレームの行数と列数をタプルで返すので，結果の x1 の列に (153, 6) が書かれている。列名の x1 はデフォルトである。quantil() はベクトルを返すが，tuple に変換しないと結果が見にくくなる。
  • nrow => :名前　のように書くと，列名を定義できる。:nrow と書いたので，結果のデータフレームで，列名が n になっている。
  • 最初の => の左にあるものはデータフレームの列を示す。シンボル(: がついている）や文字列（" " で囲まれている）や数字（列の番号）などの場合もある。
  • 最初の => の右にあるのものは関数（自作の関数も可）または無名関数である。
  • 関数を何段階か順番に適用するときには無名関数で定義する（=> で関数を数珠つなぎにすることはできない）。
  • 関数/無名関数の後に => に続くものがある場合は，それは結果のデータフレームの列名になる（重複する名前は指定できない）。ない場合にはデフォルトでもっともらしい名前が付けられる。
In [28]:
function my_std(x)
    mean = sum(x) / length(x)
    var = sum((x .- mean).^2) / (length(x) - 1)
    sqrt(var)
end

results = combine(airquality,
    nrow,
    size,
    nrow => :n,
    :Wind => std,
    "Wind" => var => :var_wind,
    3 => var => "var-wind",
    3 => (x -> sqrt(var(x))) => "std-wind",
    :Wind => my_std => "std-wind-2",
    :Ozone => (x -> length(filter(!ismissing, x))) => "valid cases",
    :Ozone => (x -> mean(skipmissing(x))) => :Mean,
    :Temp => (x -> tuple(quantile(x, [0.25, 0.5, 0.75]))) => "%-tiles(25, 50, 75)"
)
typeof(results) |> println; println()
results |> println

DataFrame

1×11 DataFrame
 Row │ nrow   x1        n      Wind_std  var_wind  var-wind  std-wind  std-wind-2  valid cases  Mean     %-tiles(25, 50, 75)   
     │ Int64  Tuple…    Int64  Float64   Float64   Float64   Float64   Float64     Int64        Float64  Tuple…                
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │   153  (153, 6)    153     3.523   12.4115   12.4115     3.523       3.523          116  42.1293  ([72.0, 79.0, 85.0],)

グループごとの記述統計量¶

データフレームをグループを表す変数でグループ化し，それぞれのサブグループごとに統計量を求める。

airquality データセットの :Month でグループ化する。

select()/select!() はデータフレームの列を選択する。Not(シンボル)，Not(シンボルベクトル) はそれらを含まないことを表す。

groupby(データフレーム，グループ化に使う列の指定) でグループ化する。列の指定はベクトルでもよい。

In [29]:
using RCall
airquality = rcopy(R"airquality");
gdf = groupby(select(airquality, Not(:Day)), :Month);
println("number of subgroup = ", length(gdf))

number of subgroup = 5

gdf[1] は :Month が 5 の場合のサブグループである。

In [30]:
first(gdf[1], 5)
Out[30]:
5×5 DataFrame
RowOzoneSolar_RWindTempMonth
Int64?Int64?Float64Int64Int64
1411907.4675
2361188.0725
31214912.6745
41831311.5625
5missingmissing14.3565

describe() を使う¶

groupby() により分割されたそれぞれの gdf[i] は，データフレームなので，前節に述べた方法をそのまま適用できる。

In [31]:
describe(select(gdf[1], Not(:Month)))
Out[31]:
4×7 DataFrame
Rowvariablemeanminmedianmaxnmissingeltype
SymbolFloat64RealFloat64RealInt64Type
1Ozone23.6154118.01155Union{Missing, Int64}
2Solar_R181.2968194.03344Union{Missing, Int64}
3Wind11.62265.711.520.10Float64
4Temp65.54845666.0810Int64
In [32]:
describe2(select(gdf[1], Not(:Month)))
Out[32]:
4×11 DataFrame
Rowvariablesummeanvarstdminimumq25medianq75maximumnmissing
SymbolFloat64Float64Float64Float64Float64Float64Float64Float64Float64Int64
1Ozone614.023.6154493.92622.22441.011.018.031.5115.05
2Solar_R4895.0181.29613242.4115.0758.072.0194.0284.5334.04
3Wind360.311.622612.47113.531455.78.911.514.0520.10
4Temp2032.065.548446.98926.8548756.060.066.069.081.00
In [33]:
quantile(select(gdf[end], Not(:Month)).Wind)
Out[33]:
5-element Vector{Float64}:
  2.8
  7.550000000000001
 10.3
 12.325
 16.6
In [34]:
quantile2(select(gdf[end], Not(:Month)))
Out[34]:
4×7 DataFrame
Rowvariablenminq25medianq75max
StringInt64Float64Float64Float64Float64Float64
1Ozone297.016.023.036.096.0
2Solar_R3014.0116.75192.0234.5259.0
3Wind302.87.5510.312.32516.6
4Temp3063.071.076.081.093.0

for ループで，すべてのサブグループについて分析できる。

In [35]:
for df in gdf
    println("\nMonth: ", df[1, :Month], "\n", describe(select(df, Not(:Month))))
end

Month: 5
4×7 DataFrame
 Row │ variable  mean      min   median   max    nmissing  eltype                
     │ Symbol    Float64   Real  Float64  Real   Int64     Type                  
─────┼───────────────────────────────────────────────────────────────────────────
   1 │ Ozone      23.6154   1       18.0  115           5  Union{Missing, Int64}
   2 │ Solar_R   181.296    8      194.0  334           4  Union{Missing, Int64}
   3 │ Wind       11.6226   5.7     11.5   20.1         0  Float64
   4 │ Temp       65.5484  56       66.0   81           0  Int64

Month: 6
4×7 DataFrame
 Row │ variable  mean      min   median   max    nmissing  eltype                
     │ Symbol    Float64   Real  Float64  Real   Int64     Type                  
─────┼───────────────────────────────────────────────────────────────────────────
   1 │ Ozone      29.4444  12       23.0   71          21  Union{Missing, Int64}
   2 │ Solar_R   190.167   31      188.5  332           0  Union{Missing, Int64}
   3 │ Wind       10.2667   1.7      9.7   20.7         0  Float64
   4 │ Temp       79.1     65       78.0   93           0  Int64

Month: 7
4×7 DataFrame
 Row │ variable  mean       min   median   max    nmissing  eltype                
     │ Symbol    Float64    Real  Float64  Real   Int64     Type                  
─────┼────────────────────────────────────────────────────────────────────────────
   1 │ Ozone      59.1154    7       60.0  135           5  Union{Missing, Int64}
   2 │ Solar_R   216.484     7      253.0  314           0  Union{Missing, Int64}
   3 │ Wind        8.94194   4.1      8.6   14.9         0  Float64
   4 │ Temp       83.9032   73       84.0   92           0  Int64

Month: 8
4×7 DataFrame
 Row │ variable  mean       min   median   max    nmissing  eltype                
     │ Symbol    Float64    Real  Float64  Real   Int64     Type                  
─────┼────────────────────────────────────────────────────────────────────────────
   1 │ Ozone      59.9615    9       52.0  168           5  Union{Missing, Int64}
   2 │ Solar_R   171.857    24      197.5  273           3  Union{Missing, Int64}
   3 │ Wind        8.79355   2.3      8.6   15.5         0  Float64
   4 │ Temp       83.9677   72       82.0   97           0  Int64

Month: 9
4×7 DataFrame
 Row │ variable  mean      min   median   max    nmissing  eltype                
     │ Symbol    Float64   Real  Float64  Real   Int64     Type                  
─────┼───────────────────────────────────────────────────────────────────────────
   1 │ Ozone      31.4483   7       23.0   96           1  Union{Missing, Int64}
   2 │ Solar_R   167.433   14      192.0  259           0  Union{Missing, Int64}
   3 │ Wind       10.18     2.8     10.3   16.6         0  Float64
   4 │ Temp       76.9     63       76.0   93           0  Int64

変数ごとに統計量を一覧表示¶

変数ごとに一覧表示するには，いわゆるロングフォーマットにして処理をする必要がある。

In [36]:
long = stack(select(airquality, Not(:Day)), [:Ozone, :Solar_R, :Wind, :Temp])
gdf2 = groupby(long, [:variable, :Month]);
first(gdf2[1], 5)
Out[36]:
5×3 DataFrame
RowMonthvariablevalue
Int64StringFloat64?
15Ozone41.0
25Ozone36.0
35Ozone12.0
45Ozone18.0
55Ozonemissing
In [37]:
result = combine(gdf2, nrow,
    :value => (x -> length(filter(!ismissing, x))) => :n,
    :value => (x -> mean(skipmissing(x))) => :Mean,
    :value => (x -> std(skipmissing(x))) => :SD,
    );
result |> println

20×6 DataFrame
 Row │ variable  Month  nrow   n      Mean       SD        
     │ String    Int64  Int64  Int64  Float64    Float64   
─────┼─────────────────────────────────────────────────────
   1 │ Ozone         5     31     26   23.6154    22.2244
   2 │ Ozone         6     30      9   29.4444    18.2079
   3 │ Ozone         7     31     26   59.1154    31.6358
   4 │ Ozone         8     31     26   59.9615    39.6812
   5 │ Ozone         9     30     29   31.4483    24.1418
   6 │ Solar_R       5     31     27  181.296    115.075
   7 │ Solar_R       6     30     30  190.167     92.883
   8 │ Solar_R       7     31     31  216.484     80.5683
   9 │ Solar_R       8     31     28  171.857     76.8349
  10 │ Solar_R       9     30     30  167.433     79.1183
  11 │ Wind          5     31     31   11.6226     3.53145
  12 │ Wind          6     30     30   10.2667     3.76923
  13 │ Wind          7     31     31    8.94194    3.03598
  14 │ Wind          8     31     31    8.79355    3.22593
  15 │ Wind          9     30     30   10.18       3.46125
  16 │ Temp          5     31     31   65.5484     6.85487
  17 │ Temp          6     30     30   79.1        6.59859
  18 │ Temp          7     31     31   83.9032     4.31551
  19 │ Temp          8     31     31   83.9677     6.58526
  20 │ Temp          9     30     30   76.9        8.35567

必要なら小分けにして表示すればよい。

result.variable の型は String なので，Wind の結果だけを抽出するには result.variable .== "Wind" とする（:Wind ではないので注意）。

In [38]:
result[result.variable .== "Wind", [:n, :Mean, :SD]]
Out[38]:
5×3 DataFrame
RownMeanSD
Int64Float64Float64
13111.62263.53145
23010.26673.76923
3318.941943.03598
4318.793553.22593
53010.183.46125
In [39]:
function temp(i, sym1, df)
    combine(filter(x -> x[sym1] == i, df),
    :Wind => (x -> length(filter(!ismissing, x))) => :n,
    :Wind => mean,
    :Wind => std)
end
result2 = temp(5, :Month, airquality);
for i in 6:9
    append!(result2, temp(i, :Month, airquality))
end

result2
Out[39]:
5×3 DataFrame
RownWind_meanWind_std
Int64Float64Float64
13111.62263.53145
23010.26673.76923
3318.941943.03598
4318.793553.22593
53010.183.46125

月ごとの結果が必要なら Month でソートして表示する。

In [40]:
sort(result, :Month)
Out[40]:
20×6 DataFrame
RowvariableMonthnrownMeanSD
StringInt64Int64Int64Float64Float64
1Ozone5312623.615422.2244
2Solar_R53127181.296115.075
3Wind5313111.62263.53145
4Temp5313165.54846.85487
5Ozone630929.444418.2079
6Solar_R63030190.16792.883
7Wind6303010.26673.76923
8Temp6303079.16.59859
9Ozone7312659.115431.6358
10Solar_R73131216.48480.5683
11Wind731318.941943.03598
12Temp7313183.90324.31551
13Ozone8312659.961539.6812
14Solar_R83128171.85776.8349
15Wind831318.793553.22593
16Temp8313183.96776.58526
17Ozone9302931.448324.1418
18Solar_R93030167.43379.1183
19Wind9303010.183.46125
20Temp9303076.98.35567

以下のようなプログラムを書くこともできる。

In [41]:
a = [describe(df) for df in gdf]
variable = [:Ozone, :Solar_R, :Wind, :Temp]
for i in variable
    df = DataFrame[]
    for j in 1:length(a)
        pos = indexin([Symbol(i)], a[j][:, 1])[1]
        ismissing(pos) && continue
        row = DataFrame(a[j][pos, :])
        if df == [] # first
            df = row
        else
            append!(df, row)
        end
    end
    println("\n", df[1,1])
    rename!(df, 1 => :GroupKey)
    nam = [replace("$key", "GroupKey: " => "") for key in keys(gdf)]
    df[!, 1] = nam
    println(df)
end

Ozone
5×7 DataFrame
 Row │ GroupKey      mean     min   median   max   nmissing  eltype                
     │ String        Float64  Real  Float64  Real  Int64     Type                  
─────┼─────────────────────────────────────────────────────────────────────────────
   1 │ (Month = 5,)  23.6154     1     18.0   115         5  Union{Missing, Int64}
   2 │ (Month = 6,)  29.4444    12     23.0    71        21  Union{Missing, Int64}
   3 │ (Month = 7,)  59.1154     7     60.0   135         5  Union{Missing, Int64}
   4 │ (Month = 8,)  59.9615     9     52.0   168         5  Union{Missing, Int64}
   5 │ (Month = 9,)  31.4483     7     23.0    96         1  Union{Missing, Int64}

Solar_R
5×7 DataFrame
 Row │ GroupKey      mean     min   median   max   nmissing  eltype                
     │ String        Float64  Real  Float64  Real  Int64     Type                  
─────┼─────────────────────────────────────────────────────────────────────────────
   1 │ (Month = 5,)  181.296     8    194.0   334         4  Union{Missing, Int64}
   2 │ (Month = 6,)  190.167    31    188.5   332         0  Union{Missing, Int64}
   3 │ (Month = 7,)  216.484     7    253.0   314         0  Union{Missing, Int64}
   4 │ (Month = 8,)  171.857    24    197.5   273         3  Union{Missing, Int64}
   5 │ (Month = 9,)  167.433    14    192.0   259         0  Union{Missing, Int64}

Wind
5×7 DataFrame
 Row │ GroupKey      mean      min   median   max   nmissing  eltype  
     │ String        Float64   Real  Float64  Real  Int64     Type    
─────┼────────────────────────────────────────────────────────────────
   1 │ (Month = 5,)  11.6226    5.7     11.5  20.1         0  Float64
   2 │ (Month = 6,)  10.2667    1.7      9.7  20.7         0  Float64
   3 │ (Month = 7,)   8.94194   4.1      8.6  14.9         0  Float64
   4 │ (Month = 8,)   8.79355   2.3      8.6  15.5         0  Float64
   5 │ (Month = 9,)  10.18      2.8     10.3  16.6         0  Float64

Temp
5×7 DataFrame
 Row │ GroupKey      mean     min   median   max   nmissing  eltype 
     │ String        Float64  Real  Float64  Real  Int64     Type   
─────┼──────────────────────────────────────────────────────────────
   1 │ (Month = 5,)  65.5484    56     66.0    81         0  Int64
   2 │ (Month = 6,)  79.1       65     78.0    93         0  Int64
   3 │ (Month = 7,)  83.9032    73     84.0    92         0  Int64
   4 │ (Month = 8,)  83.9677    72     82.0    97         0  Int64
   5 │ (Month = 9,)  76.9       63     76.0    93         0  Int64
In [42]:
using DataFrames, Statistics
for variable in [:Ozone, :Solar_R, :Wind, :Temp]
    println("\n$variable")
    result = DataFrame[]
    for i in 1:length(gdf)
        each = dropmissing(select(gdf[i], [:Month, variable]))
        res = combine(each, :Month => first => :Month, nrow => :n, variable => mean => :Mean, variable => std => :SD)
        if result == []
            result = res
        else
            result = vcat(result, res)
        end
    end
    println(result)
end

Ozone
5×4 DataFrame
 Row │ Month  n      Mean     SD      
     │ Int64  Int64  Float64  Float64 
─────┼────────────────────────────────
   1 │     5     26  23.6154  22.2244
   2 │     6      9  29.4444  18.2079
   3 │     7     26  59.1154  31.6358
   4 │     8     26  59.9615  39.6812
   5 │     9     29  31.4483  24.1418

Solar_R
5×4 DataFrame
 Row │ Month  n      Mean     SD       
     │ Int64  Int64  Float64  Float64  
─────┼─────────────────────────────────
   1 │     5     27  181.296  115.075
   2 │     6     30  190.167   92.883
   3 │     7     31  216.484   80.5683
   4 │     8     28  171.857   76.8349
   5 │     9     30  167.433   79.1183

Wind
5×4 DataFrame
 Row │ Month  n      Mean      SD      
     │ Int64  Int64  Float64   Float64 
─────┼─────────────────────────────────
   1 │     5     31  11.6226   3.53145
   2 │     6     30  10.2667   3.76923
   3 │     7     31   8.94194  3.03598
   4 │     8     31   8.79355  3.22593
   5 │     9     30  10.18     3.46125

Temp
5×4 DataFrame
 Row │ Month  n      Mean     SD      
     │ Int64  Int64  Float64  Float64 
─────┼────────────────────────────────
   1 │     5     31  65.5484  6.85487
   2 │     6     30  79.1     6.59859
   3 │     7     31  83.9032  4.31551
   4 │     8     31  83.9677  6.58526
   5 │     9     30  76.9     8.35567

二変量統計¶

2 変量の統計量は，カテゴリー変数の場合はクロス集計表，間隔尺度変数の場合はピアソンの積率相関係数 cor()，共分散 var()，順序尺度変数の場合はスピアマンの順位相関係数 corspearman()，ケンドールの順位相関係数 corkendall() がある。

二重クロス集計表¶

2 変量の少なくとも 1 つがカテゴリー変数（名義尺度変数，順序尺度変数）の場合に用いる。対になる変数が間隔尺度変数，非尺度変数の場合は，階級分けして集計しなければならない。

用いるのは 1 変量の場合の度数分布を作成するのと同じ freqtable() である。

esoph はフランスにおける喫煙・飲酒と食道癌のケース・コントロール研究の結果である。集計データなので，元データに展開したものを用いる。

In [43]:
esoph = rcopy(R"esoph")
insertcols!(esoph, 6, :which => " ")
size(esoph) |> println
first(esoph, 2) |> println
(sum(esoph.ncases), sum(esoph.ncontrols)) |> println

(88, 6)
2×6 DataFrame
 Row │ agegp  alcgp      tobgp     ncases   ncontrols  which  
     │ Cat…   Cat…       Cat…      Float64  Float64    String 
─────┼────────────────────────────────────────────────────────
   1 │ 25-34  0-39g/day  0-9g/day      0.0       40.0
   2 │ 25-34  0-39g/day  10-19         0.0       10.0
(200.0, 775.0)
In [44]:
ESOPH = DataFrame()
for row in eachrow(esoph)
    row.which = "Case"
    for j in 1:row.ncases
        push!(ESOPH, row)
    end
    row.which = "Control"
    for j in 1:row.ncontrols
        push!(ESOPH, row)
    end
end
In [45]:
freqtable(ESOPH.which)
Out[45]:
2-element Named Vector{Int64}
Dim1    │ 
────────┼────
Case    │ 200
Control │ 775
In [46]:
select!(ESOPH, Not([:ncases, :ncontrols])); # 症例数の 2 列を除く
In [47]:
size(ESOPH) |> println
first(ESOPH, 5) |> println

(975, 4)
5×4 DataFrame
 Row │ agegp  alcgp      tobgp     which   
     │ Cat…   Cat…       Cat…      String  
─────┼─────────────────────────────────────
   1 │ 25-34  0-39g/day  0-9g/day  Control
   2 │ 25-34  0-39g/day  0-9g/day  Control
   3 │ 25-34  0-39g/day  0-9g/day  Control
   4 │ 25-34  0-39g/day  0-9g/day  Control
   5 │ 25-34  0-39g/day  0-9g/day  Control
In [48]:
freqtable(ESOPH.agegp, ESOPH.alcgp)
Out[48]:
6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-39g/day      40-79     80-119       120+
────────────┼───────────────────────────────────────────
25-34       │        61         45          5          5
35-44       │        89         80         20         10
45-54       │        78         81         39         15
55-64       │        89         84         43         26
65-74       │        71         53         29          8
75+         │        27         12          2          3
In [49]:
freqtable(ESOPH.agegp, ESOPH.tobgp)
Out[49]:
6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-9g/day     10-19     20-29       30+
────────────┼───────────────────────────────────────
25-34       │       70        19        11        16
35-44       │      109        46        27        17
45-54       │      104        57        33        19
55-64       │      117        65        38        22
65-74       │       99        38        20         4
75+         │       26        11         3         4
In [50]:
using RCall
airquality = rcopy(R"airquality");
airquality2 = dropmissing(select(airquality, Not([:Month, :Day])));

相関係数，共分散¶

欠損値を含まない場合¶

それぞれの関数を使う¶

変数の対ごとに，関数を適用する。変数の指定法は色々ある。

In [51]:
cor(airquality2.Ozone, airquality2.Temp)         |> println
corspearman(airquality2.Ozone, airquality2.Temp) |> println
corkendall(airquality2.Ozone, airquality2.Temp)  |> println
cov(airquality2.Ozone, airquality2.Temp)         |> println
cov(airquality2[!, :Ozone], airquality2.Temp)    |> println
cov(airquality2[!, "Ozone"], airquality2.Temp)   |> println
cov(airquality2[!, 1], airquality2.Temp)         |> println

0.6985414096486394
0.7729319330689582
0.5861471249834473
221.52072072072067
221.52072072072067
221.52072072072067
221.52072072072067
combine() を使う¶
In [52]:
combine(airquality2,
    [:Ozone, :Temp] => cor,
    [:Ozone, :Wind] => cor,
    [:Ozone, :Solar_R] => cor,
    [:Ozone, :Solar_R] => corspearman,
    )
Out[52]:
1×4 DataFrame
RowOzone_Temp_corOzone_Wind_corOzone_Solar_R_corOzone_Solar_R_corspearman
Float64Float64Float64Float64
10.698541-0.6124970.3483420.348186

欠損値を含む場合¶

それぞれの関数を使う¶

欠損値が含まれるデータの場合は，両方の変数がともに missing でないデータ対のみが分析に使われる。

方法はいくつかあるが，airquality2 を作ったときと同じように，select() で 2 変数のみを含むデータフレームを作り，dropmissing() で欠損値を含むデータを除くのがもっとも簡単であろう。

In [53]:
df = dropmissing(airquality[:, [:Ozone, :Temp]]);
first(df, 3)
Out[53]:
3×2 DataFrame
RowOzoneTemp
Int64Int64
14167
23672
31274
In [54]:
cor(df.Ozone, df.Temp)         |> println
corspearman(df.Ozone, df.Temp) |> println
corkendall(df.Ozone, df.Temp)  |> println
cov(df.Ozone, df.Temp)         |> println

0.6983603421509317
0.7740429554613014
0.5862988215264407
218.5212143928036

combine() を使う¶

欠損値処理を行う無名関数の記述が長いので，何回も使う場合は関数定義をしておくとよい。

In [55]:
cor2(x, y) = cor(Matrix(dropmissing(DataFrame(x=x, y=y))))[1,2]
cov2(x, y) = cov(Matrix(dropmissing(DataFrame(x=x, y=y))))[1,2]
combine(airquality,
    [:Ozone, :Temp] => ((x, y) -> cor(Matrix(dropmissing(DataFrame(x=x, y=y))))[1,2]),
    [:Ozone, :Temp] => cor2,
    [:Ozone, :Wind] => cor2,
    [:Ozone, :Wind] => ((x, y) -> cov(Matrix(dropmissing(DataFrame(x=x, y=y))))[1,2]),
    [:Ozone, :Wind] => cov2,
    )
Out[55]:
1×5 DataFrame
RowOzone_Temp_functionOzone_Temp_cor2Ozone_Wind_cor2Ozone_Wind_functionOzone_Wind_cov2
Float64Float64Float64Float64Float64
10.698360.69836-0.601547-70.9385-70.9385
In [56]:
corspearman(airquality2.Ozone, airquality2.Temp) |> println
corspearman(df.Ozone, df.Temp) |> println

0.7729319330689582
0.7740429554613014

相関係数行列，分散・共分散行列¶

3 個以上の変数を含むデータにおいて，すべての2変数の組み合わせにおける統計量は，相関係数行列，分散共分散行列として得られる。

欠損値を含まない場合¶

関連する変数すべてに欠損値がない場合はそれぞれの関数に行列変換したデータフレームを与えればよい。

In [57]:
cor(Matrix(airquality2)) |> display

4×4 Matrix{Float64}:
  1.0        0.348342  -0.612497   0.698541
  0.348342   1.0       -0.127183   0.294088
 -0.612497  -0.127183   1.0       -0.49719
  0.698541   0.294088  -0.49719    1.0
In [58]:
corspearman(Matrix(airquality2)) |> display

4×4 Matrix{Float64}:
  1.0        0.348186   -0.605136    0.772932
  0.348186   1.0        -0.0616964   0.209537
 -0.605136  -0.0616964   1.0        -0.499323
  0.772932   0.209537   -0.499323    1.0
In [59]:
corkendall(Matrix(airquality2)) |> display

4×4 Matrix{Float64}:
  1.0        0.240319   -0.440459    0.586147
  0.240319   1.0        -0.0430135   0.142902
 -0.440459  -0.0430135   1.0        -0.362387
  0.586147   0.142902   -0.362387    1.0
In [60]:
cov(Matrix(airquality2)) |> display

4×4 Matrix{Float64}:
 1107.29    1056.58    -72.5112  221.521
 1056.58    8308.74    -41.2448  255.468
  -72.5112   -41.2448   12.6573  -16.8572
  221.521    255.468   -16.8572   90.8203

欠損値を含む場合¶

対処法は 2 つある。1 つは，dropmissing() によって欠損値を含まない（リストワイズ除去）データフレームにして 「2.2.1.　欠損値を含まない場合」 を適用する。

もう 1 つは，すべての 2 変数の組み合わせにおいて，2 変数の場合に示したペアワイズ除去により統計量を計算する。

そのための関数を以下に示しておく。

In [61]:
function multi(df, FUNC; samplesize=false)
    n = ncol(df)
    A = ones(n, n)
    N = zeros(Int, n, n)
    for i in 1:n
        x = filter(!ismissing, df[:, i])
        FUNC == cov && (A[i, i] = var(x))
        N[i, i] = length(x)
    end
    for i in 1:n-1
        for j in i+1:n
            df2 = Matrix(dropmissing(df[:, [i, j]]))
            A[j, i] = A[i, j] = FUNC(df2[:, 1], df2[:, 2])
            N[j, i] = N[i, j] = size(df2, 1)
        end
    end
    return (samplesize ? (A, N) : A)
end;

samplesize=false を指定すれば，2変数の組み合わせにおける有効サンプルサイズの行列を返す。そのときには，戻り値を 2 個受け取るので，左辺には 2 個の変数を指定しなければならない。1 つ目は統計量の行列，2 つ目がサンプルサイズ行列である。

In [62]:
df = select(airquality, Not([:Month, :Day]))
A, N = multi(df, cor, samplesize=true)
A|> display
N|> display

4×4 Matrix{Float64}:
  1.0        0.348342   -0.601547    0.69836
  0.348342   1.0        -0.0567917   0.27584
 -0.601547  -0.0567917   1.0        -0.457988
  0.69836    0.27584    -0.457988    1.0
4×4 Matrix{Int64}:
 116  111  116  116
 111  146  146  146
 116  146  153  153
 116  146  153  153
In [63]:
multi(df, corspearman) |> display

4×4 Matrix{Float64}:
  1.0        0.348186     -0.590155      0.774043
  0.348186   1.0          -0.000977333   0.207428
 -0.590155  -0.000977333   1.0          -0.446541
  0.774043   0.207428     -0.446541      1.0
In [64]:
multi(df, corkendall) |> display

4×4 Matrix{Float64}:
  1.0       0.240319    -0.42836      0.586299
  0.240319  1.0          0.00067856   0.144234
 -0.42836   0.00067856   1.0         -0.322242
  0.586299  0.144234    -0.322242     1.0
In [65]:
multi(df, cov) |> display

4×4 Matrix{Float64}:
 1088.2     1056.58   -70.9385  218.521
 1056.58    8110.52   -17.946   229.16
  -70.9385   -17.946   12.4115  -15.2721
  218.521    229.16   -15.2721   89.5913

多変量統計¶

多重クロス集計表¶

3 重，4 重...のクロス集計も，2.1. で用いた freqtable() で行う。

REPL など，環境によっては集計表すべてを表示するのに工夫が必要なことがある。

In [66]:
x = freqtable(ESOPH.agegp, ESOPH.alcgp, ESOPH.tobgp)
for i in 1:size(x, 3)
    println("\nESOPH.tobgp: ", names(x)[3][i])
    display(x[:, :, i])
end

ESOPH.tobgp: 0-9g/day

6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-39g/day      40-79     80-119       120+
────────────┼───────────────────────────────────────────
25-34       │        40         27          2          1
35-44       │        60         35         11          3
45-54       │        46         38         16          4
55-64       │        49         40         18         10
65-74       │        48         34         13          4
75+         │        18          5          1          2
ESOPH.tobgp: 10-19

6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-39g/day      40-79     80-119       120+
────────────┼───────────────────────────────────────────
25-34       │        10          7          1          1
35-44       │        14         23          6          3
45-54       │        18         21         14          4
55-64       │        22         21         15          7
65-74       │        14         10         12          2
75+         │         6          3          1          1
ESOPH.tobgp: 20-29

6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-39g/day      40-79     80-119       120+
────────────┼───────────────────────────────────────────
25-34       │         6          4          0          1
35-44       │         7         14          2          4
45-54       │        10         15          5          3
55-64       │        12         17          6          3
65-74       │         7          9          3          1
75+         │         0          3          0          0
ESOPH.tobgp: 30+

6×4 Named Matrix{Int64}
Dim1 ╲ Dim2 │ 0-39g/day      40-79     80-119       120+
────────────┼───────────────────────────────────────────
25-34       │         5          7          2          2
35-44       │         8          8          1          0
45-54       │         4          7          4          4
55-64       │         6          6          4          6
65-74       │         2          0          1          1
75+         │         3          1          0          0