4 
set.seed(1234) 
mydata <- data.frame (
 individual = factor(1:10), 
 M1a = factor (sample (c(1,2),10, replace = T)),
 M1b = factor (sample (c(1,2),10, replace = T)), 
  pop = factor (c(rep(1, 5), rep (2, 5))), 
 yld = rnorm(10, 10, 2))

这里 M1a，M1b 是固定的，但个体是随机的。

    require(lme4)
    model1 <- lmer(yld ~  M1a + M1b + pop + (1|individual), data = mydata)
    model1
     Error in function (fr, FL, start, REML, verbose)  : 
      Number of levels of a grouping factor for the random effects
    must be less than the number of observations

我们可以在 lme4 中做到这一点吗？这些被称为动物模型，asrmel 可以做一些这样的事情 (l ink )。

编辑：我忘了提到关系矩阵是必需的。以下是这样做的谱系结构。为了使示例适合大小，我将样本大小减少到 10。

 peddf <- data.frame (individual = factor(1:10), 
   mother = c(NA, NA, NA, 1, 1, 1, 1,3, 3,3), 
    father = c(NA, NA, NA, 2, 2, 2, 2, 2, 2, 2))

  individual mother father
1           1     NA     NA
2           2     NA     NA
3           3     NA     NA
4           4      1      2
5           5      1      2
6           6      1      2
7           7      1      2
8           8      3      2
9           9      3      2
10         10      3      2

在矩阵方面如下（仅显示下半三角形加对角线）：

1   NA  NA  NA  NA  NA  NA  NA  NA  NA
0   1   NA  NA  NA  NA  NA  NA  NA  NA
0   0   1   NA  NA  NA  NA  NA  NA  NA
0.25    0.25    0   1   NA  NA  NA  NA  NA  NA
0.25    0.25    0   0.25    1   NA  NA  NA  NA  NA
0.25    0.25    0   0.25    0.25    1   NA  NA  NA  NA
0.25    0.25    0   0.25    0.25    0.25    1   NA  NA  NA
0   0.25    0.25    0.125   0.125   0.125   0.125   1   NA  NA
0   0.25    0.25    0.125   0.125   0.125   0.125   0.25    1   NA
0   0.25    0.25    0.125   0.125   0.125   0.125   0.25    0.25    1

图片形式：

在此处输入图像描述

4

2 回答 2

9

我正在扩展亚伦所说的话，所以所有的功劳都应该归功于亚伦的回答。

 kmat <- kinship(peddf$individual , peddf$father,peddf$mother)
    kmat 
              1    2    3     4     5     6     7     8     9    10
        1  0.50 0.00 0.00 0.250 0.250 0.250 0.250 0.000 0.000 0.000
        2  0.00 0.50 0.00 0.250 0.250 0.250 0.250 0.250 0.250 0.250
        3  0.00 0.00 0.50 0.000 0.000 0.000 0.000 0.250 0.250 0.250
        4  0.25 0.25 0.00 0.500 0.250 0.250 0.250 0.125 0.125 0.125
        5  0.25 0.25 0.00 0.250 0.500 0.250 0.250 0.125 0.125 0.125
        6  0.25 0.25 0.00 0.250 0.250 0.500 0.250 0.125 0.125 0.125
        7  0.25 0.25 0.00 0.250 0.250 0.250 0.500 0.125 0.125 0.125
        8  0.00 0.25 0.25 0.125 0.125 0.125 0.125 0.500 0.250 0.250
        9  0.00 0.25 0.25 0.125 0.125 0.125 0.125 0.250 0.500 0.250
        10 0.00 0.25 0.25 0.125 0.125 0.125 0.125 0.250 0.250 0.500

# 没有关联结构 

model1 <- lmekin(yld ~  M1a + M1b + pop , random = ~ 1|individual, data = mydata)
Linear mixed-effects kinship model fit by maximum likelihood
  Data: mydata 
  Log-likelihood = -20.23546 
  n= 10 

Fixed effects: yld ~ M1a + M1b + pop 
              Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  8.6473627   1.977203  4.3735334 0.004701001
M1a2         1.6722908   1.671041  1.0007477 0.355584122
M1b2        -0.7939123   1.671041 -0.4751003 0.651516161
pop2         0.5265145   1.671041  0.3150817 0.763369802

Wald test of fixed effects =  1.343476 df =  3 p =  0.718836

Random effects: ~1 | individual 
              individual     resid
Standard Dev:  0.9493070 1.5651426
% Variance:    0.2689414 0.7310586

# 用A矩阵 

model2 <- lmekin(yld ~  M1a + M1b + pop , random = ~ 1|individual, varlist=list(kmat), data = mydata) 
Linear mixed-effects kinship model fit by maximum likelihood
  Data: mydata 
  Log-likelihood = -20.23548 
  n= 10 

Fixed effects: yld ~ M1a + M1b + pop 
              Estimate Std. Error    t value    Pr(>|t|)
(Intercept)  8.6473583   1.977206  4.3735251 0.004701044
M1a2         1.6722972   1.671042  1.0007511 0.355582600
M1b2        -0.7939228   1.671044 -0.4751057 0.651512529
pop2         0.5265200   1.671040  0.3150851 0.763367298

Wald test of fixed effects =  1.343489 df =  3 p =  0.7188331

Random effects: ~1 | individual 
 Variance list: list(kmat) 
               individual    resid
Standard Dev: 5.78864e-03 1.830529
% Variance:   9.99990e-06 0.999990
于 2011-11-23T18:56:23.303 回答
6

尝试kinship基于nlme. 有关详细信息，请参阅r-sig-mixed-models 上的此线程。

对于非正常响应，您需要修改lme4和pedigreemm包；有关详细信息，请参阅此问题。

于 2011-11-23T18:40:03.487 回答