Hi
Pravin,
There
are numerous answers to this question. A (not extensive) set of possible
answers include:-
·
Use a loop
·
Apply over columns of a matrix
containing your dependant variables
·
Lapply your function over a vector
of dependant column names
·
Stack your dependant variables
into a column (called "Y", say with a vector of dependant names
called "myVar", say) then do a by/tapply "by" rep AND “myVar”
**
The
slowest bit of your code may involve your call to "by". The “by”
function simply tapply's an index to a dataset, which means the dataset is read
in over and over again. I tend avoid using “by” for this reason
(+ the fact I don’t like the by class anyway).
The
choice of solution will depend very much on the amount of data you are using + what
you want back from the model fit. For example, if you’re just
looking to grab a log likelihood, you could tapply by factors rep & “myVar”
(**) and return just one value into a neat matrix of numbers. However, if
you want to return a large number of models and interrogate them afterwards,
you’ll need to think about the way in which you would reference the output
list / by object.
“An”
answer to this could look something like …..
# Setup dummy data
myDf <-
data.frame(matrix(rnorm(1000), 100)) #
Get data frame of X data
mySample <-
factor(sample(letters[1:3], 100, T)) #
Get sample "by" vector
xVar <-
rnorm(100) #
Get X data (time in your example)
# Set up function to use
myFun <- function(df) try(lme(y~x,
data="" #
Function to apply and fit model
# Do the work
myFits <- apply(myDf, 2, function(y,
x, z) by(cbind(y, x), z, myFun, simplify=F), x=xVar, z=mySample)
# Example of using output = Getting a
matrix of log likelihood values
sapply(myFits, function(byObj) sapply(byObj,
function(x) if (class(x) == "Error") return(NA) else
return(x$logLik)))
In
the above example, “myDf” contains my dependant values (a data
frame or matrix containing resp1 … resp10 in your example). The “xVar”
vector is effectively your “time” variable. The “mySample”
vector is acting as your “rep” column.
I
have fitted a basic model within a random term, and used the “try”
function to capture any models that blow up …
IMPORTANT:
The above code has been quickly put together to illustrate the way an “apply-based”
answer be implemented. It has not really been put together in the most “efficient”
way, or the “most robust” way (eg. I’ve not thought about
factor levels, column names, etc).
I
can, however, recommend a very good book which should help you create some
efficient code: “S Programming” by Venables & Ripley.
Good
luck,
Rich.
S-PLUS & R Consulting & Training
mangosolutions
Tel
+44 118 902 6617
Fax
+44 118 902 6401
-----Original
Message-----
From: s-news-owner@lists.biostat.wustl.edu
[mailto:s-news-owner@lists.biostat.wustl.edu] On Behalf Of Pravin Jadhav
Sent: 17 February 2005 22:20
To: s-news@lists.biostat.wustl.edu
Subject: [S] help with a function
Hello,
I
am writing a a small function to fit an LME model.
fit.data<-function(data)
{
lme.fit<-by(data, data$rep, function (x) lme(resp~time, data="">
random=list(id=pdDiag(~1+time))))
return(lme.fit)
}
It
WORKS for fitting response (resp) versus time for a given dataset
by
replications (rep).
PROBLEM:
Now, I wanted to fit other variables in the same dataset. One
option
(inefficient) was to rewrite the function for resp1,
resp2,.......respN
and do the fitting.
However,
I would like to pass on the vector information to the
function
used in "by". I know the following would not work..because it
is
not supposed to...but I tried.
fit.dataN<-function(data,vec)
{
lme.fit<-by(data, data$rep, function (x) lme(vec~time, data="">
random=list(id=pdDiag(~1+time))))
}
So
that I can do the following
fit1<-fit.dataN(data,vec1)
fit2<-fit.dataN(data,vec2)
.............
fitN<-fit.dataN(data,vecN)
Any
suggestions? How to pass on the dataset as well as vector information.
Thanks.
Pravin
--
Pravin
Jadhav
Graduate
Student
Department
of Pharmaceutics
MCV/Virginia
Commonwealth University
DPE1/CDER/OCPB/Food
and Drug Administration
Phone:
(301) 594-5652
Fax:
(301) 480-3212
--------------------------------------------------------------------
This
message was distributed by s-news@lists.biostat.wustl.edu. To
unsubscribe
send e-mail to s-news-request@lists.biostat.wustl.edu with
the
BODY of the message: unsubscribe s-news