R help - Apr 2011 - lme4 problem: model defining and effect estimation ------ question from new bird to R community from SAS community

Hi R community,

I am new bird to R and moved recently from SAS. I am no means expert on
either but very curious learner. So your help crucial for me to learn R.
I have already got positive expression.

I was trying to fit a mixed model in animal experiment but stuck at simple
point. The following similar example is from SAS mixed model pp 212.

# data

genetic_evaluation <- read.table(textConnection("
sire dam adg
1  1  2.24
1  1  1.85
1  2  2.05
1  2  2.41
2  1  1.99
2  1  1.93
2  2  2.72
2  2  2.32
3  1  2.33
3  1  2.68
3  2  2.69
3  2  2.71
4  1  2.42
4  1  2.01
4  2  1.86
4  2  1.79
5  1  2.82
5  1  2.64
5  2  2.58
5  2  2.56"), header = TRUE)

# my R practice codes
require (lme4)
 lmer(adg ~ 1 + (1|sire) + (1|dam/sire), data=genetic_evaluation)

 ****error message********************************************88
 Error: length(f1) == length(f2) is not TRUE
In addition: Warning messages:
1: In sire:dam : numerical expression has 20 elements: only the first used
2: In sire:dam : numerical expression has 20 elements: only the first used

**********************how can I estimate the BLUP effects?*************
#equavalent code in SAS
proc mixed data=genetic_evaluation;
   class sire dam;
   model adg= / ddfm=kr;
   random sire dam(sire);
   estimate 'sire 1 BLUP "broad" '
            intercept 1 | sire 1 0;
   estimate 'sire 1 BLUP "narrow" '
            intercept 2 | sire 2 0
            dam(sire) 1 1  0 0  0 0  0 0  0 0  / divisor=2;
   estimate 'sire 1 BLUP with dam 1'
            intercept 1 | sire 1 0
            dam(sire) 1 0;
   ods select CovParms Estimates;
run;

# Estimate statement define predictable functions. All fixed effect
cofficient must appear first and then random effect coefficients. The fixed
and random
#effect cofficient are seperated by |

****************Expected outputs according to SAS
*************************************
                      Estimate
sire 1 BLUP "broad    2.2037
sire 1 BLUP "narrow"  2.1609
sire 1 BLUP with dam1 2.1002

 Data details:
The data is animal science data in which five sires were randomly sampled
from the population and were randomly mated with two dams.
Two offspring per sire dam combination were measured. Average daily gain was
recorded. We are interested in breeding value of ith sire(means that which
gives offsprings with higher gain

NIL

	[[alternative HTML version deleted]]

My question was how can we estimate effects and define correct
model equivalent to SAS code provided.

On Fri, Apr 15, 2011 at 10:21 PM, Nilaya Sharma
<nilaya.sharma@gmail.com>wrote:
> Hi R community,
>
> I am new bird to R and moved recently from SAS. I am no means expert on
> either but very curious learner. So your help crucial for me to learn R.
> I have already got positive expression.
>
> I was trying to fit a mixed model in animal experiment but stuck at simple
> point. The following similar example is from SAS mixed model pp 212.
>
> # data
>
> genetic_evaluation <- read.table(textConnection("
> sire dam adg
> 1  1  2.24
> 1  1  1.85
> 1  2  2.05
> 1  2  2.41
> 2  1  1.99
> 2  1  1.93
> 2  2  2.72
> 2  2  2.32
> 3  1  2.33
> 3  1  2.68
> 3  2  2.69
> 3  2  2.71
> 4  1  2.42
> 4  1  2.01
> 4  2  1.86
> 4  2  1.79
> 5  1  2.82
> 5  1  2.64
> 5  2  2.58
> 5  2  2.56"), header = TRUE)
>
> # my R practice codes
> require (lme4)
>  lmer(adg ~ 1 + (1|sire) + (1|dam/sire), data=genetic_evaluation)
>
>  ****error message********************************************88
>  Error: length(f1) == length(f2) is not TRUE
> In addition: Warning messages:
> 1: In sire:dam : numerical expression has 20 elements: only the first used
> 2: In sire:dam : numerical expression has 20 elements: only the first used
>
> **********************how can I estimate the BLUP effects?*************
> #equavalent code in SAS
> proc mixed data=genetic_evaluation;
>    class sire dam;
>    model adg= / ddfm=kr;
>    random sire dam(sire);
>    estimate 'sire 1 BLUP "broad" '
>             intercept 1 | sire 1 0;
>    estimate 'sire 1 BLUP "narrow" '
>             intercept 2 | sire 2 0
>             dam(sire) 1 1  0 0  0 0  0 0  0 0  / divisor=2;
>    estimate 'sire 1 BLUP with dam 1'
>             intercept 1 | sire 1 0
>             dam(sire) 1 0;
>    ods select CovParms Estimates;
> run;
>
> # Estimate statement define predictable functions. All fixed effect
> cofficient must appear first and then random effect coefficients. The fixed
> and random
> #effect cofficient are seperated by |
>
> ****************Expected outputs according to SAS
> *************************************
>                       Estimate
> sire 1 BLUP "broad    2.2037
> sire 1 BLUP "narrow"  2.1609
> sire 1 BLUP with dam1 2.1002
>
>  Data details:
> The data is animal science data in which five sires were randomly sampled
> from the population and were randomly mated with two dams.
> Two offspring per sire dam combination were measured. Average daily gain
> was recorded. We are interested in breeding value of ith sire(means that
> which
> gives offsprings with higher gain
>
> NIL
>
	[[alternative HTML version deleted]]

On Apr 16, 2011, at 04:21 , Nilaya Sharma wrote:
> genetic_evaluation <- read.table(textConnection("
> sire dam adg
> 1  1  2.24
> 1  1  1.85
> 1  2  2.05
> 1  2  2.41
....
> 4  2  1.86
> 4  2  1.79
> 5  1  2.82
> 5  1  2.64
> 5  2  2.58
> 5  2  2.56"), header = TRUE)
> 
> # my R practice codes
> require (lme4)
> lmer(adg ~ 1 + (1|sire) + (1|dam/sire), data=genetic_evaluation)

You're missing the equivalent of the SAS class statement. Grouping variables
need to be factors:

genetic_evaluation<-transform(genetic_evaluation,
 sire=factor(sire),
 dam=factor(dam))

Also, you probably don't want a random main effect of dam, so

lmer(adg ~ 1 + (1|sire) + (1|dam:sire), data=genetic_evaluation)

or even

lmer(adg ~ 1 + (1|sire/dam), data=genetic_evaluation)


-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

Nilaya Sharma wrote:
> 
> I was trying to fit a mixed model in animal experiment but stuck at simple
> point. The following similar example is from SAS mixed model pp 212.
> 
> 
> genetic_evaluation <- read.table(textConnection("
> sire dam adg
> 1  1  2.24
> 1  1  1.85
> ...
> 5  2  2.58
> 5  2  2.56"), header = TRUE)
> 
> # my R practice codes
> require (lme4)
>  lmer(adg ~ 1 + (1|sire) + (1|dam/sire), data=genetic_evaluation)
> 
>  ****error message********************************************88
>  Error: length(f1) == length(f2) is not TRUE
> In addition: Warning messages:
> 
Thanks for providing a self-contained example. The error message is really a
bit confusing (anybody around who understands what lme thinks here?), but
the solution is simple. Just make sure that dam and sire are factors:

genetic_evaluation$dam = as.factor(genetic_evaluation$dam)
genetic_evaluation$sire = as.factor(genetic_evaluation$sire)

I recommend the excellent course notes in the MCMCglmm package as a
complementary reading.

Dieter




--
View this message in context:
http://r.789695.n4.nabble.com/lme4-problem-model-defining-and-effect-estimation-question-from-new-bird-to-R-community-from-SAS-comy-tp3453530p3453690.html
Sent from the R help mailing list archive at Nabble.com.