[R] Variance Components in R

Iuri Gavronski iuri at proxima.adm.br
Sun Aug 20 15:21:13 CEST 2006


Harold,

I have tried the following syntax:

> fm <- lmer(RATING ~ CHAIN*SECTOR*RESP +(1|CHAIN*SECTOR*RESP), gt)
> summary(fm)
Linear mixed-effects model fit by REML
Formula: RATING ~ CHAIN * SECTOR * RESP + (1 | CHAIN * SECTOR * RESP)
   Data: gt
      AIC      BIC    logLik MLdeviance REMLdeviance
 2767.466 2807.717 -1374.733   2710.253     2749.466
Random effects:
 Groups                Name        Variance Std.Dev.
 CHAIN * SECTOR * RESP (Intercept) 5.7119   2.3900
 Residual                          2.8247   1.6807
number of obs: 647, groups: CHAIN * SECTOR * RESP, 71

Fixed effects:
                    Estimate Std. Error  t value
(Intercept)        4.5760000  2.6193950  1.74697
CHAIN             -0.2014603  0.7984752 -0.25231
SECTOR            -0.1093434  2.3516722 -0.04650
RESP               0.0184237  0.0276326  0.66674
CHAIN:SECTOR       0.1423668  0.3005919  0.47362
CHAIN:RESP         0.0024786  0.0083782  0.29584
SECTOR:RESP       -0.0046001  0.0240517 -0.19126
CHAIN:SECTOR:RESP -0.0011219  0.0030762 -0.36470

Correlation of Fixed Effects:
              (Intr) CHAIN  SECTOR RESP   CHAIN:SECTOR CHAIN:R SECTOR:
CHAIN         -0.435
SECTOR        -0.845 -0.050
RESP          -0.778  0.345  0.645
CHAIN:SECTOR   0.886 -0.732 -0.635 -0.680
CHAIN:RESP     0.351 -0.782  0.038 -0.466  0.566
SECTOR:RESP    0.666  0.038 -0.786 -0.822  0.500       -0.046
CHAIN:SECTOR: -0.709  0.586  0.500  0.879 -0.789       -0.729  -0.635
>

Again, my problem is: there are no fixed effects...
The same dataset, when running at SPSS (I have a subset with 647
records), using the syntax I showed somewhere before, gives me the
following output:

Variance Components Estimation
Variance Estimates
Component      Estimate
Var(CHAIN)     ,530
Var(SECTOR)    ,000(a)
Var(RESP)      2,734
Var(ASPECT)    ,788
Var(ITEM)      ,000(a)
Var(SECTOR *   ,061
RESP)
Var(SECTOR *   ,000(a)
ASPECT)
Var(SECTOR *   ,031
ITEM)
Var(CHAIN *    2,183
RESP)
Var(CHAIN *    ,038
ASPECT)
Var(CHAIN *    ,003
ITEM)
Var(RESP *     ,467
ASPECT)
Var(RESP *     ,279
ITEM)
Var(SECTOR *   ,000(a)
RESP * ASPECT)
Var(SECTOR *   ,077
RESP * ITEM)
Var(CHAIN *    ,773
RESP * ASPECT)
Var(Error)     ,882
Dependent Variable: RATING
 Method: Restricted Maximum Likelihood Estimation
a This estimate is set to zero because it is redundant.

That's what I would like to get from R.

Any help would be appreciated.

Best regards,

Iuri

On 8/20/06, Iuri Gavronski <iuri at ufrgs.br> wrote:
>
> Harold, I have tried to adapt your syntax and got some problems. Some responses from lmer:
>
> On this one, I have tried to use "1" as a grouping variable. As I understood from Bates (2005), grouping variables are like nested design, which is not the case.
> > fm <- lmer(RATING ~ CHAIN*SECTOR*RESP +(CHAIN*SECTOR*RESP|1), gt)
> Erro em lmer(RATING ~ CHAIN * SECTOR * RESP + (CHAIN * SECTOR * RESP |  :
>         Ztl[[1]] must have 1 columns
>
> Then I have tried to ommit the fixed effects...
> > fm <- lmer(RATING ~ (CHAIN*SECTOR*RESP|1), gt)
> Erro em x[[3]] : não é possível dividir o objeto em subconjuntos
> (the error message would be something like "not possible to divide the object in subsets"... I don't know the original wording of message because my R is in Portuguese...)
>
> Then... I have tried to specify RESP (the persons) as the grouping variable (which doesn't make any sense to me, but...)
> > fm <- lmer(RATING ~ CHAIN*SECTOR*RESP +(CHAIN*SECTOR|RESP), gt)
> Warning message:
> nlminb returned message false convergence (8)
>  in: "LMEoptimize<-"(`*tmp*`, value = list(maxIter = 200, tolerance = 1.49011611938477e-08,
> >
>
> Any idea?
>
>
> Regards,
>
> Iuri.
>
>
> On 8/17/06, Doran, Harold <HDoran at air.org> wrote:
> >
> >
> >
> > Iuri:
> >
> > Here is an example of how a model would be specified using  lmer using a couple of your factors:
> >
> > fm <- lmer(response.variable ~ chain*sector*resp  +(chain*sector*resp|GroupingID), data)
> >
> > This will give you a main effect for each factor and all  possible interactions. However, do you have a grouping variable? I wonder if aov  might be the better tool for your G-study?
> >
> > As a side note, I don't see that you have a factor for  persons. Isn't this also a variance component of interest for your  study?
> >
> >
> >    ________________________________
   From: prof.iuri at gmail.com    [mailto:prof.iuri at gmail.com] On Behalf
Of Iuri    Gavronski
> > Sent: Thursday, August 17, 2006 1:26 PM
> > To:    Doran, Harold
> >
> > Cc: r-help at stat.math.ethz.ch
> >
> > Subject: Re:    [R] Variance Components in R
> >
> >
> >
> >
> > I am trying to replicate Finn and Kayandé (1997) study on G-theory    application on Marketing. The idea is to have people evaluate some aspects of    service quality for chains on different economy sectors. Then, conduct a    G-study to identify the generalizability coefficient estimates for different    D-study designs.
> > I have persons rating 3 different items on 3 different    aspects of service quality on 3 chains on 3 sectors. It is normally assumed on    G-studies that the factors are random. So I have to specify a model to    estimate the variance components of CHAIN SECTOR RESP ASPECT ITEM, and the interaction of    SECTOR*RESP SECTOR*ASPECT SECTOR*ITEM CHAIN*RESP CHAIN*ASPECT CHAIN*ITEM    RESP*ASPECT RESP*ITEM SECTOR*RESP*ASPECT SECTOR*RESP*ITEM CHAIN*RESP*ASPECT.    '*' in VARCOMP means a crossed design.
> > Evaluating only the two dimensions    interactions (x*y) ran in few minutes with the full database. Including three    interactions (x*y*z) didn't complete the execution at all. I have the data and    script sent to a professor of the department of Statistics on my university    and he could not run it on either SPSS or SAS (we don't have SAS licenses here    at the business school, only SPSS). Nobody here at the business school has any    experience with R, so I don't have anyone to ask for help.
> > Ì am not    sure if I have answered you question, but feel free to ask it again, and I    will try to restate the problem.
> >
> > Best regards,
> >
> > Iuri
> >
> >
> >
> >
> > On 8/17/06, Doran,    Harold <HDoran at air.org>    wrote:
> >
> > >
> > >
> > >
> > >
> > > This      will (should) be a piece of cake for lmer. But, I don't speak SPSS. Can      you write your model out as a linear model and give a brief description of      the data and your problem?
> > >
> > > In      addition to what Spencer noted as help below, you should also check out the      vignette in the mlmRev package. This will give you many      examples.
> > >
> > > vignette('MlmSoftRev')
> > >
> > >
> > >
> > >
> > >
> > >
> > >
> > >        ________________________________
       From: prof.iuri at gmail.com        [mailto:prof.iuri at gmail.com]
     On Behalf Of Iuri Gavronski
> > > Sent: Thursday, August 17,        2006 11:16 AM
> > > To: Doran, Harold
> > >
> > >
> > > Subject: Re: [R] Variance Components in        R
> > >
> > >
> > >
> > >
> >
> >
> >
> > 9500 records. It didn`t run in SPSS or SAS on Windows machines,      so I am trying to convert the SPSS script to R to run in a RISC station at      the university.
> >
> >
> >
> >
> > On 8/17/06, Doran,      Harold <HDoran at air.org>      wrote:
> >
> > >
> >
> > Iuri:
> >
> > The lmer function is optimal for large data with crossed random        effects.
> > How large are your data?
> >
> > > -----Original        Message-----
> > > From: r-help-bounces at stat.math.ethz.ch
> >
> > > [mailto: r-help-bounces at stat.math.ethz.ch] On Behalf Of Iuri        Gavronski
> >
> > > Sent: Thursday, August 17, 2006 11:08 AM
> > > To:        Spencer Graves
> > > Cc: r-help at stat.math.ethz.ch
> >
> > > Subject: Re: [R]        Variance Components in R
> > >
> > > Thank you for your reply.
> > >        VARCOMP is available at SPSS advanced models, I'm not sure
> > > for        how long it exists... I only work with SPSS for the last
> > > 4        years...
> > > My model only has crossed random effects, what perhaps        would
> > > drive me to lmer().
> > > However, as I have unbalanced        data (why it is normally called
> > > 'unbalanced design'? the data was        not intended to be
> > > unbalanced, only I could not get responses for        all cells...),
> > > I'm afraid that REML would take too much CPU,        memory and time
> > > to execute, and MINQUE would be faster and provide        similar
> > > variance estimates (please, correct me if I'm wrong on        that point).
> > > I only found MINQUE on the maanova package, but as my        study
> > > is very far from genetics, I'm not sure I can use this        package.
> > > Any comment would be appreciated.
> > >        Iuri
> > >
> >
> > > On 8/16/06, Spencer Graves <spencer.graves at pdf.com > wrote:
> > > >
> > >        >       I used SPSS over 25 years ago,        but I don't recall
> > > ever fitting a
> > > > variance components        model with it.  Are all your random
> > > effects        nested?
> > > > If they were, I would recommend you use 'lme' in the        'nlme' package.
> > > > However, if you have crossed random effects,        I suggest you
> > > try 'lmer'
> > > > associated with the 'lme4'        package.
> > > >
> > > >       For        'lmer', documentation is available in Douglas
> > > Bates.        Fitting
> > > > linear mixed models in R. /R News/, 5(1):27-30, May        2005
> >
> > > > (www.r-project.org ->        newsletter).  I also recommend you try the
> >
> > > > vignette        available with the 'mlmRev' package (see, e.g.,
> >
> > > >  http://finzi.psych.upenn.edu/R/Rhelp02a/archive/81375.html        ).
> >
> > > >
> > >        >        Excellent        documentation for both 'lme' (and indirectly for
> > > > 'lmer') is        available in Pinheiro and Bates (2000)
> > > Mixed-Effects        Models
> > > > in S and S-Plus (Springer).  I have        personally recommended
> > > this book
> > > > so many times on        this listserve that I just now got 234 hits for
> > > >        RSiteSearch("graves pinheiro").  Please don't hesitate to pass        this
> > > > recommendation to your university        library.  This book is
> > > the primary
> > > >        documentation for the 'nlme' package, which is part of the
> > >        standard R
> > > > distribution.  A subdirectory        "~library\nlme\scripts" of your R
> > > > installation includes        files named "ch01.R", "ch02.R", ...,
> > > "ch06.R",
> > > >        "ch08.R", containing the R scripts described in the book.  These        R
> > > > script files make it much easier and more enjoyable to        study that
> > > > book, because they make it much easier to try the        commands
> > > described
> > > > in the book, one line at a time,        testing modifications to check you
> > > > comprehension,        etc.  In addition to avoiding problems with
> > > >        typographical errors, it also automatically overcomes a few
> > > minor        but
> > > > substantive changes in the notation between S-Plus and        R.
> > > >
> > > >       Also, the        "MINQUE" method has been obsolete for over
> > > 25 years.
> > > >        I recommend you use method = "REML" except for when you want to
> > >        > compare two nested models with different fixed        effects;  in
> > > that case,
> > > > you should use        method = "ML", as explained in Pinheiro and
> > > Bates (2000).
> > >        >
> > > >       Hope this        helps.
> > > >       Spencer        Graves
> > > >
> > > > Iuri Gavronski wrote:
> > > > >        Hi,
> > > > >
> > > > > I'm trying to fit a model using        variance components in R, but if
> > > > > very new on it, so I'm        asking for your help.
> > > > >
> > > > > I have imported        the SPSS database onto R, but I don't know how to
> > > > >        convert the commands... the SPSS commands I'm trying to
> > > convert        are:
> > > > > VARCOMP
> > > >        >    RATING BY CHAIN SECTOR RESP ASPECT        ITEM
> > > > >    /RANDOM = CHAIN SECTOR RESP        ASPECT ITEM
> > > > >    /METHOD = MINQUE        (1)
> > > > >    /DESIGN = CHAIN SECTOR RESP        ASPECT ITEM
> > > >        >                SECTOR*RESP        SECTOR*ASPECT SECTOR*ITEM CHAIN*RESP
> > > > > CHAIN*ASPECT        CHAIN*ITEM RESP*ASPECT RESP*ITEM
> > > >        >                SECTOR*RESP*ASPECT        SECTOR*RESP*ITEM
> > > CHAIN*RESP*ASPECT
> > > >        >    /INTERCEPT = INCLUDE.
> > > >        >
> > > > > VARCOMP
> > > >        >    RATING BY CHAIN SECTOR RESP ASPECT        ITEM
> > > > >    /RANDOM = CHAIN SECTOR RESP        ASPECT ITEM
> > > > >    /METHOD = REML
> > > > >    /DESIGN = CHAIN SECTOR RESP        ASPECT ITEM
> > > >        >                SECTOR*RESP        SECTOR*ASPECT SECTOR*ITEM CHAIN*RESP
> > > > > CHAIN*ASPECT        CHAIN*ITEM RESP*ASPECT RESP*ITEM
> > > >        >                SECTOR*RESP*ASPECT        SECTOR*RESP*ITEM
> > > CHAIN*RESP*ASPECT
> > > >        >    /INTERCEPT = INCLUDE.
> > > >        >
> > > > > Thank you for your help.
> > > > >
> > >        > > Best regards,
> > > > >
> > > > > Iuri.
> > >        > >
> > > > >        _______________________________________
> > > > > Iuri Gavronski -        iuri at ufrgs.br
> >
> > > >        > doutorando
> > > > > UFRGS/PPGA/NITEC - www.ppga.ufrgs.br        Brazil
> > > > >
> > > > >        ______________________________________________
> > > > > R-help at stat.math.ethz.ch mailing list
> >
> > > > >        https://stat.ethz.ch/mailman/listinfo/r-help
> > >        > > PLEASE do read the posting guide
> > > > http://www.R-project.org/posting-guide.html
> > >        > > and provide commented, minimal, self-contained, reproducible        code.
> > > > >
> > >        >
> > >
> > >       [[alternative        HTML version deleted]]
> > >
> > >        ______________________________________________
> >
> > > R-help at stat.math.ethz.ch  mailing list
> >
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > >        PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > and        provide commented, minimal, self-contained, reproducible code.
> > >
> >
> >
> >
> >
> >
>
>



More information about the R-help mailing list