StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

MacのKarabinerが削除できずターミナルでも"su: Sorry"が出るとき

MacでUSキーボードを使っている人が、漢字かなの変換をしやすくするのに、Karabiner-Elementsというソフトがよく使われている。で、そのソフトのアップデートが上手くいかず、アンインストールしてやりなおそうとして、うまく行かなかったという報告がネット上にけっこうある。
自分も今日、その現象にハマってしまったのだが、大事なことを忘れていただけだった。


Macでアプリを削除するときは、AppCleanerとかを使うことが多いけど、Karabinerはその方法では削除できない。権限がどうのこうのという警告がでるので、ターミナルから

su rm -fr Karabiner-Elements.app


とかやるのだが、

su: Sorry


って出て先に進めない。
Mac OSでデフォルトではrootのパスワードが設定されていないために、こうなることはあるらしいのだが、
「su: Sorry」Macでrootユーザーに変更できない - いづいづブログ


Karabinerに関してはそんな問題でなく、単に、GUIで「Misc」のタブからアンインストールすればいいだけだったw
以前も同じことがあって学習したはずだったのだが、完全に忘れていて、30分ぐらい悩んでしまった……。
将来、検索で思い出せるように、ブログに書いておこうと思った。

「シャブ漬け」という比喩は世間でそれなりに使われている

吉野家の役員のくだんの発言については、特に非難したいとも擁護したいとも思わず、「生娘をシャブ漬けにしてる暇があったらうちの職場の近くに店をつくってくれ」ぐらいの感想なのだが(飲食店がほとんどないので)、「シャブ漬け」という喩えはそれなりに耳にする割に、えらい勢いで世間が怒っているなと感じた。まぁ、「生娘」のほうに怒っている人が多いのかもしれないが。
 

マーケティングの先生が、
 


とおっしゃっていて、それはそうだろうと思う一方で、「しかし『シャブ漬け』っていう比喩はそれなりに聞くけど、どんな場面だっけ?」としばらく悩んで思い出したのだが、「補助金」だな。


沖縄では、地元経済界が基地利権に依存しており、それが基地問題の政治的・社会的な解決を阻んでいることを指して、「シャブ漬け」と呼ぶことがよくあるらしい。確かに、何度か聞いたことがある気がする。

俯瞰的に見れば、基地を押し付け、基地と向き合おうとしない本土の人間がいて、沖縄からいくら声を投げ掛けても無視されてしまう状況がある。とりわけ復帰後は沖縄を黙らせるため、徹底的に各種補助金をつぎ込み既得権益層を引き込むという植民地主義的な戦略があった。「シャブ漬け」という言い方をよくするが、悪いのは明らかにシャブを打った方だ。打たれた側がそこから逃れられなくなるという構造が結果として沖縄にある。竜一君は、そのことへの苛立ちを時々口にしていた。
 
琉球新報(2015/12/17)「対談・石川竜一×初沢亜利 沖縄を撮る意味とは」より


社会学者の上野千鶴子氏は以下のインタビューで、沖縄の基地のみならず、国策プロジェクトを進めるに際して多額の補助金で地元の協力を取り付けることを、「シャブ漬け」と呼んでいた。

——沖縄、特に名護は新基地への賛否をめぐり二つに割れることを強いられてきた。地元としては自分たちが招いたわけでもない新基地に振り回されている。
 
「それは国策だから。巨額の金が動く。石原伸晃前環境相が福島について『最後は金目でしょ』と言ったが、沖縄も同じ構図だ。原発だけじゃない、三里塚(成田闘争)も八ッ場ダムも同じ。国策による補助金漬けの構図。シャブ漬け、だ」
 
——政府からすると、シャブ漬けにしているのに言うことを聞かない沖縄にいら立つのだろう。
 
「住民に権利意識が生まれ、受忍限度が下がったから。良いことだ。強姦だって昔は女たちは泣き寝入りしていた。訴えたら、落ち度があるんじゃないかと女の方が逆に責められた。今は女がガマンしないし、ガマンしなくてもよくなった。それに復帰後、沖縄県民も他府県の県民と同じだという平等意識が生まれた。だから『どうして沖縄だけがガマンしなければいけないのか』と思うようになった。当然の変化だと思う」
 
琉球新報(2014/10/29 )「<憲法・女性・沖縄 上野千鶴子さんに聞く>自らの未来決める時/金で魂買う構図変えよう」より


東海道新幹線のグリーン車に置いてある雑誌『WEDGE』では、ふるさと納税がやり玉に挙がっていた。

この2年で激増している「ふるさと納税」は、地方自治体と地方の生産者を”シャブ漬け”に追い込んでいる。
「下りると損」のチキンレースであるため、自治体や個人に自制を求めるのは不可能だ。国が適切な規制をかけるしかない。
(中略)
一方、自治体も”シャブ漬け”のようになっている。そもそも地元から得られる税収が少ない地方の自治体にとって、ふるさと納税による税収増は魅力的だ。通常の税収を超えるふるさと納税金額が集まる自治体さえ出てきており、後に続けと躍起になっている。
 
WEDGE(2016/1/20)「OPINION」より


JALが破綻した当時の前原誠司国土交通大臣と自民党の三ッ矢憲生議員のやり取りのなかで、親方日の丸の上にあぐらを書いた経営を行ってきたことや、公的資金が注ぎ込まれようとしていることが、「シャブ漬け」と呼ばれていた。
(議事録の表現を事後に修正しているかもしれないが、三ッ矢議員が「薬物中毒のどら息子みたいなもの」と言ったのを受けて、前原大臣が「シャブ漬けの息子」と言い換えている。)

○三ッ矢委員 (略)日本航空というのは、私に言わせると、薬物中毒のどら息子みたいなものですよ。生活に困りましたといって、じゃ、生活費の支援をしてあげましょうと。
(略)
○前原国務大臣 これは控えようと思ったんですが、どら息子とおっしゃいましたね。それはだれがどら息子にしたんですかと私は申し上げたいです。シャブ漬けの息子とおっしゃいましたね。だれがシャブ漬けにしたんだと。それを私はあえてここでは申し上げておきたいと思います。
 
第174回国会 国土交通委員会 第14号(平成22年4月21日(水曜日))


経済学者の金子勝教授は2021年のインタビューで、アベノミクス期の日銀の異次元緩和を指して「シャブ漬け」と呼んでいる。

「異次元策の麻薬効果、僕は『シャブ漬け』と呼ぶが、その最大の被害は産業衰退だ。産業の新陳代謝、必要な構造転換を阻んできた。特に製造業だ。その労働生産性は2000年代初めまでは世界一だったが、今や16位(OECD加盟主要31カ国中=19年データ)。あらゆる指標が落ちている。本来は潰れなければいけない電力会社がまだ生きている。原発部門が大赤字の東芝、三菱重工、日立など重電メーカーを異次元緩和で政策的に支えてきた。軽電のパナソニックも身売りの連続。自動車ですら電気自動車、自動運転になったら世界についていけない」
 
サンデー毎日(2021/8/22)「安倍・菅政治の罪と罰/7 金子勝・慶應大名誉教授が五輪後経済と産業衰退の悪夢を語る」より


古い話だが、テリー伊藤氏が小沢一郎氏について、与党の利権に「シャブ漬け」になっていると批判していた。これは補助金とは違うが。

小沢氏の態度について、演出家のテリー伊藤氏は「てっきりテーブルをひっくり返してくれると思ったのに、裏切られた思いだ。『妻とは別れる』とか言いながらセックスして、不倫関係を清算できずにいるダメな男のノリ。与党として甘い汁を吸ったらもう冷や飯食いには戻れないという一種のシャブ漬け状態」とバッサリ。
 
日刊スポーツ(1999/08/14)「自由党 小沢一郎党首 小渕恵三首相と会談し「離脱」の意思を引っ込める」


要するに、国からもらう補助金をはじめとして、「既得権益」「不労所得」「濡れ手に粟」のようなリソースに依存している状態を指して「シャブ漬け」ということがよくあるわけである。
これは「楽な方法に依存してしまうこと」を意味するので、吉野家の役員が単に「牛丼大好き」の意味で「シャブ漬け」と言ったならそれは不適切な用法だろう。仮に、「牛丼という安易な選択肢に頼ってしまい自炊しなくなる」みたいな意味であれば適切であり、大学生時代の私は生娘ではなかったが、シャブ漬けではあったことになる。


ところで、私が「シャブ漬け」という表現を「それなりに耳にする」と思ったのは、仕事で公共事業の関連の話題に触れることが多いためかも知れないので、どれだけ一般的かはわからない。どうなんでしょうね。

ステップワイズ法で特定の変数を残したい場合

私が住んでいる町では、「ステップワイズ法で、ただし特定の変数は残るようにした上で、最良のモデルを選択してほしい」と言われ、手法の妥当性を説明するのが難しいので気が進まないものの、やらざるを得ないことがあります。
ここにやり方が載っていて、


以下、ロジスティック回帰の例でやってみます。

library(dplyr)
d <- iris 

# 従属変数を2カテゴリにして、値を0,1に変換しておきます(しなくてもいいけど)
d <- d %>%
  mutate(Species = as.numeric(Species)-2) %>% 
  filter(Species != -1)

# まず、すべての変数を入れたロジスティック回帰
full.model <- glm(Species ~., family = binomial(logit), data = d, )
summary(full.model)

# AIC最小となるベストのモデルを探索
best.model.1 <- step(full.model)
summary(best.model.1)

# 特定の変数を残したい場合は、scopeの中でlowerに最小限モデルを設定する
#  (今回は結果的にフルモデルと同じものが採択される)
best.model.2 <- step(full.model, scope=list(lower=Species ~ Sepal.Length, upper=full.model))
summary(best.model.2)

Rでパワーアナリシス(検出力・検定力分析)を行う

たとえば重回帰分析とかを行って、効果が「有意ではない」ということを積極的に主張に取り入れたいときは、検出力が十分であったかを確認しなければならない。というのも、わざと出来の悪いモデルを使ったり、サンプルサイズを少なくしたりすれば、有意でないような結果を得ることは簡単だからだ。まあ検出力も、必要な効果量をどれぐらいに設定するかとか、その効果量に含まれる残差の分散の仮定の置き方などによって、ある程度恣意的に引き上げることができてしまうが、やらないよりはいいでしょう。


G*Powerというソフトが便利なのだが(MacでもWindowsでも使える)、シンプルなケースならRのpwrパッケージでできる。
解説はこのあたりを見ればいいかと。
Power Analysis in R
Quick-R: Power Analysis


G*Powerの説明は、ちょっとわかりにくいけど、たとえばこの解説。
https://www.mizumot.com/method/mizumoto-takeuchi.pdf


検出力の分析は「有意水準」「効果量」「サンプルサイズ」「検出力」のどれかを求める手続きで、ネット上や教科書の解説のほとんどは「サンプルサイズの決定」に関するものだが、べつにそれに限らず、これらのうち3つについて仮定を与えれば残りの1つを決められるので、用途は意外と広い。
効果量については、Cohenが提案している、
f^2 > 0.02 ・・・小さめの効果
f^2 > 0.15 ・・・中くらいの効果
f^2 > 0.35 ・・・大き目の効果
という基準を採用する場合は難しく考えなくてよいのだが、たとえば「回帰係数が0.1」というように具体的な値を検出したいと考えている場合は、効果量の表現に平方和が必要になって、データのばらつき具合についても仮定を置かなくてはならなくなる。


ここでは回帰分析で「有意でない」という結果が得られたときに、それを積極的に主張してよいのかどうか判断する方法を考えたいのだが、正直良くわからないので、以下のような方法で問題がある場合はご指摘いただけると助かります。


データはもう手元にあるわけで、サンプルサイズは決まっている。有意水準はたいがい0.05とかだろう。そのうえで、「検出力0.8以上になる効果量を知りたい」あるいは「◯◯という効果量を有意水準0.05で検出できる確率(検出力)を知りたい」ということになるわけである。
で、効果量について上述のCohenのような基準を考えておけば悩まくていいのだが、こんなものでは含意がよくわからないので、現実的な値で語りたくなる。


たとえば回帰分析を行った結果、x1の係数が凄く小さくて有意でなかったとする。しかしこれだけでは、「効果がない」のか「検出力が低い」のか区別がつかない。
で、「変数x1の回帰係数が少なくとも0.1ぐらいあれば、現実社会へのインプリケーションとして意味あるよね」みたいなことが言える時がある。その場合は、「真の係数が0.1だった場合に、今回と同じサンプルサイズ等でその効果が検出できる確率」を計算して、それが0.8とか0.9ぐらいあれば、「十分検出可能なのに、分析の結果有意でなかったのだから、たぶん意味があるレベルの効果は無いのだろう」と積極的に言えるようになる。「本当は0.1ぐらいの効果があるのに、検出力不足で有意にならなかった」わけではなさそうだからだ。


しかしこの計算をするには、効果量を求めるのに、回帰係数だけでなく残差についても仮定を置かないといけない。そこが面倒なのだが、残差(期待値の周りのばらつき)はもとの回帰分析と同じということにしておけば、話がわかりやすい。要するに、予測値からの乖離は同じぐらいになるという前提にして、「意味のある回帰係数」と「元のx」から「yの予測値」を出し、それに「最初の回帰分析をやったときの残差」を足して、「社会的に意味のある効果があった場合の、目的変数y」のデータを仮につくって、そのデータに対して同じように回帰分析(と分散分析)を行って、効果量を求めるわけである。期待値のまわりのばらつきが同じでも、「意味のある回帰係数」が大きければ、yの変動に占める残差の割合は小さくなるので、検出力が上がることになる。
数学が得意な人が考えればわざわざ仮のデータを作ったりしなくても求められそうだけど、自分は「説明変数同士に相関があるかもしれないからタイプIII平方和を使う」というあたりでもう考える気がなくなったので↓のようにやってみた。


以下のような関数を作っておいて、

library(car)
library(pwr)

calculate_power <- function(model, val.number, sufficient.effect, sig.level=0.05) {
  # library{car}{pwr}が必要。
  # modelにはlmの推定結果を与える。推定するときにlm(..., x=TRUE, y=TRUE, model=TRUE)としておく。
  # val.numberには、切片から数えて何個目の変数について計算するか番号を入れる。
  # sufficient.effectには、検出したい大きさの回帰係数を入れる。
  # val.number, sufficient.effectはベクトルで複数指定してもよい。
  
  coef.sufficient <- model$coefficients
  coef.sufficient[val.number] <- sufficient.effect
  
  # 回帰係数が指定した「社会的に意味ある値」だった場合のy
  y.sufficient <- t(coef.sufficient) %*% t(model$x) + model$residuals  # xに切片も入ってる
  
  # データセットの形にする
  df.sufficient <- data.frame(y.sufficient=t(y.sufficient),
                              as.data.frame(model$model[,-1]))  # modelの1列目はyになってる
  
  # 再び回帰して分散分析結果を得る
  lm.sufficient <- lm(y.sufficient ~ ., data=df.sufficient)
  anova.sufficient <- Anova(lm.sufficient, type='III')
  ss.factor.sufficient <- anova.sufficient$`Sum Sq`[val.number]
  ss.resid.sufficient <- anova.sufficient$`Sum Sq`[length(anova.sufficient$`Sum Sq`)]
  df.factor.sufficient <- anova.sufficient$Df[val.number]
  df.resid.sufficient <- anova.sufficient$Df[length(anova.sufficient$Df)]
  
  # 効果量を計算
  eta.sq <- ss.factor.sufficient / (ss.factor.sufficient + ss.resid.sufficient)
  f.sq <- eta.sq/(1-eta.sq)
  
  # この効果量を検出できる確率
  pow <- pwr.f2.test(u = df.factor.sufficient, v = df.resid.sufficient, f2 = f.sq, sig.level = sig.level)
  
  output <- list(pow$power, f.sq)
  names(output) <- c('Power', 'f2')
  return(output)
}
}


まず、目的変数y、説明変数x1〜x3を作って、回帰分析をやってみる。

> set.seed(1)
> x1 <- rnorm(100,0,1)
> x2 <- rnorm(100,5,1)
> x3 <- rnorm(100,10,1)
> y <- 2*x1 + 3*x2 + 3*x3 +rnorm(100, 0, 10)
> 
> model.1 <- lm(y~x1+x2+x3, x=T, y=T, model=T)
> summary(model.1)

Call:
lm(formula = y ~ x1 + x2 + x3, model = T, x = T, y = T)

Residuals:
     Min       1Q   Median       3Q      Max 
-25.8201  -5.2738   0.0281   6.5058  18.2364 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  -7.1879    11.3108  -0.635    0.5266    
x1            1.4204     1.1169   1.272    0.2065    
x2            2.4506     1.0485   2.337    0.0215 *  
x3            4.0462     0.9712   4.166 0.0000677 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.98 on 96 degrees of freedom
Multiple R-squared:  0.1981,	Adjusted R-squared:  0.1731 
F-statistic: 7.908 on 3 and 96 DF,  p-value: 0.00009087


x1の真の係数は2だけど、1.4204と推定され、pは0.2065だから有意ではない。これが、サンプルサイズ(もしくはばらつきが大きい)のせいなのか、本当に効果がないのかを考えたいなら、「x1がどの程度大きければ意味のある効果なのか」を決める必要がある。
仮にそれが1ぐらいだとすると……

> calculate_power(model=model.1, val.number=2, sufficient.effect=1)$Power
[1] 0.1457216


今回のようなサンプルサイズやばらつき具合のデータでは、「x1の係数が1」という効果が仮にあったとしても0.146ぐらいの確率でしか検出できないみたいだから、「x1に効果はない」と言うのは憚られる。
特に、あってほしい効果が検出されないなら単に「残念です、サンプルサイズを増やすか、実験計画を工夫して誤差を減らしましょう」で済むが、自分の説を主張する上で「あってほしくない効果」だった場合は、ついうっかり「どうだ、俺の言ったとおり効果はないじゃないか!」と言いたくなってしまい危険である。


ちなみに、真の値である「係数は2」は検出できるのかというと……

> calculate_power(model=model.1, val.number=2, sufficient.effect=2)$Power
[1] 0.4330289


4割ぐらいしか検出できないようだ。

Rでファイルをダウンロードするかどうか確認させる関数

ネット上にあるcsvファイルをダウンロードしてきて使う場合、read.table()にurlを与えて直接データフレームをつくってしまう場合もあれば、ファイルとしてダウンロードして置いておきたい場合もあります(実行のたびにダウンロードしたくない等の理由で)。
で、実行するときに、すでにダウンロード済みなのであれば再ダウンロード(上書き)はしたくないという場合があるので、確認しながらダウンロードできると便利です。

簡単な処理ではありますが、コードの使い回しのためにここに書いておきます。
そのファイルがすでに存在するかどうかをメッセージで表示した上で、ダウンロードするかどうか訊いてくるので、コンソールでyを入力したらダウンロード、nを入力したら中止ということにしておきます。
ファイルがすでに存在する場合、そのファイルの更新日時を確認してから(たとえば古かったら)ダウンロードしたいという場合もあるので、更新日時も取得して表示するようにしておきます。
処理確認ダイアログを出す関数askYesNo()の使い方は、askYesNoのドキュメントをみてください。

confirm_file_dl <- function(url,dst){
  # dstには保存するときのパス&ファイル名を与える
  
  if(file.exists(dst)) {
    message(paste('\n\"', dst, '\"', ' already exists!', sep=''))
    mtime <- file.info(dst)$mtime  # ファイルの更新日時の取得
    message(paste('(Last modified at ', mtime, ')\n', sep=''))
  } else {
    message(paste('\n\"', dst, '\"', ' doesn\'t exist!\n', sep=''))
  }

    dl_or_not <- askYesNo(msg='Do you want to download the file?', default = FALSE, prompts = 'y/n/na')
  if(dl_or_not==TRUE) {
    download.file(url, destfile = dst)
  }
}


たとえばGoogleのモビリティ・レポートのファイルは、国別データは1ファイルにまとまってるのですが、現時点で600MBもあって、外出中にテザリングで作業してる時なんかは、必要ないならダウンロードしたくないですね。


ファイルがすでに存在する場合。

> confirm_file_dl(url='https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv',
+                 dst='source/Global_Mobility_Report.csv')

"source/Global_Mobility_Report.csv" already exists!
(Last modified at 2021-08-17 18:13:17)

Do you want to download the file? (y/n/na) 


ファイルが存在しない場合。(nを入力してダウンロードを実行)

> confirm_file_dl(url='https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv',
+                 dst='source/Global_Mobility_Report.csv')

"source/Global_Mobility_Report.csv" doesn't exist!

Do you want to download the file? (y/n/na) y
trying URL 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
Content type 'text/csv' length 637979115 bytes (608.4 MB)
==================================================
downloaded 608.4 MB

Rのsource()関数で呼び出すスクリプトに引数を渡すとき

学生に説明する必要が発生したためエントリ起こしておく。
コマンドライン引数みたいな感じで、source()で呼び出されるスクリプトに呼ぶ側のスクリプトから引数を渡すときは、以下のようにやればよい。
source1.Rからsource3.Rまでのスクリプトを準備してあり、それをreader.Rからsource()関数で呼び出す。
その時に、commandArgs()という関数を使う。

#呼ばれる側のスクリプト1(source1.R)
# 文字列が1つ渡される想定

print(commandArgs())
#呼ばれる側のスクリプト2(source2.R)
# ベクトルが一つ渡される想定

mean(commandArgs())
#呼ばれる側のスクリプト3(source3.R)
# 文字列と数値がリストで渡される想定

args <- commandArgs()
print(rep(args[[1]], args[[2]]))
# 呼ぶ側のスクリプト(reader.R)

# 渡す引数が1個なら単にそれを書けばいい
commandArgs <- function(...) {'ばか'}
source('source1.R')

# ベクトルを渡すこともできる
commandArgs <- function(...) {c(1,2,5,7,8)}
source('source2.R')

# 複数の引数をリストで渡す
commandArgs <- function(...) {list('あほ', 5)}
source('source3.R')


なおここで、function(...)の三点ドットは、定義されていない引数を表す記号。
以下、実行結果。

> # 渡す引数が1個なら単にそれを書けばいい
> commandArgs <- function(...) {'ばか'}
> source('source1.R')
[1] "ばか"
> 
> # ベクトルを渡すこともできる
> commandArgs <- function(...) {c(1,2,5,7,8)}
> source('source2.R')
[1] 4.6
> 
> # 複数の引数をリストで渡す
> commandArgs <- function(...) {list('あほ', 5)}
> source('source3.R')
[1] "あほ" "あほ" "あほ" "あほ" "あほ"