今日の演習で使う関数の説明

データ読み込み

はじめに第５回で使用したデータを読み込み、その中でtunaのデータを使用。

Ex <- read.csv("Ex.csv")

このデータセットにはtunaとmedakaの身長が格納されているので、tunaのデータだけをsubset関数を利用して抽出する。

td <- subset(Ex, species=="tuna")　# species列にmedakaと書かれている行だけを抽出

ヒストグラム

ヒストグラムを描画。

hist(td$bl, xlim=c(0,500), breaks=seq(0,500,by=5), xlab="Body length (cm)", main="tuna (population)")

算術平均

avg_total <- round(mean(td$bl), 2) # 全員の身長の算術平均

cat("算術平均(マグロ)：", avg_total, "cm\n")

## 算術平均(マグロ)： 90.19 cm

中央平均

中央平均はmean関数のオプションをtrim=0.25とすることで算出できる。

tavg_total <- round(mean(td$bl, trim=0.25), 1) # 全員の身長の算術平均

cat("中央平均（マグロ）：", tavg_total, "cm\n")

## 中央平均（マグロ）： 81.7 cm

母標準偏差

R標準のsd関数は不偏標準偏差を算出してしまうので、前回と同様、母標準偏差を求めるmsd関数を作成して計算した。

msd <- function(x){
  N <- length(x)
  freedom <- N -1
  mvar_x <-(freedom/N)*var(x)
  msd_x <- sqrt(mvar_x)
  return(msd_x)
}

msd_total <- round(msd(td$bl), 2) # 全員の身長の母標準偏差 

cat("母標準偏差（マグロ）：", msd_total, "\n")

## 母標準偏差（マグロ）： 51.91

変動係数

変動係数は母標準偏差÷平均値。従って上で求めたmsd_totalをつかえば以下のようにシンプルに読みやすく出来る（これを毎回全て最初から式を書いていると，括弧などでタイポミスをしやすくなる。自分で読みやすくなるように書くのが大事）。

cv_total <- round(msd_total / mean(td$bl), 2) # 全員身長の変動係数

cat("変動係数（マグロ）：", cv_total, "\n")

## 変動係数（マグロ）： 0.58

標準偏差も変動係数もどちらもデータのバラツキを数量化して表してくれる指標ですが，変動係数のほうは単位の異なるデータのばらつきや、平均値に対するデータとばらつきの関係を相対的に評価する際に用いると便利です。具体的にどのような時に便利なのかは，以下のurl https://bellcurve.jp/statistics/course/5929.html を参照のこと。

無作為抽出とヒストグラムの描画

以下のコードでは２つのヒストグラムが描画される。上は母集団、下は300回の無作為抽出結果のヒストグラム。

total_sample <- sample(td$bl, 30)
par(mfrow=c(2,1))
hist(td$bl, xlim=c(0,500), breaks=seq(0, 500, 20), xlab="Body length (cm)", main="population")
hist(total_sample, xlim=c(0,500), breaks=seq(0, 500, 20), xlab="Body length (cm)", main="sample (n =30)")

結果をみると多くの場合，30標本もあれば多くの場合、母集団の性質をほぼ反映じたかたちで標本抽出ができることができる（Rチャンクの実行を繰り返してみるとときに、母集団の分布とは異なる分布のヒストグラムが見える時があるでしょう）。

標本の算術平均・標準誤差

# 算術平均
savg_total <- round(mean(total_sample), 1) 

#標準誤差
se_total <- round(sd(total_sample)/sqrt(length(total_sample)),2)

#計算結果の表示
cat ("マグロ標準体長　\n")

## マグロ標準体長

cat ("母平均 = ", round(mean(td$bl), 2), "cm\n")

## 母平均 =  90.19 cm

cat ("標本平均 ± 標準誤差 ＝", savg_total, " ± ", se_total, "cm\n")

## 標本平均 ± 標準誤差 ＝ 90.7  ±  10.52 cm

となり，ほぼ母集団の体長の平均値と近い値が得られる。そこでこの標本平均が母平均値をどの程度反映しているのかを区間推定によってしらべる。

母平均の区間推定

信頼度に応じて信頼区間の上限下限を計算する関数(get_confintと命名）を定義する。

演習内の説明で，ｋの値は標本のデータ数によって変わるがデータ数が30以上の場合は正規分布に従うと考えてよい。kの値は正規分布の確率密度（確率変数がある範囲の値をとる確率を積分した物）を求めるqnorm関数によって計算することが出来る。

係数の算出

例えば95％信頼区間の係数ｋの値は，

Q_95 <- qnorm((1 - 0.95) / 2, lower.tail=FALSE)
round(Q_95, 2)

## [1] 1.96

と求めることが出来る。

信頼区間を算出する関数get_confintを作成

上述のqnorm関数を利用して標本平均の信頼区間を計算してくれる関数get_confintを作ってみると以下の様になる。

get_confint <- function(x, level = 0.95) {
  ## 標準正規分布を利用して信頼区間を求める
  ## 匹数：x = 推定に使う標本（サンプル）
  ##       level = 信頼度。何も入力しない（デフォルト）と0.95で計算する
  ## 返り値：信頼区間の下限値，点推定による平均値，信頼区間の上限値の３つを含むベクトル
  N <- length(x)
  mn_x <- mean(x)
  SE <- sd(x) / sqrt(N)
  
  # k値を計算
  k <- qnorm((1 - level) / 2, lower.tail = FALSE)
  
  # 信頼区間の下限値・上限値
  lb <- mn_x - k * SE
  ub <- mn_x + k * SE
  
  # 算出した値をベクトルestimatesに格納
  estimates <- c(round(lb, 2), round(mn_x, 1), round(ub, 2)) 
  
  # estimatesに入れた要素に名前をつけている
  names(estimates) <- c("下限", "平均（cm）", "上限")
  
  return(estimates)
}

この関数が意図したとおりに動くか試してみると、

全体

cat ("全体の90%,95%,99%信頼区間\n")

## 全体の90%,95%,99%信頼区間

get_confint(total_sample, 0.9) #total=sampleの90%信頼区間を求める

##       下限 平均（cm）       上限 
##      73.39      90.70     107.98

get_confint(total_sample, 0.95)

##       下限 平均（cm）       上限 
##      70.07      90.70     111.29

get_confint(total_sample, 0.99)

##       下限 平均（cm）       上限 
##      63.60      90.70     117.77

今日の課題

同時期に捕獲されたマサバScomber japonicusとマアジTrachurus japonicusの体長データを用いて，以下の課題１・２を実施してください。

csvファイルを使ったデータの読み込み

read.csvコマンドを使用してcsvファイル(Lesson6.csv)からマサバとマアジの体長データをFishという名称のデータセットにして読み込めるようにしてあります。これにより以降の操作ではマアジの体長は”Fish＄maaji”，マサバの体長は”Fish＄masaba”と指定することで解析に使用することが出来るようになるようになっています。

Fish <- read.csv("Lesson6.csv")

課題１

マアジ，マサバの体長データを母集団として扱った場合、

ヒストグラムの描画

下のRチャンク内にマアジとマサバの体長分布を比較しやすいようにX軸幅・階級幅を揃えたヒストグラムを作成し、左右に並べて表示する命令文を記述せよ。

hist(Fish$masaba,breaks=seq(0,50,2), xlab="Body Length(cm)",main="masaba Body length(n=181)")

hist(Fish$maaji,breaks=seq(0,50,2),xlab="Body Length(cm)",main="maaji Body length(n=181)")

算術平均

算術平均を計算・表示する命令文をこの下のRチャンク内に記述せよ。

cat("masaba算術平均=",mean(Fish$masaba))

## masaba算術平均= 27.27845

cat(",maaji算術平均=",mean(Fish$maaji))

## ,maaji算術平均= 20.56409

中央平均

中央平均を計算・表示する命令文をこの下のRチャンク内に記述せよ。

cat("masaba中央平均=",mean(Fish$masaba, trim=0.25))

## masaba中央平均= 27.06813

cat(",maaji中央平均=",mean(Fish$maaji, trim=0.25))

## ,maaji中央平均= 20.51319

母標準偏差

母標準偏差を計算・表示するする命令文をこの下のRチャンク内に記述せよ。

N<- length(Fish$masaba)
N<- length(Fish$maaji)
freedom<- N -1
seldom<- N -1
msd<- sqrt((freedom/N)*var(Fish$masaba))
msd2<- sqrt((seldom/N)*var(Fish$maaji))
cat("masaba母標準偏差=",round(msd,1))

## masaba母標準偏差= 7.5

cat(", maaji母標準偏差=",round(msd2,1))

## , maaji母標準偏差= 4.3

変動係数

変動係数（母標準偏差/平均値)を計算・表示する命令文をこの下のRチャンク内に記述せよ。

cat("mssaba変動変数=",mean(Fish$masaba)/round(msd,1))

## mssaba変動変数= 3.637127

cat(",maaji変動変数=",mean(Fish$maaji)/round(msd2,1))

## ,maaji変動変数= 4.782346

課題１の考察

ここまでで求めた代表値を用いてマアジとマサバの体長分布について比較・考察せよ。母標準偏差の値がマアジと比べて、マサバの方が大きく、体長のばらつきが多いことが考察できる。また、変動変数に注目して、マアジの方がマサバよりも変動変数が大きいため、平均値に対して、データのばらつきが多いことが分かり、マサバはヒストグラム上で山の裾が広く、マアジは中央値の付近に平均値が来るが、極端に体長が短かったり、長かったりする個体がいることが分かった。

課題２

標本抽出とヒストグラム作成

マサバの体長データについて，30個の無作為抽出を行いベクトルに格納。その上でマサバの母集団および無作為抽出標本のヒストグラムを、左右に並べて比較できるように(par関数もしくはlayout関数を使用して）描画せよ。

masaba30 <- sample(Fish$masaba, 30)
par(mfrow=c(2,1))
hist(Fish$masaba, xlim=c(0,50), breaks=seq(0, 50, 2), xlab="Body length (cm)", main="population")
hist(masaba30, xlim=c(0,50), breaks=seq(0, 50, 2), xlab="Body length (cm)", main="masaba30 (n =30)")

標本の算術平均・標準誤差

標本の算術平均・標準誤差を計算するためのRチャンクをこの下に挿入せよ。

なお、算術平均と標準誤差の値はcatコマンドを利用して「算術平均　±　標準誤差　（cm)」という標記で出力すること。

cat("masaba30標本算術平均=", mean(masaba30))

## masaba30標本算術平均= 27.61333

cat(", masaba30標本算術誤差=", sd(masaba30)/sqrt(30))

## , masaba30標本算術誤差= 1.201357

母平均の区間推定

注：計算には今日の課題に必要な関数の説明部分で定義したget_confint関数を使用することができる。

信頼度９０％の場合

信頼度90％の場合の信頼区間の下限・上限を算出するRチャンクを記入せよ。

get_confint(masaba30,0.9)

##       下限 平均（cm）       上限 
##      25.64      27.60      29.59

信頼度９５％の場合

信頼度95％の場合の信頼区間の下限・上限を算出するRチャンクを記入せよ。

get_confint(masaba30,0.95)

##       下限 平均（cm）       上限 
##      25.26      27.60      29.97

信頼度９９％の場合

信頼度99％の場合の信頼区間の下限・上限を算出するRチャンクを記入せよ。

get_confint(masaba30,0.99)

##       下限 平均（cm）       上限 
##      24.52      27.60      30.71

課題２の考察

課題１で求めた標準偏差と，課題２で求めた標準誤差の意味の違いを説明せよ。また、区間推定の信頼度を変化させたことで何が生じたかを説明する共に，その意味を考察せよ。

標準偏差はマサバの体長データそのものに注目し、データのばらつきの大きさを示した値であるが、標準誤差はデータの平均値に注目して算出し、平均値とそれぞれのデータのズレの大きさを示している。また、信頼度を上げることで、データの平均値は変化しないが、下限、上限の幅が大きくなった。信頼度を上げることで、母平均が区間内に入る確率が上がり、推測の正確さを上げることができる。

ここまでの内容に関して感想や質問があったらこの行の下に記述して下さい（阿部）

Ｒ操作が難しすぎて、困難を極めました。

生物情報処理演習　第６回　課題（グラフで日本語利用不可バージョン）

092320095・北中裕登

2024-11-14

今日の演習で使う関数の説明

データ読み込み

ヒストグラム

算術平均

中央平均

母標準偏差

変動係数

無作為抽出とヒストグラムの描画

標本の算術平均・標準誤差

母平均の区間推定

係数の算出

信頼区間を算出する関数get_confintを作成

全体

今日の課題

csvファイルを使ったデータの読み込み

課題１

ヒストグラムの描画

算術平均

中央平均

母標準偏差

変動係数

課題１の考察

課題２

標本抽出とヒストグラム作成

標本の算術平均・標準誤差

母平均の区間推定

信頼度９０％の場合

信頼度９５％の場合

信頼度９９％の場合

課題２の考察

ここまでの内容に関して感想や質問があったらこの行の下に記述して下さい（阿部）

生物情報処理演習 第６回 課題（グラフで日本語利用不可バージョン）

092320095・北中裕登

2024-11-14

今日の演習で使う関数の説明

データ読み込み

ヒストグラム

算術平均

中央平均

母標準偏差

変動係数

無作為抽出とヒストグラムの描画

標本の算術平均・標準誤差

母平均の区間推定

係数の算出

信頼区間を算出する関数get_confintを作成

全体

今日の課題

csvファイルを使ったデータの読み込み

課題１

ヒストグラムの描画

算術平均

中央平均

母標準偏差

変動係数

課題１の考察

課題２

標本抽出とヒストグラム作成

標本の算術平均・標準誤差

母平均の区間推定

信頼度９０％の場合

信頼度９５％の場合

信頼度９９％の場合

課題２の考察

ここまでの内容に関して感想や質問があったらこの行の下に記述して下さい（阿部）

生物情報処理演習　第６回　課題（グラフで日本語利用不可バージョン）