1

我有一些问题，我很感激一些帮助。

 head(new.data)
  WSZ_Code Treatment_Code Year Month TTHM CL2_FREE      BrO3 Colour  PH  TURB seasons
1        2              3 1996     1 30.7     0.35 0.5000750   0.75 7.4 0.055  winter
2        6              1 1996     2 24.8     0.25 0.5001375   0.75 6.9 0.200  winter
3        7              4 1996     2 60.4     0.05 0.5001375   0.75 7.1 0.055  winter
4        7              4 1996     2 58.1     0.15 0.5001570   0.75 7.5 0.055  winter
5        7              4 1996     3 62.2     0.20 0.5003881   2.00 7.6 0.055  spring
6        5              2 1996     3 40.3     0.15 0.5003500   2.00 7.7 0.055  spring

        library(nlme)
    > mod3 <- lme(TTHM ~ CL2_FREE, random= ~ 1| Treatment_Code/WSZ_Code, data=new.data, method ="ML")
    > mod3
    Linear mixed-effects model fit by maximum likelihood
      Data: new.data 
      Log-likelihood: -1401.529
      Fixed: TTHM ~ CL2_FREE 
    (Intercept)    CL2_FREE 
       54.45240   -40.15033 

    Random effects:
     Formula: ~1 | Treatment_Code
            (Intercept)
    StdDev: 0.004156934

     Formula: ~1 | WSZ_Code %in% Treatment_Code
            (Intercept) Residual
    StdDev:    10.90637 13.52372

    Number of Observations: 345
    Number of Groups: 
                  Treatment_Code WSZ_Code %in% Treatment_Code 
                               4                            8 
    > plot(augPred(mod3))
    Error in plot(augPred(mod3)) : 
      error in evaluating the argument 'x' in selecting a method for function 'plot': Error in sprintf(gettext(fmt, domain = domain), ...) : 
      invalid type of argument[1]: 'symbol'

我不确定为什么会收到此错误。ranef 情节似乎还可以

plot(ranef(mod3))

但这仅给出了随机截距的值，没有给出 TTHM 预测。我正在寻找一种方法来绘制典型的 augPred 中的预测，这将显示每个区域的所有随机效应。希望这是有道理的。

4

2 回答 2

2

您需要一个 groupedData 对象才能使用 augPred。我希望这有帮助。

最良好的祝愿@CSJCampbell

con <- textConnection("
WSZ_Code Treatment_Code Year Month TTHM CL2_FREE      BrO3 Colour  PH  TURB seasons
        2              3 1996     1 30.7     0.35 0.5000750   0.75 7.4 0.055  winter
        6              1 1996     2 24.8     0.25 0.5001375   0.75 6.9 0.200  winter
        7              4 1996     2 60.4     0.05 0.5001375   0.75 7.1 0.055  winter
        7              4 1996     2 58.1     0.15 0.5001570   0.75 7.5 0.055  winter
        7              4 1996     3 62.2     0.20 0.5003881   2.00 7.6 0.055  spring
        5              2 1996     3 40.3     0.15 0.5003500   2.00 7.7 0.055  spring
")
new.data <- read.table(con, header = TRUE)

library(nlme)

new.data.grp <- groupedData(TTHM ~ CL2_FREE | Treatment_Code/WSZ_Code, data = new.data)
mod3 <- lme(TTHM ~ CL2_FREE, random= ~ 1| Treatment_Code/WSZ_Code, data=new.data.grp, method ="ML")
mod3
ap3 <- augPred(mod3)
plot(ap3)
于 2013-05-03T13:34:26.797 回答
0

我意识到大多数人可能正在使用ggplot2andlme4在这一点上，但我有点笨拙。

以下是我发现使用适合使用的响应变量列表的几件事lme()。

所以，我一直在处理一些我想适应一组特定输入的响应变量。简而言之，我的代码看起来像

mymodels = list()
for(resp in my_response_vars){
    f = as.formula(paste(resp,paste(my_input_vars,collapse='+'),sep='~'))
    mymodels[[resp]] = lme(fixed=f,random=~wave|group,method="ML",
        data=mydata, na.action=na.exclude)
}

我已经成功地将结果列表中的条目视为普通lme()对象。当我想通过augPred(). 具体来说，我收到以下错误，

Error in tapply(object[[nm]], groups, FUN[["numeric"]], ...) : 
  arguments must have same length

所以，经过大量搜索，我决定看看augPred()via的底层debug()。以下是我得出的一些见解……我不确定这些是否属于错误或是否需要补丁，但我希望它们可以帮助其他有类似问题的人。

  1. 调用augPred()该函数时，会查找原始调用中使用的数据的名称，然后通过调用lme()从 继承此对象。我不确定这是默认为对象框架还是全局，但是，当我在调试中将其更改为时，一切正常。因此，从表面上看，如果您在模型中使用了这些数据的子集，它将调用完整的数据集。parent.frame()eval()data = object$data
  2. 如果一个响应缺少值而您对没有的响应感兴趣，则上述问题会导致问题。由于它包含 data.frame 中的所有内容，作为最终调用gsummary()非响应变量中缺失值的一部分，这将给事情带来麻烦。

  3. 所以，缺失值会把事情搞砸。我默认使用感兴趣的列制作一个临时 data.frame，然后在拟合模型complete.cases()之前运行它。lme()

    mymods = list()
for(resp in my_response_vars){
  f = as.formula(paste(resp,paste(my_input_vars,collapse='+'),sep='~'))
  v2keep = all.vars(f) # grab terms
  smdat = mydata[,c(v2keep,'group')] # include group
  smdat=smdat[complete.cases(smdat),] # scrub missing
  tmpmod = lme(fixed=f, random=~wave|group, 
                method='ML', data=smdat)
  mymods[[resp]] = tmpmod
  # include augPred() call here
}

  4. 如果您没有在调用中包含主要参数，augPred()则需要您的data.frameis 一个groupedData()对象。

因此，如果您遇到arguments must have the same length错误，请尝试：首先以不同的名称对数据进行子集化，确保在拟合模型之前明确清除缺失的行。

于 2018-05-20T00:38:07.710 回答