Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? Thanks a lot! Wen ################# library(nlme) data(Machines) new.data = Machines[c(1:6, 25:30, 49:54), ] fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) fm1 fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) fm2 intervals(fm2)
Hello Wen: On 2009.09.06 10:49:03, Wen Huang wrote:> Hello, > > I wanted to fit a linear mixed model to a data that is similar in > terms of design to the 'Machines' data in 'nlme' package except that > each worker (with triplicates) only operates one machine. I created a > subset of observations from 'Machines' data such that it looks the > same as the data I wanted to fit the model with (see code below). > > I fitted a model in which 'Machine' was a fixed effect and 'Worker' > was random (intercept), which ran perfectly. Then I decided to > complicate the model a little bit by fitting 'Worker' within > 'Machine', which was saying variation among workers was nested within > each machine. The model could be fitted by 'lme', but when I tried to > get > confidence intervals by 'intervals(fm2)' it gave me an error: > > Error in intervals.lme(fm2) : > Cannot get confidence intervals on var-cov components: Non-positive > definite approximate variance-covariance > > I am wondering if this is because it is impossible to fit a model like > 'fm2' or there is some other reasons?The problem doesn't seem to be the model specification but is most likely the result of estimating a more complicated model with very little data. Using the complete Machines dataset with the same model specification seems to work fine: # ----------------------------------------------------------------------------- #> fm3 <- lme(score ~ Machine, random = ~ Machine - 1 | Worker, data = Machines) > intervals(fm3)Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 48.972459 52.355556 55.73865 MachineB 3.093747 7.966667 12.83959 MachineC 10.816607 13.916667 17.01673 attr(,"label") [1] "Fixed effects:" Random Effects: Level: Worker lower est. upper sd(MachineA) 2.1702468 4.0792807 7.6675752 sd(MachineB) 4.6301082 8.6252908 16.0677975 sd(MachineC) 2.3387870 4.3894795 8.2382579 cor(MachineA,MachineB) 0.1992744 0.8027499 0.9647702 cor(MachineA,MachineC) -0.1702480 0.6225047 0.9260744 cor(MachineB,MachineC) 0.1235115 0.7708309 0.9579666 Within-group standard error: lower est. upper 0.7629124 0.9615766 1.2119736 # ----------------------------------------------------------------------------- # With the restricted dataset, there are only 18 observations in 6 groups. This is probably too little data for the (restricted) maximum likelihood technique used by lme(). Hope that helps, ~Jason -- Jason W. Morgan Graduate Student Department of Political Science *The Ohio State University* 154 North Oval Mall Columbus, Ohio 43210
Dear Wen, Since each worker only works on one machine, your model fm2 does not make sense. Your random effects tries to model how the effect of each worker differs between machines. But you don't have that kind of information if a work only works one machine. HTH, Thierry ------------------------------------------------------------------------ ---- ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 Thierry.Onkelinx at inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -----Oorspronkelijk bericht----- Van: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] Namens Wen Huang Verzonden: zondag 6 september 2009 17:49 Aan: r-help at r-project.org Onderwerp: [R] linear mixed model question Hello, I wanted to fit a linear mixed model to a data that is similar in terms of design to the 'Machines' data in 'nlme' package except that each worker (with triplicates) only operates one machine. I created a subset of observations from 'Machines' data such that it looks the same as the data I wanted to fit the model with (see code below). I fitted a model in which 'Machine' was a fixed effect and 'Worker' was random (intercept), which ran perfectly. Then I decided to complicate the model a little bit by fitting 'Worker' within 'Machine', which was saying variation among workers was nested within each machine. The model could be fitted by 'lme', but when I tried to get confidence intervals by 'intervals(fm2)' it gave me an error: Error in intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance I am wondering if this is because it is impossible to fit a model like 'fm2' or there is some other reasons? Thanks a lot! Wen ################# library(nlme) data(Machines) new.data = Machines[c(1:6, 25:30, 49:54), ] fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) fm1 fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) fm2 intervals(fm2) ______________________________________________ R-help at r-project.org 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. Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.
Wen, to follow up on Thierry, your workers are nested in machines (since each worker only works one machine). Consider fitting a nested model. Though, with few observations, you might run into the same problem. Further, if you have observation triplets, and you expect systematic differences between each trial, then you would have to include a trial effect in some way. Have you plotted the data? This can be very informative. Finally, you may want to consider a fixed-effects approach. Throw in worker dummies only (without intercept) and test the linear hypothesis that the sum of the worker dummy coefficients is equal between two machines. The following example does that for two machines (if you have n machines you would have binomial coefficient (n,2) linear hypotheses): #simulate data machines=rep(c(0,1),each=9) workers=rep(c(1:6),each=3) workers.re=rep(rnorm(6),each=3) error=rnorm(18,0.3) output=2*machines+workers.re+error #code machines and workers as factors/dummies machines=as.factor(machines) #run OLS (fixed effects) of output on worker dummies without intercept fm4=lm(output~-1+workers) summary(fm4) #test the linear hypothesis that coefficients on the worker dummies are equal for both machines #the equivalent formulation used is: is the sum of the coefficients across machines equal to zero? library(car) linear.hypothesis(fm4,hypothesis.matrix=c(1,1,1,-1,-1,-1),rhs=c(0)) hope that helps, daniel workers=as.factor(workers) Wen Huang-3 wrote:> > Hello, > > I wanted to fit a linear mixed model to a data that is similar in > terms of design to the 'Machines' data in 'nlme' package except that > each worker (with triplicates) only operates one machine. I created a > subset of observations from 'Machines' data such that it looks the > same as the data I wanted to fit the model with (see code below). > > I fitted a model in which 'Machine' was a fixed effect and 'Worker' > was random (intercept), which ran perfectly. Then I decided to > complicate the model a little bit by fitting 'Worker' within > 'Machine', which was saying variation among workers was nested within > each machine. The model could be fitted by 'lme', but when I tried to > get > confidence intervals by 'intervals(fm2)' it gave me an error: > > Error in intervals.lme(fm2) : > Cannot get confidence intervals on var-cov components: Non-positive > definite approximate variance-covariance > > I am wondering if this is because it is impossible to fit a model like > 'fm2' or there is some other reasons? > > Thanks a lot! > > Wen > > ################# > > library(nlme) > data(Machines) > new.data = Machines[c(1:6, 25:30, 49:54), ] > fm1 = lme(score ~ Machine, random = ~1|Worker, data = new.data) > fm1 > fm2 = lme(score ~ Machine, random = ~Machine-1|Worker, data = new.data) > fm2 > intervals(fm2) > > ______________________________________________ > R-help at r-project.org 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. > >-- View this message in context: http://www.nabble.com/linear-mixed-model-question-tp25318961p25333956.html Sent from the R help mailing list archive at Nabble.com.