使用 nmle 包的 R 中的多级模型

Question

我正在使用nlme 包来学习多级模型，并在它发生时遵循教科书“使用 R 发现统计”中的示例。

Mixed Models Code

数据集是 Honeymoon 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))

此时发生错误，当我尝试更新“timeRS”模型时。 错误信息如下：

as.character.factor(X[[i]], ...) 中的错误：因子格式错误

这里有知道这意味着什么的统计人员/程序员吗？

Answer 1

我看过这本书。看来编码是错误的。你得到的错误是因为你的时间变量是一个因素，你需要它是数字的。在书中作者的第一个图中，他将时间表示为数字 (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)

模型比较读出的方差分析现在与书中显示的完全一样。