[R] linear mixed model question

Daniel Malter daniel at umd.edu
Mon Sep 7 19:39:53 CEST 2009


The workers=as.factor(workers) codeline in my post dropped below my name. It
should be in the code before the command line for the linear model.



Daniel Malter wrote:
> 
> 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-tp25318961p25333981.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list