作者ching0629 (Syameroke)
看板Statistics
標題Re: [程式] 在R裡面用多組資料作線性迴歸
時間Wed Jan 9 19:53:23 2013
※ 引述《swwf (十字路口)》之銘言：
: ------------------------------------------------------------------------
: [軟體程式類別]:R
: [程式問題]: 資料處理、迴歸
: [軟體熟悉度]:
: 低(1~3個月)
: [問題敘述]:
:
: 是這樣，我有多組數據需要跑線性迴歸，y~ ax+b
: 如果只是單純跑一組的話我還會弄，但如果是很多分組分層的話就不知道怎麼寫迴圈
: 原始數據是X(時間)，和Y(存活比率)作簡單直線迴歸，
: 現在如果我數據有多項不同的處理，整理結果如下樣式
: 處理A 處理B 處理C 重複 時間 存活比率
: 7個層級 2個層級 6個 4個 X1~Xn Y1~ Yn
: 要注意的是，每個組合的時間n是不一樣的觀察數量，有的組合觀察數少(早早全死光)，
: 有的高(等很久才死完)
: 現在想要跑這7*2*6*4個直線迴歸，把作迴歸得來的值給輸出，
: 並把這些直線迴歸作圖，請問我的迴圈程式要怎麼寫?
: 搞不清楚要用apply還是tapply...
: 請詳盡敘述遭遇到的問題，可能的話，分點敘述你要處理的流程
: [程式範例]:
: 雖然張貼程式很可怕，但基本上有些程式還是要張貼才能解決
: for (i in treat1)
: {for (j in treat2)
: {for (k in treat3)
: {for (l in rep)
: {lm.reg=lm(x~y)
: }
: }
: }
: }
: summary[i,j,k,l](lm.reg)
: plot[i,j,k,l](x,y)
: -----------------------------------------------------------------------------
: 不好意思，我是R和寫程式的新手，真的不太知道要怎麼樣跑這種迴圈，另外想要一頁
: 多圖把結果的圖都畫一下，想請教各位大大指導我該如何寫這迴圈，謝謝。
: 打擾各位不好意思 >_<
問題實在問的很不清楚,我就按照我假設你在想的做吧
你想要做7*2*4*6個迴歸,而輸出的項目是每個迴歸的"斜率"及"標準誤"
如果你還想輸出更多東西就自己寫在迴圈內吧
Output=matrix(NA,nrow=7*2*4*6,ncol=6)
data$treat1=factor(data$treat1)
data$treat2=factor(data$treat2)
data$treat3=factor(data$treat3)
data$rep=factor(data$rep)
for (i in 1:length(levels(treat1))) {
for (j in 1:length(levels(treat2))) {
for (k in 1:length(levels(treat3))) {
for (l in 1:length(levels(rep))) {
dataset=data[data$treat1==levels(treat1)[i]&data$treat2==levels(treat2)[j]&data
$treat3==levels(treat3)[k]&data$rep==levels(rep)[l],]
n=2*4*6*(i-1)+4*6*(j-1)+6*(k-1)+l
Output[n,1]=levels(treat1)[i]
Output[n,2]=levels(treat2)[j]
Output[n,3]=levels(treat3)[k]
Output[n,4]=levels(rep)[l]
if (sum(dataset$x)==0) {
Output[n,5]=NA
Output[n,6]=NA
}
else {
Output[n,5]=lm(y~x,data=dataset)$coefficients[2]
Output[n,6]=summary(lm(y~x,data=dataset))$coefficients[2,2]
}
}
}
}
}
colnames(Output)=c("treat1","treat2","treat3","rep","coefficients","se")
最後Output這項物件就有你要的結果
共有7*2*4*6列
第一行是treat1
第二行是treat2
第三行是treat3
第四行是rep
第五行是符合前4項條件的資料的斜率
第六行是符合前4項條件的資料的斜率的標準誤
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 61.64.175.188
※ 編輯: ching0629 來自: 61.64.175.188 (01/09 19:54)
1F：推 swwf:謝謝!!! 01/10 00:14
2F：→ swwf:圖的話就無所謂一定要一頁全擺完，但想要批次畫系列圖不會 01/10 00:15
※ 編輯: ching0629 來自: 27.105.2.161 (01/10 00:19)
3F：推 swwf:發現一個問題是，有的數據是存活率直接0，結果就沒X值，因為 01/10 00:39
4F：→ swwf:只有一個點，我就沒把它放進資料中，因為無存活時間@@ 這樣造 01/10 00:40
5F：→ swwf:成各個數據不一樣長@@，不知道要怎麼設定跳過那些n/a的 @@ 01/10 00:40
6F：→ ching0629:不太知道你說的是什麼,假設有些列X是0吧,幫你跳過了 01/10 07:59
7F：→ ching0629:如果是NA什麼的,先用x[is.na(x)==1]<-0,把NA設定為0 01/10 08:01
※ 編輯: ching0629 來自: 210.60.122.8 (01/10 08:02)
8F：推 swwf:結果跟我說 'closure' 類型的物件無法具有子集合<--請教何解? 01/10 13:19
9F：→ gsuper:as.factor(X) 01/11 20:14