★新サイト完成しました!
３秒後に自動的に移動します
変わらない方は こちらからどうぞ
http://logics-of-blue.com/r%E3%81%AA%E4%BA%88%E6%B8%AC/
統計ゼミで実演したプログラムが載ってます。
『一般線形モデルによる生物科学のための現代統計学―あなたの実験をどのように解析するか』の内容と、Rの使い方に関するものです。サンプルデータもこの本の物です。
とりあえずRのスクリプトに張り付けてF5を押せば実行できます。
#わたしのプレゼンの内容
#第2部 計算例 の計算
1-pf(5,2,12)
data<-data.frame(Tl=c(13,11,6,11,10,7,7,5,8,7,5,5,4,3,3),Factor=c(rep("A",3),rep("B",5),rep("C",7)))
data
aov.model<-aov(Tl~Factor,data=data)
summary(aov.model)
lm.model<-lm(Tl~Factor,data=data)
summary(lm.model)
anova(lm.model)
# 第1章 分散分析への招待
######################################
# データの読み込み クリップボードより #
######################################
# まずは、Excelなりなんなりのデータをコピーして、クリップボードに入れてやる。
# その作業をしてから以下のプログラムを実行（カーソルを合わせてF5を押す）する。
# ちゃんとCSVファイルに入れてあるんだったら、まっとうな読み込み方法で読み込めるんだけど、
# 今回はExcelで保存されてしまっているので、こんな方法をとった。
# CSVから読み込む方法は、第1回のゼミで紹介されていたので略。
data<-read.delim("clipboard")
# これで、Rという「環境」の中にdataという「変数」が読み込まれたことになる。
# 注意 基本的に、実行されたプログラムは勝手に表示されるので、「表示」させるというコマンドは普 #通いらない。
# 対話型プログラムなので、C++とかとは根本的に扱いが違う。C++より100倍簡単。
# 「定義」をしたときには「表示」コマンドがいる。けれども、「定義しましたよ」という結果はちゃんとコン #ソールに表示されている。
# 「表示」コマンドは、単に変数名を入力して実行（F5を押す）するだけ
data
# このデータのままだと、本来は「Factor」であるところのFERTILが「数値」になってしまっているので、#「要因」として変換する
d1<-data.frame(YIELD=data$YIELD,FERTIL=as.factor(data$FERTIL))
d1 # 結果の表示。見た目は変わらないのだけれども･････
class(data$YIELD) # これは数値
class(data$FERTIL) # 元データは整数として扱われている。
class(d1$FERTIL) # 要因に変換されている
data$FERTIL
d1$FERTIL # 返還前と後でちょっと違っているのが分かる
# 変数 ＄ 要素名 で、「データ内のこういう名前の要素」を表示することができるので、$とかいう変 #な記号が上にはついている
d1
d1$YIELD
d1$FERTIL
# たんに要素名だけを表示させようとしてもダメ
YIELD # 単にこれを実行すると、「 エラー： オブジェクト 'YIELD' がありません 」と怒られる
FERTIL # 上に同じ
# attach すると、「変数」内の「要素」が、Rという環境内に存在するようになる。
# だから、$ をつかう必要はなくなる
attach(d1)
YIELD # attach すると YIELR が環境内に存在するようになる。
d1$YIELD==YIELD # さっきやった$を使った結果と一致している
FERTIL # もちろんこれも、ちゃんと環境中に存在するようになってる
##########################
# 作図 #
##########################
# 2通りの図を作る
# まずは普通の散布図
# FERTILが整数として読み込まれている「data」という変数を用いる
plot(YIELD~FERTIL,data=data,xlab="肥料",ylab="収穫量",ylim=c(0,8),xlim=c(0,4)
,main="3肥料が与えられた30区画での収穫量のグラフ")
win.graph() #新たにグラフ描画画面を作れと言う指示
# 今度は、FERTILが「要素」として扱われたグラフを書く
# 勝手にボックスプロットになってくれている
plot(YIELD~FERTIL,data=d1,xlab="肥料",ylab="収穫量",ylim=c(0,8),xlim=c(0,4)
,main="3肥料が与えられた30区画での収穫量のグラフ",sub="ボックスプロット FERTILを水準にしたので")
############################
# 分散分析モデルを作る #
############################
# ANOVA1 という変数は「分散分析モデルの結果」だと定義させる
ANOVA1<-aov(YIELD~FERTIL,data=d1) # 分散分析するモデルとして「定義」されただけなので･････
summary(ANOVA1) # この関数で結果の要約を「表示」させる。
names(ANOVA1) # 今回作ったモデルの結果の「要素の名前」を「表示」させる
# この名前と、さっきやったように$を使えば、中身を表示できる
# 分かりにくいとは思うが、以下の式で、教科書p１４に載ってる「総計された標準偏差」なるものが計 #できる
# 教科書の定義どおりに計算してる
# 最初のうちは、数式を言葉で置き換えて計算ができる R の”雰囲気”が味わえれば、それだけでい#いと思う
# 数式の展開が自分でできなくても、数式の意味さえ理解できれば、簡単に計算できるんですね
# この直感的な操作方法は、すごいことだと思う
sqrt(sum(ANOVA1$residual^2)/ANOVA1$df.residual)
# 上で計算された標準偏差をつかえば、表1.4の結果も簡単に出せる
SD<-sqrt(sum(ANOVA1$residual^2)/ANOVA1$df.residual) # SDという変数に、さっきの計算結果を#閉じ込める
mean.of.YIELD<-tapply(YIELD,FERTIL,mean) # 各肥料ごとに、収穫量の平均をを計算した
mean.of.YIELD
# こんな結果になってる
t.value<-qt(p=0.975,df=ANOVA1$df.residual) # これで、t値が計算できる
t.value
# こんな結果
SE<-SD/sqrt(10) # 標準誤差の計算
SE
mean.of.YIELD-t.value*SE # 95%信頼区間 下側
mean.of.YIELD+t.value*SE # 95%信頼区間 上側
# 線形モデルを推定するための関数 lm でも、モデリングできる
lm.model<-lm(YIELD~FERTIL,data=d1)
anova(lm.model)
summary(lm.model)
############################
# 練習問題 メロン #
############################
data<-read.delim("clipboard")
d2<-data.frame(YIELDM=data$YIELDM,VARIETY=as.factor(data$VARIETY))
tapply(d2$YIELDM,d2$VARIETY,mean)
ANOVA2<-aov(YIELDM~VARIETY,data=d2)
summary(ANOVA2)
############################
# 練習問題 樹木 #
############################
data<-read.delim("clipboard")
d3<-data.frame(FLOWERS=data$FLOWERS,SEX=as.factor(data$SEX))
d3
tapply(d3$FLOWERS,d3$SEX,mean)
#ANOVA
ANOVA3<-aov(FLOWERS~SEX,data=d3)
summary(ANOVA3)
# t検定でも、結果は一緒。 デフォルトの t.test は、分散が異なる検定をかってにしちゃうので、var.equar=T にした。
t.test(FLOWERS~SEX,var.equal = T,data=d3)
plot(FLOWERS~SEX,data=d3)
