×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

HomeRな予測分散分析

分散分析


★新サイト完成しました!
3秒後に自動的に移動します

変わらない方は こちらからどうぞ
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)                # 今回作ったモデルの結果の「要素の名前」を「表示」させる
                        # この名前と、さっきやったように$を使えば、中身を表示できる

# 分かりにくいとは思うが、以下の式で、教科書p14に載ってる「総計された標準偏差」なるものが計  #できる
# 教科書の定義どおりに計算してる
# 最初のうちは、数式を言葉で置き換えて計算ができる 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)

ページトップへ
Copyright (C) 2011- 海と魚と統計解析 All Rights Reserved.