4
$\begingroup$

I have nested structure design with the following model:

Yijklt = m + ai + bij + rijk + sikl + (bs)ijkl + eijkl 

where i = 1, 2, ..................., s; number of sets 
j = 1, 2, ........................, r; number of replications 
k = 1, 2, ........................., m; males 
l = 1,2,............................, f; females 
 t = 1, 2, ................, n; progeny 

where m is the general mean
      ai is the effect of the ith set
      bij is effect of jth replication in the ith set
      rik is the effect of kth male in the ith set
      sikl is the the effect of the lth female to kth male in the ith set
      (bs)ijkl is the interaction term

The "non-mixed" model formula in R is as follows:

 model <- lm(y ~ set + replication %in% set + male %in% 
            set + female %in% male %in% set + replication %in% 
            female %in% male %in% set)
 model <- lm(y ~ set + replication:set + male:set + female:male:set + replication:
            female:male:set)

What would be its corresponding model in mixed model using lme4 or nlme package in R?

I want to do it right.

Please note that all terms in the new mixed model are random, I just want to estimate variance components associated.

Edits: work out example based on the answers below by @Aaron

library(lme4)
set.seed(5)
d <- expand.grid(set=factor(1:3), male=factor(1:3),female=factor(1:3), 
progeny = factor(1:3), replication=factor(1:4))
d$y <- rnorm(nrow(d))
Improve this question
$\endgroup$
2
  • $\begingroup$ As I am not getting any responses in this statistical forum, I am not sure if this is proper forum to post this question, or somebody (suggest) migrate this question to stackoverflow $\endgroup$
    – John
    Nov 15, 2011 at 11:00
  • $\begingroup$ Could you provide some sample data (preferably generated with just a few lines of code)? It's not clear to me how the (bs)ijkl term is different from the eijkl term, or what role the t variable plays. (I'll post an example as an answer.) $\endgroup$
    – Aaron left Stack Overflow
    Nov 15, 2011 at 15:41

2 Answers 2

Reset to default
3
$\begingroup$

How about something like this?

library(lme4)
model <- lmer(y ~ (1|set) + (1|replication:set) + (1|male:set) +  
                  (1|female:male:set) +
                  (1|replication:female:male:set), data=d)
Improve this answer
$\endgroup$
7
  • $\begingroup$ you need to add another factor progeny, which can be 1:3, this will add term to calculate error df (not 0). thanks; $\endgroup$
    – John
    Nov 15, 2011 at 16:08
  • $\begingroup$ Done; data creation code now moved to the question. $\endgroup$
    – Aaron left Stack Overflow
    Nov 15, 2011 at 16:37
  • $\begingroup$ thank you so much, I still wonder difference between ":", "/", and "|" notation in modelling in lme4 $\endgroup$
    – John
    Nov 15, 2011 at 17:36
  • $\begingroup$ proably this question can be different question in itself, too much in the same post $\endgroup$
    – John
    Nov 15, 2011 at 17:42
  • 1
    $\begingroup$ @Wayne: No, it doesn't seem to be in the internal documentation. But see section 11.1 of An Introduction to R. $\endgroup$
    – Aaron left Stack Overflow
    Nov 15, 2011 at 20:16
1
$\begingroup$

I would recommend that you look at the text book "Linear Mixed Models: a practical guide using statistical software" by West, Welch and Galecki. It goes through step-wise the process of model specification, testing, intrepretation for real world datasets for a variety of software packages. the R component is using nlme. I found it very helpful at developing my code. I haven't seen variance component estimates within the book, but have seen this discussed in "the R book", by Crawley.
Good luck Natasha

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.