今天处理数据，用predict()获得SED，但预测值、标准误及SED都缺失。估计状态是aliased，即不可估计，上论坛查找原因。

View topic - Predicting over three seasons for many individuals | Forum | VSN International

kiwi Posted: Wed Mar 24, 2010 2:02 am

数据示例：

#                Season      Vine Block     Flower    DM 
# 2007-52-04-14g   2007 52-04-14g    52      first 18.87 
# 2008-52-04-14g   2008 52-04-14g    52 previously 18.64 
# 2009-52-04-14g   2009 52-04-14g    52 previously 17.70 
# 2007-52-04-15a   2007 52-04-15a    52      first 18.08 
# 2008-52-04-15a   2008 52-04-15a    52 previously 17.98 
# 2009-52-04-15a   2009 52-04-15a    52 previously 17.20 
# 2007-52-04-15b   2007 52-04-15b    52      first 17.27 
# 2008-52-04-15b   2008 52-04-15b    52 previously 18.52 
# 2009-52-04-15b   2009 52-04-15b    52 previously 18.50 
# 2007-52-04-15d   2007 52-04-15d    52      first    NA 
# 2008-52-04-15d   2008 52-04-15d    52 previously 18.20
# 
# with(model.df, levels(Block)) 
# [1] "3"  "39" "45" "52" "54" 
# with(model.df, levels(Season)) 
# [1] "2007" "2008" "2009" 
# with(model.df, levels(Flower)) 
# [1] "first"      "previously" 
# length(levels(model.df$Vine)) 
# [1] 2345

模型：

dm.asr <- try(asreml(fixed = DM ~ Season*Flower + Block, 
                     random = ~ped(Vine):corv(Block) + 
                       ped(Vine):idv(Season), 
                     rcov = ~id(Vine):idv(Season), 
                     na.method.X = "include", 
                     ginverse = list(Vine = G.inv$ginv),
                     data = model.df, 
                     workspace = 50e6, maxiter = 50)) 

预测

dict <- predict(dm.asr, classify = "Vine:Season", maxiter = 1) 
use.preds <- dict$predictions$pvals[,1:4]

出现2个问题：

1. Season和Vine的顺序和原始数据中不一致

tail(use.preds) 
#      Season      Vine predicted.value standard.error 
# 7030   2009 54-13-01b        19.17404       1.100054 
# 7031   2009 54-13-04e        20.88258       1.100054 
# 7032   2009 54-13-12b        19.23295       1.100054 
# 7033   2009 54-14-04b        18.05464       1.100054 
# 7034   2009 54-14-04e        20.76475       1.100054 
# 7035   2009 54-14-04g        19.35078       1.100054

tail(model.df) 
#                Season      Vine Block Flower   DM 
# 2007-54-14-04e   2007 54-14-04e  <NA>   <NA>   NA 
# 2008-54-14-04e   2008 54-14-04e  <NA>   <NA>   NA 
# 2009-54-14-04e   2009 54-14-04e    54  first 23.1 
# 2007-54-14-04g   2007 54-14-04g  <NA>   <NA>   NA 
# 2008-54-14-04g   2008 54-14-04g  <NA>   <NA>   NA 
# 2009-54-14-04g   2009 54-14-04g    54  first 20.7

有点费解，但不是大问题。但是，2007年的每一个Vine的预测都是NA，而其他Season的每个Vine都预测，无论它们是否有读数（根据需要）。 如果我去掉Season 2007，那么2008年的每个预测就都会变成NA。 难道这两种现象是相关的，得到的预测值是可信的吗？ 我也想知道如何在模型的随机部分中进行三向交互。 我知道可以在空间模型中完成，但我能够收敛的唯一模型具有“singular”约束，而且z ratios极高。根据我的经验，这个模型所得到的预测值并不可信（即使残差图看起来很好）。

david.butler Posted: Wed Mar 31, 2010 4:15 am

你的模型在固定效应中有Flower x Season，然而，就你的预测而言，“问题”是这是一个不平衡的组合，Season 1 (2007)没有Flower的水平2。 这些项在平均集中并将导致aliasing，在这种情况下，所有2007年预测都将aliased【这几句是解决问题的关键】。 有几个方法可以处理，但一个例子是使用“present”平均，即对数据中存在的组合进行“平均”预测。 在asreml-R中，可以使用'present'参数指定它以进行预测。 如果你添加present=c('Season','Flower')，然后，预测将由Flower x Season列表中的成员组成。
不清楚你是如何去掉2007的，但利用levels=list(Season=c("2008","2009"))你可以只获得2008和2009的预测值。
当前版本的asreml-R不按照分类集中指定项的顺序，目前这由R遇到项的顺序决定。预测值应该是正确的。

---

具体到我的问题

查看预测结果，可以得到如下信息：

asr_pv$predictions$pvals

Notes:
- The predictions are obtained by averaging across the hypertable
  calculated from model terms constructed solely from factors in
  the averaging and classify sets.
- Use "average" to move ignored factors into the averaging set.

- The SIMPLE averaging set:  Rep 
- Aliased cells:  79.000000
- The overall SED includes non-estimable predictions

知道，平均集是Rep，套到上面那个例子中，就是预测因子和Rep的组合不平衡，即有些水平中有，有些没有，所以加一个present=c()即可解决，尽管我的Rep是固定因子。