2016-12-20 68 views
0

我正在使用nlme軟件包來學習多層次模型,並遵循教科書「Discovering Statistics Using R」的示例。錯誤使用nlme;畸形因子

Mixed Models Code

的數據集是蜜月Period.dat,在他們的同伴網站也可下載。

Data Set - Multilevel Models

require(nlme) 
require(reshape2) 
satisfactionData = read.delim("Honeymoon Period.dat", header = TRUE) 

restructuredData<-melt(satisfactionData, id = c("Person", "Gender"), measured = c("Satisfaction_Base", "Satisfaction_6_Months", "Satisfaction_12_Months", "Satisfaction_18_Months")) 
names(restructuredData)<-c("Person", "Gender", "Time", "Life_Satisfaction") 


#print(restructuredData) 
#restructuredData.sorted<-restructuredData[order(Person),] 

intercept <-gls(Life_Satisfaction~1, data = restructuredData, method = "ML", na.action = na.exclude) 
randomIntercept <-lme(Life_Satisfaction ~1, data = restructuredData, random = ~1|Person, method = "ML", na.action = na.exclude, control = list(opt="optim")) 
anova(intercept, randomIntercept) 

timeRI<-update(randomIntercept, .~. + Time) 
timeRS<-update(timeRI, random = ~Time|Person) 
ARModel<-update(timeRS, correlation = corAR1(0, form = ~Time|Person)) 

的錯誤發生在這一刻,當我試圖更新「定時器」模式。錯誤如下:

Error in as.character.factor(X[[i]], ...) : malformed factor 

任何統計人員/程序員在這裏誰知道這是什麼意思?

回答

3

我看過這本書。看起來編碼是錯誤的。你得到的錯誤是因爲你的時間變量是一個因素,你需要它是數字。在作者的書中第一個數字中,他將時間表示爲數字(0 - 3),但他的模型代碼不正確。我已經爲你重新編寫了它:

## First, Time needs to be recoded as a numeric 

restructuredData$Time.Numeric <- with(restructuredData, ifelse(Time == "Satisfaction_Base", 0, 
     ifelse(Time == "Satisfaction_6_Months", 1, 
     ifelse(Time == "Satisfaction_12_Months", 2, 
     ifelse(Time == "Satisfaction_18_Months", 3, NA))))) 

## Baseline Model 

intercept <-gls(Life_Satisfaction~1, data = restructuredData, method = "ML", na.action = na.exclude) 

summary(intercept) 

## Model where intercept can vary for Individuals 

randomIntercept <- lme(Life_Satisfaction ~ 1, data = restructuredData, random = ~1|Person, method = "ML", na.action = na.exclude, control = list(opt = "optim")) 

summary(randomIntercept) 

## Add time as a fixed effect 

timeRI <- lme(Life_Satisfaction ~ Time.Numeric, data = restructuredData, random = ~1|Person, method = "ML", na.action = na.exclude, control = list(opt = "optim")) 

summary(timeRI) 

## Add a random slope to the model by nesting the Individual within the test time 

timeRS <- lme(Life_Satisfaction ~ Time.Numeric, data = restructuredData, random = ~Time.Numeric|Person, method = "ML", na.action = na.exclude, control = list(opt = "optim")) 

summary(timeRS) 


## Modeling the covariance structure structure of the errors with a first-order autoregressive covariance structure 

ARModel <- update(timeRS, correlation = corAR1(0, form = ~Time.Numeric|Person)) 

summary(ARModel) 

anova(intercept, randomIntercept, timeRI, timeRS, ARModel) 

讀出的模型比較的anova現在完全如書中所示。

+1

很好的建議!我簡化了一下:'restructuredData [,'Time'] < - as.numeric(restructuredData [,'Time「]) - 1'這可以在第一個模型定義之前運行:'攔截<-gls(Life_Satisfaction〜1 ,data = restructuredData,method = 「ML」,na.action = na.exclude)'而不需要改變其餘的代碼。 –