2013-03-12 98 views
2

我試圖將回歸函數應用於因子(主題)的每個單獨的級別。這個想法是,對於每個主題,我可以根據他們的實際閱讀時間(RT)和相應打印字符串(WordLen)的長度來獲得預測閱讀時間。一位同事幫助我解決了一些基於(Subject)中另一個函數(Region)的每個級別應用函數的代碼。但是,無論是原始代碼還是我的嘗試修改(在單個因素間使用跨功能的功能)都可以使用。應用迴歸,同時循環R中的因子水平

下面是一些樣本數據的嘗試:

test0<-structure(list(Subject = c(101L, 101L, 101L, 101L, 101L, 101L, 
101L, 101L, 101L, 101L, 102L, 102L, 102L, 102L, 102L, 102L, 102L, 
102L, 102L, 102L, 103L, 103L, 103L, 103L, 103L, 103L, 103L, 103L, 
103L, 103L), Region = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L), RT = c(294L, 241L, 346L, 339L, 332L, NA, 399L, 
377L, 400L, 439L, 905L, 819L, 600L, 520L, 811L, 1021L, 508L, 
550L, 1048L, 1246L, 470L, NA, 385L, 347L, 592L, 507L, 472L, 396L, 
761L, 430L), WordLen = c(3L, 3L, 3L, 3L, 3L, 3L, 5L, 7L, 3L, 
9L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 5L, 7L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 5L, 7L, 3L)), .Names = c("Subject", "Region", "RT", "WordLen" 
), class = "data.frame", row.names = c(NA, -30L)) 

不幸的是,這個數據正在恢復,我不跟我的完整數據集得到了一個問題:

"Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
    0 (non-NA) cases" 

也許這是因爲樣本數據太小?

無論如何,我希望有人看到這個問題的代碼,儘管我提供工作數據的能力......

這是原來的代碼(不工作):

for(i in 1:length(levels(test0$Subject))) 
    for(j in 1:length(levels(test0$Region))) 
    {tmp=predict(lm(RT~WordLen,test0[test0$Subject==levels(test0$Subject)[i] & test0$Region==levels(test0$Region)[j],],na.action="na.exclude")) 
    test0[names(tmp),"rt.predicted"]=tmp 
    } 

,這是修改後的代碼(這並不奇怪,也不起作用):

for(i in 1:length(levels(test0$Subject))) 
    {tmp=predict(lm(RT~WordLen,test0[test0$Subject==levels(test0$Subject)[i],],na.action="na.exclude")) 
    test0[names(tmp),"rt.predicted"]=tmp 
    } 

我將非常感謝任何建議。

+2

也看到'? 'nlme'包中的lmList'。 – 2013-03-12 12:44:18

回答

3

您可以使用庫plyr中的函數ddply()獲得結果。 這將根據Subject拆分數據幀,計算迴歸模型的預測,然後作爲新列添加到數據幀。

ddply(test0,.(Subject),transform, 
    pred=predict(lm(RT~WordLen,na.action="na.exclude"))) 

    Subject Region RT WordLen  pred 
1  101  1 294  3 327.9778 
...... 
4  101  1 339  3 327.9778 
5  101  1 332  3 327.9778 
6  101  2 NA  3  NA 
7  101  2 399  5 363.8444 
....... 
13  102  1 600  3 785.4146 

要通過Subject和拆分數據Region你應該把兩個變量中.()

ddply(test0,.(Subject,Region),transform, 
    pred=predict(lm(RT~WordLen,na.action="na.exclude"))) 
+0

這很好用,謝謝。我如何修改這個也是按區域分割的(對每個主題的每個區域進行迴歸)? – 2013-03-12 12:31:49

+0

@DT更新了我的答案。 – 2013-03-12 12:35:46

+0

非常好。我仍然好奇原始循環方法爲什麼不起作用。我意識到循環不應該成爲我與R的第一線攻擊,但它是很好的知道。 – 2013-03-12 12:41:34

2

在測試數據的唯一問題是,SubjectRegion不是因素。

test0$Subject <- factor(test0$Subject) 
test0$Region <- factor(test0$Region) 

for(i in 1:length(levels(test0$Subject))) 
    for(j in 1:length(levels(test0$Region))) 
    {tmp=predict(lm(RT~WordLen,test0[test0$Subject==levels(test0$Subject)[i] & test0$Region==levels(test0$Region)[j],],na.action="na.exclude")) 
    test0[names(tmp),"rt.predicted"]=tmp 
    } 
# 26  27  28  29  30 
# 442.25 442.25 560.50 678.75 442.25 

原因你讓你的錯誤(0 non-NA cases)是當你子集,你在做它是不是因素的變量水平。在你原始數據集,嘗試:

test0[test0$Subject==levels(test0$Subject)[1],] 

你得到:

# [1] Subject Region RT  WordLen 
# <0 rows> (or 0-length row.names) 

這是什麼lm()試圖用

+0

謝謝你收到這個錯誤。在我的原始數據中,它們是因素,但是在裁減數據時我錯過了這一點。 – 2013-03-12 12:30:23

0

工作,我會認爲這是由以下事實引起的兩個分類變量的組合不存在數據。你可以做的是首先提取子集,檢查它是否不等於NULL,並且只有在有數據時才執行lm。

2

雖然你的問題好像是問錯誤的解釋,這人已經回答(數據不被因素在所有),這裏是一個辦法做到這一點只用base

test0$rt.predicted <- unlist(by(test0[, c("RT", "WordLen")], list(test0$Subject, test0$Region), FUN = function(x) predict(lm(RT ~ 
    WordLen, x, na.action = "na.exclude")))) 

test0 
## Subject Region RT WordLen rt.predicted 
## 1  101  1 294  3  310.4000 
## 2  101  1 241  3  310.4000 
## 3  101  1 346  3  310.4000 
## 4  101  1 339  3  310.4000 
## 5  101  1 332  3  310.4000 
## 6  101  2 NA  3  731.0000 
## 7  101  2 399  5  731.0000 
## 8  101  2 377  7  731.0000 
## 9  101  2 400  3  731.0000 
## 10  101  2 439  9  731.0000 
## 11  102  1 905  3  448.5000 
## 12  102  1 819  3   NA 
## 13  102  1 600  3  448.5000 
## 14  102  1 520  3  448.5000 
## 15  102  1 811  3  448.5000 
## 16  102  2 1021  3   NA 
## 17  102  2 508  3  399.0000 
## 18  102  2 550  5  408.5000 
## 19  102  2 1048  7  389.5000 
## 20  102  2 1246  3  418.0000 
## 21  103  1 470  3  870.4375 
## 22  103  1 NA  3  870.4375 
## 23  103  1 385  3  877.3750 
## 24  103  1 347  3  884.3125 
## 25  103  1 592  3  870.4375 
## 26  103  2 507  3  442.2500 
## 27  103  2 472  3  442.2500 
## 28  103  2 396  5  560.5000 
## 29  103  2 761  7  678.7500 
## 30  103  2 430  3  442.2500 
+0

謝謝你的替代 - 。因素水平問題只是次要的。真正的問題是我的代碼不適用於真正的數據集(正確編碼的因子水平)。或者我錯了,你是說我的原代碼應該已經工作? – 2013-03-12 21:14:34