```Date: Wed, 23 Jul 2003 12:17:07 -0700 Reply-To: Dale McLerran Sender: "SAS(r) Discussion" From: Dale McLerran Subject: Re: GEE with multiple treatments per subject: how to model? Comments: To: Margarita In-Reply-To: <3d935f86.0307220914.35274940@posting.google.com> Content-Type: text/plain; charset=us-ascii Margarita, There are a couple of possible ways to analyze these data. One possible approach would be to employ a bootstrap process, although I am not certain right off-hand what the resample units should be. Given that you have correlations due to both visit and treatment and you have an ordinal response, I would suggest that you try to fit an ordinal response mixed model using the procedure NLMIXED. Note that other correspondents suggested GLIMMIX (which is not a procedure, but a macro) and MIXED. Neither of these would be suitable for an ordinal response variable. However, the NLMIXED procedure has the flexibility to fit a proportional odds mixed model. The only problem that you might face with NLMIXED is the number of random effects that you must estimate. NLMIXED is limited in this regard, at least in comparison with the MIXED procedure. To use the NLMIXED procedure, you will have to write down the likelihood function for an ordinal response subject to random effects. For an ordinal response with k+1 levels of response, the probability that y=i is given by p(y=0) = exp(eta0) / (1 + exp(eta0)) p(y=1) = exp(eta1) / (1 + exp(eta1)) - exp(eta0) / (1 + exp(eta0)) p(y=2) = exp(eta2) / (1 + exp(eta2)) - exp(eta1) / (1 + exp(eta1)) ... p(y=k) = 1 - exp(eta[k-2]) / (1 + exp(eta[k-2)) where eta0, eta1, ..., eta[k-2] are linear combinations of the predictor variables, including any random effects. Given the complexity of your model, I would suggest using a common random intercept for all k-1 linear combinations for the construction of p(y=i) by treatment and visit. That is, I would suggest eta0 = trt1_0*(trt=1) + trt2_0*(trt=2) + trt3_0*(trt=3) + trt1_vis2*(trt=1)*(visit=2) + trt1_vis3*(trt=1)*(visit=3) + trt2_vis2*(trt=2)*(visit=2) + trt2_vis3*(trt=2)*(visit=3) + trt3_vis2*(trt=3)*(visit=2) + trt3_vis3*(trt=3)*(visit=3) + gtrt1_visit1*(trt=1)*(visit=1) + gtrt1_visit2*(trt=1)*(visit=2) + gtrt1_visit3*(trt=1)*(visit=3) + gtrt2_visit1*(trt=2)*(visit=1) + gtrt2_visit2*(trt=2)*(visit=2) + gtrt2_visit3*(trt=2)*(visit=3) + gtrt3_visit1*(trt=3)*(visit=1) + gtrt3_visit2*(trt=3)*(visit=2) + gtrt3_visit3*(trt=3)*(visit=3); eta1 = trt1_1*(trt=1) + trt2_1*(trt=2) + trt3_1*(trt=3) + trt1_vis2*(trt=1)*(visit=2) + trt1_vis3*(trt=1)*(visit=3) + trt2_vis2*(trt=2)*(visit=2) + trt2_vis3*(trt=2)*(visit=3) + trt3_vis2*(trt=3)*(visit=2) + trt3_vis3*(trt=3)*(visit=3) + gtrt1_visit1*(trt=1)*(visit=1) + gtrt1_visit2*(trt=1)*(visit=2) + gtrt1_visit3*(trt=1)*(visit=3) + gtrt2_visit1*(trt=2)*(visit=1) + gtrt2_visit2*(trt=2)*(visit=2) + gtrt2_visit3*(trt=2)*(visit=3) + gtrt3_visit1*(trt=3)*(visit=1) + gtrt3_visit2*(trt=3)*(visit=2) + gtrt3_visit3*(trt=3)*(visit=3); eta2 = trt1_2*(trt=1) + trt2_2*(trt=2) + trt3_2*(trt=3) + trt1_vis2*(trt=1)*(visit=2) + trt1_vis3*(trt=1)*(visit=3) + trt2_vis2*(trt=2)*(visit=2) + trt2_vis3*(trt=2)*(visit=3) + trt3_vis2*(trt=3)*(visit=2) + trt3_vis3*(trt=3)*(visit=3) + gtrt1_visit1*(trt=1)*(visit=1) + gtrt1_visit2*(trt=1)*(visit=2) + gtrt1_visit3*(trt=1)*(visit=3) + gtrt2_visit1*(trt=2)*(visit=1) + gtrt2_visit2*(trt=2)*(visit=2) + gtrt2_visit3*(trt=2)*(visit=3) + gtrt3_visit1*(trt=3)*(visit=1) + gtrt3_visit2*(trt=3)*(visit=2) + gtrt3_visit3*(trt=3)*(visit=3); eta3 = trt1_3*(trt=1) + trt2_3*(trt=2) + trt3_3*(trt=3) + trt1_vis2*(trt=1)*(visit=2) + trt1_vis3*(trt=1)*(visit=3) + trt2_vis2*(trt=2)*(visit=2) + trt2_vis3*(trt=2)*(visit=3) + trt3_vis2*(trt=3)*(visit=2) + trt3_vis3*(trt=3)*(visit=3) + gtrt1_visit1*(trt=1)*(visit=1) + gtrt1_visit2*(trt=1)*(visit=2) + gtrt1_visit3*(trt=1)*(visit=3) + gtrt2_visit1*(trt=2)*(visit=1) + gtrt2_visit2*(trt=2)*(visit=2) + gtrt2_visit3*(trt=2)*(visit=3) + gtrt3_visit1*(trt=3)*(visit=1) + gtrt3_visit2*(trt=3)*(visit=2) + gtrt3_visit3*(trt=3)*(visit=3); You will observe that the effects trt1_i, trt2_i, and trt3_i are intercepts modelling the probability of response level i for each of the three outcomes at the first visit. The terms trt1_vis2, trt1_vis3, ..., trt3_vis3 are deviations from the intercepts at subsequent visits. We are fitting a proportional odds model, so that the same deviations from the intercept apply to all levels of the response. Terms gtrt1_visit1, gtrt1_visit2, ..., gtrt3_visit3 are visit by treatment specific random intercept terms that apply to all levels of the response variable. Having constructed eta[i] and p[i], the log likelihood contribution can be constructed as log(p[i]). We can now put all this together as NLMIXED code. I will construct a model in which the random effects have an unstructured covariance for the multivariate response, and have an AR(1) structure over time. You may wish to modify this to have a compound symmetric error structure over time or an unstructured error structure for the repeated measures component. proc nlmixed data=mydata; parms b0=0 b1=0 b2=0 ... b[k-1]=0 trt1=0 trt2=0 visit1=0 visit2=0 trt1_visit1=0 trt1_visit2=0 trt2_visit1=0 trt2_visit2=0 logsig_v1=0 logsig_v2=0 logsig_v3=0 rho12=0 rho13=0 rho23=0 rho_t=0; bounds -1 wrote: > Jay Weedon wrote in message > news:... > > On 21 Jul 2003 15:54:40 -0700, rlopatin@yahoo.com (Margarita) > wrote: > > > > >-Hello all, > > > > > >I am trying to fit a GEE model for the following study design: > every > > >subject in a study was exposed to all study treatments (the > treatments > > >are local in nature and do not interact) at baseline. The > response, an > > >irritation score on the ordinal scale with more than two levels, > is > > >recoded at each treatment site at several time points. Some > people > > >are more likely to have irritation than others (no matter what > > >treatment is applied) so one would expect a correlation between > > >treatment responses within a subject, as well as between time > points. > > > The objective is to see if some treatments are more likely to > lead to > > >higher irritation scores than others. > > > > > >I would appreciate any advice on what would be the best way to > model > > >this either using SAS version 8 procedures, Splus or custom > macros. > > > > > >I've tried proc GENMOD so far but am not convinced it's doing the > > >right thing. It allows only independent working correlation > matrix > > >with multinomial response, which it seems would ignore the > correlation > > >between time measurements & treatment responses. > > > > Post details of your problem and SAS code to comp.soft-sys.sas. > > > > JW > > -Jay, > > Here is an example of the data: > subject trt visit y > 1 A 1 0 > 1 B 1 1 > 1 C 1 1 > 2 A 1 1 > 2 B 1 1 > 2 C 1 0 > .... > 1 A 2 0 > 1 B 2 0 > 1 C 2 1 > .... > Each subject was exposed to, say, 3 treatments at the start of the > study. The irritation levels are recorded on ordinal scale (0 to 4) > at several visits. I used the following code for analysis: > > proc genmod data=irritation descending; > class trt visit subject; > model y= trt visit /dist=mult link=clogit type3; > repeated subject=subject / type=ind withinsubject=visit*trt; > run; > > This what I received from the SAS technical support: > "Only TYPE=IND is available for multinomial responses. Once there is > ample statistical literature supporting the other correlation > matrices, they will be added to GENMOD. Currently this is an area of > research for many statisticians. > > Something to keep in mind is that the Emperical based standard errors > make no assumptions regarding the correlation matrix. They require > only that your mean model be correctly specified. This means that > they are robust no matter what your choice of correlation matrix." > > So, it seems there is no way to model ordinal data accounting for > within subject correlation in SAS. If anyone heard of a program that > can do that, please let me know. > > I have received a reference to GLIMMIX from Iyue Sung. My assumption > has been that GLM deals with continuos responses but I'll check it > out, just in case. > > In terms of analysis, I also tried dichotomizing the response and > analyze it with AR(1) correlation matrix. The results seem to be > consistent. > > Thanks > Margarita ===== --------------------------------------- Dale McLerran Fred Hutchinson Cancer Research Center mailto: dmclerra@fhcrc.org Ph: (206) 667-2926 Fax: (206) 667-5977 --------------------------------------- __________________________________ Do you Yahoo!? Yahoo! SiteBuilder - Free, easy-to-use web site design software http://sitebuilder.yahoo.com ```

Back to: Top of message | Previous page | Main SAS-L page