Question

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

我们能够做到这一点。这些产品被称为动物模型,并且作为熔炉可以做一些(link )。

EDITS:我忘记提到关系矩阵。下面是挖掘结构。为了与规模相适应,我将抽样规模减少到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

在矩阵中,矩阵如下(仅显示低温三角和 di):

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

图片:

Answer 1

Try the kinship Pack, which is based on nlme. 见 ,关于细节的线上-sig-mixed-models。

对于非常规答复,请修改<代码>lme4和pedigreemm。页: 1

Answer 2

我正在扩大Aaron所说的话,因此,所有信贷都应用于Aaron的回答。

 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

<>0>无相关结构>

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

# with Amap

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

友情链接