No.15800 Re: 重回帰分析について 【青木繁伸】 2011/11/28(Mon) 09:51
文字通り，一次従属な独立変数が含まれていたので除いたということです。
その独立変数を除外すればよいだけのことです。
No.15806 Re: 重回帰分析について 【統計学初心者】 2011/11/29(Tue) 12:43
青木先生
早速のご返答をありがとうございます。
なるほど，説明変数どうしでは一次独立でなければならなかったですね。
だから，一次従属な説明変数を除外すれば解決できるのですが，今36個のダミー変数を使っていますので，37個の回帰係数を求めたいのです。（1つを基準としてそれの回帰係数は0です。）
例えば，東京都と大阪府と愛知県の3つのダミー変数を作るときに，東京都を基準にしてダミー変数を東京都を（0，0），大阪府を（1，0），愛知県を（0，1）とするということです。
一次従属の独立変数を一次独立の独立変数にすることは，不可能ですよね？
ダミー変数の数を減らすしか方法はないのでしょうか？
No.15809 Re: 重回帰分析について 【青木繁伸】 2011/11/29(Tue) 15:09
「今36個のダミー変数」って，それら同士は独立でしょう?しかし，なんかへん。
例えば，東京のデータは数個しかないということでは?
データの個数が97個しかないのに，36個もダミー変数を作ったら，たまたま一次従属になってしまうこともあるでしょう。それらの中のいくつかと他のダミー変数（使っているのかどうか分からないけど）が従属になっているのでは?
> ダミー変数の数を減らすしか方法はないのでしょうか
それはそうでしょう。R が自動的に一次従属な変数を除去してくれているのに，どんなことをやったって，独立な変数になるわけがありません。
No.15814 Re: 重回帰分析について 【統計愛好家】 2011/11/30(Wed) 04:02
青木先生
ご返答ありがとうございます。
青木先生のおっしゃる通り，ある地域のデータは1個のところもあれば，2個，3個あるところもあります。
ただ，地域性を調べようと思っているため，出来る限り地域を多くしようとした結果，このように地域間のデータ数が異なってしまいました。
ですから，一度一次従属になってNAが出てしまったものについて，除去してもう一度出力してみようと思います。
また，何かありましたら，その時は宜しくお願い致します。
この度は本当にありがとうございました。
No.15860 Re: 重回帰分析について 【統計愛好家】 2011/12/03(Sat) 16:43
お久しぶりです。
あれからNAの出たダミー変数を説明変数からはずして重回帰分析を行うと結果が出たのですが，少ししてからこのようなことに気付きました。
例えば，
ダミー1 ダミー2 ダミー3 ダミー4
東京 0 0 0 0
名古屋 1 0 0 0
神戸 0 1 0 0
京都 0 0 1 0
大阪 0 0 0 0
とダミー変数を置いて，ダミー4がNAになったとします。
そこで，説明変数からダミー4を除くと上手くはいくのですが，大阪が東京と一緒になってしまうので，おかしいことに気付きました。
そして，それを回避するために，大阪の行データそのものをなくして，東京と名古屋と神戸と京都に減らしました。
しかし，これで重回帰分析を行っても，NAが出たので同じことを繰り返すと説明変数がほとんどなくなってしまい，おかしいので，どのように対処すればよいかを教えていただきたいです。
お願い致します。
No.15861 Re: 重回帰分析について 【青木繁伸】 2011/12/03(Sat) 17:14
> それを回避するために，大阪の行データそのものをなくして
そのようにするのは間違って います。ケース（データ）を削除するのは無意味なことです。そして，そのようにせざるを得ないのは，ケースを表すダミー変数を使うということが目的にか なっていないということです。全体の中で，ダミー変数が1になるケースが数パーセントというようなダミー変数はあまり予測効率が高くない。ダミー変数が特 殊な状況を反映する（例えば，「東京」であること）ような場合，予測値の特殊な変動だけを説明する事になってしまうでしょう。
ではどうするか？地域が予測値に影響を及ぼすと言うことが，「東京である」とか[大阪である」ということではなく，「大都市を多く含むか」とか，「関西か関東か」というようなダミー変数を創るべきだということでしょう。
例えば，ダミー変数1 ダミー変数2 従属変数のようなデータがあって，予測式が
東京 0 0 3
名古屋 1 0 2
大阪 0 1 5(Intercept) ダミー変数1 ダミー変数2だったら，ガッカリでしょう？
3 -1 2
No.15865 Re: 重回帰分析について 【統計愛好家】 2011/12/04(Sun) 13:34
青木先生
ご返信ありがとうございます。
全体の中で，ダミー変数が1になるデータが数パーセントというダミー変数はよくないことは初めて知りました。
これから，「大都市を多く含むかどうか」や「関西か関東かどうか」などの，大枠のダミー変数をどのように作るかを考えていこうと思います。
この度はありがとうございました。
No.15950 重回帰分析について 【統計大好き】 2011/12/12(Mon) 18:27
お久しぶりです。
あれから色々と行い，なんとか地域を分けることに成功しました。
そして，地域ダ ミーを（地域数-1）個作り，それを説明変数の中に入れて，総当たり法を使って変数選択を行いたいのですが，全ての説明変数での総当たりは青木先生の作成 されたプログラムを用いてRで出来るのですが，今回は地域ダミーは固定させて，それ以外の説明変数+地域ダミー全てで総当たりをさせたいのですが，プログ ラム（方法）が全く分かりません。
また，日本だけではなく，世界についても調べているのですが，世界の何十カ国かのデータを白地図を使って色分けしたり（目的変数），その上に説明変数のデータを棒グラフとして載せたいのですが，これも全く分かりません。
青木先生のHPには日本の白地図で色分けすることができるのですが，それの世界バージョンで，しかも，色分け＋棒グラフと重ね合わせすることはできないでしょうか？
ご教授お願い致します。
No.15953 Re: 重回帰分析について 【青木繁伸】 2011/12/12(Mon) 20:25
> 今回は地域ダミーは固定させて，それ以外の説明変数+地域ダミー全てで総当たりをさせたいのですが，プログラム（方法）が全く分かりません。
プログラムを改編して，固定すべき変数は固定し，総当たりすべき変数部分を総当たりにして，というようにすればよいのです。
str <- name[bincomb[i,]] # どの独立変数が使われるかを割り出す
の次にでも，固定すべき変数名を"+"で繋いだ文字列を付け加えるだけ。例えば，
str <- paste(str, "+foo+bar+baz", sep="")
などと。プログラミングの基礎がないと，難しいかな。
> 青木先生のHPには日本の白地図で色分けすることができるのですが，それの世界バージョンで，しかも，色分け＋棒グラフと重ね合わせすることはできないでしょうか？
イージー・オーダーの関数はありません。ないものは（必要と思う人が）作るしかないのが現実です。
No.15962 Re: 重回帰分析について 【統計大好き】 2011/12/13(Tue) 08:56
青木先生
ご返信ありがとうございます。
総当たりの方に関してですが，青木先生のおっしゃる通りに
str <- paste(str, "+foo+bar+baz", sep="")
を付け加えて，出力してみると，地域ダミーが2個も3個も結果として出てきてしまったので，何かがおかしいように思うのですが，分かりますか？
地図に関してですが，R以外にはそのようなものがあるソフトはご存知ですか？
No.15963 Re: 重回帰分析について 【青木繁伸】 2011/12/13(Tue) 09:30
> 何かがおかしいように思うのですが，分かりますか？
まさか，str <- paste(str, "+foo+bar+baz", sep="") をそのまま挿入したのですか?
ほかにも，いろいろ調整しないと無理ですよ。
> 地図に関してですが，R以外にはそのようなものがあるソフトはご存知ですか？
知りません。
No.15967 Re: 重回帰分析について 【統計大好き】 2011/12/13(Tue) 12:33
青木先生
ご回答ありがとうございます。
他にも色々と調整しなければいけないのですね。分かりました。
色々と試してみます。
地図に関してはパソコンではできないようなので手作業で作るなりすることにします。
ありがとうございました。
No.16009 Re: 重回帰分析について 【統計大好き】 2011/12/16(Fri) 13:43
青木先生
あれからstr <- paste(str, "+foo+bar+baz", sep="")を付け加えて，色々と調整しようと試みたのですが，ベクトルの長さが定数倍になってないだとか，色々とエラーがでてしまい，できていません。
どこを調整するべきなのでしょうか？
No.16010 Re: 重回帰分析について 【青木繁伸】 2011/12/16(Fri) 14:59
> どこを調整するべきなのでしょうか？
それが分からないなら，無理ですね。
No.16011 Re: 重回帰分析について 【青木繁伸】 2011/12/16(Fri) 15:49
以下のようにすればよいでしょう。
関数の第二引数に，必ずモデルに取り入れる変数名のベクトルを指定します。
最後に使用例を付けておきます。# 総当たり法による重回帰分析を行う
# データフレームには，分析に使用する独立変数と従属変数のみを含むこと。
# また，従属変数は最終列に置くこと。
#
All.possible.subset.selection <- function( df, # データフレーム（独立変数，従属変数）
fixed, ### 追加：必ず取り入れる変数の名前（ベクトル）
limit=10) # 独立変数の個数の上限（数が多いと計算時間が指数的に増える）
{
df <- subset(df, complete.cases(df)) # 欠損値を持つケースを除く
nv <- ncol(df)-1-length(fixed) ### 変更 独立変数の個数
if (nv > limit) { # limit より多いと分析を中断する
stop(paste("独立変数が", limit,
"個以上である（多すぎる）。\n",
"limit 引数で変更できる", paste=""))
}
n <- 2^nv # 独立変数を取り出す取り出し方
bincomb <- matrix(FALSE, nrow=n, ncol=nv) # e1071 パッケージの bincombinations より
for (j in 1:nv) {
bincomb[, j] <- rep(c(rep(FALSE, n/2^j), rep(TRUE, n/2^j)), length = n)
}
bincomb <- bincomb[-1,]
n <- n-1
name <- names(df) # 変数名を取り出す
depname <- name[ncol(df)] ### 変更
name <- setdiff(name[1:(ncol(df)-1)], fixed) ### 変更
result4 <- character(n) # 数値型ベクトル確保
result1 <- result2 <- result3 <- numeric(n) # 数値型ベクトル確保
for (i in 1:n) { # 独立変数の全ての組み合わせについて，
str <- c(name[bincomb[i,]], fixed) ### 変更 どの独立変数が使われるかを割り出す
form <- reformulate(str, depname) # モデル式を作る（"formula" クラス）
ans <- lm(form, df) # 重回帰分析の結果
result <- summary(ans)
result1[i] <- result$r.square # 重相関係数の二乗（決定係数）
result2[i] <- result$adj.r.square # 自由度調整済み重相関係数の二乗
result3[i] <- AIC(ans) # AIC
temp <- as.character(form) # モデル式を文字列に変換
result4[i] <- paste(temp[2], "~", temp[3]) # モデル式を記録
}
return(structure(list(rsq=result1, adj=result2, AIC=result3, form=result4),
class="all.possible.subset.selection"))
}
print.all.possible.subset.selection <- function( obj, # "all.possible.subset.selection" クラスのオブジェクトをプリント
sort.by=c("adj", "rsq", "AIC"), # 結果を何で並べ替えるかを指示
models=20) # 良い方から何番目まで出力するか
{
result <- data.frame(obj$rsq, obj$adj, obj$AIC, obj$form)
sort.by <- match.arg(sort.by)
o <- order(switch (sort.by, "rsq"=result[,1], "adj"=result[,2], "AIC"=result[,3]), decreasing = sort.by != "AIC")
result <- result[o,]
cat("\nR square Adjusted R square AIC Formula\n") # 表頭
models <- min(models, nrow(result))
apply(result[1:models,], 1, function(x) # 各行の出力
cat(sprintf("%8.5f %13.5f %13.5f %s\n",
as.double(x[1]), as.double(x[2]), as.double(x[3]), x[4])))
invisible(result)
}
# 使用例
data(fat, package="UsingR")
df <- fat[, c(7,10:19, 6)]
All.possible.subset.selection(df, fixed=c("height", "hip"), limit=20)
# 実行結果
R square Adjusted R square AIC Formula
0.96875 0.96745 1568.61700 weight ~ neck + chest + abdomen + thigh + knee + ankle + bicep + wrist + height + hip
0.96882 0.96739 1570.07300 weight ~ neck + chest + abdomen + thigh + knee + ankle + bicep + forearm + wrist + height + hip
0.96852 0.96735 1568.46900 weight ~ neck + chest + abdomen + thigh + knee + ankle + bicep + height + hip
0.96862 0.96732 1569.63700 weight ~ neck + chest + abdomen + thigh + knee + ankle + bicep + forearm + height + hip
0.96830 0.96712 1570.23300 weight ~ neck + chest + thigh + knee + ankle + bicep + wrist + height + hip
0.96807 0.96702 1570.03900 weight ~ neck + chest + thigh + knee + ankle + bicep + height + hip
0.96832 0.96701 1572.02000 weight ~ neck + chest + thigh + knee + ankle + bicep + forearm + wrist + height + hip
0.96812 0.96694 1571.63500 weight ~ neck + chest + thigh + knee + ankle + bicep + forearm + height + hip
0.96807 0.96688 1572.05800 weight ~ neck + chest + abdomen + knee + ankle + bicep + forearm + height + hip
0.96793 0.96688 1571.10800 weight ~ neck + chest + abdomen + knee + ankle + bicep + height + hip
0.96804 0.96685 1572.30300 weight ~ neck + chest + abdomen + knee + ankle + bicep + wrist + height + hip
0.96814 0.96682 1573.44900 weight ~ neck + chest + abdomen + knee + ankle + bicep + forearm + wrist + height + hip
0.96794 0.96661 1575.07800 weight ~ neck + chest + abdomen + thigh + knee + ankle + forearm + wrist + height + hip
0.96751 0.96658 1572.40100 weight ~ neck + chest + knee + ankle + bicep + height + hip
0.96762 0.96655 1573.60100 weight ~ neck + chest + knee + ankle + bicep + wrist + height + hip
0.96759 0.96652 1573.83100 weight ~ neck + chest + knee + ankle + bicep + forearm + height + hip
0.96767 0.96647 1575.17400 weight ~ neck + chest + knee + ankle + bicep + forearm + wrist + height + hip
0.96767 0.96647 1575.18800 weight ~ neck + chest + abdomen + thigh + knee + ankle + forearm + height + hip
0.96762 0.96641 1575.59500 weight ~ neck + chest + abdomen + thigh + knee + ankle + wrist + height + hip
0.96754 0.96633 1576.18700 weight ~ neck + chest + thigh + knee + ankle + forearm + wrist + height + hip
No.16012 Re: 重回帰分析について 【青木繁伸】 2011/12/17(Sat) 20:38
応答が途絶えたのでしょうか。
自分で何らかのプログラムを書いた経験のない人にとっては，既存のプログラムを自分の用途に応じて書き換えるということは難しいのでしょうか。
No.16017 Re: 重回帰分析について 【統計大好き】 2011/12/18(Sun) 14:34
青木先生
ご返信ありがとうございます。
最初の返信を見てから，もう一度色々と試してから，コメントをしようと思っていて，あれから色々としましたが全くできなかったということを報告しようと思って，今ここに来て，青木先生のプログラムを見させていただきました。
かなり難しいのですね，変数を固定するというのは。
わざわざプログラムを載せていただき，ありがとうございます。
このプログラムで行ってみます。
No.16020 Re: 重回帰分析について 【青木繁伸】 2011/12/18(Sun) 19:34
> かなり難しいのですね，変数を固定するというのは。
難しいことだとは思いませんけどねえ。
ま あ，人により，どの程度のことが難しいかは程度問題ですけど，この程度が難しいレベルだとすれば，プログラムの応用というのはとてつもなくむずかしいとい うことになりますね。この程度のことが難しいと言われると，「提示したプログラムを必要に応じて改変してください」というのは無意味なプロポーザルになっ てしまいますね。
No.16092 Re: 重回帰分析について 【統計大好き】 2011/12/27(Tue) 02:45
お久しぶりです。
あれから青木先生のプログラムで，説明変数をいくつか固定させてRに出力してみました。
すると，結果を出すことができたので，説明変数を増やして，25個で行ってみると，メモリー不足がでたので，それを解決するために次のようなプログラムで，2進数表示で，25次元のベクトル表示にして，総当たりをすることにしました。
こうすると，メモリー不足のエラーが出ることはないのですが，かなり時間がかかるようです。
ですので時間をもっと速くする方法があるのかどうかと，このプログラムでは説明変数の固定の仕方が分からないので，もしお分かりであるのであれば，お手数おかけしますが，教えていただけないでしょうか？
お願い致します。data1<-read.table("data.txt")
nisinsu<-NULL
n<-25
bit<-rep(0,n)
for (i in 0:(2^n-1))
{
k<-i
for (j in 1:n)
{
bit[j]<- k %% 2
k<-k %/% 2
}
c(paste(i,":"),bit)
nisinsu<-rbind(nisinsu,c(bit))
}
w5<-numeric(2^25)
set<-data1[,-c(1)]
set_new<-cbind(c(1),set)
mok<-t(t(data1[,1]))
for (i in 2:2^25)
{
p<-sum(nisinsu[i,])+2
b<-nisinsu[c(i),]
b_new<-t(t(rbind(c(1),t(t(b)))))
s<-b_new*c(1:26)
z<-set_new[,t(s)]
data2<-cbind(mok,z)
result1<-lm("mok~.",data=data.frame(data2))
w5[i]<-AIC(result1)
if(i%%1000==0)
{
print(i)
}
}
order(w5)[c(1,2,3,4,5)]
nisinsu[order(w5)[2],]
No.16094 Re: 重回帰分析について 【青木繁伸】 2011/12/27(Tue) 07:38
最後の手段は，コンパイル。
プログラム自体をコンパイラ言語で書き直すとまではいかなくても，R のバイトコンパイラ（compiler パッケージ），include パッケージによる関数のコンパイルなどを検討する。
いつもモデルに含める変数からなるデータ行列をfix とすると
data2 <- cbind(mok, z)
のところで
data2 <- cbind(mok, z, fix)
とするだけでしょう。
プログラムを書くときは，せめて，カンマの後や <- の前後などに空白を置くのがお作法です。
No.16099 Re: 重回帰分析について 【統計大好き】 2011/12/27(Tue) 14:32
青木先生
早速のご返信ありがとうございます。
プログラムの件は以後気をつけます。
compilerパッケージやincludeパッケージは全く知らないので，これから調べてみようと思います。
プログラムについても青木先生の助言を元にもう一度行ってみようと思います。
No.16103 Re: 重回帰分析について 【青木繁伸】 2011/12/27(Tue) 14:52
間違いました include ではなく inline でした。
● 「統計学関連なんでもあり」の過去ログ--- 045 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る