相関係数の差の検定と、回帰係数の差の検定を、Rでやる

たまに、２つの相関係数が有意に異なるのかや、１つの重回帰モデル中の２つの回帰係数が有意に異なるかを示せると、主張が通りやすいという場面がある。
まぁ、あまり必要になることはないのだが、相関係数の差の検定や回帰係数の差の検定について、日本語でググるとRでパッと使える方法がまとまっているわけではなさそうだったので、ここに書いておこう。

相関係数の差の検定

相関係数の差の検定は、Rの{psych}パッケージのr.test()関数で簡単にできる。
下記に関数の解説があるが、
https://www.rdocumentation.org/packages/psych/versions/2.0.12/topics/r.testr.test function - RDocumentation
Test4のtable1が崩れていたり、「Case A: where Masculinity at time 1 (M1) correlates with Verbal Ability .5 (r12), femininity at time 1 (F1) correlates with Verbal ability r13 =.4」とあるところはM1とF1が逆になっちゃってたり、出力の説明が不足していたりするので要注意。そこから参照されている論文Steiger(1980)をみると、いくつかの疑問は解決する。

1) この関数は1つの相関係数が有意かどうか検定したい場合にも使え、たとえば変数1・2ともにサンプルサイズが100で相関係数が0.66だとすると、

library(psych)
r.test(n=100, r12=0.66)

というふうにする。多くの場合両側で検定すると思うが、片側で検定したい場合は
twotailed = F
という引数を付ければよい。

2) ここから2つの相関係数の差の話に移る。まずは、2つの独立した（サンプルを共有していない）相関係数の差を検定する場合。サンプルサイズは異なっても良い。
変数1と2はサンプルサイズが100で相関係数が0.3、変数3と4はサンプルサイズが120で相関係数が0.25だった場合、

r.test(n=100, r12=0.30, r34=0.25, n2=120)

サンプルサイズが等しいのであれば、n2は付けなくてもよい（デフォルトでnと同じになる）。

3) サンプルを共有している変数が3つ（1・2・3）あって、1と2、1と3というふうに1つの変数を共有する形で相関係数同士の差を検定したい場合。

r.test(n, r12, r23, r13)

というふうにする。この場合、r12とr13を比べた検定結果が返される。このとき、r23も引数として与える点に注意。（使い方の説明が不親切だが、冒頭の解説リンクのtest4のCaseAと同じことで、そちらについては論文Steiger(1980)と照合すると分かる。）

4) サンプルは共有してるけど変数は共有していないような、2組の相関係数を比べる場合。（冒頭の解説リンクでいうとtest4のCaseBに相当）

r.test(n, r12, r34, r23, r13, r14, r24)

冒頭の解説リンクは使い方の説明が不親切なのだが、論文Steiger(1980)と照合すると分かる。差を検定したい関心のある相関係数がr12とr34で、この2つの相関係数の差の検定結果を返している。

他に参考になるページ
http://www.snap-tck.com/room04/c01/stat/stat05/stat0501.html統計学入門−第5章
母相関係数の差の検定の概要と結果 :: 【公式】株式会社アイスタット｜統計分析研究所

回帰係数の差の検定

重回帰分析の回帰係数の差の検定については、たとえば複数の回帰係数の信頼区間を出して、大きい方の下限値と小さい方の上限値が重なるかどうかという基準で検定すると、（たとえば5%水準の検定をしたいときに95%信頼区間を用いると）厳しすぎる検定になるらしい。*1

で、使えるケースが多少限られるものの、心理統計の分野で分散分析の延長で共分散分析を習うときに出てくる「平行性の検定」の考え方で、「交互作用が有意になるかどうか」という観点で検定すると簡単にできる。*2

d <- airquality  # 練習用データセットの読み込み
d <- d[complete.cases(d),]  # 欠損値削除
summary(lm(Ozone~Wind+Temp, data = d))

このようにすると以下のような結果が得られる。

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept) -67.3220    23.6210  -2.850         0.00524 ** 
Wind         -3.2948     0.6711  -4.909 0.0000032617283 ***
Temp          1.8276     0.2506   7.294 0.0000000000529 ***

この-3.29と1.83が有意に異なるかを知りたいという話で、まあ標準誤差とかをみても余裕で異なりそうではあるが、以下のように（従属変数であるOzoneは残して）ロング型にデータを変更し、交互作用の検定をすればいい。これは要するに、2つの変数を1つの変数にまとめて、変数名を新たな変数（ダミー変数）として設け、いわば「値」と「変数名」の交互作用をみるような感じ。

d2 <- d %>%
  select(Ozone, Wind, Temp) %>%  # 変数を限定
  pivot_longer(-Ozone, names_to = 'Variable', values_to = 'Value')   # ロング型に変換

summary(lm(Ozone~Variable*Value, data = d2))

以下のような結果が得られる。

Coefficients:
                    Estimate Std. Error t value             Pr(>|t|)    
(Intercept)        -147.6461    19.7613  -7.471      0.0000000000019 ***
VariableWind        246.6873    21.0072  11.743 < 0.0000000000000002 ***
Value                 2.4391     0.2522   9.673 < 0.0000000000000002 ***
VariableWind:Value   -8.1679     0.7210 -11.329 < 0.0000000000000002 ***

「VariableWind:Value」のところが交互作用（Tempが基準=0になって要するにWindダミーになってるという意味）で、有意になってるので、さっきの回帰係数の差は有意ってことになります。

後で知ったのですが、回帰係数の差の検定については、{car}パッケージのlinearHypothesis関数を使って、簡単に調べることができる。これは、2つの説明変数の回帰係数が等しいという帰無仮説を検定するもの。
あまりよく考えてないんですが、従属変数を複数もつ多変量回帰（Multivariate Regression。carパッケージのManovaで検定したりする）の場合、「すべてのモデルにおいて両変数の回帰係数に差がない」という帰無仮説を検定してるような気がする。モデルの数が増えると、帰無仮説がほとんど棄却されるようになります。上のように交互作用項を入れた上でManovaする方法の場合、交互作用はそれほど有意になりやすくはないです。

# 必要なパッケージの読み込み
library(car)

# airqualityデータセットを使用
data("airquality")

# NA値の除去
airquality_complete <- na.omit(airquality)

cor(airquality_complete)

# 線形回帰モデルの作成
model <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality_complete)

# Solar.R と Wind の回帰係数が等しいという仮説を検定
linearHypothesis(model, "Solar.R = Wind")

# Multivariate Regression（多変量回帰：従属変数が複数ある）でもできる
model_multi <- lm(cbind(Ozone, Wind) ~ Solar.R + Temp + Month + Day, data = airquality_complete)

# Manovaで説明変数の総合的な有意性を確認
Manova(model_multi)

# Solar.R と Tempの回帰係数が等しいという仮説を検定（棄却できる）
linearHypothesis(model_multi, "Solar.R = Temp")

# Month と Dayの回帰係数が等しいという仮説を検定（棄却できない）
linearHypothesis(model_multi, "Month = Day")

*1:何かで読んだけど出典を忘れてしまった。

*2:回帰係数の差は交互作用で検定すればいいよっていうやり方自体は、以前、SPSSのマニュアルかなにかで読んだのだが、出典は忘れてしまった。でもまあ心理統計では共分散分析とかの解説でよく出てくる内容だと思います。

「不偏分散の平方根」は不偏標準偏差ではない件（メモ）

不偏分散の平方根を取っても不偏標準偏差にはならないという話があり、私は不真面目な研究者なのでそもそもそんなこと考えたこともなかったですが、知り合いが「数学的な導出はみれば分かるが、“平方根を取っては駄目な理由”が直観的に理解できなくて気持ち悪い」と言っていました。

ちなみに↓の記事によると、「不偏分散の平方根」を「不偏標準偏差」と記述しているケースはけっこうあるらしく、統計ソフトの計算もそうなってるらしいです（自分では確認してないですが）。
不偏分散の平方根は不偏標準偏差じゃなかった - 静粛に、只今統計勉強中
不偏標準偏差とは?：統計検定を理解せずに使っている人のために

で、ググってみたら下記のような解説があり、
Econometrics Beat: Dave Giles' Blog: Unbiased Estimation of a Standard Deviation
要するに直観的には、

• 不偏推定というのは「期待値」に一致する値を求める手続きである
• 期待値には線形性がある（たとえばE(X+Y)=E(X) + E(Y)）
• 平方根を取るのは非線形な変換である
• だから「不偏分散の平方根」は「母分散の平方根の期待値」にはならない

というような話みたいです。
15分ぐらいしか考えてないのでまちがってたらすいません。

それと合わせて思ったのは、不偏推定値というのはそもそも、母集団からランダムサンプリングして推定値を得るという手続きを何回も繰り返して平均を取ると、だんだん母集団の真の値に近づいていくという性質の推定値で、１回１回得られた「不偏標準偏差」や「不偏分散」というものは、「何かの標準偏差」「何かの分散」ではないと考えておくべきな気もする。95%信頼区間というのが、１回１回の推定で得られる「その区間」に何かの95%が収まるという意味での区間ではないのと同様。

Rのループ中に進捗率を表示するプログレスバーを作る

Rのループで使えるプログレスバーは、いくつかのパッケージで提供されているみたいなのですが、自分で書くのも簡単なので、単純な関数でつくってみた。
進捗が知りたいのは時間のかかる処理をするときであり、時間がかかるなら無駄な計算は省きたいので、ループ1回ごとには出さずに100回とか1000回ごとに表示したほうがいいかもしれない。

### 自作プログレスバーをつくる

# バーは40文字からなり、-が#に置き換わっていくようにした。
# 右端に％表示を出しておく。
# 処理は単に、ループ全体の何番目をやってるのかをwhichとlengthで
# 取得して、それを文字列でのバーに置き換えてmessageとして吐くだけ。
# '\r'は行頭復帰を意味し、appendLF=FALSEは改行しないことを意味する。
# これらの組み合わせで、「表示をまるごと更新する」の意味になる。
# '\r'は、message()内の引数としては、prg.barの前でも後ろでもいいし、
# paster時にくっつけといてもよい。
# （messageは、複数の文字列を任意個数並べられる。）

show.progress <- function(i, x){
  # iとxにはforのループ変数とリストを与える
  prg <- round(which(x==i)/length(x)*100)
  done   <- paste(rep('#', round(prg/2.5)),    collapse='')
  remain <- paste(rep('-', 40-round(prg/2.5)), collapse='')
  prg.bar <- paste('|', done, remain, '|  ', as.character(prg), '%', sep='')
  message('\r', prg.bar, appendLF=FALSE)
}

### 試用してみる

# 素因数分解できるパッケージ
library(gmp)
factorize(as.bigz("5656"))

# 1億〜1億3万までの素因数分解を進捗バーつきでやってみる
x<- 100000000:100030000
primes <- list()
for (xi in x){
  show.progress(xi,x)
  xi <- as.character(xi)
  p <- factorize(as.bigz(xi))
  p <- as.numeric(p)
  primes <- c(primes, list(p))
}

使うと、こんな感じで昔のパソコンみたいな動きになる。

ちなみに、環境によっては、messageの行のあとにflush.console()を挟まないと、100%まで行ってからいきなり表示されるような場合もあるかもしれない。自分の場合は不要だったが。

for、apply、ベクトル演算の処理速度の比較

ツイッターのライムラインで、forループをapplyに置き換えた場合の高速化の話が流れていて（こちら）、気になって検索したところ、applyよりむしろwithを使えと言っている人がいた。
r - apply() is slow - how to make it faster or what are my alternatives? - Stack Overflow

ちなみにこれは、with関数が早いというより、要するに「データフレームの行ごとの計算」を「ベクトル同士の計算」に置き換えることによって速くなっていると考えたほうがよい。withのメリットは、データフレームの列を指定するときに'd$'を書かなくでよくなるという点にある。

以下では、1万行×3列のデータで、横向きに標本分散を計算して1万個取り出す処理を、

• forループで空のベクトルにアペンドしていく
• forループで長さを指定したベクトルを更新していく（Pythonのリストの処理でこうすると速くなる場合があったのでやってみる）
• applyを使う
• applyを使い、かつ関数を事前にコンパイルしておく
• ベクトル同士の計算にする（$で列を取り出す）
• ベクトル同士の計算にする（indexで列を取り出す）
• ベクトル同士の計算にする（with関数をつかう）

の7通りで試してみた。

# パッケージよみこみ
library(rbenchmark)  # 速度を測る作業のラッパー
library(compiler)    # 関数をコンパイルする

# データフレームさくせい
d <- data.frame(
  x1 = runif(10000),
  x2 = runif(10000),
  x3 = runif(10000)
)

# 標本分散を出す関数を定義しておく
svar <- function(x){
  return(sum((x-mean(x))^2)/(length(x)))
}


# 関数をコンパイルする
c.svar <- cmpfun(svar)

# 処理時間を比較する
bm1 <- benchmark(
  'for_append' = {
    v1 <- c()
    for(i in 1:nrow(d)){
      v1 <- c(v1, svar(unlist(d[i,])))
    }
  },
  'for_replace' = {
    v2 <- rep(NA, nrow(d))
    for(i in 1:nrow(d)){
      v2[i] <- svar(unlist(d[i,]))
    }
  },
  'apply' = {
    v3 <- apply(d, 1, svar)
  },
  'apply+compile' = {
    v4 <- apply(d, 1, c.svar)
  },
  'vectorize_$' = {
    v5 <- ((d$x1-(d$x1+d$x2+d$x3)/3)^2 + (d$x2-(d$x1+d$x2+d$x3)/3)^2 + (d$x3-(d$x1+d$x2+d$x3)/3)^2)/3
  },
  'vectorize_index' = {
    v6 <- ((d[,1]-(d[,1]+d[,2]+d[,3])/3)^2 + (d[,2]-(d[,1]+d[,2]+d[,3])/3)^2 + (d[,3]-(d[,1]+d[,2]+d[,3])/3)^2)/3
  },
  'vectorize_with' = {
    v7 <- with(d, ((x1-(x1+x2+x3)/3)^2 + (x2-(x1+x2+x3)/3)^2 + (x3-(x1+x2+x3)/3)^2)/3)
  },
  replications = 100)

# いらん列を消して経過時間でソートしておく
bm1 <- bm1[1:4] %>%
  arrange(elapsed)

replications = 100というのは100回処理を繰り返して比較してるということ。
以下のような結果が得られた。

             test replications elapsed relative
1  vectorize_with          100   0.032    1.000
2     vectorize_$          100   0.034    1.063
3 vectorize_index          100   0.042    1.313
4   apply+compile          100   6.618  206.812
5           apply          100   7.565  236.406
6     for_replace          100  52.767 1648.969
7      for_append          100  93.641 2926.281

elapsedが経過時間を表し、relativeは最速のものとの比を示したもの。
とにかく、ベクトル同士の演算に置き換えられる場合はめちゃめちゃ速くなることが分かる。アクセスの仕方で差はほぼ無く、withを使うと記述が簡単になるのでwithを使っておけばよさそうだ。
applyはforループより1桁速いが、それでもベクトル演算に比べると2桁倍の時間がかかっている。コンパイルしても凄く速くはならない。
forループは、上書き方式にするとアペンド方式より速いが、いずれにしてもベクトル演算の数千倍の時間がかかっている。

↓のように値を置き換えるような処理だと、ベクトル演算にはできないので、とりあえずforよりはapplyを使っておくのが良さそうな気がする。コンパイルしたほうが遅くて草。【追記】Twitterで指摘をもらって気づいたけど、↓の例だと縦方向にも処理できるのでベクトル化容易やなｗｗ 横方向にしか処理できない内容であとで試そう。

> ### 値を置換する処理をforとapplyで比較する
> 
> # ベクトル中の0.5より小さい値を0に置換する関数
> rpl <- function(x){
+   replace(x, which(x < 0.5), 0)
+ }
> 
> c.rpl <- cmpfun(rpl)
> 
> bm2 <- benchmark(
+   'for' = {
+     d.new1 <- d
+     for(i in 1:nrow(d)){
+       d[i,] <- rpl(d[i,])
+     }
+   },
+   'for+compile' = {
+     d.new2 <- d
+     for(i in 1:nrow(d)){
+       d[i,] <- c.rpl(d[i,])
+     }
+   },
+   'apply' = {
+     d.new3 <- d
+     apply(d.new2, 1, rpl)
+   },
+   'apply+compile' = {
+     d.new4 <- d
+     apply(d.new2, 1, c.rpl)
+   },
+   replications = 5
+ )
> 
> bm2[1:4] %>% arrange(elapsed)
           test replications elapsed relative
1         apply            5   0.246    1.000
2 apply+compile            5   0.247    1.004
3           for            5  20.957   85.191
4   for+compile            5  22.890   93.049

ggplot2で2軸グラフを描く時の軸スケーリングの作業

ggplot2で2軸のグラフを描くときは、先日のエントリでも書いたように、ggplot2自身は左軸（第1軸）と右軸（第2軸）を別々の情報として持つことはできないので、左軸と右軸の尺度の違いを自分で設定して変換しなければならない。
あとで使いまわすので、このスケーリングの作業を行うスケーラを以下のように定義しておいた。
最初は1行で書いてggplotの描画の中に埋め込んでたけどあとで混乱しないように分けて書いておいた。

library(ggplot2)
data(airquality)  # 練習用データのよみこみ

### 変数を追加
# x軸用に月と日の列を日付の変数にしておく（2行に分けてかいてる）
airquality <- airquality %>%
  mutate(Date = paste(as.character(Month), '/', as.character(Day), sep='')) %>%
  mutate(Date = as.Date(Date, format='%m/%d'))

### y1, y2の目盛り範囲を決めておく
# ここでは恣意的に決めてるが、最大と最小を取るとかでもいいと思う
y1.lim <- c(0, 25)
y2.lim <- c(50, 100)

### スケーラの関数を書いておく
# 上でつくった、y1とy2の目盛りの範囲を定めた要素数2のベクトルを使って、
# いい感じにスケールを合わせる。

# 変数のスケーラ。
# pにy2の値ベクトルを与えると、y1の尺に合わせた数字に変換。
# y2をまずゼロ基準に戻し、y2とy1のlimの幅の比でスケーリング
# した後で、y1のゼロ基準からの乖離分を足す。

variable_scaler <- function(p, lim1, lim2){
  to_zero <- p-lim2[1]
  y1_range <- lim1[2]-lim1[1]
  y2_range <- lim2[2]-lim2[1]
  scaled <- to_zero*y1_range/y2_range
  from_zero <- scaled + lim1[1]
  return(from_zero)
}

# 第2軸の目盛りのスケーラ。
# pは、sec_axis()の'.'になる。y1の目盛りをy2の目盛りに読み替えるもの。
# y1の目盛りをまずゼロ基準に戻し、y1とy2のlimの幅の比スケーリング
# した後で、y2の目盛りのゼロ基準からの乖離分を足す。

axis_scaler <- function(p, lim1, lim2){
  to_zero <- p-lim1[1]
  y1_range <- lim1[2]-lim1[1]
  y2_range <- lim2[2]-lim2[1]
  scaled <- to_zero*y2_range/y1_range
  from_zero <- scaled + lim2[1]
  return(from_zero)
}

### 描画してみる
airquality %>%
  ggplot(aes(x=Date)) +
  geom_line(aes(y=Wind, colour='Wind')) + 
  geom_line(aes(y=variable_scaler(Temp, y1.lim, y2.lim), colour='Temp')) + 
  scale_y_continuous(limit=y1.lim,  # 第1軸の範囲
                     breaks=c(0, 5, 10, 15, 20, 25),  # 第1軸の目盛り
                     sec.axis=sec_axis(
                       ~(axis_scaler(., y1.lim, y2.lim)), # 軸スケーリング
                       breaks=c(50, 60, 70, 80, 90,100), # 第2軸の目盛り
                       name='Temp')  # y2のラベルはここで設定する
                     ) + 
  labs(title = 'DUAL AXIS CHART', 
       x='Date',
       y='Wind',
       colour = 'Variable')

もちろんy2のラベルもscale_y_continuousの中でnameをつかって指定してもよい。
2つ目の折れ線のyを指定するときに変数スケーラを使い、sec.axisを指定するときにformulaの中で軸スケーラを使う。

折れ線グラフの端っこにラベルを付けるやつ（ggplot2）

最近になって遅ればせながらggplot2を頻繁に使うようになってきました。
で、↑こういうふうに、折れ線の端っこにラベルを置きたいと思いました。白黒の記事原稿で4本もの折れ線を重ねるのは見づらいのでそもそもやめたほうがいいですが、人生いろいろあるわけです。
そして、ggplot2の線種は見分けがつきにくいので、凡例だけで示すのは難しい。

ググると解説ブログがいくつか見つかりますが、パット見では何をやっているのか意味がわかりやすくない気もしたので、メモしておきます。
考え方としては、

• ggplot2の他にggrepelパッケージを入れておく。これは互いに重ならない「いい感じのラベル」を書くときに使うもの。
• x軸の範囲を左右に少し広げておく。
• 折れ線の末端を「点」とみなして、そこに新たに散布図を描くようなイメージで、geom_text_repelでテキストを置く。
• 散布図にnudge（点とラベルの距離）を設定することで、強制的にラベルを左右方向にズラす。するとggrepelの機能で、「点」とラベルが線で結ばれるようになる。（点とラベルというか、正確に言うと、geom_text_repelは座標に文字を配置する散布図を描くものだが、座標で指定した位置が点として認識されている。）

という手順になります。「折れ線の端を点とみなした散布図を描いて、そこに指示線つきのラベルを添える」という考え方がポイントだと思います。
イメージとしては、

こういうタイプのグラフをまず思い浮かべればよい。この、ラベルと点が線で繋がれてる感じのものを、折れ線グラフの端っこに重ねてやって、派手な赤点を目立たないようにしてやればいいわけです。

冒頭のグラフは以下のようなコードで描きましたが、

library(ggplot2)
library(ggrepel)

dat %>%
  ggplot(aes(x=YEAR, y=FREQ)) +
  geom_line(aes(size=CODE, linetype=CODE, colour=CODE)) +   # ここ書いとかないとmanual設定も動かない
  theme_classic() +
  geom_text_repel(
    data = dat %>% filter(YEAR==max(YEAR)),    # 折れ線の右端にうつ点
    aes(x=YEAR, y=FREQ, label = CODE),
    nudge_x = 4,
    segment.alpha = 0.5,
    size = 6,
    family="MS Gothic") + 
  geom_text_repel(
    data = dat %>% filter(YEAR==min(YEAR)),
    aes(x=YEAR, y=FREQ, label = CODE),
    nudge_x = -4,                 #  点とラベルの距離
    segment.alpha = 0.5,     # 間に引かれる線をすこし薄くする
    size = 6,                          # ラベルの大きさ
    family="MS Gothic") +  # 日本語を表示するときはフォント指定しておく
  scale_linetype_manual(values = c("solid", "dashed", "solid", "twodash")) +         # 線のタイプ
  scale_size_manual(values = c(1.1,0.9,1.3,0.9)) +                                                      # 線の太さ
  scale_color_manual(values = c('#000000','#000000','#BBBBBB','#000000')) +  # 線の色
  scale_x_continuous(limits=c(1960, 2023),
                     breaks=c(1970, 1980, 1990, 2000, 2010, 2020)) + 
  theme(text = element_text(family="MS Gothic"),
        axis.text=element_text(size=15,color="black"),
        axis.title.y = element_text(size=15,color="black"),
        legend.position = 'none') + 
  labs(x="", y="\n頻度/掲載論文数\n", color = "")

• geom_lineの行で、線種、色、太さを分けるグループを指定しておく。
• theme_classicはシンプルなスタイルを選択してる設定。
• geom_text_repelを2回やってますが、1つ目はxが最大の点を利用して右側にラベルを書く作業、2つ目はxが最小の点を利用して左側にラベルを書く作業です。要するにこのグラフは、1つの折れ線と2つの散布図が重なったものってことです。
• nudgeしない場合、点の座標にそのままラベルが書かれることになり、ラベルのテキスト同士が重なる場合だけggrepelの機能でずらされて、場合によっては指示線が出るのですが、正負方向に4年分だけnudgeする設定にしておくことで、強制的にこの指示線を出すようにします。
• scale_linetype_manual、scale_size_manual、scale_color_manualは、線の種類・太さ・色を手動設定するもので、これを設定しなければ、グループ変数（ここではCODE）ごとに適当に割当られます。
• x軸のデータは1965年から2017年までしかないのですが、ラベルを表示する領域を確保するために、scale_x_continuousのlimitsで左右の範囲を少し広げています。
• legend.position = 'none'で凡例は無しにしている。

Rで棒グラフと折れ線グラフを重ねた2軸グラフを描く

さっき、Rで棒グラフと折れ線グラフを重ねたものを作ろうとして、けっこう手間取りました。最終的に描いたのは↓のようなものなのですが。

「2軸グラフの書き方」「種類の異なるグラフの重ね方」についていろいろ調べたところ、barplot()とplot()を組み合わせるやり方もあるんですが、ggplot2でやるほうがやりやすかったです。以下、まず単なる2軸グラフの書き方をおさらいした後で、棒と折れ線を組み合わせる方法をメモしておきます。
　
　

2軸グラフの作り方

左右の軸をつかってスケールの異なるグラフを重ねたいだけなら、さほど難しくはなく、plot()ですぐできます。
例えば折れ線グラフを２つ重ねたいのだとしたら、まず1つめの変数のグラフを軸無しで書いたあとに、x軸と左軸を描く。その後、par(new=T)で、2つめの変数のグラフを軸なしで重ね描きした上で、右軸を描けばいいです。
で、最後にbox()で枠線を入れ、凡例を付けたければ付けます。

以下、適当に乱数でつくった変数で実行例を書いておきますが、冒頭の成果物にあわせてx軸の変数を日付型にしてるので、そこだけ多少ややこしくなっています。ただの数字であれば、x軸を描くときにaxis.Dateを使う必要はないです。

# 練習用の変数を適当に乱数でこしらえる
x.date <- seq(from=as.Date('2020-04-01'), to=as.Date('2020-05-31'), by=1)
y1 <- sin(seq(length(x))/7)+rnorm(n=length(x), mean=0, sd=0.2) + 2
y2 <- sin(seq(length(x))/5)*100 + rnorm(n=length(x), mean=0, sd=20) + 200

# y軸の範囲をそれぞれ決めておく
y1.lim <- c(min(y1), max(y1))
y2.lim <- c(min(y2), max(y2))

# グラフの左右の余白を少し多めにするためparを設定しとく（特に右）
par(oma = c(0, 1, 0, 3))

# 1枚目のグラフをかく（y1）
# いったん軸なしにするためaxes=Fにしてる
plot(x=x.date, y=y1, ylim=y1.lim, type='l', lwd=1.5,
     xlab='Date', ylab='y1', 
     axes = F,
     main='plot関数での二軸グラフ')

# x軸を追加（日付データなのでaxis.Dateを使う）
axis.Date(1,at=seq(min(x.date), max(x.date),"week"),format="%m/%d")

# 左の軸をかく
axis(2)

# 2枚目のグラフを重ねる（y2）
par(new=T)
plot(x=x, y=y2, ylim=y2.lim, type='l', lty='dotted', lwd=1.5,
     xlab='', ylab='', axes = F)

# 右側の軸の名前をかく
mtext('y2', side = 4, line = 3)

# 右側の軸を表示
axis(4)

# 枠をかく
box()

# 凡例
legend("bottomleft", legend = c("y1", "y2"), lty = c('solid', 'dotted'), lwd=1.5)

　
　

ggplot2で棒グラフを折れ線グラフを重ねる

さて、今度は種類の異なるグラフを重ねるやり方ですが、barplot()とplot()を重ねるやり方だと、今回は日付のデータを使ってることもあって、x軸のコントロールが難しかったので、ggplot2でやることにしました。
あと、今日気づきましたが、MacのRStudioだとplotするときに余白が足りませんというエラーが出まくるのが、ggplot2だと出ないんですね。いままで基本的にplot派でしたが、ggplot派に改宗しようかなと思いました……。

さて作図ですが、半分ぐらいは、西浦博（8割おじさん）氏が5月ぐらいに出していた、新型コロナの実行再生算数を計算するプログラムの作図のところを参考にさせて頂きました。
最大のポイントは、左右の軸のスケールの調整です。
上述の「plotの2軸化」の場合、右側の軸の幅は、2つめの変数の値の幅がそのまま反映されていました。plotはそもそもスケールの異なるグラフを重ねることができるようになってて*1、今回の場合でいうと、後で描いたグラフの主軸を単に右側表示にしただけというわけです。

一方、ggplot2の場合は、2軸グラフを描くときも、y軸のスケール（軸の最小値と最大値）はあくまで共通になります。今回の場合、y軸の縦幅はあくまで、1つめの変数y1の最小〜最大の幅に合わせた尺度で固定される感じになります。
じゃあ、そこにどうやって2つ目の変数を重ねるのかというと、両変数の縮尺を先に計算しておいて、２つ目の変数を1つ目の変数にあわせて縮めたり伸ばしたりして収めるわけです。
で、その後で、右側の軸のところに好きなように（ただし左軸からの変換という形で）目盛りを打つことができるので、この目盛りを、もとの第2変数に対応するものにしておけばよい。

少し分かりづらいですが、たとえば左の軸で表現したい第1変数が0〜10ぐらいのレンジで分布していて、右の軸で表現したい第2変数が0〜100ぐらいのレンジで分布してるとすると、縮尺は1:10になるので、

1. グラフ領域はまず第1変数にあわせて、y軸が0〜10になるように描く。目盛りは左軸に表示される。
2. 第1変数をプロットする。
3. 第2変数を0.1倍して縮め、同じ領域に第2変数のグラフを重ねる。
4. sec.axisというオプションをつかって、左軸を好きなように変換した軸を右側に設定することができるので、ここで「左軸を10倍する」という変換設定をする。
5. breaksで目盛りも適当な間隔で設定する。

という手順で、2軸グラフが描かれるわけです。

以下の実行例では、さっき乱数でつくった変数をもっかい使ってるのですが、y1（左軸）y2（右軸）の縮尺を先に計算してscalerという変数に入れてあります。

# データフレームにまとめる
d <- data.frame(Date=x.date, Y1=y1, Y2=y2)

# 各軸の範囲をきめる
# yも最小値〜最大値の形で設定してよいが明示的に与えたい場合が多い気がする
x.lim  <- c(min(d$Date), max(d$Date))
y1.lim <- c(0, 4)
y2.lim <- c(0, 400)


# 左軸と右軸の関係を表すスケーラをつくる
# 各軸の最大最小差の比をとっている
scaler <- (y1.lim[2] - y1.lim[1])/(y2.lim[2] - y2.lim[1])

d %>% 
  ggplot() + 
  geom_bar(aes(x=Date, y=Y1), stat='identity', width=0.7) +
  geom_line(aes(x=Date,y=Y2*scaler, colour = "Y2"), size=1) +
  scale_x_date(date_labels="%m/%d",date_breaks="7 day", 
               limits=x.lim, expand=c(0, 0)) +
  scale_y_continuous(limit=y1.lim, expand = c(0, 0), 
                     sec.axis=sec_axis(trans = ~ ./scaler, 
                        breaks=seq(from=y2.lim[1], to=y2.lim[2], by=50), 
                        name="\nY2\n")) +
  theme(text=element_text(size=12, family="MS Gothic",color="black"),
        axis.text=element_text(size=10, family="MS Gothic",color="black"),
        legend.position="top",
        plot.subtitle=element_text(size=10, color="#666666")) + 
  labs(x="\nDate\n", y="\nY1\n", color = "",
       title='\nggplot2で重ねたグラフ', 
       subtitle='(折れ線を複数追加することもできます)')

sec.axisの中のtransというところには、左軸と右軸の対応関係をformula形式で書くのですが、「.」はデータ全体を表してて、これをスケーラで割るという変換を設定してあります。
x軸が日付なので、scale_x_dateをつかって日付表示の設定をしています。
breaksってところで、右2軸の目盛りを設定しています。水平のグリッドがある場合、左軸の目盛りと右軸の目盛りが噛み合ってたほうが綺麗なので、最初にy1とy2の幅を設定する時に、いい感じの公約数がある値を選ぶのがいいと思います。
themeでフォントを指定してるのは、日本語を文字化けなく表示させるためです。

冒頭に貼った成果物のように、折れ線を2本引きたいときは、geom_line()をもう1個プラスすればいいですね。
グラフのタイトルや軸のタイトルの前後に改行（\n）を入れているのはなんとなく隙間を開けるためです。

*1:だから逆に、重ねる時にスケールが揃ってないことを忘れたりすることがありますねｗ