3

作为 MCMCglmm 模型选择的一个选项（另请参阅此相关问题），我正在尝试使用包 MuMIn 进行模型平均。它似乎不起作用- 请参阅下面的输出。任何想法为什么？输出看起来很废话。特别是，z 值有一堆 NA 值，当这些不是 NA 时，它们都正好是 1。这可能源于除一个模型之外的所有模型都被分配了 0 的权重，这又看起来不切实际.

请注意，在 MuMIn 的文档中，它被列为与 MCMCglmm 对象兼容。

可重现的例子：

set.seed(1234)
library(MCMCglmm)

data(bird.families) 
n <- Ntip(bird.families)

# Create some dummy variables
d <- data.frame(taxon = bird.families$tip.label, 
                X1 = rnorm(n),
                X2 = rnorm(n),
                X3 = sample(c("A", "B", "C"), n, replace = T),
                X4 = sample(c("A", "B", "C"), n, replace = T))

# Simulate a phenotype composed of phylogenetic, fixed and residual effects
d$phenotype <- rbv(bird.families, 1, nodes="TIPS") +
               d$X1*0.7 + 
               ifelse(d$X3 == "B", 0.5, 0) + 
               ifelse(d$X3 == "C", 0.8, 0) +
               rnorm(n, 0, 1)  

# Inverse matrix of shared phyloegnetic history
Ainv <- inverseA(bird.families)$Ainv

# Set priors
prior <- list(R = list(V = 1, nu = 0.002), 
              G = list(G1 = list(V = 1, nu = 0.002)))

uMCMCglmm <- updateable(MCMCglmm)
model <- uMCMCglmm(phenotype ~ X1 + X2 + X3 + X4, 
                  random = ~taxon, 
                  ginverse = list(taxon=Ainv),
                  data = d, 
                  prior = prior, 
                  verbose = FALSE)

# Explore possible simplified models
options(na.action = "na.fail")
dred <- dredge(model)

# Calculate a model average
avg <- model.avg(dred)
summary(avg)

输出：

Call:
model.avg(object = dred)

Component model call: 
uMCMCglmm(fixed = phenotype ~ <16 unique rhs>, random = ~taxon, data = d, 
     prior = prior, verbose = FALSE, ginverse = list(taxon = Ainv))

Component models: 
       df  logLik   AICc  delta weight
3       5  -49.24 108.93   0.00      1
4       5  -71.18 152.82  43.89      0
(Null)  3  -76.98 160.13  51.20      0
34      7  -90.35 195.56  86.63      0
23      6  -95.03 202.71  93.78      0
24      6 -105.79 224.22 115.29      0
1       4 -134.87 278.04 169.11      0
123     7 -137.36 289.59 180.66      0
2       4 -154.82 317.93 209.00      0
234     8 -162.69 342.51 233.58      0
13      6 -167.74 348.12 239.19      0
124     7 -171.06 356.99 248.05      0
14      6 -172.53 357.70 248.77      0
134     8 -171.60 360.33 251.40      0
12      5 -181.16 372.78 263.84      0
1234    9 -189.33 398.07 289.14      0

Term codes: 
X1 X2 X3 X4 
 1  2  3  4 

Model-averaged coefficients:  
(full average) 
              Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.642e-01         NA      NA       NA
X3B          6.708e-01  6.708e-01       1    0.317
X3C          9.802e-01  9.802e-01       1    0.317
X4B         -9.505e-11  9.505e-11       1    0.317
X4C         -7.822e-11  7.822e-11       1    0.317
X2          -3.259e-22  3.259e-22       1    0.317
X1           1.378e-37  1.378e-37       1    0.317

(conditional average) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.76421         NA      NA       NA
X3B          0.67078         NA      NA       NA
X3C          0.98025         NA      NA       NA
X4B         -0.32229         NA      NA       NA
X4C         -0.26522         NA      NA       NA
X2          -0.07528         NA      NA       NA
X1           0.72300         NA      NA       NA

Relative variable importance: 
                     X3    X4    X2    X1   
Importance:              1 <0.01 <0.01 <0.01
N containing models:     8     8     8     8
4

0 回答 0