初心者 (による) セッション

SappoRo.R #8

中西大輔 (広島修道大学健康科学部)

24th Feb, 2018

自己紹介

  • 社会心理学者 (社会的影響、協力行動)
  • 広島修道大学健康科学部心理学科 (修道院とは関係ありません)
  • 経歴: 北教大函館 (学部) → (学歴ロンダ) → 北大 (博士前期・後期) → 広島修道大
  • @daihiko

Rとは

Rの利点

  • フリーで (SPSSとかSASはバカ高い)
  • 美しく (Excelのグラフは吐き気がする)
  • CUI (GUIはクソ)

Rを使うためにすること

  • SPSSを物理的に棄てる。
    • 物理的に、永久に。
  • Windowsを棄てる。
    • HADに頼ってしまうから。
    • しかし、JASPという強敵が……。

ken

Rを使っているとガンになっても大丈夫

  • 人文学部から健康科学部に配置転換になったらガンになりました。
  • 何が健康だ。

hospital01 hospital02

メラノーマ データセットを使って気を紛らわす

  • MASSパッケージのMelanomaデータセット
##   time status sex age year thickness ulcer
## 1   10      3   1  76 1972      6.76     1
## 2   30      3   1  56 1968      0.65     0
## 3   35      2   1  41 1977      1.34     0
## 4   99      3   0  71 1968      2.90     0
## 5  185      1   1  52 1965     12.08     1
## 6  204      1   1  28 1971      4.84     1

## 
##  Pearson's product-moment correlation
## 
## data:  MASS::Melanoma$thickness and MASS::Melanoma$time
## t = -3.451, df = 203, p-value = 0.0006793
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3608044 -0.1016529
## sample estimates:
##        cor 
## -0.2354087
## 
##  Welch Two Sample t-test
## 
## data:  MASS::Melanoma$time by MASS::Melanoma$ulcer
## t = 3.8629, df = 181.1, p-value = 0.0001559
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  292.1275 902.1807
## sample estimates:
## mean in group 0 mean in group 1 
##        2414.965        1817.811

Rの基本

  • RはRStudioとセットで使うと便利なので、最初からRとRStudio両方インストールする。
  • データはExcelで入力しても良いが、csv形式にする (csvが分からない人はググれください)。
    • データに日本語が入っている場合にはWindowsとMacとで対応が異なるが、you can write in English。

Rで読み込むデータ

  • ふつうのcsv形式のデータ
  • いわゆるコンマで区切られたデータ
  • 1オブザベーション1行のデータ (でなくてもいいけれど、とりあえず)
  • 欠損値には「NA」を入れる。

Rを使うときには

  • Rをインストールする (当然)
    • cran.r-project.org
  • RStudioをインストールする (これも当然)
    • www.rstudio.com

RStudioの使い方

プロジェクトを作る

  • プロジェクトを作る
    • なぜプロジェクトを作った方がよいかはまた後で自分で調べておいてください (ヒントはWorking Directory)。とにかく、プロジェクトを作る。

Working Directoryの指定

  • 要するに「ここにデータがあって、ここで作業するよ」というフォルダ (= directory) のこと。
    • Session -> Set Working Directory -> Choose Directory
    • コンソールウィンドウに「setwd(“~/Dropbox/nakanisi/2017/SappoRoR”)」とか表示されたら成功。

ソースファイルの作成

  • ソースファイルとは、データのハンドリングや分析の命令を書くための手順書みたいなもの。
    • File -> New File -> R Script
      • とりあえず、名前を付けて保存しておく。
        • File -> Save As…
        • ここでは「test」という名前を付ける。

データファイルの作成について

  • データファイルはコンマ区切りのcsv形式
  • 基本的に1オブザベーションにつき1行で入力 (そうじゃないほうがいい場合もあるがとりあえず)
  • Excelで入力した場合には
    • ファイル -> 名前を付けて保存… -> フォーマットから「CSV (コンマ区切り) (.csv)」を選んで保存する。

データをいじる

サンプルデータの構造

  • data.csv
  • 1行目は変数名 (英数字が無難)
  • 欠損値はNA (本データにはない)
  • 880人分の福島県産品についてのWeb調査データ

どんなデータか (1)

  • ID……連番
  • SEX……1が男、2が女
  • AGE……年齢
  • AREAdis……居住県の県庁所在地から福島市への距離
  • MARRIED……1が未婚、2が既婚
  • CHILD……1が子どもなし、2が子どもあり

どんなデータか (2)

  • 以下6段階
    • Undesirable 福島第一原子力発電所の事故で起こった風評被害は望ましくない現象である
    • ConservativeBuying 国の基準が守られていれば、福島第一原子力発電所の近辺で採れた作物を食べることに抵抗はない

どんなデータか (3)

  • 以下5段階
    • FoodLiteracy01 食品の安全性のことで知らないことがあると気になる
    • FoodLiteracy02 テレビの健康情報番組を見て、安易にまねしないようにしている
    • FoodLiteracy03 食品の安全性のための情報を日頃から積極的に集めている
    • FoodLiteracy04 同じ食品を毎日食べることはリスクがあるため、多くの種類の食品をバランスよく食べるようにしている
    • FoodLiteracy05 食品を買うときや食べるときは遺伝子組み換え食品であるかどうかを気にする
    • FoodLiteracy06 野菜や果物を買うときや食べるときは無農薬であるかどうかを気にする
    • FoodLiteracy07 加工食品を買うときや食べるときは食品添加物が入っているかを気にする
    • FoodLiteracy08 野菜や果物を買うときや食べるときは中国産であるかどうかを気にする
    • FoodLiteracy09 牛肉を買うときや食べるときはアメリカ産であるかどうかを気にする

どんなデータか (4)

Quiz01

  • 人工放射線と自然放射線とで、どちらが人体への悪い影響を及ぼしますか。最も適切なものを1つ選んでください。
    • 人工放射線の方が悪影響を及ぼす
    • 自然放射線の方が悪影響を及ぼす
    • 自然か人工かで違いはない
    • どれも正しくない
    • 分からない

どんなデータか (5)

Quiz02

  • 現在、原発事故による大気中や食品中の放射性物質が健康に及ぼす悪影響として何が一番おきやすいですか。最も適切なものを1つ選んでください。
    • 白血球が減少する
    • 鼻血が出る
    • 妊娠しにくくなる
    • これらのことは起こらない
    • 分からない

どんなデータか (6)

Quiz03

  • 農産物の放射能のレベルが「ただちに健康に影響があるレベルではない」という発表はどのような意味ですか? 最も適切なものを1つ選んでください。
    • 将来的には影響がある
    • 将来的には影響があるかもしれない
    • 毎日一定量、長期間食べ続けないかぎり健康に影響するレベルに達しない
    • 毎日一定量、長期間食べ続けた場合の健康に影響するレベルはわかっていない
    • 分からない

どんなデータか (7)

Quiz04

  • 人間が放射線を浴びたときの影響について「吸収線量に放射線の危険度に応じた値を掛けた」数値の単位をなんと言いますか。最も適切なものを1つ選んでください。
    • ベクレル
    • シーベルト
    • キュリー
    • グレイ
    • 分からない

どんなデータか (8)

Quiz05

  • 1マイクロ・グラムは1ミリ・グラムの何分の1ですか?最も適切なものを1つ選んでください。
    • 10分の1
    • 100分の1
    • 1000分の1
    • 1万分の1
    • 分からない

どんなデータか (9)

Quiz06

  • 原発事故によって、農産物に付着したり水道水に溶けたりした放射性物質の中で、乳幼児への健康被害が心配された、半減期が8日程度と比較的短いものは何ですか。
    • セシウム
    • ヨウ素
    • プルトニウム
    • ウラン
    • 分からない

分析編

データの読み込み

d<-read.csv("data.csv", head=TRUE)
  • data.csvをさきほど指定したWorking Directoryに移動またはコピーしよう。
  • さきほど作成したソースファイルの一行目に上の文字列を入力する。
    • (data.csvの内容を読み込み、dというデータフレームに入れなさいという意味)
    • 「データフレームとは何か」など考えてはいけない (知りたいひとは、教えてぞうさん)。

csvが嫌いな方

  • SPSSとかExcelのファイルをそのまま読み込むこともできる。やり方はググれ。

データをいじろう

  • まずはデータの全体像を知りたい。
  • summary関数でdデータフレームの全ての変数の最小値、第1四分位、中央値、平均値、第3四分位、最大値が表示される。
summary (d)
##        ID             SEX           AGE           AREAdis      
##  Min.   :  1.0   Min.   :1.0   Min.   :20.00   Min.   :   0.0  
##  1st Qu.:220.8   1st Qu.:1.0   1st Qu.:29.75   1st Qu.: 156.3  
##  Median :440.5   Median :1.5   Median :39.50   Median : 342.1  
##  Mean   :440.5   Mean   :1.5   Mean   :39.68   Mean   : 411.2  
##  3rd Qu.:660.2   3rd Qu.:2.0   3rd Qu.:49.25   3rd Qu.: 594.8  
##  Max.   :880.0   Max.   :2.0   Max.   :59.00   Max.   :1756.4  
##     MARRIED          CHILD        Undesirable   ConservatibeBuying
##  Min.   :1.000   Min.   :1.000   Min.   :1.00   Min.   :1.000     
##  1st Qu.:1.000   1st Qu.:1.000   1st Qu.:4.00   1st Qu.:2.000     
##  Median :2.000   Median :2.000   Median :5.00   Median :3.000     
##  Mean   :1.583   Mean   :1.503   Mean   :4.58   Mean   :3.422     
##  3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:6.00   3rd Qu.:4.000     
##  Max.   :2.000   Max.   :2.000   Max.   :6.00   Max.   :6.000     
##  FoodLiteracy01  FoodLiteracy02  FoodLiteracy03 FoodLiteracy04 
##  Min.   :1.000   Min.   :1.000   Min.   :1.00   Min.   :1.000  
##  1st Qu.:3.000   1st Qu.:3.000   1st Qu.:2.00   1st Qu.:3.000  
##  Median :3.000   Median :3.000   Median :3.00   Median :3.000  
##  Mean   :3.407   Mean   :3.481   Mean   :2.89   Mean   :3.307  
##  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:3.00   3rd Qu.:4.000  
##  Max.   :5.000   Max.   :5.000   Max.   :5.00   Max.   :5.000  
##  FoodLiteracy05  FoodLiteracy06  FoodLiteracy07  FoodLiteracy08 
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:2.000   1st Qu.:2.000   1st Qu.:3.000  
##  Median :3.000   Median :3.000   Median :3.000   Median :4.000  
##  Mean   :3.067   Mean   :2.893   Mean   :3.034   Mean   :3.689  
##  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:5.000  
##  Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.000  
##  FoodLiteracy09      Quiz01          Quiz02          Quiz03     
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:2.000  
##  Median :3.000   Median :3.000   Median :3.000   Median :3.000  
##  Mean   :3.076   Mean   :3.073   Mean   :3.073   Mean   :2.843  
##  3rd Qu.:4.000   3rd Qu.:5.000   3rd Qu.:5.000   3rd Qu.:4.000  
##  Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.000  
##      Quiz04         Quiz05         Quiz06     
##  Min.   :1.00   Min.   :1.00   Min.   :1.000  
##  1st Qu.:2.00   1st Qu.:3.00   1st Qu.:1.000  
##  Median :2.00   Median :4.00   Median :2.000  
##  Mean   :2.65   Mean   :3.78   Mean   :2.717  
##  3rd Qu.:5.00   3rd Qu.:5.00   3rd Qu.:5.000  
##  Max.   :5.00   Max.   :5.00   Max.   :5.000

特定の変数の平均値だけ出す

  • 「d$AGE」とは、「dというデータフレームのAGEという変数」という意味。
    • なんでそんな面倒くさいことをしなきゃいけないのかと思われるかもしれないが、複数のデータフレームを扱うこともあるのだから、しかたがない。大人になれば分かる。
mean(d$AGE)
## [1] 39.67614

男女別に年齢の平均値を出す (汎用的な方法)

  • na.rm=TRUE……欠損値は無視しろという命令。
  • function(data)……おまじない。気にしない。
  • meanをsdにすれば標準偏差が出る
by (d, d$SEX, function(d) mean(d$AGE, na.rm=TRUE))
## d$SEX: 1
## [1] 39.85455
## -------------------------------------------------------- 
## d$SEX: 2
## [1] 39.49773

男女別に年齢の平均値を出す (簡単な方法: 複数の値を返す関数には使えない)

tapply(d$AGE, d$SEX, mean, na.rm=TRUE)
##        1        2 
## 39.85455 39.49773

データの構造

  • データ (変数) にはいろいろな種類がある。よく使うのは、
    • numeric: 数値
    • integer: 整数型
    • factor: 因子型 (カテゴリカルな変数)

データ構造を確認する

str(d)
## 'data.frame':    880 obs. of  23 variables:
##  $ ID                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ SEX               : int  1 1 1 1 1 1 2 1 2 1 ...
##  $ AGE               : int  53 54 48 59 52 39 55 58 41 43 ...
##  $ AREAdis           : num  196 946 0 0 141 ...
##  $ MARRIED           : int  2 2 2 2 1 2 2 2 1 2 ...
##  $ CHILD             : int  2 2 1 2 2 2 2 2 1 2 ...
##  $ Undesirable       : int  4 6 6 5 4 5 4 4 4 4 ...
##  $ ConservatibeBuying: int  3 2 6 3 3 5 4 3 1 3 ...
##  $ FoodLiteracy01    : int  3 4 5 1 4 3 4 3 4 4 ...
##  $ FoodLiteracy02    : int  3 4 5 4 3 3 3 2 3 3 ...
##  $ FoodLiteracy03    : int  3 3 5 2 3 3 3 2 3 3 ...
##  $ FoodLiteracy04    : int  3 4 5 3 3 5 3 3 4 4 ...
##  $ FoodLiteracy05    : int  3 3 5 2 3 3 3 2 4 4 ...
##  $ FoodLiteracy06    : int  3 2 5 2 3 3 3 2 4 4 ...
##  $ FoodLiteracy07    : int  3 2 5 2 3 3 3 2 4 3 ...
##  $ FoodLiteracy08    : int  3 2 5 2 5 4 3 2 4 4 ...
##  $ FoodLiteracy09    : int  3 2 5 2 1 3 3 2 3 4 ...
##  $ Quiz01            : int  4 3 3 4 3 5 5 3 2 3 ...
##  $ Quiz02            : int  1 1 3 1 1 5 5 5 1 1 ...
##  $ Quiz03            : int  1 2 2 2 2 2 5 5 2 2 ...
##  $ Quiz04            : int  2 1 1 2 1 2 5 2 5 5 ...
##  $ Quiz05            : int  3 3 4 3 3 3 3 4 2 5 ...
##  $ Quiz06            : int  2 3 2 1 1 5 5 2 1 2 ...
  • ほとんどint型になる。小数はnum型。
  • でも、SEXとかCHILDはfactor型が正しい (なんで正しい型じゃなければいけないのかと思うかもしれませんが、分散分析のときにfactor型じゃないものを独立変数にするとうまく動かないとかそういうことがあるので、とにかく問答無用でやっておく)。

変数一覧を出力

names(d)
##  [1] "ID"                 "SEX"                "AGE"               
##  [4] "AREAdis"            "MARRIED"            "CHILD"             
##  [7] "Undesirable"        "ConservatibeBuying" "FoodLiteracy01"    
## [10] "FoodLiteracy02"     "FoodLiteracy03"     "FoodLiteracy04"    
## [13] "FoodLiteracy05"     "FoodLiteracy06"     "FoodLiteracy07"    
## [16] "FoodLiteracy08"     "FoodLiteracy09"     "Quiz01"            
## [19] "Quiz02"             "Quiz03"             "Quiz04"            
## [22] "Quiz05"             "Quiz06"

データの型の変換

  • factor関数でカテゴリカル変数に
    • sex=1が男性、sex=2が女性。でも分かりにくい。どっちがどっちだったか分からなくなるかもしれない (ふつうはコード表作成するから忘れたりしない)。忘れる前に、以下のコードを実行しよう。
    • <- は右辺を左辺に代入するという意味
  • 終わったら、「str(d)」で、かくにん!
d$ID<-factor(d$ID)
d$SEX<-factor(d$SEX, labels=list("male", "female"))
d$MARRIED<-factor(d$MARRIED,labels=c("dontmarried","married"))
tapply(d$AGE, d$SEX, mean, na.rm=TRUE)
##     male   female 
## 39.85455 39.49773
tapply(d$AGE, d$SEX, sd, na.rm=TRUE)
##     male   female 
## 11.31418 11.06576

変数の加工

  • 複数の変数の合計
d$FoodLiteracy_all<-d$FoodLiteracy01+d$FoodLiteracy02+d$FoodLiteracy03+d$FoodLiteracy04
+d$FoodLiteracy05+d$FoodLiteracy06+d$FoodLiteracy07+d$FoodLiteracy08+d$FoodLiteracy09
hist(d$FoodLiteracy_all)

条件処理

d$Quiz01A<-ifelse(d$Quiz01==3,1,0)
d$Quiz02A<-ifelse(d$Quiz02==4,1,0)
d$Quiz03A<-ifelse(d$Quiz03==3,1,0)
d$Quiz04A<-ifelse(d$Quiz04==1,1,0)
d$Quiz05A<-ifelse(d$Quiz05==3,1,0)
d$Quiz06A<-ifelse(d$Quiz06==2,1,0)
d$Quiz<-d$Quiz01A+d$Quiz02A+d$Quiz03A+d$Quiz04A+d$Quiz05A+d$Quiz06A
hist(d$Quiz)

散布図 (1)

plot(d$Quiz,d$FoodLiteracy_all)

散布図 (2)

  • ggplot2パッケージを使う (インターネット接続必要)。
    • Tools -> Install Packages…
    • Packagesの欄にggplot2と入力しinstall
library("ggplot2")
ggplot(d,aes(x=Quiz, y=FoodLiteracy_all, colour=SEX)) + geom_jitter(size=2)

相関係数

cor.test(d$Quiz,d$FoodLiteracy_all)
## 
##  Pearson's product-moment correlation
## 
## data:  d$Quiz and d$FoodLiteracy_all
## t = -0.057895, df = 878, p-value = 0.9538
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06803188  0.06414125
## sample estimates:
##          cor 
## -0.001953845
cor.test(d$Quiz,d$FoodLiteracy_all,method="spearman")
## Warning in cor.test.default(d$Quiz, d$FoodLiteracy_all, method =
## "spearman"): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  d$Quiz and d$FoodLiteracy_all
## S = 111130000, p-value = 0.523
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02156133

総当たりの相関を求めたい

  • psychパッケージのcorr.test関数を使う (インターネット接続必要)。
    • Tools -> Install Packages…
    • Packagesの欄にpsychと入力しinstall
library(psych)
corrvar<-data.frame(d$AGE,d$AREAdis,d$ConservatibeBuying,d$Undesirable,d$FoodLiteracy_all,d$Quiz)
corr.test(corrvar)
## Call:corr.test(x = corrvar)
## Correlation matrix 
##                      d.AGE d.AREAdis d.ConservatibeBuying d.Undesirable
## d.AGE                 1.00      0.01                 0.01          0.01
## d.AREAdis             0.01      1.00                -0.07         -0.08
## d.ConservatibeBuying  0.01     -0.07                 1.00          0.29
## d.Undesirable         0.01     -0.08                 0.29          1.00
## d.FoodLiteracy_all    0.07     -0.02                -0.10          0.12
## d.Quiz                0.08     -0.10                 0.14          0.14
##                      d.FoodLiteracy_all d.Quiz
## d.AGE                              0.07   0.08
## d.AREAdis                         -0.02  -0.10
## d.ConservatibeBuying              -0.10   0.14
## d.Undesirable                      0.12   0.14
## d.FoodLiteracy_all                 1.00   0.00
## d.Quiz                             0.00   1.00
## Sample Size 
## [1] 880
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##                      d.AGE d.AREAdis d.ConservatibeBuying d.Undesirable
## d.AGE                 0.00      1.00                  1.0          1.00
## d.AREAdis             0.82      0.00                  0.3          0.11
## d.ConservatibeBuying  0.82      0.05                  0.0          0.00
## d.Undesirable         0.69      0.01                  0.0          0.00
## d.FoodLiteracy_all    0.03      0.49                  0.0          0.00
## d.Quiz                0.02      0.00                  0.0          0.00
##                      d.FoodLiteracy_all d.Quiz
## d.AGE                              0.21   0.14
## d.AREAdis                          1.00   0.03
## d.ConservatibeBuying               0.03   0.00
## d.Undesirable                      0.00   0.00
## d.FoodLiteracy_all                 0.00   1.00
## d.Quiz                             0.95   0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

複数の相関係数を眺める

library(corrplot)
corplot1<-cor(corrvar)
corrplot(corplot1)

t検定

  • Rの不親切なt.test関数で消耗してみる (1)
t.test(d$Quiz~d$SEX)
## 
##  Welch Two Sample t-test
## 
## data:  d$Quiz by d$SEX
## t = 3.6787, df = 862.28, p-value = 0.0002489
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1378203 0.4530888
## sample estimates:
##   mean in group male mean in group female 
##             1.454545             1.159091
tapply(d$Quiz, d$SEX, sd, na.rm=TRUE)
##     male   female 
## 1.269128 1.107917
table(d$SEX)
## 
##   male female 
##    440    440

t検定

  • Rの不親切なt.test関数で消耗してみる (2)
library(compute.es)
result_t<-t.test(d$Quiz~d$SEX)
result_sd<-tapply(d$Quiz, d$SEX, sd, na.rm=TRUE)
result_table<-table(d$SEX)
mes(result_t$estimate[1],result_t$estimate[2],result_sd[1],result_sd[2],result_table[1],result_table[2])
## Mean Differences ES: 
##  
##  d [ 95 %CI] = 0.25 [ 0.12 , 0.38 ] 
##   var(d) = 0 
##   p-value(d) = 0 
##   U3(d) = 59.79 % 
##   CLES(d) = 56.96 % 
##   Cliff's Delta = 0.14 
##  
##  g [ 95 %CI] = 0.25 [ 0.12 , 0.38 ] 
##   var(g) = 0 
##   p-value(g) = 0 
##   U3(g) = 59.79 % 
##   CLES(g) = 56.95 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.12 [ 0.06 , 0.19 ] 
##   var(r) = 0 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = 0.12 [ 0.06 , 0.19 ] 
##   var(z) = 0 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 1.57 [ 1.23 , 2 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = 0.45 [ 0.21 , 0.69 ] 
##   var(lOR) = 0.02 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = 13.09 
##  Total N = 880

関数の作成

便利な (でも不完全な) t検定の関数を作る

  • t.test関数が不親切だって? だったら関数作ればいいべや。
daihiko_t<-function(x,y,paired=TRUE) {
  if (paired) {
    d<-data.frame(x,y)
    d<-subset(d,complete.cases(d))
    rt<-t.test(d$x,d$y,paired=TRUE)
    rsd<-c(sd(d$x,na.rm=T),sd(d$y,na.rm=T))
    cohend<-abs(mean(d$x-d$y)/(sd(d$x-d$y)/sqrt(2*(1-cor(d$x,d$y)))))
    return(list(meanx=mean(d$x),meany=mean(d$y),SD=rsd,t=rt,cohend=cohend))
  }
  else {  
    rt<-t.test(x~y)
    rsd<-tapply(x,y,sd,na.rm=TRUE)
    rt1<-table(y,x)
    rt2<-apply(rt1,1,sum)
    cohend<-(rt$estimate[1]-rt$estimate[2])/sqrt(((rt2[1]-1)*rsd[1]^2+(rt2[2]-1)*rsd[2]^2)/(rt2[1]+rt2[2]-2))
    cohend<-unname(cohend)
    return(list(N=rt2,SD=rsd,t=rt,cohend=cohend))
  } 
}

便利な (でも不完全な) t検定の関数を使う (1)

daihiko_t(d$Quiz,d$SEX,paired=FALSE)
## $N
##   male female 
##    440    440 
## 
## $SD
##     male   female 
## 1.269128 1.107917 
## 
## $t
## 
##  Welch Two Sample t-test
## 
## data:  x by y
## t = 3.6787, df = 862.28, p-value = 0.0002489
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1378203 0.4530888
## sample estimates:
##   mean in group male mean in group female 
##             1.454545             1.159091 
## 
## 
## $cohend
## [1] 0.2480201

便利な (でも不完全な) t検定の関数を使う (2)

daihiko_t(d$FoodLiteracy01,d$FoodLiteracy02,paired=TRUE)
## $meanx
## [1] 3.406818
## 
## $meany
## [1] 3.480682
## 
## $SD
## [1] 1.0154434 0.9597542
## 
## $t
## 
##  Paired t-test
## 
## data:  d$x and d$y
## t = -1.8344, df = 879, p-value = 0.06694
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.152893991  0.005166718
## sample estimates:
## mean of the differences 
##             -0.07386364 
## 
## 
## $cohend
## [1] 0.07473953