args.ellipsis
From: Steven Paul Millard
To: 'Nicole DePriest Demers'
Cc: S-news (E-mail)
Subject: RE: [S] Integration difficulties
Date: Wednesday, March 11, 1998 2:06 PM
Hello,
Nicole DePriest Demers wrote about difficulty using the integrate()
function with a function that takes an argument called "k". The
problem with integrate(), nlmin(), nlminb(), and a few other functions
like these is that they have their own arguments as well as the
ellipsis argument (...), and S-PLUS does partial matching on argument
names.
For the integrate() function, the ellipsis argument passes argument
values from the call to integrate() to the function f that is being
integrated. Unfortunately, Nicole's function moment() takes an
argument called "k", but integrate() has an argument called "keep.xy".
When you call integrate() and set k=2, integrate thinks you are trying
to set the value of keep.xy, not the value of k in the call to the
function moment(). You can see this in the value of ans$call:
> ans$call
integrate(f = moment, lower = - Inf, upper = -1, keep.xy = 2,
theta = 0.1)
For this reason, I always use unusual argument names for functions that
I know I will be using with integrate(), nlmin(), nlminb(), etc. For
example, Nicole could rename the argument k to k.weird, and then set
k.weird=2 within the call to integrate:
> moment_function(x,theta,k.weird=1) x^k.weird*p.mix(x,theta)
> p.mix_function(x,theta) .5*dnorm(x,theta) + .5*dnorm(x,-theta)
> ans_integrate(moment,-Inf,-1,theta=0.1,k.weird=2)
> ans$call
integrate(f = moment, lower = - Inf, upper = -1, theta = 0.1,
k.weird = 2)
Sincerely,
--Steve M.
_____________
| *** | Steven P. Millard, Ph.D.
| * |
| * * * | P robability, TEL: 206-528-4877
| * * * | S tatistics & FAX: 206-528-4802
| * | I nformation E-mail: SMillard@ProbStatInfo.com
| * | Web: www.ProbStatInfo.com
| *** | 7723 44th Avenue NE
|___________| Seattle, WA 98115-5117 USA
-----Original Message-----
From: Nicole DePriest Demers [SMTP:ntd@stat.ohio-state.edu]
Sent: Tuesday, March 10, 1998 3:11 PM
To: s-news@wubios.wustl.edu
Subject: [S] Integration difficulties
I am trying to integrate a function and a few problems have arisen.
#Function to be integrated is
moment_function(x,theta,k=1) x^k*p.mix(x,theta)
#where p.mix is a mixture of two normal distributions
p.mix_function(x,theta) .5*dnorm(x,theta) + .5*dnorm(x,-theta)
Note what happens in the following cases:
1) > ans_integrate(moment,-Inf,-1,theta=0.1,k=2)
> ans$integral
[1] -0.2443884
> ans$call
integrate(f = moment, lower = - Inf, upper = -1, keep.xy = 2,
theta = 0.1)
#The integrate function ignored my k=2 argument in the call to
#integrate. The result is for k=1.
2) #I removed the default value for k
> moment_function(x,theta,k) x^k*p.mix(x,theta)
> ans_integrate(moment,-Inf,-1,theta=0.1,k=2)
Error in qsub(f, bound, sgn, 0, 1, aux = optargs(l..: singularity
encountered
Dumped
3) #I switched the order of k and theta.
> moment_function(x,k,theta=.1) x^k*p.mix(x,theta)
> ans_integrate(moment,-Inf,-1,k=2)
Error in qsub(f, bound, sgn, 0, 1, aux = optargs(l..: singularity
encountered
Dumped
4) #I supply defaults for k and theta and send a different value for
#theta when I call integrate.
> moment_function(x,k=2,theta=.1) x^k*p.mix(x,theta)
> ans_integrate(moment,-Inf,-1,theta=0.2)
> ans$call
integrate(f = moment, lower = - Inf, upper = -1, theta = 0.2)
#For some reason, I need a default value for k in my integrand function
but not for theta. When k is supplied as an argument to integrate it
is
ignored but theta is not. I want to vary both k and theta within the
integrate command. All suggestions/explanations will be appreciated.
Nicole Demers
____________________________________________________________________
__________
Nicole Demers
Statistics Department
The Ohio State University
1958 Neil Ave., #404
Columbus, OH 43210, USA
ntd@stat.ohio-state.edu
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To
unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
bootstrap.memory
This message is in response to Carl Galopin's query on memory problems
in bootstrapping, but the response may interested anyone doing
simulation in S-Plus.
We have tested a number of looping methods for bootstrapping, and
found one that performs well on both time and memory. I'll include a
version of it below. We're doing more testing, and working on a more
complete set of bootstrapping functions for inclusion in our next
version, 4.0, and looking at further improvements after that.
The basic idea of the looping procedures below is to do blocks of
bootstrap samples, within a for() loop. The block size involves
a tradeoff:
-- short block sizes require that .Random.seed be saved
many times. Because .Random.seed is slated for storage on
a permanent database, every call to any sampling function
uses 81 bytes of memory that are not reclaimed until all functions
and loops are done.
-- large block sizes require storage of a large matrix of indices.
The default block size below may need to be reduced if n is large,
e.g. if (n * block * 6 + 32 > options("object.size"))
The speed isn't very sensitive to the block size, except that
large n and large block require too much of memory, and
very small block sizes cause sample() to be called too often.
In addition to the basic idea, there are a number of subtle points
that affect memory use; these are noted below.
The function is written so that a user interrupt will save as many
blocks of bootstrap results as have been completed. It also saves
results after exits due to some memory problems.
bootstrap <- function(X, FUN, ..., B = 100, seed = 0, block = 50){
# demonstration program for nonparametric bootstrapping
# X is a matrix or data frame, rows are observations
# FUN(X, ...)
# This version of bootstrap assumes that FUN() returns a scalar
# B bootstrap replications
# seed is an integer between 0 and 1000
# block is the block size, number of bootstrap values computed simultaneously
set.seed(seed) # So results are reproducible
if(is.null(dim(X))) X <- as.matrix(X)
n <- nrow(X)
call.stat <- function(i, X, FUN, indices, ...)
FUN(X[indices[,i], ], ...)
# The call.stat() function will be called by lapply() to
# do the actual bootstrapping
nblocks <- ceiling(B/block) # number of blocks
replicates <- numeric(B) # Create space for results
indices <- matrix( integer(n*block), nrow=n)
temp <- 1:block
names(temp) <- rep("",block) # This prevents lapply from coercing temp
# to a list internally, which uses more memory.
B2 <- block # Number of bootstrap replications in a block
on.exit({
cat("Saving replications 1:",(i-1)*block," to .bootstrap.results\n")
assign(".bootstrap.results",replicates,where=1,immediate=T)
}) # In case function is interrupted
for(i in 1:nblocks){
if(i==nblocks & B %% block){ # The last block may have a different B2
B2 <- B %% block
temp <- temp[1:B2]
indices <- indices[,temp]
}
indices[] <- sample(1:n, n*B2, T) #Sample the indices
replicates[temp+block*(i-1)] <-
unlist(lapply(temp, call.stat, X, FUN, indices, ...))
NULL # A loop "returns a value"; this statement
# at the end of the loop reduces memory buildup
}
on.exit()
replicates
}
>Bootstrap and Memory Management Gurus (both Windows and UNIX):
>I would appreciate your help with the following problem.
>
>I'm running S-Plus 3.3 using Windows 95 on a Pentium 133Mhz with
>32MB RAM and 2GB free hard drive space solely dedicated to S-Plus.
>I can run under a Windows 3.1 machine, if needed.
>
>In running Tibshirani's bcanon bootstrap program with a dataset of
>1000 observations for 5,000 iterations, I run into the error
>"Unable to obtain requested dynamic memory...."
>
>In viewing past discussions of the bootstrap, it appears that previous
>versions of S-Plus had a problem of freeing memory while in a loop.
Different methods of doing a loop use significantly different
amounts of memory. Typically lapply() is better than for(),
calling lapply() directly rather than calling apply() is better,
for() loops use less memory if most of the work is done inside
a function and if the last statement is a null. I.e. replace this
for(i in 1:B){
bunch of lines of code
results[i] _ final.result
}
with this
f _ function() { bunch of lines of code; final.result }
for(i in 1:B){
results[i] <- f()
NULL
}
The "Programmer's Manual" discusses memory use in loops in more detail.
You may find it useful to call the allocated() function inside your
loop, either to compare looping methods, or to monitor memory use to
stop gracefully if memory overload appears imminent.
>1. My question to S-Plus users on UNIX workstations are as follows:
> a. How well does S-Plus free up memory while in a loop?
> b. What type of computer configuration are you running
> (Machine, OS, RAM, Hard Drive)?
>
>2. My question to S-Plus users on Windows (either Windows 95 or 3.1):
> a. S-Plus has the following settings for memory:
> Compact, Free, Memory, Object Size.
> If I don't care how long the computer runs as long it doesn't
> dump the program, what's the best configuration of the memory
> settings?
> b. If I run bootstraps with datasets of 1,000 to 5,000 observations
> and iterations between 5,000 and 10,000,
> will I have to write intermediate results to a file in order to
> avoid getting a dump?
The bootstrap function above saves results automatically if you
interrupt the function, or if it fails due to trying to allocate
too large a dataset for indices. I can't guarantee that it will save
results in all crashes.
>3. My question to all S-Plus boostrap users are as follows:
> a. With datasets of 1000, 5000, 10000 observations, how practical
> is it running the bootstrap in S-Plus?
> b. What are your experiences with bootstrapping large datasets?
The above function ran for 10000 observations, except that I had to
reduce the block size.
If you hear your hard disk working a lot then you know you are using
more memory than you have RAM, and things will slow down a lot.
In that case you should try to reduce the block size.
>Thanks in advance,
>Carl Galopin
You're welcome.
==========================================================
| Tim Hesterberg MathSoft/Statistical Sciences |
| timh@statsci.com 1700 Westlake Ave. N, Suite 500 |
| (206)283-8802x319 Seattle, WA 98109-3044, U.S.A. |
| (206)283-0347 (fax) |
==========================================================
Conditional.logistic
From s-sender@utstat.toronto.edu Mon Jul 5 11:06:40 1993
Return-Path:
Received: from utstat.toronto.edu by biostat.mc.duke.edu (4.1/SMI-4.1)
id AA15004; Mon, 5 Jul 93 11:06:39 EDT
Received: by utstat; Mon Jul 5 10:45 EDT 1993
Received: by hsph.harvard.edu (4.1/1.34)
id AA06702; Mon, 5 Jul 93 10:45:37 EDT
Date: Mon, 5 Jul 93 10:45:37 EDT
From: rossini@hsph.harvard.edu (Anthony Rossini)
Message-Id: <9307051445.AA06702@hsph.harvard.edu>
To: s-news@utstat.toronto.edu
Subject: Conditional Logistic Regression (summary)
Status: R
Well, here is the list of what's being done and what has been done about
implementing conditional logistic regression in Splus.
First method is to take advantage of the links between partial likelihood
estimation for Cox's proportional hazards model and conditional likelihood
estimation. For information on this, you might check:
1 - Breslow and Day's nice IARC monograph on case-control studies
(vol 1 of the statistical methods in cancer research series)
2 - the SAS manual for proc phglm (proprotional hazards model proc)
(There is a discussion in there about how to use the code to to
condiditional logistic regression) (this information by Phil
Smith, pjs1@apus.ddt.cdc.gov ) (also can find the information
in the docs for proc phreg, the old version, which is the only
version I can find here at my school)
3 - One of the initial articles is by Prentice and Breslow,
Biometrika 1978, v65, pp153--8. It's very interesting to see
how perspectives changed in 15 years!
Basically, you want a version that does exact adjustment for ties.
The suggestions for implementation include coxreg, agreg, and whatever comes
with survival3 from statlib (coxph?). coxph will definitely do it, but the
others might be using asymptotic methods, which might not be the proper
thing to do. See my reference (about 3 paragraphs further on) about a paper
by Gail, Lubin, and Rubinstein.
Second, a student of Rob Tibshirani's is expecting to complete a version
implementing it, including additive models as well, sometime soon. This
sounds like a good addition to the current software, though it might be
"fakable" in a backwards way using the above methods and the present
assortment of GAM/smoothing stuff (warning -- last remark was a conjecture
based on locking myself up in a room for 3 days with lots of literature on
GAM's... there are experts, including Dr. Tibshirani, who probably know
better!)
Third, one could implement it using another program as a driver (such as
AS196 : Conditional multivariate logistic analysis of stratified
case-control studies. (33/1984), in the apstat algorithms section of applied
statistics (JRSS-C).), or pipe through stata, sas, etc, using some sort of
awful kludge, with shell scripts, filters, etc... The suggestion for AS196
also came with the hint that the suggestor was unable to make it work
properly.
Fourth, someone suggested using multinomial logit regression, and faking
by using log-linear models. I think this is a similar but different
problem, though I'm not sure (and will probably look it up sometime before I
write up my final results). Come to think of it, this is probably wrong, as
I'm pretty sure that only the matched pairs design can be run in GLIM, and
one of the references given was to a GLIM manual (the nice one, Statistical
Modelling in GLIM by Aitkin et al. , relevant info in (sec 5.4)). But you
can decide for yourself.
The crucial step here is deciding on the amount of approximation to do. A
nice discussion which will help in evaluating the methods from a theoretical
standpoint is in Gail, Lubin and Rubinstein, Biometrika 1981 vol68,
pp703--7. Some of the procedures do more approximation than others.
Anyhow, for posterity, I've included relevant replies as a p.s. and
apologize in advance for any misinterpretations that I put into the above
summary. Thanks!
-tony
--
Anthony Rossini - grad student/statistician/hacker
rossini@biostat.harvard.edu
Department of Biostatistics, Harvard School of Public Health
677 Huntington Ave, Boston MA 02115 617-432-1056
===============================================
From: tibs@utstat.toronto.edu
Date: Thu, 1 Jul 93 12:28 EDT
To: rossini@hsph.harvard.edu
By coincidence, I have a grad student who is implementing the matched
case control model right now, including the additive version.
he should be done in a month. Please let me know if other versions
exist.
Rob Tibshirani
===============================================
Date: Thu, 1 Jul 93 12:44:40 EDT
From: pjs1@apus.ddt.cdc.gov (Phil)
To: rossini@hsph.harvard.edu
Subject: Re: Conditional logistic regression in SHi,
Hi, Toni:
The likelihood for conditional logistic regression and Cox's partial
(nonparametric) likelihood for survival data are identical.
This means that you can use the survival functions (coxreg) in S to do the
computations you need for matched case control studies.
Check the SAS manual for proc phglm (proprotional hazards model proc).
There is a discussion in there about how to use the code to to
condiditional logistic regression.
Phil Smith
Centers for Disease Control
================================================
Return-Path:
From: stvjc@gauss.med.harvard.edu (Vincent Carey 2-3265)
Subject: Re: Conditional logistic regression in S
To: rossini@hsph.harvard.edu (Anthony Rossini)
Date: Thu, 1 Jul 93 13:19:55 EDT
In-Reply-To: <9307011626.AA20708@hsph.harvard.edu>; from "Anthony Rossini" at Jul 1, 93 12:26 pm
X-Mailer: ELM [version 2.3 PL11]
tony, i would think it can be done using agreg. i know
people used to "trick" SAS PHGLM (the old proportional
hazards regression proc of the SAS SUGI lib) into
doing CLR, but i don't recall the details. i have notes
somewhere if you need more help.
also, i seem to remember an ApStat algorithm that is
on statlib that does the CLR, but i was never able to
get it to run.
here it is:
196 Conditional multivariate logistic analysis of stratified case-control
studies. (33/1984)
an S interface to that would probably be very easy to concoct,
but my efforts in that direction, some time ago, ran aground.
--
stvjc@gauss.med.harvard.edu
====================================================
Date: Thu, 1 Jul 1993 11:16:47 -0700
From: David Eaves
To: rossini@hsph.harvard.edu
Subject: conditional logistic regression in S
Is that the same as multinomial logit regression? If so,
yes, it can be "faked", by conditioning a log-linear model
on the total count across the outcome categories. There is a good
description of this in Statistical Modelling in GLIM by
Aitkin et al. (sec 5.4) However it becomes cumbersome with more than a
very small number of continuous explanatory variables, in which
case I use BMDPLR.
===================================================
Date: Thu, 1 Jul 93 13:31:24 PDT
From: cberry@tajo.UCSD.EDU (Charles C. Berry)
To: rossini@hsph.harvard.edu
Subject: Re: Conditional logistic regression in S
I think all the necessary machinery is there, but is only practical for
relatively small numbers of events. (I think 150 is pushing it, but this depends on memory size).
By using the Poisson family with a bit of restructuring of the data and inclusion of suitable terms for incidental parameters (which do not appear in the conditional logistic), you can perform a proper likelihood-based conditional logistic regression.
See McCullagh and Nelder's GLIM book, Section 6.4 and/or Palmgren (1981) Biometrika 68,563-6, for a more lucid rendition than I give here.
BTW, this approach is nice if you have continous varaibles for which the functional form is up-for-discussion. You can use all of the diagnostic machinery in Stat Models in S as well as smoothing splines.
____________________________________________________________________________
| | |
|Charles C. Berry | Department of Family and Preventive Medicine |
| | University of California, San Diego |
|ccberry@ucsd.edu | La Jolla, San Diego 92093-0622 |
|cberry@tajo.ucsd.edu | (619) 534-2098 (voice) (619) 534-4642 (FAX) |
|_______________________|____________________________________________________|
======================================================
Date: Thu, 1 Jul 93 16:36:28 -0500
From: therneau@mayo.edu (Terry Therneau)
To: rossini@hsph.harvard.edu
Subject: Conditional logistic
This can be done using coxph and the "exact partial likelihood". You will
need survival3 from statlib (S-plus has survival2, but will upgrade to version
3 in the fall). If n and m are large it won't be particularly fast, but I
intend to improve it.
Details: For a cox regression model with
time=1 (or any other constant)
status= 1:case, 0:control
strata= group (1 stratum per case/control group)
x = covariates
the exact partial likelihood proposed by Cox is exactly the same as the
conditional logistic likelihood. Version 2 of the survival package only
impliments the Breslow approximation, however.
Terry Therneau
Mayo Clinic
====================================================
DataInput
Last week I posted a notice about reading in some
non-standard format data. There are lots of helpful
suggestions - basically preprocess the data in UNIX
before reading it into Splus - and so here is a
full summary of the responses, in order of arrival.
A big thank you to everyone who replied.
Best Regards, Anwar
>>>> Original request
Whats the best/easiest way to read data into Splus
where are in the following format. The files
- are all of type ascii
- there are 20 fields (or columns) in each file
- most, but not all, fields are surrounded by double quotes of the
type "X"
- each field is separated by a comma
- some fields are sometimes blank, ie. "", just a pair of double quotes
A typical file (4 lines shown) might look like,
"SH","199409","199504","CANCEL","25","102","01","HU6","199409",000015000,
000733,"",000000000,000000,"",000000000,000000,"",000000000,000000
"SH","199409","199504","CANCEL","12","087","01","SA1","199409",000034000,
001141,"",000000000,000000,"",000000000,000000,"",000000000,000000
"SH","199409","","","03","063","02","UB4","199610",000130000,003553,"199509",
000125300,003438,"199409",000119000,003500,"",000000000,000000
"SH","199408","","","03","042","01","SW6","199604",000094800,000302,"199504",
000091800,000293,"199408",000090000,000287,"",000000000,000000
In Splus fields with text (e.g. field 1 = "SH", field 4 = "CANCEL",
field 8 = "HU6") should remain text. However, some fields in double
quotes (e.g. field 2 = "199409", field 3 = 199504", field 5 = "25",
field 6 = "102", field 7 = "01", field 9 = "199409") should in fact
be (converted to) numeric, so that I can analyse (mean, variance, etc.)
them. Fields which appear as only double quotes (e.g. field 12 = "",
field 15 = "", field 18 = "") should be treated as numeric but I would
like to have the option of setting it to a default value (numeric (e.g. 0)
or text (e.g. NA, X)) when only a double quote occurs; when a value is
present there is obviously no problem.
I expect I'll label the columns (e.g. F1, F2, ...) before reading them in.
Whats the best approach??
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 1 : Z. Todd Taylor
You don't say what OS you're on, but if it's Unix, the natural
thing to do is preprocess your file with 'awk' or 'sed' or
'perl' or something like that and then slurp the result into S
with scan() or read.table().
For your example above, the procedure would look something like
this ($ is the Unix prompt, S> is the S prompt):
$ # Eliminate quotes...
$ sed -e 's/"//g' original.file > new.file
S> # Read it in...
S> data <- scan("new.file", sep=",",
what=list(field1="", field2=0, )
)
S> # --- OR ---
S> data <- read.table("new.file", sep=",",
col.names=c("field1", "field2", )
)
(The above will break if there are any commas within the quoted
strings.)
I do that kind of thing so often I wrote modified versions of scan
and read.table that will read directly from a Unix pipe, so I can
do things like:
S> data <- scan('unix_cmd |', ...)
where unix_cmd is a call to sed, perl, awk or whatever to do the
preprocessing. Email me if you're interested.
If you're not running under Unix, your options are more limited.
- --Todd
- --
Z. Todd Taylor
Pacific Northwest National Laboratory
zt_taylor@pnl.gov
>>>>>>>>>>>>>
Below are my modified versions of scan and read.table (which
I've named Scan and Read.Table, respectively). I've also
included standard S help files for these. You'll have to break
these up into separate files, source the *.sfun files into
your .Data directory (or another library if you prefer) and
copy the .d files into .Data/.Help *without* the .d on the
filenames.
i.e., break this message into four files:
Scan.sfun Scan.d
Read.Table.sfun Read.Table.d
Then do something like ($ is the Unix prompt):
$ Splus < Scan.sfun
$ cp Scan.d .Data/.Help/Scan
$ Splus < Read.Table.sfun
$ cp Read.Table.d .Data/.Help/Read.Table
I hope you find these as useful as I do.
- --Todd
P.S. Note that the Scan function is a little tricky in that it
makes use of a named pipe (FIFO) to avoid creating a
temporary file. This speeds it up considerably--about a
factor of 2 or 3 for very large files. However, if FIFOs
don't work on your Unix you might have to modify the Scan
function a bit. I doubt you'll have any problem, though a
few (older) Unixes didn't do FIFOs quite right on
NFS-mounted filesystems.
- --
Z. Todd Taylor
Pacific Northwest National Laboratory
zt_taylor@pnl.gov
- ------------------------ Scan.sfun --------------------------
# Function to replicate 'scan', but to allow specification of
# a unix command in lieu of a filename.
# $Id: Scan.sfun,v 0.1 1994/06/09 21:06:35 d3a061 Exp $
# $Log: Scan.sfun,v $
# Revision 0.1 1994/06/09 21:06:35 d3a061
# Initial revision.
#
"Scan" <- function(file, ...) {
# If an ordinary filename was given, we're done...
len <- nchar(file)
if ( substring(file, len, len) != "|" ) return(scan(file=file, ...))
# Otherwise, get busy...
# Strip off the "|"...
cmd <- substring(file, 1, len-1)
# Set up a named pipe...
fifo <- tempfile("Scan")
on.exit(unlink(fifo))
unix.shell(paste("mknod", fifo, "p"), output.to.S=F, shell="/bin/sh")
# Fire off the Unix cmd...
unix.shell(paste(cmd, ">", fifo, "&"), output.to.S=F, shell="/bin/sh")
# Let scan() do the rest...
scan(file=fifo, ...)
}
- ------------------------ Scan.d --------------------------
.BG
.FN Scan
.TL
Input data from a file or a Unix pipe.
.DN
The function is just like scan() except the "file" argument can
specify either a filename or a Unix command. If the latter is
given, it scans from the command's standard output.
.CS
Scan(file, ...)
.RA
.AG file
A character string describing either a filename (see documentation
for the scan function) or a Unix command. If it is a Unix command,
its last character MUST be a pipeline symbol (|).
.OA
.AG ...
Other arguments to be passed on to the the scan function.
.RT
The data from the file (or Unix pipe) as read by the scan function.
.SE
If 'file' describes a Unix command, temporarily creates a FIFO
(named pipe) and spawns a shell to execute
the command and write its output into the FIFO. The FIFO is removed
upon completion of the Scan call.
.DT
Because it relies on Unix FIFO files (named pipes), this function
might not work on Unix systems that don't support FIFOs. Also, some
Unices don't support FIFOs on NFS-mounted filesystems. (SunOS supports
both.)
.SH REFERENCES
Written by Todd Taylor, Pacific Northwest National Laboratory.
.SA
scan
.EX
# Read data from the output of a program...
data <- Scan("myprog |",
what=list(rowname="", var1=0, var2=0)
)
.WR
- ------------------------ Read.Table.sfun --------------------------
# $Log: Read.Table.sfun,v $
# Revision 1.1 1995/02/27 18:15:42 d3a061
# Initial revision
#
#
"Read.Table" <- function(file, ...) {
Version <- "$Id: Read.Table.sfun,v 1.1 1995/02/27 18:15:42 d3a061 Exp $"
# If this is just a plain file, we're done...
len <- nchar(file)
if(substring(file, len, len) != "|") return(read.table(file = file, ...))
# Otherwise, get busy...
# Strip off the "|"...
cmd <- substring(file, 1, len - 1)
# Set up a temporary file. Can't use a named pipe because read.table
# scans the thing more than once...
tfile <- tempfile("Read.Table")
on.exit(unlink(tfile))
if ( unix.shell(paste(cmd, ">", tfile), output.to.S=F, shell="/bin/sh") ) {
stop("Creation of temporary file failed")
}
# Let read.table() do the rest...
read.table(file = tfile, ...)
}
- ------------------------ Read.Table.d --------------------------
.BG
.FN Read.Table
.TL
Input data from a file or a Unix pipe.
.DN
The function is just like read.table() except the "file" argument can
specify either a filename or a Unix command. If the latter is
given, it reads from the command's standard output.
.CS
Read.Table(file, ...)
.RA
.AG file
A character string describing either a filename (see documentation
for the read.table function) or a Unix command. If it is a Unix command,
its last character MUST be a pipeline symbol (|).
.OA
.AG ...
Other arguments to be passed on to the the read.table function.
.RT
The data from the file (or Unix pipe) as read by the read.table function.
.SE
If 'file' describes a Unix command, creates a temporary file
and spawns a shell to execute
the command and write its output into the temporary file. The temporary
file is removed upon completion of the Read.Table call.
.SH REFERENCES
Written by Todd Taylor, Pacific Northwest National Laboratory.
.SA
read.table, Scan
.EX
# Read data from the output of a program...
data <- Read.Table("myprog |", header=T)
.WR
>>>>>>>>>>>>>
> Hi Todd,
>
> Nearly there. After getting rid of the quotes using
>
> csh> sed -e 's/"//g' file1.txt > file1
>
> I replaced all the periods (.) with NA using
>
> csh> sed -e 's/\./NA/g' file1 > file2
>
> Then I loaded it into Splus using
>
> > a<-read.table("file2",sep=",")
>
> This worked fine, except that when I did a summary
> some missing variables showed up as empty fields. For
> example
>
> F4
> :5700
> CANCEL:1167
> LAPSED:1038
> NTU : 776
> REPLAC: 79
>
> The 5700 is related to missing values. I tried to use the
> is.na function to replace with something like MT (missing
> text) using
>
> > a$F4[is.na(a$F4)[<- MT
That looks like the summary for a factor object. The blank
field is not representing NA's, but rather the null string "",
which is what came in from the empty fields between commas in
your file (,,). In making the factor object, the read.table()
function considered only the string "NA" to be a missing value,
assuming the empty string "" was one of the legal values in the
column. You can transform those empty values into true missing
values in your data column with this non-obvious incantation:
a$F4 <- factor(as.character(a$F4), exclude="")
[ snip ]
> By the way your previous answers and codes are too good for
> just me. Would you mind if I post them onto the news group?
It's considered normal practice to summarize private email
responses to the mailing list, and that's fine with me. My
Read.Table and Scan functions were posted to the list in August,
so it may not be necessary to post them again unless someone
asks.
- --Todd
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 2 : Rod Tjoelker
Hello,
You asked about how to read in ASCII data where some of the
fields have values in quotes, and others do not and fields
are separated by commas.
>> Whats the best/easiest way to read data into Splus
>> where are in the following format. The files
>> "SH","199409","199504","CANCEL","25","102","01","HU6","199409",000015000,
>> 000733,"",000000000,000000,"",000000000,000000,"",000000000,000000
>> "SH","199409","199504","CANCEL","12","087","01","SA1","199409",000034000,
>> 001141,"",000000000,000000,"",000000000,000000,"",000000000,000000
>> "SH","199409","","","03","063","02","UB4","199610",000130000,003553,"199509",
>> 000125300,003438,"199409",000119000,003500,"",000000000,000000
>> "SH","199408","","","03","042","01","SW6","199604",000094800,000302,"199504",
>> 000091800,000293,"199408",000090000,000287,"",000000000,000000
I've got a little function that will strip quotes from
character strings inside Splus. It will also convert the
vector to character, numeric, integer, or factor as specified
by the user. Specify ouput = "c", "n", "i", or "f" respectively.
strip.quote <- function(x, output = "c") {
out <- x
if(is.character(x)) {
nc <- nchar(x)
out <- substring(x, 2, nc - 1)
}
out <- switch(output,
c = as.character(out),
n = as.numeric(out),
f = factor(out),
i = as.integer(out))
return(out)
}
Here's one way to do this:
input <- read.table("file",as.is=T,row.names=NULL, sep=",")
for( i in c(1,4,8)) input[[i]] <- strip.quote(input[[i]],"c")
for( i in c(2,3,5,6,7,9,12,15,18)) input[[i]] <- strip.quote(input[[i]],"i")
Hope this helps,
==================
Rod Tjoelker
Applied Statistics
The Boeing Company
==================
tjoelker@redwood.rt.cs.boeing.com
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 3 :Don MacQueen
Try this
-outside splus (on a copy of a file) with a text editor.
1) convert commas to tabs
2) remove all quotes
-inside splus-
3) read.table('filename',sep='\t',head=F)
If you have lots of files, and are on a unix platform, then sed can make it
easy to edit all of the files. On DOS or Windows, I don't know.
Oops, simpler still, just get rid of all the quotes, then use
read.table('filename',sep=',')
Check your documentation for more details on the read.table() function, for
instance, how to assign field names. If the files are very, very large, you
might have to use scan() instead of read.table().
- -Don
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 4 : Clint Bowman
I'd be tempted to pass the file (especially if you have many files) by perl and
rewrite the data and probably add the headers at the same time. If you have
just the one file and it's not too long, you could make all the changes in any
text editor.
- --
Clint Bowman INTERNET: clint@ecy.wa.gov
Air Quality Modeler INTERNET: clint@math.utah.edu
Department of Ecology VOICE: (360) 407-6815
PO Box 47600 FAX: (360) 407-6802
Olympia, WA 98504-7600
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 5 : Don McKenzie
You can go into a text editor and globally replace the
"" with "", then eliminate all the quotes.
Now scan the file into a list object, for which you've
specified some elements as numeric and some as chars. See
the blue book or V&R for how.
dmck
dmck@silvae.cfr.washington.edu
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 6 : Andrea Spano
Stupid way but it works....
Open your file in MS word
Convert your file into a table Comma as a separator)
remove all the " by a simple replace
import your file into MS excell
save it as an ASCII file
then use read.table
Believe me, Splus is a real pain when reading data.
Not very "Scientific" but it works and it is quick
Regards
Andrea
andrea@spano.demon.co.uk
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 7 : Nick Ellis
Anwar,
if you're on Unix I'd strip out all the double quotes thus
% sed -e 's/"//g' infile > outfile
then I'd read them in using scan with the "what" argument to tell
whether it's character ot numeric
>
res<-data.frame(scan("c:/outfile",what=list("",0,0,"",0,0,0,"",0,0,0,0,0,0,0
,0,0,0,0,0),sep=","))
> names(res)<-paste("F",1:20,sep="")
> res
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 F13
F14 F15 F16 F17 F18 F19 F20
1 SH 199409 199504 CANCEL 25 102 1 HU6 199409 15000 733 NA 0
0 NA 0 0 NA 0 0
2 SH 199409 199504 CANCEL 12 87 1 SA1 199409 34000 1141 NA 0
0 NA 0 0 NA 0 0
3 SH 199409 NA 3 63 2 UB4 199610 130000 3553 199509 125300
3438 199409 119000 3500 NA 0 0
4 SH 199408 NA 3 42 1 SW6 199604 94800 302 199504 91800
293 199408 90000 287 NA 0 0
To set defaults for missings you do something like this
> res$F3[is.na(res$F3)] <- 0
Good Luck
Nick
Nick Ellis
CSIRO Marine Research
PO Box 120 email Nick.Ellis@marine.csiro.au
Cleveland QLD 4163 ph +61 (07) 3826 7260
Australia fax +61 (07) 3826 7222
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 8 : Zhaorong Jiao
I am not sure if you have known a little bit of vi, a text editor in
unix. If you do then the research and replace command can easily
delete the double quotes in the files and once the double quotes are
deleted you can use the usual splus functions to import your data sets
into Splus . You do not need to learn vi in order to delete the
double quotes, say the first of your data files is called first.dat
then you follow these instructions:
(1) vi first.dat
(2) hold down the 'shift' key and the key comma, and that will
get you to the bottom of the file and you will see a line start with
a comma, then
(3) type in after the comma the following codes g/"/s// /g, note there is
a empty space between s// and /g.
(4) once you finish typing the codes you push the enter key
all the double quotes will be disappeared.
(5) to save the file you first push the "Esc" key once or twic and
hold down the "Shift" and ":/;"key the same time and that will
get you to the end of the file again and you will a line begin will a ':'
(6) once you see the ":" type x! and the file is save and you exist vi.
Using the vi editor is a very quicker way to delete the quotes
particularly if you have a very big file. Some people have also
written splus functions to delete quotes within splus. I guess it
will take much longer time to get the data set into splus and delete
the quotes the same time. If you like to have functions for delete the
quotes I will be happy to forward those functions to you.
Zhaorong Jiao
Zhaorong.Jiao@anu.edu.au
>>>>>>>>>>>>>
>>>>>>>>>>>>> RESPONSE 9 : Alejandro Munoz Del Rio
Hi Mohammed,
By now you may have received several replies to your query on the S+
list. Here is my approach. Maybe clever use of scan() may read your
input file and then you can clean it inside S+, but I'd rather do the
filtering / formatting outside S+ since I think working in the shell
is faster and more reliable. (S+ for Windows 3.1 can be clunky.)
You don't mention what platform you're in; I use sed and awk which
are available for unix and have been ported to MSDOS. The steps are:
1 - use sed to remove quotation marks, commas, extra zeros,
2 - use awk to have one line per case (all 20 variables in one line),
3 - insert column headings (variable names),
4 - read in S+
1. Here is my sed script (I'll assume it's in a file called "sedscript"):
- ----- begin sedscript
# replace commas by spaces
s/,/ /g
# replace double quotes by missing data indicator, NA
s/\"\"/NA/g
# replace two or more zeros at end of record by one 0
s/ [0]*$/ 0/g
# replace multiple zeros by single zero
s/ 0\{2,\} / 0 /g
# (needed twice to take care of consecutive cases)
s/ 0\{2,\} / 0 /g
# replace multiple zeros at beginning of line by one 0
s/^0\{2,\}//g
# delete leading zeros in nonzero entries
s/ 0\{2,\}/ /g
# strip double quotes
s/\"//g
- ----- end sedscript
2. Run the sed script and collapse cases to one record (line) per case
$ sed -f sedscript yourdata | awk '{ a=$0 ; getline ; print a $0 }' > data2
3. Add column names:
$ (echo var1label var2label ... var20label ; cat data2) > data3
4. Read into S object:
> yourframe <- read.table("data3", header=T)
Hope this helps.
alejandro
munoz@stat.wisc.edu
>>>>>>>>>>>>>
deparse
Your message could not be delivered to
'feh3k@avery.med.virginia.edu (host: avery.med.virginia.edu) (queue: servers)' for the following
reason: ' collect: Cannot write tfXAA55386'
Your message follows:
Received: from biostat.mc.duke.edu by mail.virginia.edu id aa18810;
9 May 97 21:27 EDT
Received: from utstat.toronto.edu (utstat.toronto.edu [128.100.73.1]) by biostat.mc.duke.edu (8.6.10/8.6.6) with SMTP id VAA08806 for ; Fri, 9 May 1997 21:27:57 -0400
Received: by utstat; Fri May 9 20:37 EDT 1997
Received: from localhost by grover.stat.washington.edu
(5.65/UW-NDC Revision: 2.23 ) id AA08116; Fri, 9 May 97 17:37:24 -0700
Date: Fri, 9 May 1997 17:37:23 -0700 (PDT)
From: "Jacob A. WEGELIN"
To: S-News MailingList
Subject: Re: Can't pass data of mode "call" ... summary of responses
In-Reply-To:
Message-Id:
Mime-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII
Thanks to all of you who responded to my post. Below is a summary
of some of the answers, followed by a copy of my original post.
###############################################################
Brian Ripley's response entire:
---------------------------------------------------------------
Three comments:
(1) Using match.call rather than sys.call will give a more meaningful
record in many cases.
(2) A character version of a call seems tricky: I too have used
dput. But this seems to work:
call2char <- function(x)
{
args <- x[-1]
argnames <- names(args)
havenames <- argnames != ""
args[havenames] <- paste(argnames[havenames], args[havenames],
sep = " = ")
paste(x[[1]], "(", paste(unlist(args), collapse=", "), ")", sep="")
}
(3) In actual.text I would use tempfile to get a name of the scratch file,
and unlink it afterwards.
###############################################################
Martin Maechler pointed out that one can strip unwanted attributes
from an object before passing it to one of the fussy old-S functions:
---------------------------------------------------------------
When having to use an "old-S" function, I use
oldSfunction(unattr(myobj))
instead of
oldSfunction(myobj)
unattr <- function(obj, leave.attrib = NULL)
{
## Purpose: Remove non-standard attributes from 'obj'
## To work with old-S ( .S and .Internal ..)
## Author: Martin Maechler ETH Zurich, Date: 17 Jul 92
## ----------------------------------------------------------------
leave.attrib <- c(leave.attrib, "dim","dimnames","names")
a <- attributes(obj)
attributes(obj) <- a[!is.na(match( names(a), leave.attrib))]
obj
}
###############################################################
Several people pointed out that a character string can be
turned back into an expression and then evaluated, as follows:
> eval(parse(text="print(1:5)"))
[1] 1 2 3 4 5
On the other hand, trying to use deparse to go in the other
direction doesn't work quite right when the function call is long:
> junkfn
function(...)
deparse(sys.call())
>junkfn(dog,cat, mouse, rat,
hyena,fox,wolf=itsaveryniceday,coyote=wellIamhowlingloudly,snake=IliketoeatmiceifIamallowedtocrushthemfirst)
[1] "junkfn(dog, cat, mouse, rat, hyena, fox, wolf = itsaveryniceday, coyote = "
[2] "\twellIamhowlingloudly, snake = "
[3] "\tIliketoeatmiceifIamallowedtocrushthemfirst)"
Similarly, using deparse(as.name()) only works if the call is fairly
short. The function as.name() truncates its argument, so that the following
approach does not for for complicated function calls:
> junkfn
function(...)
deparse(as.name(sys.call()))
> junkfn(dog,cat, mouse, rat, hyena, fox,wolf=itsaveryniceday,coyote=wellIamhowl>
[1] "junkfn(dog, cat, mouse, rat, hyena, fox, wolf = itsaveryniceday, coyote = "
###############################################################
S.D. BYERS held out a ray of hope for the future:
---------------------------------------------------------------
You may want to read the technical report on S 4.0 that is available on
the web. This may not be immediately applicable now but in this version
the object oriented system will be fully exploited, for instance the help
file for each function object will be part of the object itself, not some
separate .d file. This seems to be very much the sort of thing that youi
are trying to achieve and there may be suggestions in there you would
find useful.
I think you can get it off the web from Chambers home page at AT&T.
###############################################################
The original post:
---------------------------------------------------------------
>From wegelin@stat.washington.edu Fri May 9 17:33:34 1997
Date: Wed, 7 May 1997 16:08:53 -0700 (PDT)
From: Jake WEGELIN
To: S-News MailingList
Subject: Can't pass data of mode "call" to old-S function: "call" vs "character"
Often, when returning a value from a function, I include the call to the
function as an attribute. This is a way to keep track of where an object
came from, how it was created. Later, one can pass this attribute to the
"eval" function and reproduce the original results.
This causes problems, however. Certain otherwise very useful S functions
vomit when you try to pass an object to them that has an attribute of mode
"call".
I have one rather kludgey solution, but am curious whether others have better
ideas. My solution replaces the attribute of mode "call" with one of mode
"character". One disadvantage of this is that making this replacement is a
messy operation; another is that you can no longer pass the attribute to the
"eval" function.
QUESTIONS.
Has anyone else been tempted to include a "call" attribute in the return
value of a function? Other than this awkward behavior of old-S functions, is
there some reason not to do it? Is there a better way to keep track of how
an object was created?
Does anyone know of a better way to change an object of mode "call" to
one of mode "character" than the function listed below?
DETAILS.
For instance,
> myfn
function(x, y)
{
# do a bunch of stuff
answer <- c(x, y, x * y) # or whatever
attr(answer, "myinfo") <- sys.call()
answer
}
> junk_myfn(1,44)
> # much time passes; I forget where junk came from ...
> junk
[1] 1 44 44
attr(, "myinfo"):
myfn(1, 44)
I can reproduce the results:
> eval(attr(junk,"myinfo"))
[1] 1 44 44
attr(, "myinfo"):
myfn(1, 44)
But:
> dist(junk)
Error in dist(junk): Can't pass data of mode "call" to old-S function
Dumped
Remove the attribute of mode "call", and it works fine:
> attr(junk,"myinfo")_NULL
> junk
[1] 1 44 44
> dist(junk)
[1] 43 43 0
attr(, "Size"):
[1] 3
To change the object of mode "call" to an object of mode "character", I
use this awkward function:
> actual.text
function(thename)
{
# WARNING!! This overwrites any file named "deleteMe" in your working
# directory.
dput(thename, "deleteMe")
paste(scan("deleteMe", what = " ", sep = "\n"), collapse = "")
}
For instance,
> myfn
function(x, y)
{
# do a bunch of stuff
answer <- c(x, y, x * y) # or whatever
attr(answer, "myinfo") <- actual.text(sys.call())
answer
}
> junk_myfn(1,44)
> junk
[1] 1 44 44
attr(, "myinfo"):
[1] "myfn(1, 44)"
> dist(junk)
[1] 43 43 0
attr(, "Size"):
[1] 3
ecdf
From:
4:20 PM
Subject: RE: [S] empirical
CDFTo: "'Jane_2_Xu@sbphrd.com'" , s-news@wubios.wustl.edu
Try this instead. It uses the actual data points, not some interpolation.
Plus it is simpler.
plot(sort(x),ppoints(x),type="l")
Sam A. Martin
Statistician
Innovative Emergency Management
http://www.ieminc.com
8555 United Plaza Blvd. Suite 100 Baton Rouge, La 70809
(225) 952-8230 [voice]
(225) 952-8122 [fax]
> -----Original Message-----
> From: Martin, Sam [SMTP:sam.martin@ieminc.com]
> Sent: Wednesday, April 19, 2000 1:33 PM
> To: 'Jane_2_Xu@sbphrd.com'; s-news@wubios.wustl.edu
> Subject: RE: [S] empirical CDF
>
>
> Try
>
> x<- rnorm(1000)
> plot(quantile(x,seq(.01,1,,length(x))),seq(.01,1,,length(x)),type="l")
>
>
> Sam A. Martin
> Statistician
> Innovative Emergency Management
> http://www.ieminc.com
> 8555 United Plaza Blvd. Suite 100 Baton Rouge, La 70809
> (225) 952-8230 [voice]
> (225) 952-8122 [fax]
>
>
> > -----Original Message-----
> > From: Jane_2_Xu@sbphrd.com [SMTP:Jane_2_Xu@sbphrd.com]
> > Sent: Wednesday, April 19, 2000 1:19 PM
> > To: s-news@wubios.wustl.edu
> > Subject: [S] empirical CDF
> >
> >
> > Hi everyone,
> >
> > Has anyone written a function to plot the empirical cumulative
> > distribution
> > function? Splus has cdf.compare, but what I need is to just look at the
> > Empirical CDF.
> >
> > Any information would be greatly appreciated.
> >
> > Thanks,
> >
> > Jane
The empirical CDF should not be diagonal between observed points (unless
you want to interpolate a density), it should be staircase-shaped.
(Hence type="l" does not do quite the right thing.)
The following will do it for you:
delta_diff(range(x))*.05
n_length(x)
plot(c(min(x)-delta,sort(x),max(x)+delta),c(0,0:n/n),type="S")
the "delta" extension at each end is just for looks; of course in principal
the empirical CDF can be extended to plus and minus infinity. If you
plot the minimum and maximum, you can see that the steps start and end
in the right place
abline(v=range(x),lty=3)
(Now if you want to be really pedantic you can argue that the vertical
connecting lines and the last pixel of each horizontal should be deleted,
but that is too damn ugly!)
-A Z.
From:
Jane_2_Xu@sbphrd.com
9:29 AM
Subject:
[S] Summary of responses about empirical CDF
To:
s-news@wubios.wustl.edu
Hi All,
Below is a summary of the responses I received regarding plotting Empirical
CDF (and reference quantile lines). I presonally like Alan's ECDF line and
the way Gerrit Eichner plots the quantile lines, so I kinda miss and match
and
have my "own" function. Their codes are included below. There are a couple
of people pointed out that Hmisc library in Statlab has this function too.
Since
Alan's and Gerrit's codes are good enough for me, I didn't check out the
Hmisc
library. Thanks to many more people who also responded.
#################### Alan Zaslavsk's Codes #######################
delta_diff(range(x))*.05
n_length(x)
plot(c(min(x)-delta,sort(x),max(x)+delta),c(0,0:n/n),type="S")
abline(h=0:4/4,lty=3)
abline(v=c(approx(x=(1:length(x)-1)/length(x),y=sort(x),xout=c(.25,.50,.75),
f=0,method="constant")$y,max(x)),lty=3)
##################### Gerrit Eichner's Function ####################
my.ecdf <- function(x, probs=c(0.25, 0.5, 0.75), extfak=0.05, ...) {
x <- sort(x); n <- length(x)
xext <- (x[n] - x[1])*extfak; y <- (1:n)/n
plot(x, y, xlim=c(x[1] - xext, x[n] + xext), ylim=c(0, 1), ...)
segments(c(x[1] - xext, x), c(0, y), c(x, x[n] + xext), c(0, y))
abline(h=0:1, lty=2)
quants <- quantile(x, probs); nq <- length(quants)
segments(rep(par("usr")[1], nq), probs, quants, probs, lty=3)
segments(quants, probs, quants, rep(0, nq), lty=3)
names(quants) <- probs
invisible(quants)
}
From:
Bill Venables
Fri 7:54 PM
Subject:
Re: [S] Summary of responses about empirical CDF
To:
Jane_2_Xu@sbphrd.com, s-news@wubios.wustl.edu
At 09:29 AM 21/4/00 -0400, Jane_2_Xu@sbphrd.com wrote:
>
>Hi All,
>
>Below is a summary of the responses I received regarding plotting Empirical
>CDF (and reference quantile lines). I presonally like Alan's ECDF line and
>the way Gerrit Eichner plots the quantile lines, so I kinda miss and match
and
>have my "own" function. Their codes are included below. ...
Right now I am thinking about teaching again (very briefly!) so I thought
this
might make a nice little example of a generic function and a couple of
methods.
The following suite uses essentially the same ideas as in Alan Zaslavsky's
few
lines of code, but presents it in a slightly more flexible way. You can
call
the function with one sample, several samples (as separate arguments) or a
single
list of samples and it places the ECDFs on the same plot. The argument q
may be
used to request quantile lines, e.g. q=1:3/4 asks for the quartiles. In
addition
given any object inheriting from "lm" results in an ECDF plot of the
residuals.
(Note that statistically normal scores plots are better for looking at
residuals
ECDF plots, though.)
Note: this is essentially S version 3 OOP. Converting to S version 4
without
using backward compatibility would require a comprehensive overhaul.
As usual, there are few little tricky bits, Njoi. Bill Venables.
---(cut here)---
ecdf <- function(x, ...) {
if(is.null(class(x)))
class(x) <- data.class(x)
UseMethod("ecdf")
}
ecdf.default <- function(x, ..., xlab = "x", q = NULL, ylab = "Empirical
CDF") {
x <- range(unlist(m <- list(x, ...)), na.rm = T)
plot(x, 0:1, xlab = xlab, ylab = ylab, type = "n")
i <- 0
pu <- par("usr")
for(x in m) {
sx <- sort(as.numeric(x))
n <- length(x)
lines(c(pu[1], sx, pu[2]), c(0, 0:n/n), type = "S",
col = i <- i + 1)
if(!is.null(q)) {
qv <- quantile(sx, q)
segments(pu[1], q, qv, q, col = i)
segments(qv, q, qv, pu[3], col = i)
}
}
invisible(m)
}
ecdf.list <- function(x, ...)
invisible(do.call("ecdf", c(x, list(...))))
ecdf.lm <- function(x, ..., xlab = "Residuals") {
x <- resid(x)
invisible(NextMethod("ecdf", x = x, ..., xlab = xlab))
}
factor.problems
I do believe that >90% of my Splus misery episodes are
spent trying to cope with factors where I don't want them.
Surely all of you intense
S-plus users have internalized factor-lore, but for the novices
and occasionals out there, here are a few of
my (exaggerated) misadventures (and a few solutions):
1. You've read in a "numeric" variable into a data frame,
& you go ahead to summarize it. Unbeknownst to you,
among the 1000 observations,
there was a character in there somewhere (e.g. "Lab error"),
so the whole thing was converted to a factor. You almost give
your colleagues a boxplot of the factor codes instead of the
actual numbers. Luckily, you catch the error: the minimum blood
pressure is 1 (a little low!) and the maximum is 1000 (a little
high!). But next time you might not be so lucky.
SOLUTION: Always assume every column of a data frame is
a factor until proven otherwise. Since you'll forget (if you're me),
paste a note to your CRT: "it's a factor, stupid".
-------------------------------------
2. You try to add a new observation to a data frame. Unfortunately,
you're introducing a new value to one of the factors, & it
bombs with the message
"replacement values not all in levels(x): NA's generated"
None of the obvious things to try will work.
SOLUTION: Put this definition somewhere above the standard definition
of "[<-.factor":
>get("[<-.factor")
function(x, ..., value)
{
cl <- class(x)
class(x) <- NULL
l <- levels(x)
co <- attr(x, "contrasts")
levels(x) <- NULL
is.na.value <- is.na(value)
if(inherits(value, "factor")) {
decode <- function(value)
{
return(levels(value)[value])
}
value <- decode(value)
class(value) <- NULL
}
jj <- match(value, l)
ind.newvalues <- (is.na(jj) & !is.na.value)
if(any(ind.newvalues)) {
l <- c(l, unique(value[ind.newvalues]))
jj <- match(value, l)
}
x[...] <- jj
levels(x) <- l
attr(x, "contrasts") <- co
class(x) <- cl
x
}
This lets you combine factors whose levels are not
identical. For example,
> tmp1<-data.frame ( z="A", w="B")
> tmp1[3,]<-list("z","w")
> tmp1
z w
1 A B
2 z w
Maybe someone knows why, in the standard version of "[<-.factor",
new level values are not handled??
Also needed are versions of c() and rbind() that handle
factors correctly. Just replacing "[<-.factor" isn't enough.
-------------------------------------
3. You've got a dataset of clinical laboratory results, including
columns "result" and "code", where result is the numerical value
produced by the lab test, and code represents the type of
lab test: "WBC" for white blood count, "K" for potassium, "NA"
for blood sodium etc.. Code became a factor when you read in
the dataset (say, using sas.get()).
Strangely, there are no results for sodium, even though
you know there should be. Furthermore, when you try to use
"natural.join()", it fails with huge dynamic memory requirements.
"natural.join" thinks there are lots of missing values for the key,
which includes "code", even though you know there aren't any
such missing values.
Well, it's obvious, isn't it. (It wasn't to me.)
The string "NA" was converted by factor() into the missing NA.
SOLUTION:
Change "NA"'s to "SO"'s or something in the input data set.
You could also change levels(mydataframe$code), but watch
out; don't turn genuine missing values into sodiums.
-------------------------------------
4. You want to merge data frames, but with an "outer join"
(as SAS merge does by default), not an inner join
(as natural.join and equijoin do),
so that if any key value appears in one data frame and
not the other, you'll see it in the resulting data frame.
SOLUTION: A major pain! But the following works pretty well:
(It relies on "<-.factor" shown above. It is a minor modification of
natural.join() and equijoin().) Not too efficient, but it works.
> get("outer.join")
function(xarg, yarg, debug = F)
{
x <- xarg
y <- yarg
colnames <- function(x)
{
if(!length(val <- dimnames(x)[[2]]))
stop("No column names")
val
}
by <- intersect(colnames(x), colnames(y))
by.x <- by
by.y <- by
for(i in 1:ncol(x))
if(data.class(x[[i]]) == "factor") x[[i]] <- decode(x[[i]])
for(i in 1:ncol(y))
if(data.class(y[[i]]) == "factor") y[[i]] <- decode(y[[i]])
# database equijoin based on columns given in by's
by.to.int <- function(by, mat)
{
# convert column name or number 'by' to column number
# in 'mat' and check that it is valid. 'by' might refer to
# several columns.
if(is.character(by)) {
cn <- dimnames(mat)[[2]]
if(any(duplicated(cn)))
stop("duplicate column names")
by <- match(by, cn)
if(any(is.na(by)))
stop("by is not a column name")
}
else if(is.numeric(by)) {
nc <- dim(mat)[2]
if(any(by > nc) || any(by < 1))
stop(
"by is not between 1 and number of columns")
if(any(trunc(by) != by))
stop("by is not an integer")
}
else {
stop("by must be column name or number")
}
by
}
join.index <- function(x, y)
{
lev <- sort(unique(c(x <- as.vector(x), y <- as.vector(y))))
lev.split <- function(x, lev)
{
split(seq(along = x), factor(x, lev))
}
xs <- lev.split(x, lev)
ys <- lev.split(y, lev) #if(debug)cat("xs\n")
list(x = unlist(lapply(seq(along = lev), function(i, xs, ys)
{
rep(xs[[i]], length(ys[[i]]))
}
, xs, ys)), y = rep(unlist(ys), rep(unlist(lapply(xs, length)),
unlist(lapply(ys, length)))))
}
by.y <- by.to.int(by.y, y)
by.x <- by.to.int(by.x, x)
if(length(by.x) != length(by.y))
stop("length(by.x) != length(by.y)")
if(length(by.x) == 0) {
# should return cartesian product of x and y ...
stop("length(by.x) == 0 -- cannot handle this case yet")
}
else if(length(by.x) > 1) {
# paste together 'by' columns to make keys
x.by.x <- as.vector(x[, by.x[1]])
y.by.y <- as.vector(y[, by.y[1]])
for(i in 2:length(by.x)) {
x.by.x <- paste(x.by.x, as.vector(x[, by.x[i]]), sep =
"\001")
y.by.y <- paste(y.by.y, as.vector(y[, by.y[i]]), sep =
"\001")
}
ij <- join.index(x.by.x, y.by.y)
}
else {
ij <- join.index(x[, by.x], y[, by.y])
}
xmissing <- indices(is.na(match(1:nrow(x), ij[[1]])))
ymissing <- indices(is.na(match(1:nrow(y), ij[[2]])))
result <- cbind(x[ij[[1]], , drop = F], y[ij[[2]], - by.y, drop = F])
x.xmissing <- (x[xmissing, , drop = F])
#dimnames(x.xmissing)[[1]] <- NULL
y.ymissing <- (y[ymissing, , drop = F])
#dimnames(y.ymissing)[[1]] <- NULL
n.xmissing <- nrow(x.xmissing)
n.ymissing <- nrow(y.ymissing)
if(n.xmissing + n.ymissing > 0) {
nullrow <- makeNullRow(result)
nlastrow <- nrow(result)
result[nlastrow + (1:(n.xmissing + n.ymissing)), ] <- nullrow
xcols <- 1:ncol(x)
ycols <- (ncol(x) + 1):ncol(result) #not including y keys.
result[nlastrow + (1:n.xmissing), xcols] <- x.xmissing
y.ymissing.keys <- y.ymissing[, by.y]
result[nlastrow + n.xmissing + (1:n.ymissing), by.x] <-
y.ymissing.keys
y.ymissing.not.keys <- y.ymissing[, - by.y]
result[nlastrow + n.xmissing + (1:n.ymissing), ycols] <-
y.ymissing.not.keys
}
return(result)
}
Example:
> outer.join(tmp1,tmp2)
dtypes factor factor numeric numeric
z w X1 X2
1 A B 1 2
3 z w 1 2
3 1 1 1 NA
4 2 2 NA 2
:~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:
: Roger Day _/ _/ _/_/_/ _/_/ _/_/_/ :
: Director _/ _/ _/ _/ _/ _/ _/ :
: Biostatistics Center _/ _/ _/ _/ _/ _/ :
: U.Pittsburgh Cancer Institute _/ _/ _/_/_/ _/ _/ :
: 3600 Forbes Ave Suite 400A _/ _/ _/ _/ _/ _/ :
: Pittsburgh, PA 15213 _/_/ _/ _/_/ _/_/_/ :
: (412)-647-8008 fax:647-5687 :
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
find.closest.element
From:
9:07 AM
Subject: Re: [S] Finding most closely matching
elementTo: s-news
Thanks to Alan Zaslavsky for a very fast
solution.
Here is what Alan was so kind to send:
y.order_order(y)
y.sort_y[y.order]
y.order[round(approx(y.sort, 1:length(y), xout=z, rule=2,
method='linear')$y)]
or equivalently
order(y)[round(approx(sort(y), 1:length(y), xout=z, rule=2,
method='linear')$y)]
The idea is to interpolate (after arranging in order) and then round the
index, since the fractional part of the index represents the relative
distance from the two adjacent values.
Thanks Alan. Thanks also to Mark Bravington for sending another solution.
-Frank
Frank E Harrell Jr wrote:
> Suppose we have a vector y and z (of differing lengths) and we wish
> to find, for each z, which element of (subscript) of y is closest to each
value of z.
> The folllowing works:
>
> apply(outer(z, y), function(a,b)abs(a-b)), 1, order)[1,]
>
> But this can require a lot of memory for some cases. The following
almost
> works, is blazing fast, and requires minimal memory:
>
> approx(y, 1:length(y), xout=z, rule=2, method='constant')$y
>
> For example:
>
> y <- c( -1.5, 1.22, 0.705, 0.0653, -0.436, -0.66)
> z <- c(1.09, -0.974, -0.149, -0.253, -0.084, 1.33, 1.27, 0.368)
> apply(outer(z, y, function(a,b)abs(a-b)), 1, order)[1,]
> [1] 2 6 4 5 4 2 2 4 # The right answer
>
> > approx(y, 1:length(y), xout=z, rule=2, method='constant')$y
> [1] 3 1 5 5 5 2 2 4
> > approx(y, 1:length(y), xout=z, rule=2, method='constant', f=1)$y
> [1] 2 6 4 4 4 2 2 3
>
> Does anyone see a very fast way to get the desired result?
>
> Thanks -Frank
You can compute the index if the z just above and just below each entry
in y and take the index that gives the smallest distance. E.g.
f1<-function(y,z)
approx(y, 1:length(y), xout=z, rule=2, method='constant')$y
f2<-function(y,z)
approx(y, 1:length(y), xout=z, rule=2, method='constant', f=1)$y
f3<-function(y,z){
i1<-f1(y,z)
i2<-f2(y,z)
ifelse(abs(y[i1]-z) length(reference)] <- NA
# i is length of x and has values in range 1:length(reference) or NA
# following needed if reference is not sorted (used to have 3 orders!)
i <- order(reference)[i] #
# following needed if x is not sorted
i[order(order(x))]
}
which.just.below <-
function(x, reference, strict = T)
{
# output[k] will be index of largest value in reference vector
# less than x[k]. If strict=F, replace 'less than' by
# 'less than or equal to'. Neither x nor reference need be
# sorted, although they should not have NA's (in theory, NA's
# in x are ok, but not in reference).
which.just.above( - x, - reference, strict = strict)
}
----------------------------------------------------------------------------
Bill Dunlap 22461 Mt Vernon-Big Lake Rd
Data Analysis Products Div. of MathSoft, Inc. Mount Vernon, WA 98274
bill@statsci.com 360-428-8146