1Stat 511 HW#8 Spring 2004 This assignment consists of a number of problems on mixed linear models. Some of the problems are repeats of ones from HW 10 Spring 2003. You are responsible for the content of all the problems. However, only problems 2,3,4 and 6 should be turned in. 1. Below is a very small set of fake unbalanced 2 level nested data. Consider the analysis of these under the mixed effects model ijk i ij ijkyµαβε=++ + where the only fixed effect isµ, the random effects are all independent with the()2 iid N 0,iαασ, the()2 iid N 0,ijββσ, and the()2 iid N 0,ijkεσ. Level of A Level of B within A Response 1 1 6.0, 6.1 2 8.6, 7.1, 6.5, 7.4 2 1 9.4, 9.9 2 9.5, 7.5 3 6.4, 9.1, 8.7 Again after adding the MASS and nlme libraries to your R session (and specifying the use of the sum restrictions), do the following. a) Enter your data by typing > y<-c(6.0,6.1,8.6,7.1,6.5,7.4,9.4,9.9,9.5,7.5,6.4,9.1,8.7) > A<-c(1,1,1,1,1,1,2,2,2,2,2,2,2) > B<-c(1,1,2,2,2,2,1,1,2,2,3,3,3) Run and view some summaries of a random effects/mixed effects analysis of these data by typing > lmef.out<-lme(y~1,random=~1|A/B) > summary(lmef.out) Get approximate confidence intervals for the parameters of the mixed model by typing > intervals(lmef.out) Compute an exact confidence interval for σ based on the mean square error (or pooled variance from the 5 samples of sizes 2,4,2,2, and 3). How do these limits compare to what R provides in this analysis? b) A fixed effects analysis can be made here and some summaries viewed here by typing > AA<-as.factor(A) > BB<-as.factor(B) > lmf.out<-lm(y~1+AA/BB) > summary(lmf.out)2(Note that the estimate of σ produced by this analysis is exactly the one based on the mean square error referred to in a).) Type > predict(lmef.out2) > predict(lmf.out) Identify the predictions from the fixed effects model as simple functions of the data values. (What are these predictions?) Notice that the predictions from the mixed effects analysis are substantially different from those based on the fixed effects model. 2. The article “Variability of Sliver Weights at Different Carding Stages and a Suggested Sampling Plan for Jute Processing” by A. Lahiri (Journal of the Textile Institute, 1990) concerns the partitioning of variability in “sliver weight.” (A sliver is a continuous strand of loose, untwisted wool, cotton, etc., produced along the way to making yarn.) For a particular mill, 3 (of many) machines were studied, using 5 (10 mm) pieces of sliver cut from each of 5 rolls produced on the machines. The weights of the (75) pieces of sliver were determined and a standard hierarchical (balanced data) ANOVA table was produced as below. (The units of weight were not given in the original article.) Source SS df Machines 1966 2 Rolls 644 12 Pieces 280 60 Total 2890 74 Use the same mixed effects model as in problem 1 and do the following. a) Make 95% confidence intervals for each of the 3 standard deviations , and αβσσσ. Based on these, where do you judge the largest part of variation in measured weight to come from? (From differences between pieces for a given roll? From differences between rolls from a given machine? From differences between machines?) (Use Cochran- Satterthwaite for and αβσσ and give an exact interval for σ.) b) Suppose for sake of illustration that the grand average of all 75 weight measurements was in fact ...35.0y= . Use this and make a 95% confidence interval for the model parameterµ. 3. According to a posting to an R-help list by Doug Bates “It happens that the lme function is much better suited to nested random effects than to crossed random effects. To estimate crossed random effects you must create an awkward formulation with a grouping factor that has one level and the random effects model matrix based on based on the indicators for factor a and for factor b. These two sets of random effects each have variance-covariance matrix that are multiples of an identity and are grouped together as a block-diagonal matrix.” (The context here is a 2-way factorial and a “no-interaction” analysis.) Consider the fake unbalanced 33× factorial data of problem 4 on HW 5. Here let’s attempt a random effects analyses of those th response at level of A and level of Bijkyk i j= based on a two-way random effects model without interaction3 ijk i j ijyµαβε=++ + where µ is an unknown constant, the iα are iid ()2N0,ασ, the jβ are iid ()2N0,βσ, the ijε are iid ()2N0,σ and all the sets of random effects are all independent. Load the nlme and MASS packages and then read the data into R and prepare a first appropriate “groupedData” object using the following commands. > Y<-c(12,14,10,12,9,11,12,6,7,10,11,7) > a<-c(1,1,1,1,2,2,2,2,2,3,3,3) > b<-c(1,2,3,3,1,2,2,3,3,1,2,3) > group<-c(1,1,1,1,1,1,1,1,1,1,1,1) > A<-as.factor(a) > B<-as.factor(b) > Group<-as.factor(group) > Fake<-data.frame(Y,A,B,Group) > FakeGrouped<-groupedData(Y~1|Group,Fake) Then the following command will run the mixed effects analysis on these data. > twoway.out<-lme(Y~1,data=FakeGrouped,random=pdBlocked(list(pdIdent(~A-1),pdIdent(~B-1)))) Run the summary(), predict(), and intervals() functions on the object twoway.out. What are 95% confidence limits for the 4 model parameters , , , and αβµσσ σ? 4. The data set below is taken from page 54 of Statistical Quality Assurance Methods for Engineers by Vardeman and Jobe. It gives burst strength measurements made by 5 different technicians on small pieces of paper cut from 2 different large sheets (each technician tested 2 small pieces from both sheets). Technician (B) 1 2 3 4 5 1 13.5, 14.8 10.5, 11.7 12.9, 12.0 8.8, 13.5 12.4, 16.0Sheet (A) 2 11.3, 12.0 14.0, 12.5 13.0, 13.1 12.6, 12.7 11.0, 10.6 Vardeman and Jobe present an ANOVA table for these data based on a two-way random effects model with interaction ijk i j ij ijkyµαβαβε=+ + + + (where µ is the only fixed effect). Source SS dfEMS Sheet (A) .5445 122210 2ααβσσσ++Technician (B) 2.692 422242βαβσσσ++ Sheet*Technician (A*B) 24.498 4222αβσσ+ Error 20.955 102σ Total 48.6895 194a) Use R and an ordinary linear model analysis to verify that the sums of squares above are correct. b) In the context of a “gauge study” like this, σ is usually called the “repeatability” standard deviation. Give a 95% exact confidence interval for σ (based