1Stat 511 HW#6 Spring 2008 1. The following is a "calibration" scenario met in an engineering lab. Several resistance temperature devices (RTDs) from a large batch of the devices (that can be used to provide cheap temperature measurements in terms of resistance they register) were placed in a water bath with an extremely high quality thermometer (that we will assume produces the "true" temperature of the bath). At some regular intervals, both the "true temperature" and measured resistances were recorded. We'll let the measured resistance produced by RTD at time ijyij= and suppose that these are linear functions of the true temperatures plus random error. That is, we'll assume that for the th measured temperaturejtj= it is appropriate to model as ij i i j ijytαβε=++ where the and iiαβare an intercept and a slope specific to the particular RTD being studied. Further suppose that the and iiαβcan be described as and iii iααγ β βδ=+ =+ where and αβare unknown constants and the and iiγδare unobservable random effects. We'll assume that all the random effects have mean 0, that is that E 0 , E 0 , and E 0 ,ii ijii ijγδε=∀ =∀ =∀ We'll further assume that variances are as 22 2Var , Var , and Var ,ii ijii ijγδγσ δσ εσ=∀ =∀ =∀ and that all of the random effects are uncorrelated (have 0 covariances). We then have a model with the 5 parameters 22 2, , , , and γδαβσ σ σ The first two of these are "fixed effects" and the last three are "variance components." This scenario fits into the "Mixed Linear Model" framework introduced in class. For sake of argument, suppose that there are 3 RTDs and only 3 different measured temperatures, with respective (coded) values 12 30, 1, and 4tt t== =. a) Write out in matrix form the mixed linear model for ()11 12 13 21 22 23 31 32 33,,,,,,,, 'yyyyyyyyy=Y. (What are , , , and X β Zu?) b) What isEY ? Write out and simplify as much as you can the covariance matrix,Var Y . (Do what you can in terms of using multiplication of partitioned matrices to make the form look simple.) c) Suppose that in fact()99.8,108.1,136.0,100.3,109.5,137.7,98.3,110.1,142.2 '=Y and that I somehow know that 22 21, 1, and .25γδσσ σ== =. I can then use generalized least squares to estimate the fixed effects vector(),αβ. Do this. (Note that this is estimation of the average intercept and slope for calibration of RTDs of this type.) Does the answer change if I know only that 22 24γδσσσ==? Explain. (Indeed, can I even get a BLUE for(),αβwith only this knowledge?)2d) Suppose that I know that 22 21, 1, and 1γδσσ σ== =. Use again the data vector from part c). What is the BLUE of(),αβunder these circumstances? e) Suppose that it is your job to estimate the variance components in this model. One thing you might consider is maximum likelihood under a normal distribution assumption. This involves maximizing the likelihood as a function of the 5 parameters and it's clear how to get the profile likelihood for the variance components alone. That is, for a fixed set of variance components ()222,,γδσσσ one knows Var Y , and may use generalized least squares to estimate EY and plug that into the likelihood function in order to get the profile likelihood for ()222,,γδσσσ. Consider the two vectors of variance components used in parts c) and d) of this problem and the data vector from c). Which set has the larger value of profile likelihood (or profile loglikelihood)? f) Add the MASS and nmle packages to your R session. Then type > I3<-diag(c(1,1,1)) > J3<-matrix(c(1,1,1,1,1,1,1,1,1),3,3) > M<-matrix(c(0,0,0,0,1,4,0,4,16),3,3) > X<-matrix(c(rep(1,9),rep(c(0,1,4),3)),9,2) > Y<-matrix(c(99.8,108.1,136.0,100.3,109.5,137.7,98.3,110.1,142.2),9,1) Using these, you can create a negative profile loglikelihood function for ()222,,γδσσσ as below. (We're using a negative profile loglikelihood here because the R optimizer seeks to minimize rather than maximize a function.) > minuslLstar<-function(s2,Y) + { + temp0<-kronecker(I3,((s2[1]*J3)+(s2[2]*M)+(s2[3]*(I3)))) + temp1<-ginv(temp0) + temp2<-X%*%ginv(t(X)%*%temp1%*%X)%*%t(X)%*%temp1%*%Y + temp3<-(Y-temp2) + (.5*(log(det(temp0))))+ + (4.5*log(2*pi))+ + (.5*(t(temp3)%*%temp1%*%temp3)) + } Evaluate this function for the two sets of variance components in parts c) and d) as below and verify that your answers are consistent with what you got above. > minuslLstar(c(1,1,1),Y) > minuslLstar(c(1,1,.25),Y) I did some "hunt and peck"/by hand" optimization of this function, and came up with what we'll use as starting values of 22 2.07, .45, and .37γδσσ σ== =. Use the R function optim to find better values (MLEs) by typing > optim(c(.07,.45,.37),minuslLstar,Y=Y,hessian=TRUE)3This will set in motion an iterative optimization procedure with the starting value above. This call produces a number of interesting summaries, including a matrix of second partials of the objective function evaluated at the optimum (the Hessian). What are the MLEs? g) Now consider seeking REML estimates of the variance components. Compute the projection matrices XP and =−XNIP. Notice that since ()−=XIPX 0, every row of Ν is (after transposing) in ()C⊥′X and so any set of 7 rows of Ν that make a 79×matrix, say B , of rank 7 will serve to create the vector of "error contrasts" =rBY which one uses to find REML estimates here. Verify that the first 7 rows of Ν will work in the present situation by typing > B<-rbind(N[1,],N[2,],N[3,],N[4,],N[5,],N[6,],N[7,]) > qr(B)$rank Note that by construction, ()()()()2227~MVN , , , 'γδσσσr0BV B and "REML" is maximum likelihood for the variance components based on r . You may create a negative loglikelihood function here as > minuslLstar2<-function(s2,Y,B) + { + temp0<-kronecker(I3,((s2[1]*J3)+(s2[2]*M)+(s2[3]*(I3)))) + temp1<-B%*%temp0%*%t(B) + temp2<-ginv(temp1) + temp3<-B%*%Y + (.5*(log(det(temp1))))+ + (3.5*log(2*pi))+ + (.5*(t(temp3)%*%temp2%*%temp3)) + } Then you may find REML estimates of the variance components by optimizing this function. Use the starting values from part f) and type > optim(c(.07,.45,.37),minuslLstar2,Y=Y,B=B,hessian=TRUE) How do your estimates here compare to the MLEs you found in part f)? (Statistical folklore says that usually REML estimates are bigger than MLEs.) h) The (unavailable) BLUE of any (vector of) estimable function(s) Cβ is