ST 762, HOMEWORK 6, FALL 2009These problems are to be turned in on the due date.1. Recall the data in Homework 5, Problem 3, from a clinical trial studying the effectiven ess ofa treatment for patients with respiratory illness. See that problem for a description of thestudy and the contents of the data file. As in that p roblem, let Yijbe the respiratory statusfor patient i at week tij, where tij= 0, 1, 2, 3, 4 weeks f or j = 1, . . . , 5 for all subjects andi = 1, . . . , 111. Let δi= 0 if patient i was assigned to placebo and 1 if assigned to activetreatment. Let wij= 0 if tij= 0 and wij= 1 if tij> 0.In that problem, we adopted a population-average perspective, focusing on the modelE(Yij|δi) = P (Yij= 1|δi) =exp(β1+ β2wij+ β3wijδi)1 + exp(β1+ β2wij+ β3wijδi), j = 1, . . . , 5. (1)We also assumed a working model for the marginal covariance matrix for the Yijgiven δi, andwe fit the overall population-average mean-covariance model using “GEE-1” metho ds, withthe covariance parameters estimated different ways as implemented in SAS proc genmod,proc glimmix, and the nlinmix macro.An alternative app roach is to adopt a subject-specific perspective, and instead start with amodel of the formE(Yij|δi, bi) = P (Yij= 1|δi, bi) =exp(β1+ β2wij+ β3wijδi+ bi)1 + exp(β1+ β2wij+ β3wijδi+ bi), (2)where bi∼ N (0, D) is a scalar random effect with variance D. Following the discussion onpages 355–356, this is a different model from (1) in th e sense that it implies a different spec-ification for E(Yij|δi) upon integration over the distrib ution of the random effect. Moreover,the interpretations of the parameters β1, β2, β3in (1) are different from those for β1, β2, β3in(2).(a) Give an interpretation of the parameters β1, β2, β3in (2). In particular, what issue doesβ3address?(b) Fit model (2) to the data using proc glimmix with the method=mspl option and againwith the method=quad option. (You will have to consult th e documentation for a descriptionof what these options do.)(c) Fit the model again using proc nlmixed with the default option for carrying out theintegration (adaptive Gaussian quadrature).(d) Fit the model again using the glimmix macro with method=ml in the procopt statement.(e) C ompare the results of the fits in (b)–(d). Which ones would you expect to be similar, andare they? Also compare the results to those obtained from fitting the mod el (1) in Homework5. Are the results comparable? Should they be?2. Interferon-α-2b (IFN) is one of many experimental therapies being evaluated f or treatingpatients infected with hepatitis C virus (HCV). A key way of evaluating such treatmentsis in terms of h ow they impact the within-patient dynamics of the HCV virus. That is,investigators are interested in understanding the beh avior of viral load, wh ich is roughly ameasure of the concentration of virus in the patient’s system, following the start of IFNtherapy.1Researchers have developed mathematical models in the form of systems of differential equa-tions to characterize th e dynamics, i.e., how the virus behaves in the body and how the bodyreacts and how IFN therapy may affect these processes over time. I n a recent article in thejournal Science (Volume 282, October 2, 1998, p. 103), Neumann et al. propose such a modelto describe the dynamics over the first two days of IFN therapy. If the differential equationsare solved, they imply that, if V (t) is the viral load at time x (days) following the start ofIFN therapy at day x = 0,V (x) = V0[1 − ǫ + ǫ exp{−c(x − t0)+}], 0 ≤ x ≤ 2, V0, c, ǫ > 0, (3)where x+= x for x > 0, and 0 otherwise; V0is the viral load at day 0 (prior to start oftherapy); c is the virion clearance rate; t0is the so-called pharmacological delay such thatdecay in viral load is not seen until IFN has made sufficient progress in distributing throughthe body; and ǫ is an efficacy parameter. The efficacy parameter h as the interpretation thatthe effect of IFN therapy is to reduce the production of new virus particles, or virions, in thesystem by the fraction (1 − ǫ); if ǫ = 1 then the therapy is said to be “perfect.” Here, we willtake the pharmacological delay to be known: t0= 0.20 days.In the file hcvmix.dat on the class web page you will find data from 30 HCV-infected subjectswho began IFN therapy on day 0. The first column is subject id number, the second columnis day (measured from day 0) at which a measurement of HCV viral load (copies of HC Vviral RNA per ml) was taken, and the thir d column is the viral load itself divided by 106.For each subject, viral load measurements were taken over the firs t two days of therapy.For subject i, assume that f (zij, βi) is given as above, where zijis the jth day value forthis su bject, βi= (β1i, β2i, β3i)T, β1i= log(V0i) is the logarithm of the viral load at day0, β2i= log(ci) is the log of virion clearance rate, and β3i= log(ǫi) is the log of efficacy.Viral loads are well-known to be subject to constant co efficient of variation at the individuallevel. Moreover, the researchers were willing to assume that within-subject correlation isnegligible. They proposed the following subject-specific model. Letting Yijbe the jth viralload measurement for the ith subject,E(Yij|zij, βi) = f (zij, βi), var(Yij|zij, βi) = σ2f2θ(zij, βi),where θ is kn own. The variance parameters σ2and θ are common across subjects, reflectingthe belief of a similar pattern of intra-subject variation. The second-stage model is given asβi= Aiβ + Bibi,where Aiand Biare (3 × 3) identity matrices, β = (β1, β2, β3)T, and bi= (bi1, bi2, bi3)T.The researchers would like to fit this model and get estimates and standard errors for thetypical values of log day-0 viral load, log clearance rate, and log efficacy.You will help them by doing this three ways:(a) Use the SAS macro nlinmix to fi t this model to the data using the “refined” approximationdiscussed on p ages 449–452 of the notes. Specify method=ml in the procopt statement.Because nlinmix cannot estimate intra-individual variance parameters, set θ equal to somelikely value of your choice. (You may wish to revise your choice after seeing the answers to(b) and (c).)(b) Use the R function nlme(). Let nlme() estimate the variance power parameter θ from thedata, as in the example on the web page, so you can see if the data s upport the investigators’2choice of θ = 1. Full information on the current syntax for the nlme() function