Dear Harold and others,
I have changed the syntax for lmer() and used this one:
require(lme4)
gt <- read.table("gt5.txt")
sink("GT output.txt")
attach(gt)
system.time(fm <- lmer(RATING ~ 1
+(1|CHAIN)
+(1|SECTOR)
+(1|RESP)
+(1|ASPECT)
+(1|ITEM)
+(1|SECTOR*RESP)
+(1|SECTOR*ASPECT)
+(1|SECTOR*ITEM)
+(1|CHAIN*RESP)
+(1|CHAIN*ASPECT)
+(1|CHAIN*ITEM)
+(1|RESP*ASPECT)
+(1|RESP*ITEM)
+(1|SECTOR*RESP*ASPECT)
+(1|SECTOR*RESP*ITEM)
+(1|CHAIN*RESP*ASPECT),
gt)
)
options(digits = 4)
options(OutDec = ",")
summary(fm, digits = 4)
sink()
Then the output I got from summary(lm) was:
Linear mixed-effects model fit by REML
Formula: RATING ~ 1 + (1 | CHAIN) + (1 | SECTOR) + (1 | RESP) + (1 |
ASPECT) + (1 | ITEM) + (1 | SECTOR * RESP) + (1 | SECTOR *
ASPECT) + (1 | SECTOR * ITEM) + (1 | CHAIN * RESP) + (1 | CHAIN *
ASPECT) + (1 | CHAIN * ITEM) + (1 | RESP * ASPECT) + (1 | RESP *
ITEM) + (1 | SECTOR * RESP * ASPECT) + (1 | SECTOR * RESP * ITEM)
+ (1 | CHAIN * RESP * ASPECT)
Data: gt
AIC BIC logLik MLdeviance REMLdeviance
2386 2462 -1176 2353 2352
Random effects:
Groups Name Variance Std.Dev.
CHAIN * RESP * ASPECT (Intercept) 5,89e-01 0,7675133
SECTOR * RESP * ITEM (Intercept) 4,91e-02 0,2216137
RESP * ITEM (Intercept) 2,75e-01 0,5242572
CHAIN * RESP (Intercept) 1,98e+00 1,4068696
SECTOR * RESP * ASPECT (Intercept) 5,17e-10 0,0000227
CHAIN * ITEM (Intercept) 5,17e-10 0,0000227
RESP * ASPECT (Intercept) 4,77e-01 0,6908419
SECTOR * RESP (Intercept) 3,42e-01 0,5848027
CHAIN * ASPECT (Intercept) 1,61e-02 0,1269306
SECTOR * ITEM (Intercept) 2,24e-02 0,1495102
ITEM (Intercept) 5,17e-10 0,0000227
CHAIN (Intercept) 8,88e-01 0,9424441
RESP (Intercept) 2,80e+00 1,6747970
SECTOR * ASPECT (Intercept) 5,17e-10 0,0000227
ASPECT (Intercept) 8,07e-01 0,8984151
SECTOR (Intercept) 5,17e-10 0,0000227
Residual 1,03e+00 1,0172221
number of obs: 647, groups: CHAIN * RESP * ASPECT, 138; SECTOR * RESP
* ITEM, 138; RESP * ITEM, 70; CHAIN * RESP, 70; SECTOR * RESP *
ASPECT, 47; CHAIN * ITEM, 36; RESP * ASPECT, 24; SECTOR * RESP, 24;
CHAIN * ASPECT, 18; SECTOR * ITEM, 18; ITEM, 9; CHAIN, 9; RESP, 8;
SECTOR * ASPECT, 6; ASPECT, 3; SECTOR, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 5,797 0,891 6,51
Comparing the output I had from R and SPSS, for the same database:
Component Estimate SPSS R
Var(CHAIN) ,530 0,888
Var(SECTOR) ,000(a) 0,000
Var(RESP) 2,734 2,800
Var(ASPECT) ,788 0,807
Var(ITEM) ,000(a) 0,000
Var(SECTOR * ,061 0,342
RESP)
Var(SECTOR * ,000(a) 0,000
ASPECT)
Var(SECTOR * ,031 0,022
ITEM)
Var(CHAIN * 2,183 1,980
RESP)
Var(CHAIN * ,038 0,016
ASPECT)
Var(CHAIN * ,003 0,000
ITEM)
Var(RESP * ,467 0,477
ASPECT)
Var(RESP * ,279 0,275
ITEM)
Var(SECTOR * ,000(a) 0,000
RESP * ASPECT)
Var(SECTOR * ,077 0,049
RESP * ITEM)
Var(CHAIN * ,773 0,589
RESP * ASPECT)
Var(Error) ,882 1,030
As can be seen on the previous table, the results are different. Am I
specifing a different model on R and SPSS?
Is it possible to have the output from summary(lmer()) in #,###
format, instead of scientific format?
Best regards,
Iuri.
On 8/20/06, Iuri Gavronski <[EMAIL PROTECTED]> wrote:
> 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 <[EMAIL PROTECTED]> 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 <[EMAIL PROTECTED]> 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf
> Of Iuri Gavronski
> > > Sent: Thursday, August 17, 2006 1:26 PM
> > > To: Doran, Harold
> > >
> > > Cc: r-help@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 <[EMAIL PROTECTED]> 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> 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 <[EMAIL PROTECTED]> wrote:
> > >
> > > >
> > >
> > > Iuri:
> > >
> > > The lmer function is optimal for large data with crossed random
> > > effects.
> > > How large are your data?
> > >
> > > > -----Original Message-----
> > > > From: [EMAIL PROTECTED]
> > >
> > > > [mailto: [EMAIL PROTECTED] On Behalf Of Iuri Gavronski
> > >
> > > > Sent: Thursday, August 17, 2006 11:08 AM
> > > > To: Spencer Graves
> > > > Cc: r-help@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 <[EMAIL PROTECTED] > 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 - [EMAIL PROTECTED]
> > >
> > > > > > doutorando
> > > > > > UFRGS/PPGA/NITEC - www.ppga.ufrgs.br Brazil
> > > > > >
> > > > > > ______________________________________________
> > > > > > R-help@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@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.
> > > >
> > >
> > >
> > >
> > >
> > >
> >
> >
>
______________________________________________
R-help@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.
