Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

single level variables in mixed model (lme4) error in R

Tags:

r

mixed

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))

Here M1a, M1b are fixed however individual is random.

    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    

Can we do this in lme4. These are known as animal model and asrmel can do some of such things (link).

EDITS: I forget to mention the relationship matrix is rquired. The following is pedigree structure to do so. To fit the example to size, I reduce the sample size to 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

In term of matrix is following (only lowerhalf triangle plus diagnonal is shown):

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

In picture form:

enter image description here

like image 205
jon Avatar asked Dec 07 '22 18:12

jon


2 Answers

I am expanding what Aaron said, so all credit should go to the Aaron's answer.

 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

# without relatedness structure

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 A matrix

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
like image 78
jon Avatar answered Dec 09 '22 07:12

jon


Try the kinship package, which is based on nlme. See this thread on r-sig-mixed-models for details.

For non-normal responses, you'd need to modify lme4 and the pedigreemm package; see this question for details.

like image 22
Aaron left Stack Overflow Avatar answered Dec 09 '22 06:12

Aaron left Stack Overflow