одноуровневые переменные в смешанной модели (lme4) ошибка в R - PullRequest
Новая
       47

одноуровневые переменные в смешанной модели (lme4) ошибка в R

3 голосов
/ 23 ноября 2011 
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

В картинке:

enter image description here

Ответы [ 2 ]

6 голосов
/ 23 ноября 2011

Я расширяю то, что сказал Аарон, так что вся заслуга должна идти в ответ Аарона. 

 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
4 голосов
/ 23 ноября 2011

Попробуйте пакет kinship, основанный на nlme.Подробности смотрите в этой теме на моделях r-sig-mixed.

Для нестандартных ответов вам необходимо изменить lme4 и пакет pedigreemm;см. этот вопрос для деталей.

...