s-news
[Top] [All Lists]

Re: passing arguments to a function

To: s-news@lists.biostat.wustl.edu
Subject: Re: passing arguments to a function
From: Pravin Jadhav <pravinj@gmail.com>
Date: Mon, 19 Dec 2005 16:59:18 -0500
Domainkey-signature: a=rsa-sha1; q=dns; c=nofws; s=beta; d=gmail.com; h=received:message-id:date:from:to:subject:in-reply-to:mime-version:content-type:references; b=Shy6rrOf+gHfKlNlcy1LUQ5xnAbpAeFPHmG4rKoYNc0vwP6KCYTJil7ObkwCZ4TVvy16UV5H5Ng7+PSQdHpU284CpXW5Wt1sSgPOysq5yMenHZDLtpIGFZG2wCmgzI5sUwmpL6ko2+EQhXjNf3Q2euiixXSh8cQxMw2BqMFpvVA=
In-reply-to: <679f8b3f0512191224l1173159dwa92e0f7f574e7329@mail.gmail.com>
References: <679f8b3f0512191224l1173159dwa92e0f7f574e7329@mail.gmail.com>
Hi again,

I have made it work with a very straightforward modification (rather ugly one). Please ignore the previous message unless you have a better solution for me.

fit.F<-function (data, byvar, fullF, ranF,...)
                {
                  ......
                  fit<-lme(.........)
                  fit$call$fixed<-fullF    #ADDITION
                  pred.fit<-
                  ........ 
                 }
Thanks,

Pravin

On 12/19/05, Pravin Jadhav <pravinj@gmail.com> wrote:
Dear all,

For my previous query (obtain predictions for the missing records), Dr. Pinheiro pointed me to the "predict" function. Thanks! It is doing exactly what I needed. But I am having some trouble putting it inside a generic function. Here is an example (dummy):

I want to describe response (resp) as a function of time (time) using a linear mixed effects model for several trials with subjects (id) as a grouping variable. Of the 3 trials, the following code evaluates ONLY the 1st  trial. Then the lmeObject is used for prediction. Although, there is no missing data here, the actual dataset has some missing observations. This step works just fine.

       #CREATE DATA
       data<-data.frame(trial=rep(1:3,each=30),id=rep(rep(1:10,each=3),3),time=rep(c(1:3),30))
       data$resp<-rep(rnorm(30,5,1),each=3)+rep(rnorm(30,0.5,0.2),each=3)*time+rnorm(90,0,1)
      #FIT DATA AND PREDICT
       fit1<-lme(resp~time,data="" na.action="">        pred1<-predict(fit1, newdata=data, level=0:1)

Now, if I take this inside a generic funciton (only the relevant part is shown) to evaluate all the trials, there is an error. The first part, LME works fine but has problems with predict function. It turns out that the lmeObject, "fit", does retain the information about the linear two sided formula used to describe fixed effcts. If you take a look at the object "fit1" from the previous step and "fit2" from the next step, the difference is clear. "fit1" has explicit two sided formla, whereas, "fit2" retains only list(fullF) and thus causing problems in the predict step.

fit.F<-function (data, byvar, fullF, ranF,...)
                {
                    fit<-by(data, byvar, function(bydata, fullF, ranF,...)
                                                {
                                                    fit<-lme(fixed=fullF, data="" random=ranF, method="ML", na.action="">                                                 #    pred.fit<-predict(fit,newdata=bydata, level=0:1)
                                                    return(fit)
                                                #    retrun(pred.fit)
                                                },fullF=fullF,ranF=ranF)
                }
fit2<-fit.F(data, data$trial, resp~time,list(id=pdDiag(~1+time)))

How do I tweak the function so that the predict step is also working?

Sorry for a lengthy email (if you are still reading it). Any comments are greatly appreciated.

Thanks,

Pravin


<Prev in Thread] Current Thread [Next in Thread>