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:
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
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With