Rテクニック 幾何平均/幾何標準偏差の算出

Rによる幾何平均/幾何標準偏差の算出

幾何平均とは...
高校で相加・相乗平均の関係をやったと思うが、あれの相乗平均のほう
詳しい計算はwikiみましょか
https://ja.wikipedia.org/wiki/%E5%B9%BE%E4%BD%95%E5%B9%B3%E5%9D%87

対数変換→平均値/標準偏差を算出→指数変換
で、幾何平均や幾何標準偏差を求めることができる

SASで算出する場合は以下を参照
https://qiita.com/saspy/items/a1c300e9f0aa922ce502

今回もSQLコードとrパッケージによる方法を紹介

使用するデータ

SASHELP.CLASSのHEIGHTでためす(ついでに年齢ごとの計算にしてみる)

準備

# パッケージ
Packages <- c( "sqldf" , "data.table" , "tidyverse" , "haven" )
for( i in Packages ) {
  library( i , character.only = TRUE )
}

# 読み込み
Sas_Pass <- "パス/class.sas7bdat"
SASHELP_CLASS <- haven::read_sas( Sas_Pass )

sqldfパッケージで直接的に計算する方法

とりあえず幾何平均、幾何標準偏差、幾何CVを計算
信頼区間もめんどいが式を書けば出せるはず

CLASS_STAT <- fn$sqldf(
  " select AGE , count(*) as N , 
           power( 10 , avg( log10( WEIGHT ) ) ) as GEOMEAN1 ,
           power( exp( 1 ) , avg( log( WEIGHT ) ) ) as GEOMEAN2 ,
           power( 10 , stdev( log10( WEIGHT ) ) ) as GEOSTD , 
           sqrt( exp( power( stdev( log( WEIGHT ) ) , 2 ) ) - 1 ) * 100 as GEOCV
   from SASHELP_CLASS
   group by AGE
 " )

RSQLite(dbにSASHELP_CLASSをSQLデータとして格納してハンドリング)

conn = dbConnect( RSQLite::SQLite() , "" , synchronous = "off" )
dbWriteTable( conn , "SASHELP_CLASS_SQL" , SASHELP_CLASS , overwrite = TRUE )
query <- " select AGE , count(*) as N , 
           power( 10 , avg( log10( WEIGHT ) ) ) as GEOMEAN1 ,
           power( exp( 1 ) , avg( log( WEIGHT ) ) ) as GEOMEAN2 ,
           power( exp( 1 ) , stdev( log( WEIGHT ) ) ) as GEOSTD ,
           sqrt( exp( power( stdev( log( WEIGHT ) ) , 2 ) ) - 1 ) * 100 as GEOCV
           from SASHELP_CLASS_SQL
           group by AGE "

CLASS_STAT2 <- dbGetQuery( conn , query )

dbDisconnect( conn )

dplyr

CLASS_STAT3 <- SASHELP_CLASS %>%
               group_by( Age ) %>% 
               summarise( n = length( Weight ) , GeoMean = exp( mean( log( Weight ) ) ) , GeoStd = exp( sd( log( Weight ) ) ) , .groups = "drop" )

結果

ってな感じでね