Re: [R] Functional data analysis - problems with smoothing

2009-10-23 Thread Peter Ehlers

The error message is pretty clear: regardless of what
*you* think, R says that 'isi' is not numeric.

Are you sure that 'isi' is not a *factor* object?
I'm willing to bet that it is.

Use str() to check your data.

 -Peter Ehlers

Benjamin Cheah wrote:

Hi all,

I'm having major issues with smoothing my functional data

I'm referring to Jim Ramsay's
examples in his books. The following error message keeps appearing,
despite all my data being numeric can anyone kindly offer any suggestions?

isi - vector of argument values - i.e. the independent variable of the curves
rlz - data array
TMSfdPar - functional data parameter.


TMSfdPar = fdPar(TMSbasis, 4, 0.01)
TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar)

Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric.

I don't understand why the error message keeps popping up. I've tried playing 
around with Jim Ramsay's datasets and I think my data is organised in a similar 
manner to his, so can't understand what's going on oh, the frustration!

Thanks in advance,

Ben
(the amateur R data analyst and statistician.)



  
__


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] coxph() and survfit()

2009-10-23 Thread Mura Tamakou
Dear All,

I have a question regarding the output of survfit() when I supply a Cox model. 
Lets say for example:

library(survival)

fit - coxph(Surv(time, status == 2) ~ factor(spiders), data = pbc)
fit # HR for spiders is significant


newdata - data.frame(spiders = factor(0:1))
sf - survfit(fit, newdata = newdata)
sum.sf - summary(sfit, times = c(2000, 2500, 3000))

# survival estimates for the yes/no spiders
# and the 3 follow up times
sum.sf$surv

# corresponding lower limits of the
# 95% CI
sum.sf$low

# corresponding upper limits of the
# 95% CI
sum.sf$up


we observe that the 95% CIs overlap!! How is this possible since the HR for 
spiders is significant.

Regards,

Mura

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] cut2 once, bin twice...

2009-10-23 Thread sdanzige

Hello,

I'm using the Hmisc cut2 function to bin a set of data.  It produces bins
that I like with results like this:

[96,270]:171
[69, 96): 54
[49, 69): 40
[35, 49): 28
[28, 35): 14
[24, 28):  8
(Other) : 48

I would like to take a second set of data, and assign it to bins based on
factors defined by my call to cut 2.

Does anyone know how I can do this?

Thank you,
-S
-- 
View this message in context: 
http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26020736.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] cut2 once, bin twice...

2009-10-23 Thread Dieter Menne



sdanzige wrote:
 
 
 I'm using the Hmisc cut2 function to bin a set of data.  It produces bins
 that I like with results like this:
 
 [96,270]:171
 [69, 96): 54
 [49, 69): 40
 [35, 49): 28
 [28, 35): 14
 [24, 28):  8
 (Other) : 48
 
 I would like to take a second set of data, and assign it to bins based on
 factors defined by my call to cut 2.
 

It used to be quite tricky, but on popular request Brian Ripley has added an
example how to extract the intervals using regular expression on the bottom
of the examples for cut (note:cut in base, not cut2 in Hmisc).

If someone knows of an easier way, please correct me. How about adding this
information as attribute to the standard cut?

Dieter




-- 
View this message in context: 
http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26022244.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] twoord.plot y lab size

2009-10-23 Thread Jim Lemon

On 10/22/2009 10:26 PM, Jacob Kasper wrote:

I am using twoord.plot and my Y axis units on the left y overlap. I tried
using cex.axis in my par command but that only adjusted the x units, not the
y. cex.axis in twoord.plot did not help. How do I adjust Y units in the
twoord.plot?
example:
twoord.plot (lx=myear, ly=z, rx=myear, ry=sta,xlim=c(1985,2010),
xlab=Year,ylab=Individuals,
rylab=# Stations)
   

Hi Jacob,
While there is no way to do this in the current version, the upcoming 
version 2.7-2 will have an axislab.cex argument that allows control of 
the axis and tick labels. Thanks for the idea.


Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bayesian regression stepwise function?

2009-10-23 Thread spencerg

Charles C. Berry wrote:

On Thu, 22 Oct 2009, Ben Bolker wrote:





Allan.Y wrote:


Hi everyone,

I am wondering if there exists a stepwise regression function for the
Bayesian regression model.  I tried googling, but I couldn't find
anything.  I know step function exists for regular stepwise 
regression,

but nothing for Bayes.



Why?  That seems so ... un-Bayesian ...


If 'fools rush in where angels fear to tread', then Bayesians 'jump' 
in where frequentists fear to 'step'...


Seriously, there are Bayesian regression approaches that priorize the 
model size (sometimes only implicitly by assigning a prior for the 
inclusion of each candidate regressors). Then they 'jump' between 
models of different sizes.


On CRAN, Package qtlbim (which is specialized to a particular genetics 
problem) implements one such, I think.


Package bqtl does not implement the jumping approach, but does explore 
a model space with differing numbers of regressors for the same (qtl) 
problem.


Perhaps the closest to a general purpose 'stepwise flavored' Bayesian 
regression is implemented in Package BMA, which IIRC borrows step() 
for some of its work.


But CRAN now has more packages than my cortex has neurons, so there are
probably more packages that do something like this. Try

RSiteSearch(jump regression, restric='functions')

and start reading.

library(sos)
jr - ???'jump regression'
summary(jr)
# Downloaded 22 links in 14 packages.
# packages with 3 help pages matching 'jump regression': 
# monomvn, CalciOMatic, polydect; 
# qtl has 1 ...

bs - ???'Bayesian step'
bsw - ???'Bayesian stepwise'
bs. - bs|bsw
summary(bs.)

Total number of matches: 119
Downloaded 115 links in 63 packages.

Packages with at least 2 matches using search pattern 'Bayesian+step | 
Bayesian+stepwise':

   Package Count MaxScore TotalScoreDate
1 DPpackage 82 11 2009-03-17 06:38:42
2   BMA 61  6 2009-05-01 12:21:01
3   ensembleBMA 51  5 2009-07-01 07:16:43
4  MMIX 43  8 2009-08-16 05:47:37
5mclust 42  8 2009-07-14 05:39:10
6adlift 41  4 2009-02-01 06:31:16
7   tgp 34  8 2009-07-29 11:14:36
8   geoRglm 34  6 2009-10-19 19:21:57
9 G1DBN 32  5 2008-01-24 16:27:29
10   bayesm 31  3 2008-06-14 15:31:09
11 bqtl 31  3 2009-01-15 07:02:43
12   qtlbim 2   12 13 2009-07-29 10:50:13
13  BAYSTAR 2   11 12 2009-10-19 19:17:18
14  gap 24  5 2009-01-15 07:03:58
15   MSBVAR 24  5 2009-07-29 10:49:07
16bayesSurv 23  5 2008-12-15 10:06:16
17  bnlearn 22  3 2009-09-01 08:25:43
18  spBayes 22  3 2009-03-28 07:33:41
19  dlm 21  2 2009-07-01 07:16:28
20GOFSN 21  2 2009-07-01 07:47:10
21 MCMChybridGP 21  2 2009-10-01 05:27:43
22 mgcv 21  2 2009-08-24 16:03:12
23nlreg 21  2 2007-11-02 13:14:56
24  nlt 21  2 2009-03-02 07:30:34
25  pARtial 21  2 2006-11-01 02:59:23
26plink 21  2 2009-07-01 07:19:44
27  qtl 21  2 2009-09-17 15:46:05
28   timsac 21  2 2009-03-28 07:27:28

#  This took more longer to describe it than run. 
# Hope this helps. 
# Spencer


HTH,

Chuck






--
View this message in context: 
http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p26015081.html 


Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.



Charles C. Berry(858) 534-2098
Dept of Family/Preventive 
Medicine

E mailto:cbe...@tajo.ucsd.eduUC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 
92093-0901


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.




--
Spencer Graves, PE, PhD
President and Chief Operating Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help

Re: [R] Functional data analysis - problems with smoothing

2009-10-23 Thread spencerg
Hi, Ben: 



 Which chapter in which of Ramsay's books? 



 For his most recent book (Functional Data Analysis with R and 
Matlab, with Giles Hooker and your's truly), there is one script file 
for each chapter with names like fdarm-ch01.R, ... fdarm-ch11.R.  
These script files reproduce nearly all the examples in the book 
including all but one of the 76 figures (according to what I wrote for 
back cover).  The preface tells you how to find those script files: 


system.file('scripts', package='fda')
[1] c:/Users/sgraves/R/R-2.9.1/library/fda/scripts


 Obviously, the answer is the path on my installation;  the answer 
you get should tell you where to look on the computer (and start-up 
sequence) you are using. 



 This same scripts directory contains similar partial scripts for 
Functional Data Analysis and Applied Functional Data Analysis.  If 
you'd care to send me updates to any of those partial scripts, I can see 
that they are included in the package to make those examples easier for 
others in the future. 



 Hope this helps. 
 Spencer Graves
   


Peter Ehlers wrote:

The error message is pretty clear: regardless of what
*you* think, R says that 'isi' is not numeric.

Are you sure that 'isi' is not a *factor* object?
I'm willing to bet that it is.

Use str() to check your data.

 -Peter Ehlers

Benjamin Cheah wrote:

Hi all,

I'm having major issues with smoothing my functional data

I'm referring to Jim Ramsay's
examples in his books. The following error message keeps appearing,
despite all my data being numeric can anyone kindly offer any 
suggestions?


isi - vector of argument values - i.e. the independent variable of 
the curves

rlz - data array
TMSfdPar - functional data parameter.


TMSfdPar = fdPar(TMSbasis, 4, 0.01)
TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar)

Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric.

I don't understand why the error message keeps popping up. I've tried 
playing around with Jim Ramsay's datasets and I think my data is 
organised in a similar manner to his, so can't understand what's 
going on oh, the frustration!


Thanks in advance,

Ben
(the amateur R data analyst and statistician.)



  
__ 




[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.




--
Spencer Graves, PE, PhD
President and Chief Operating Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] PDF too large, PNG bad quality

2009-10-23 Thread Jim Lemon

On 10/23/2009 06:07 AM, Lasse Kliemann wrote:

I wish to save a scatter plot comprising approx. 2 million points
in order to include it in a LaTeX document.

Using 'pdf(...)' produces a file of size about 20 MB, which is
useless.

Using 'cairo_pdf(...)' produces a smaller file, around 3 MB. This
is still too large. Not only that the document will be too large,
but also PDF viewers choke on this. Moreover, Cairo has problems
with text: by default text looks ugly, like scaled bitmaps. After
hours of trying different settings, I discovered that choosing a
different font family can help, e.g.: 'par(family=Mono)'. This
gives good-looking text. Yet, the problem with the file size
remains.

There exists the hint to produdc EPS instead and then convert to
PDF using 'epstopdf'. The resulting PDF files are slightly
smaller, but still too large, and PDF viewers still don't like
it.

So I gave PNG a try. PNG files are much smaller and PDF viewers
have no trouble with them. However, fonts look ugly. The same
trick that worked for Cairo PDF has no effect for PNG. When I
view the PNGs with a dedicated viewer like 'qiv', even the fonts
look good. But not when included in LaTeX; I simply use
'\includegraphics{...}' and run the document through 'pdflatex'.

I tried both, creating PNG with 'png(...)' and converting from
PDF to PNG using 'convert' from ImageMagick.

So my questions are:

- Is there a way to produce sufficiently lean PDFs directly in R,
   even when the plot comprises several million points?

- How to produce a PNG that still looks nice when included in a
   LaTeX PDF document?

Any hints will be greatly appreciated.
   

Hi Lasse,
I may be right off the track, but I can't make much sense of 2 million 
points in a scatterplot. If you are interested in the density of points 
within the plot, you could compute this using something like the bkde2 
function in the KernSmooth package and then plot that using something 
like image.


Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] data file with columns of unequal length

2009-10-23 Thread William Simpson
I am running an expt that presents a point process input x and
measures a point process output y. The times of each event are
recorded. The lengths of the data records of x and y are necessarily
different, and can be different by a factor of 10. I would like to
save these data after each experiment as a file with two columns, one
for x and one for y.

However, R dataframes require columns of equal length. One solution is
to fill the empty places in y with NAs so it has the same length as
x. I view that as unsatisfactory (there are in reality no missing
values). Another possibility is to store x and y in separate files. I
also view that as unsatisfactory (it is too easy to lose track of the
y file corresponding to a given x file).

Can anyone suggest a way to deal with this situation?

Thanks very much for any help.

Bill

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread Jim Lemon

On 10/23/2009 07:58 PM, William Simpson wrote:

I am running an expt that presents a point process input x and
measures a point process output y. The times of each event are
recorded. The lengths of the data records of x and y are necessarily
different, and can be different by a factor of 10. I would like to
save these data after each experiment as a file with two columns, one
for x and one for y.

However, R dataframes require columns of equal length. One solution is
to fill the empty places in y with NAs so it has the same length as
x. I view that as unsatisfactory (there are in reality no missing
values). Another possibility is to store x and y in separate files. I
also view that as unsatisfactory (it is too easy to lose track of the
y file corresponding to a given x file).

Can anyone suggest a way to deal with this situation?

   

Hi Bill,

xy-list(x=1:10,y=1:100)

Note that this cheerfully ignores how you are going to figure out which 
x goes with which y(s).


Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] standard errors in bbmle

2009-10-23 Thread alexander russell
Hello,
Mle2 is a little unforthcoming in the matter of standard errors? Is there a
way to ask the program to supply standard errors along with estimates in
cases when it doesn't print them 'voluntarily'?
regards,
s

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread Peter Dalgaard
Jim Lemon wrote:
 On 10/23/2009 07:58 PM, William Simpson wrote:
 I am running an expt that presents a point process input x and
 measures a point process output y. The times of each event are
 recorded. The lengths of the data records of x and y are necessarily
 different, and can be different by a factor of 10. I would like to
 save these data after each experiment as a file with two columns, one
 for x and one for y.

 However, R dataframes require columns of equal length. One solution is
 to fill the empty places in y with NAs so it has the same length as
 x. I view that as unsatisfactory (there are in reality no missing
 values). Another possibility is to store x and y in separate files. I
 also view that as unsatisfactory (it is too easy to lose track of the
 y file corresponding to a given x file).

 Can anyone suggest a way to deal with this situation?


 Hi Bill,
 
 xy-list(x=1:10,y=1:100)
 
 Note that this cheerfully ignores how you are going to figure out which
 x goes with which y(s).

As I understand it, they don't come in pairs anyway. For the same reason
a data frame is just the wrong kind of data structure. If you don't want
separate data files, you can use one file with two columns where the
second column is (say) 1 for the x and 2 for the y.

(There are other options, like concatenating the x and the y with some
sort of separator inbetween, but it easily gets painful to read them
back in.)

-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] standard errors in bbmle

2009-10-23 Thread Peter Dalgaard
alexander russell wrote:
 Hello,
 Mle2 is a little unforthcoming in the matter of standard errors? Is there a
 way to ask the program to supply standard errors along with estimates in
 cases when it doesn't print them 'voluntarily'?
 regards,
 s

What did you do? Which cases? Did you look at the examples? (as in,
example(mle2) )

-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread William Simpson
The way you do it is to compute the cross-intensity function (you can
google this; a key name is David Brillinger). The general problem is
that of system identification for point processes.

Bill

On Fri, Oct 23, 2009 at 10:31 AM, Jim Lemon j...@bitwrit.com.au wrote:
 On 10/23/2009 07:58 PM, William Simpson wrote:

 I am running an expt that presents a point process input x and
 measures a point process output y. The times of each event are
 recorded. The lengths of the data records of x and y are necessarily
 different, and can be different by a factor of 10. I would like to
 save these data after each experiment as a file with two columns, one
 for x and one for y.

 However, R dataframes require columns of equal length. One solution is
 to fill the empty places in y with NAs so it has the same length as
 x. I view that as unsatisfactory (there are in reality no missing
 values). Another possibility is to store x and y in separate files. I
 also view that as unsatisfactory (it is too easy to lose track of the
 y file corresponding to a given x file).

 Can anyone suggest a way to deal with this situation?



 Hi Bill,

 xy-list(x=1:10,y=1:100)

 Note that this cheerfully ignores how you are going to figure out which x
 goes with which y(s).

 Jim



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread William Simpson
 As I understand it, they don't come in pairs anyway.
Correct.

 For the same reason
 a data frame is just the wrong kind of data structure. If you don't want
 separate data files, you can use one file with two columns where the
 second column is (say) 1 for the x and 2 for the y.
Could you explain this further? I don't really get it.
My xs and ys have different lengths. How would I read them in?

Ideally I would just read in the x column and y columns separately

x-read.file(file.dat, column1)
y-read.file(file.dat, column2)
But I know of no way to do that...

[Oh, by the way, I never mentioned this, the x and ys are event times
and are in ascending order (by time of occurrence).]

 (There are other options, like concatenating the x and the y with some
 sort of separator inbetween, but it easily gets painful to read them
 back in.)


Thanks
Bill

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread William Simpson
Thanks Jim. BTW the times in x and y are in ascending order (time of
occurrence).

If I do it this way, how do I actually read the data in and store in
the file? Toy code, please.

Bill


 Hi Bill,

 xy-list(x=1:10,y=1:100)

 Note that this cheerfully ignores how you are going to figure out which x
 goes with which y(s).

 Jim


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread Peter Dalgaard
William Simpson wrote:
 As I understand it, they don't come in pairs anyway.
 Correct.
 
 For the same reason
 a data frame is just the wrong kind of data structure. If you don't want
 separate data files, you can use one file with two columns where the
 second column is (say) 1 for the x and 2 for the y.
 Could you explain this further? I don't really get it.
 My xs and ys have different lengths. How would I read them in?

Same format as when you have data from two different groups. Don't
really know how to make it clearer than that. Like the sleep dataset
(OK, so they are actually paired...).

 
 Ideally I would just read in the x column and y columns separately
 
 x-read.file(file.dat, column1)
 y-read.file(file.dat, column2)
 But I know of no way to do that...
 
 [Oh, by the way, I never mentioned this, the x and ys are event times
 and are in ascending order (by time of occurrence).]
 
 (There are other options, like concatenating the x and the y with some
 sort of separator inbetween, but it easily gets painful to read them
 back in.)

-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread William Simpson
OK thanks, I look at sleep and get it

Bill

On Fri, Oct 23, 2009 at 12:21 PM, Peter Dalgaard
p.dalga...@biostat.ku.dk wrote:
 William Simpson wrote:
 As I understand it, they don't come in pairs anyway.
 Correct.

 For the same reason
 a data frame is just the wrong kind of data structure. If you don't want
 separate data files, you can use one file with two columns where the
 second column is (say) 1 for the x and 2 for the y.
 Could you explain this further? I don't really get it.
 My xs and ys have different lengths. How would I read them in?

 Same format as when you have data from two different groups. Don't
 really know how to make it clearer than that. Like the sleep dataset
 (OK, so they are actually paired...).


 Ideally I would just read in the x column and y columns separately

 x-read.file(file.dat, column1)
 y-read.file(file.dat, column2)
 But I know of no way to do that...

 [Oh, by the way, I never mentioned this, the x and ys are event times
 and are in ascending order (by time of occurrence).]

 (There are other options, like concatenating the x and the y with some
 sort of separator inbetween, but it easily gets painful to read them
 back in.)

 --
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
  (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
 ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread Jim Lemon

On 10/23/2009 10:07 PM, William Simpson wrote:

Thanks Jim. BTW the times in x and y are in ascending order (time of
occurrence).

If I do it this way, how do I actually read the data in and store in
the file? Toy code, please.

   

Hi Bill,
This seems a bit like some heartbeat data that I had to deal with some 
years ago. There were two output streams, one the time of each R wave 
and the other an asynchronous stimulus. In that case, I had to work out 
code to interdigitate the signals and do some processing of the 
resulting data. It sounds like you have two files, each recording times 
of input and output respectively one value to a line. If so, my first 
guess is:


x-as.vector(read.table(input_times.dat))
y-as.vector(read.table(output_times.dat))
xy-list(x,y)
# to store this list in a file
write.csv(xy,file=xy_23_10_2009.csv)

Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data file with columns of unequal length

2009-10-23 Thread Jim Lemon

On 10/23/2009 10:07 PM, William Simpson wrote:

Thanks Jim. BTW the times in x and y are in ascending order (time of
occurrence).

If I do it this way, how do I actually read the data in and store in
the file? Toy code, please.



Hi Bill,
This seems a bit like some heartbeat data that I had to deal with some 
years ago. There were two output streams, one the time of each R wave 
and the other an asynchronous stimulus. In that case, I had to work out 
code to interdigitate the signals and do some processing of the 
resulting data. It sounds like you have two files, each recording times 
of input and output respectively one value to a line. If so, my first 
guess is:


x-as.vector(read.table(input_times.dat))
y-as.vector(read.table(output_times.dat))
xy-list(x,y)
# to store this list in a file
# oops, bit too quick on the keyboard
dput(xy,file=xy_23_10_2009.R)

Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima crashes too

2009-10-23 Thread Alberto Monteiro
Barry Rowlingson wrote:

 However, arima crashes for this:

 arima(c(1.71, 1.78, 1.95, 1.59, 2.13), order=c(1,0,0))
 
  I'm not getting what I'd call 'crashes' with your arma or arima
 examples- I get an error message and a warning:
 
 arma(c(2.01, 2.22, 2.09, 2.17, 2.42), order=c(1,0))
 Error in AA %*% t(X) : requires numeric/complex matrix/vector arguments
 In addition: Warning message:
 In ar.ols(x, order.max = k, aic = FALSE, demean = FALSE, intercept =
 include.intercept) :
   model order: 2singularities in the computation of the projection
 matrixresults are only valid up to model order1
 
  You've not told us what you get, 

Because I get a message in Portuguese. But the meaning is precisely
that one: Error in AA %*% t(X).

 and the phrase 'crash' normally
 means some kind of memory error that *terminates* a running R 
 session. 

Ah, maybe then crash is not the correct word. It stops running
whatever it is running. Maybe abort is the best word?

 Are you really crashing R such that it terminates? In which 
 case, what version number/platform etc?
 
I mean that, if I run a loop, it doesn't finish. Or, more 
catastrophically, if I am running a loop and saving data to an
open file, it terminates the loop and does not close the file.

Reproducible example:

test.arima - function() {
  lets.crash.arima - c(71, 78, 95, 59) # , 113
  for (x in 90:120) {
reg - arima(c(lets.crash.arima, x), order = c(1,0,0))
cat(ok for x =, x, \n)
  }
  cat(close file and prepare a nice summary\n)
  return(arima passed the test)
}

test.arima()

As you can see, the loop aborts, the function never returns, with
potentially nasty effects (namely: I have to finish R with q() to
close the files and examine them).

Alberto 'Darth Albmont' Monteiro

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] opposite estimates from zeroinfl() and hurdle()

2009-10-23 Thread Tord Snäll

Dear all,
A question related to the following has been asked on R-help before, but 
I could not find any answer to it. Input will be much appreciated.


I got an unexpected sign of the slope parameter associated with a 
covariate (diam) using zeroinfl(). It led me to compare the estimates 
given by zeroinfl() and hurdle():


The (significant) negative estimate here is surprising, given the 
biology of the species:


 summary(zeroinfl(bnl ~ 1| diam, dist = poisson, data = valdaekar, 
EM = TRUE))

Count model coefficients (poisson with log link):
 Estimate Std. Error z value Pr(|z|)   (Intercept)  3.74604
0.02635   142.2   2e-16 ***


Zero-inflation model coefficients (binomial with logit link):
 Estimate Std. Error z value Pr(|z|)  (Intercept)  21.7510 
7.6525   2.842  0.00448 **

diam -1.1437 0.3941  -2.902  0.00371 **

Number of iterations in BFGS optimization: 1
Log-likelihood: -582.8 on 3 Df


The hurdle model gives the same estimates, but with opposite (and 
expected) signs of the parameters:


summary(hurdle(bnl ~ 1| diam, dist = poisson, data = valdaekar))
Count model coefficients (truncated poisson with log link):
 Estimate Std. Error z value Pr(|z|)   (Intercept)  3.74604
0.02635   142.2   2e-16 ***

Zero hurdle model coefficients (binomial with logit link):
 Estimate Std. Error z value Pr(|z|)  (Intercept) -21.7510 
7.6525  -2.842  0.00448 **

diam  1.1437 0.3941   2.902  0.00371 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 8
Log-likelihood: -582.8 on 3 Df

Why is this so?

thanks,
Tord
Windows NT, R 2.8.1, pcsl 1.03

--

Tord Snäll
Department of Ecology / Swedish Species Information Centre
Swedish University of Agricultural Sciences (SLU)
P.O. 7044, SE-750 07 Uppsala, Sweden
Office/Mobile/Fax
+46-18-672612/+46-76-7662612/+46-18-673537
www.ekol.slu.se/staff_tordsnall
www.artdata.slu.se/personal/fototsn.asp

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] a problem about integrate function in R .thank you .

2009-10-23 Thread andrew
I don't seem to get a problem with this.  Have you tried a Monte Carlo
approach to verify that you are getting incorrect answers?

For me, I get when the upper is 1 that

 integrate(e2, lower = 0, upper = 1)
-0.2820948 with absolute error  5e-05

 sum(e2(runif(1)))/1
[1] -0.2825667

which seems to be consistent, while for upper = 0.75 I get

 0.75*sum(e2(runif(1, min=0, max=0.75)))/1
[1] -0.2333506
 integrate(e2,lower=0,upper = 0.75)
-0.2341178 with absolute error  7.8e-05



On Oct 23, 2:52 pm, fuzuo xie xiefu...@gmail.com wrote:
 e2 - function(x) {
     out - 0*x
     for(i in 1:length(x))
     out[i] -integrate(function(y) qnorm(y),lower=0,upper=x[i])$value
     out }
 integrate(e2,lower=0, upper=a)$value

 above is my code , when a is small , say a0.45  the result is right .
 however , when a0.5
 the result is incorrect .  why ?    thank you .

         [[alternative HTML version deleted]]

 __
 r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima crashes too

2009-10-23 Thread Barry Rowlingson
On Fri, Oct 23, 2009 at 12:32 PM, Alberto Monteiro
albm...@centroin.com.br wrote:

 I mean that, if I run a loop, it doesn't finish. Or, more
 catastrophically, if I am running a loop and saving data to an
 open file, it terminates the loop and does not close the file.

 Reproducible example:

 test.arima - function() {
  lets.crash.arima - c(71, 78, 95, 59) # , 113
  for (x in 90:120) {
    reg - arima(c(lets.crash.arima, x), order = c(1,0,0))
    cat(ok for x =, x, \n)
  }
  cat(close file and prepare a nice summary\n)
  return(arima passed the test)
 }

 test.arima()

 As you can see, the loop aborts, the function never returns, with
 potentially nasty effects (namely: I have to finish R with q() to
 close the files and examine them).

 If you're doing anything in a loop that has the potential to fail
because of singularities or other conditions when your model can't be
fitted, you need to stick what you are doing in a 'try' clause. This
lets you trap errors and do something with them.

 Plenty of examples in help(try) or this from me:

 for(i in 1:10){
 print(solve(matrix(c(3,3,3,i),2,2)))
 }

 This stops the loop at i=3. Now stick it in a try() clause:

 for(i in 1:10){
 print(try(solve(matrix(c(3,3,3,i),2,2
 }

 and it gives a warning and carries on. If you want your code to do
something with the failure cases then the help for try() tells you
what to look for.

 I'm not sure why your arima produces an error, but I'm assuming the
numbers are such that the model can't be fitted. I don't really know
what arima is doing.

Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] reset par() within plot layout

2009-10-23 Thread Janke ten Holt
Dear list,

I would like to produce a matrix of plots, where par() is reset after
each plot (see below [simplified] example). When I use layout() to do
so, I seem to also reset the layout. I have not been able to figure out
how to prevent this from happening.

Any help is greatly appreciated!
Janke

Example code:
#Desired result is a layout of 2 plots: one red and one black
layout(matrix(1:2, nr=2))
par.ini - par(no.readonly=TRUE)
par(col=red)
plot(1:100)

par(par.ini)

plot(1:10)

--
Janke ten Holt
Dept. of Psychology/Sociology
University of Groningen, the Netherlands

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reset par() within plot layout

2009-10-23 Thread stephen sefick
Maybe something like this?

#Desired result is a layout of 2 plots: one red and one black
par(mfrow=c(2,1))

par(col=red)
plot(1:100)

par(col=black)
plot(1:10)

On Fri, Oct 23, 2009 at 7:26 AM, Janke ten Holt j.c.ten.h...@rug.nl wrote:
 Dear list,

 I would like to produce a matrix of plots, where par() is reset after
 each plot (see below [simplified] example). When I use layout() to do
 so, I seem to also reset the layout. I have not been able to figure out
 how to prevent this from happening.

 Any help is greatly appreciated!
 Janke

 Example code:
 #Desired result is a layout of 2 plots: one red and one black
 layout(matrix(1:2, nr=2))
 par.ini - par(no.readonly=TRUE)
 par(col=red)
 plot(1:100)

 par(par.ini)

 plot(1:10)

 --
 Janke ten Holt
 Dept. of Psychology/Sociology
 University of Groningen, the Netherlands

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Functional data analysis - problems with smoothing

2009-10-23 Thread Benjamin Cheah
Hi Peter and Spencer (and everyone else out there),

Thank you for your prompt reply to my post re. FDA - FYI - i'm very new to R 
and FDA. Currently half way through a biostats masters degree in Sydney, 
Australia! Absolutely loving it and I hope to include R computing on my future 
CV!

Just to clarify things:

When I entered what you suggested, this is what I obtained with my own 
dataset
 str(isi)
 num [1:14] 1 1.5 2 2.5 3 3.5 4 5 7 10 ...

Similarly, for Jim Ramsay's dataset (chapter five of Use R book):
 str(age)
 num [1:31] 1 1.25 1.5 1.75 2 3 4 5 6 7 ...

Doesn't this mean that my data object, 'isi' is numeric?

Also, I was looking through Jim Ramsay's datasets - I'm fairly sure that my 
vector, 'isi' was similarly organised to his vector, 'age' in chapter 5 of his 
new Use R FDA book, so I'm fairly certain that the data is numeric but it 
obviously isn't

Hope this email clarified things - apologies if it hasn't - still very new to 
the terminology and overall feel of R.

The R user network is terrific! Thank you for caring!

Ben





From: Peter Ehlers ehl...@ucalgary.ca

Cc: R-help@r-project.org
Sent: Fri, 23 October, 2009 5:07:15 PM
Subject: Re: [R] Functional data analysis - problems with smoothing

The error message is pretty clear: regardless of what
*you* think, R says that 'isi' is not numeric.

Are you sure that 'isi' is not a *factor* object?
I'm willing to bet that it is.

Use str() to check your data.

-Peter Ehlers

Benjamin Cheah wrote:
 Hi all,
 
 I'm having major issues with smoothing my functional data
 
 I'm referring to Jim Ramsay's
 examples in his books. The following error message keeps appearing,
 despite all my data being numeric can anyone kindly offer any suggestions?
 
 isi - vector of argument values - i.e. the independent variable of the curves
 rlz - data array
 TMSfdPar - functional data parameter.
 
 TMSfdPar = fdPar(TMSbasis, 4, 0.01)
 TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar)
 Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric.
 
 I don't understand why the error message keeps popping up. I've tried playing 
 around with Jim Ramsay's datasets and I think my data is organised in a 
 similar manner to his, so can't understand what's going on oh, the 
 frustration!
 
 Thanks in advance,
 
 Ben
 (the amateur R data analyst and statistician.)
 
 
 
   
 __
 
 
 [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 



  
__


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reset par() within plot layout

2009-10-23 Thread Janke ten Holt
OK, maybe I oversimplified the example, because you are of course correct...

I am working on a function, which I feed data to plot in layouts of
specified dimensions. I want to be able to set any par() variables for
each plot in the layout. Because the next plot does not know which par()
variables were changed in the previous plot, I want to reset all par()
variables after each plot, rather than simply change back one specific
par() variable.
So I really need to reset par() without 'messing with' my layout.

stephen sefick wrote:
 Maybe something like this?
 
 #Desired result is a layout of 2 plots: one red and one black
 par(mfrow=c(2,1))
 
 par(col=red)
 plot(1:100)
 
 par(col=black)
 plot(1:10)
 
 On Fri, Oct 23, 2009 at 7:26 AM, Janke ten Holt j.c.ten.h...@rug.nl wrote:
 Dear list,

 I would like to produce a matrix of plots, where par() is reset after
 each plot (see below [simplified] example). When I use layout() to do
 so, I seem to also reset the layout. I have not been able to figure out
 how to prevent this from happening.

 Any help is greatly appreciated!
 Janke

 Example code:
 #Desired result is a layout of 2 plots: one red and one black
 layout(matrix(1:2, nr=2))
 par.ini - par(no.readonly=TRUE)
 par(col=red)
 plot(1:100)

 par(par.ini)

 plot(1:10)

 --
 Janke ten Holt
 Dept. of Psychology/Sociology
 University of Groningen, the Netherlands

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

 
 


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Wavelets

2009-10-23 Thread Ashta
Hi all,

I am trying to do wavelets and I got  an error message saying The
length of data is not a power of 2
Is there a way of handing  that? or should the data length be  exactly
the power of  2?
I am using   R version 2.9.2 (2009-08-24)
The is  library(wavethresh).

wds - wd(ds$v,filter.number=1)

Thanks

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Wavelets

2009-10-23 Thread stephen sefick
I believe that you can extend it to the power of 2 by padding it with
zeros.  I don't remember if it is at the begining or the end of the
time series.
hth

stephen

On Fri, Oct 23, 2009 at 7:56 AM, Ashta sewa...@gmail.com wrote:
 Hi all,

 I am trying to do wavelets and I got  an error message saying The
 length of data is not a power of 2
 Is there a way of handing  that? or should the data length be  exactly
 the power of  2?
 I am using   R version 2.9.2 (2009-08-24)
 The is  library(wavethresh).

 wds - wd(ds$v,filter.number=1)

 Thanks

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Wavelets

2009-10-23 Thread Josué Polanco
hi all,

I don't use the package wavethresh, but I suppose that
should be possible to fill of ZEROS (zero padding) the time series,
such that its length has a 2^n elements. Another way is to
cut (remove) some elements, so that, the new length
has 2^n elements.

You can get more info. @
A practical guide to Wavelet Analysis
C. torrence and G. P.Compo
Bulleting of the American Meteorology. Society.
1998

I hope this help you,

--
Josue Polanco


On Fri, Oct 23, 2009 at 2:56 PM, Ashta sewa...@gmail.com wrote:
 Hi all,

 I am trying to do wavelets and I got  an error message saying The
 length of data is not a power of 2
 Is there a way of handing  that? or should the data length be  exactly
 the power of  2?
 I am using   R version 2.9.2 (2009-08-24)
 The is  library(wavethresh).

 wds - wd(ds$v,filter.number=1)

 Thanks

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Josué Mosés Polanco Martínez
Correo-e alternativo jom...@linuxmail.org

It is a wasted day unless you have learned something new and made
someone smile -Mark Weingartz.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] twoord.plot y lab size

2009-10-23 Thread Peter Ehlers

For now, you might be able to use par(las=1). You may have
to make enough margin room with par(mar=...) and you'll
probably want to set ylab= and add the y-label with
mtext().

 -Peter Ehlers

Jim Lemon wrote:

On 10/22/2009 10:26 PM, Jacob Kasper wrote:

I am using twoord.plot and my Y axis units on the left y overlap. I tried
using cex.axis in my par command but that only adjusted the x units, 
not the

y. cex.axis in twoord.plot did not help. How do I adjust Y units in the
twoord.plot?
example:
twoord.plot (lx=myear, ly=z, rx=myear, ry=sta,xlim=c(1985,2010),
xlab=Year,ylab=Individuals,
rylab=# Stations)
   

Hi Jacob,
While there is no way to do this in the current version, the upcoming 
version 2.7-2 will have an axislab.cex argument that allows control of 
the axis and tick labels. Thanks for the idea.


Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to find moving averages within each subgroup of a data frame

2009-10-23 Thread clair.crossup...@googlemail.com
Dear all,

If I have the following data frame:

 set.seed(21)
 df1 - data.frame(col1=c(rep('a',5), rep('b',5), rep('c',5)), 
 col4=rnorm(1:15))

   col1 col4
1 a  0.793013171
2 a  0.522251264
3 a  1.74641
4 a -1.271336123
5 a  2.197389533
6 b  0.433130777
7 b -1.570199630
8 b -0.934905667
9 b  0.063493345
10b -0.002393336
11c -2.276781240
12c  0.757412225
13c -0.548405554
14c  0.172549478
15c  0.562853068

How can i make a 2 point moving average within each group? i.e.

col1col4SMA
a   0.793013171 NA
a   0.522251264 0.657632218
a   1.74641 1.134236753
a   -1.2713361230.237443059
a   2.197389533 0.463026705
b   0.433130777 NA
b   -1.57019963 -0.568534427
b   -0.934905667-1.252552649
b   0.063493345 -0.435706161
b   -0.0023933360.030550005
c   -2.27678124 NA
c   0.757412225 -0.759684508
c   -0.5484055540.104503336
c   0.172549478 -0.187928038
c   0.562853068 0.367701273

From what i've read, i was thinking it would be something along the

lines of:

 aggregate(df1$col4, by=list(df1$col1), function(x) {filter(x, rep(1/2,2), 
 sides=1 )} )

Error in aggregate.data.frame(as.data.frame(x), ...) :
  'FUN' must always return a scalar

But this tells me (i think) that aggregate should only return a single
value per group. So what i need, i guess, is something that takes all
the values in a given group, and returns a vector of the same length.
Not sure which function to use for that.

Thanks in advance,

C xx

P.S. on a side note, is it possible to extract the values which
aggregate passes to the function(x) in my example above?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data frame is killing me! help

2009-10-23 Thread bbslover



Steve Lianoglou-6 wrote:
 
 Hi,
 
 On Oct 22, 2009, at 2:35 PM, bbslover wrote:
 
 Usage
 data(gasoline)
 Format
 A data frame with 60 observations on the following 2 variables.
 octane
 a numeric vector. The octane number.
 NIR
 a matrix with 401 columns. The NIR spectrum

 and I see the gasoline data to see below
 NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR.1694 nm NIR.1696  
 nm
 NIR.1698 nm NIR.1700 nm
 1 1.242645 1.250789 1.246626 1.250985 1.264189 1.244678 1.245913  
 1.221135
 2 1.189116 1.223242 1.253306 1.282889 1.215065 1.225211 1.227985  
 1.198851
 3 1.198287 1.237383 1.260979 1.276677 1.218871 1.223132 1.230321  
 1.208742
 4 1.201066 1.233299 1.262966 1.272709 1.211068 1.215044 1.232655  
 1.206696
 5 1.259616 1.273713 1.296524 1.299507 1.226448 1.230718 1.232864  
 1.202926
 6 1.24109 1.262138 1.288401 1.291118 1.229769 1.227615 1.22763  
 1.207576
 7 1.245143 1.265648 1.274731 1.292441 1.218317 1.218147 1.73  
 1.200446
 8 1.222581 1.245782 1.26002 1.290305 1.221264 1.220265 1.227947  
 1.188174
 9 1.234969 1.251559 1.272416 1.287405 1.211995 1.213263 1.215883  
 1.196102

 look at this NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR. 
 1694 nm
 NIR.1696 nm NIR.1698 nm NIR.1700 nm

 how can I add letters NIR to my variable, because my 600  
 independents never
 have NIR as the prefix. however, it is needed to model the plsr.   for
 example aa=plsr(y~NIR, data=data ,), the prefix NIR is  
 necessary, how
 can I do with it?
 
 I'm not really sue that I'm getting you, but if your problem is that  
 the column names of your data.frame don't match the variable names  
 you'd like to use in your formula, just change the colnames of your  
 data.frame to match your formula.
 
 BTW - I have no idea where to get this gasoline data set, so I'm just  
 imagining:
 
 eg.
 colnames(gasoline) - c('put', 'the', 'variable', 'names', 'that',  
 'you', 'want', 'here')
 
 -steve
 
 --
 Steve Lianoglou
 Graduate Student: Computational Systems Biology
|  Memorial Sloan-Kettering Cancer Center
|  Weill Medical College of Cornell University
 Contact Info: http://cbio.mskcc.org/~lianos/contact
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 

thanks for you. but the numbers of indenpendence are so many, it is not easy
to identify them one by one,  is there some better way?


-- 
View this message in context: 
http://www.nabble.com/data-frame-is-killing-me%21-help-tp26015079p26024985.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] opposite estimates from zeroinfl() and hurdle()

2009-10-23 Thread Tord Snäll

Dear all,
A question related to the following has been asked on R-help before, but 
I could not find any answer to it. Input will be much appreciated.


I got an unexpected sign of the slope parameter associated with a 
covariate (diam) using zeroinfl(). It led me to compare the estimates 
given by zeroinfl() and hurdle():


The (significant) negative estimate here is surprising, given the 
biology of the species:


 summary(zeroinfl(bnl ~ 1| diam, dist = poisson, data = valdaekar, 
EM = TRUE))

Count model coefficients (poisson with log link):
  Estimate Std. Error z value Pr(|z|)   (Intercept)  
3.746040.02635   142.2   2e-16 ***


Zero-inflation model coefficients (binomial with logit link):
  Estimate Std. Error z value Pr(|z|)  (Intercept)  
21.7510 7.6525   2.842  0.00448 **

diam -1.1437 0.3941  -2.902  0.00371 **

Number of iterations in BFGS optimization: 1
Log-likelihood: -582.8 on 3 Df


The hurdle model gives the same estimates, but with opposite (and 
expected) signs of the parameters:


summary(hurdle(bnl ~ 1| diam, dist = poisson, data = valdaekar))
Count model coefficients (truncated poisson with log link):
  Estimate Std. Error z value Pr(|z|)   (Intercept)  
3.746040.02635   142.2   2e-16 ***

Zero hurdle model coefficients (binomial with logit link):
  Estimate Std. Error z value Pr(|z|)  (Intercept) 
-21.7510 7.6525  -2.842  0.00448 **

diam  1.1437 0.3941   2.902  0.00371 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 8
Log-likelihood: -582.8 on 3 Df

Why is this so?

thanks,
Tord
Windows NT, R 2.8.1, pcsl 1.03

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Wavelets

2009-10-23 Thread Josué Polanco
hi all,

On Fri, Oct 23, 2009 at 3:05 PM, stephen sefick ssef...@gmail.com wrote:
 I believe that you can extend it to the power of 2 by padding it with
 zeros.  I don't remember if it is at the begining or the end of the
 time series.
 hth

at the end of the time series

Cheers,

--
Josue Polanco



 stephen

 On Fri, Oct 23, 2009 at 7:56 AM, Ashta sewa...@gmail.com wrote:
 Hi all,

 I am trying to do wavelets and I got  an error message saying The
 length of data is not a power of 2
 Is there a way of handing  that? or should the data length be  exactly
 the power of  2?
 I am using   R version 2.9.2 (2009-08-24)
 The is  library(wavethresh).

 wds - wd(ds$v,filter.number=1)

 Thanks

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




 --
 Stephen Sefick

 Let's not spend our time and resources thinking about things that are
 so little or so large that all they really do for us is puff us up and
 make us feel like gods.  We are mammals, and have not exhausted the
 annoying little problems of being mammals.

                                                                -K. Mullis

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Josué Mosés Polanco Martínez
Correo-e alternativo jom...@linuxmail.org

It is a wasted day unless you have learned something new and made
someone smile -Mark Weingartz.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] access elements of a named list using a factor

2009-10-23 Thread Robin Hankin

Hi

I have a factor 'f' and a named list 'jj'.

I want names(jj) to match up with levels(f).

How do I use levels(f) to access elements of jj?


 f - factor(c(pigs,pigs,slugs))
 f
[1] pigs  pigs  slugs
Levels: pigs slugs

 jj - list(pigs=1:10,slugs=1:3)



My attempts to produce jj$pigs all give errors:

 jj$levels(f)[1]
Error: attempt to apply non-function

 do.call($,jj,levels(f)[1])
Error in if (quote) { : argument is not interpretable as logical

  $(jj,levels(f)[1])
Error in jj$levels(f)[1] : invalid subscript type 'language'






--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Memory Problems with CSV and Survey Objects

2009-10-23 Thread Anthony Damico
I'm working with a 350MB CSV file on a server that has 3GB of RAM, yet I'm
hitting a memory error when I try to store the data frame into a survey
design object, the R object that stores data for complex sample survey data.

When I launch R, I execute the following line from Windows:
C:\Program Files\R\R-2.9.1\bin\Rgui.exe --max-mem-size=2047M
Anything higher, and I get an error message saying the maximum has been set
to 2047M.

Here are the commands:
 library(survey)

#this step takes more than five minutes
 data08-read.csv(data08.csv,header=TRUE,nrows=210437)

 object.size(data08)
#329877112 bytes

#Looking at Windows Task Manager, Mem Usage for Rgui.exe is already 659,632K

 brr.dsgn -svrepdesign( data = data08 , repweights = data08[, grep(
^repwgt , colnames( data08)) ], type = BRR , combined.weights = TRUE ,
weights = data08$mainwgt )
#Error: cannot allocate vector of size 254.5 Mb

#The survey design object does not get created.

#This also causes Windows Task Manager, Mem Usage to spike to 1,748,136K

#And here are some memory diagnostics
 memory.limit()
[1] 2047
 memory.size()
[1] 1449.06
 gc()
   used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells   131148   3.6 593642   15.9  15680924  418.8
Vcells 45479988 347.0  173526492 1324.0 220358611 1681.3

A description of the survey package can be found here:
http://faculty.washington.edu/tlumley/survey/

I tried creating a work-around by using the database-backed survey objects
(DB SO), included in the survey package to conserve memory on larger
datasets like this one.  Unfortunately, I don't think the survey package
supports database connections for replicate weight designs yet, since I've
only been able to get a database connection working after creating a
svydesign object and not a svrepdesign object - and also because neither the
DB SO website nor the svrepdesign help page make any mention of those
parameters.

The DB SOs are described in detail here:
http://faculty.washington.edu/tlumley/survey/svy-dbi.html

Any advice would be truly appreciated.

Thanks,
 Anthony Damico

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] access elements of a named list using a factor

2009-10-23 Thread Dimitris Rizopoulos

do you mean:

f - factor(c(pigs, pigs, slugs))
jj - list(pigs = 1:10, slugs = 1:3)

jj[levels(f)[1]]
jj[[levels(f)[1]]]


Best,
Dimitris


Robin Hankin wrote:

Hi

I have a factor 'f' and a named list 'jj'.

I want names(jj) to match up with levels(f).

How do I use levels(f) to access elements of jj?


  f - factor(c(pigs,pigs,slugs))
  f
[1] pigs  pigs  slugs
Levels: pigs slugs
 
  jj - list(pigs=1:10,slugs=1:3)



My attempts to produce jj$pigs all give errors:

  jj$levels(f)[1]
Error: attempt to apply non-function

  do.call($,jj,levels(f)[1])
Error in if (quote) { : argument is not interpretable as logical

   $(jj,levels(f)[1])
Error in jj$levels(f)[1] : invalid subscript type 'language'








--
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Inserting rows

2009-10-23 Thread Ashta
Hi all,

I have the data set  df with three varaibles,

x1 x2 x3
1  2   5
2  4   1
5  6   0
1  1   2

I want to insert more rows ( eg, 3 rows with value  filled with zeros)
1  2   5
2  4   1
5  6   6
1  1   2
0  0  0
0  0  0
0  0  0

Can any body help me out?

Thanks

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] access elements of a named list using a factor

2009-10-23 Thread Robin Hankin

Hello Dimitris

thanks for this.   It works!  I guess I was fixated on the dollar sign.

I must confess that I don't really understand any of the error
messages below.  Can anyone help me interpret them?

rksh



Dimitris Rizopoulos wrote:

do you mean:

f - factor(c(pigs, pigs, slugs))
jj - list(pigs = 1:10, slugs = 1:3)

jj[levels(f)[1]]
jj[[levels(f)[1]]]


Best,
Dimitris


Robin Hankin wrote:

Hi

I have a factor 'f' and a named list 'jj'.

I want names(jj) to match up with levels(f).

How do I use levels(f) to access elements of jj?


  f - factor(c(pigs,pigs,slugs))
  f
[1] pigs  pigs  slugs
Levels: pigs slugs
 
  jj - list(pigs=1:10,slugs=1:3)



My attempts to produce jj$pigs all give errors:

  jj$levels(f)[1]
Error: attempt to apply non-function

  do.call($,jj,levels(f)[1])
Error in if (quote) { : argument is not interpretable as logical

   $(jj,levels(f)[1])
Error in jj$levels(f)[1] : invalid subscript type 'language'











--
Robin K. S. Hankin
Uncertainty Analyst
University of Cambridge
19 Silver Street
Cambridge CB3 9EP
01223-764877

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Inserting rows

2009-10-23 Thread Henrique Dallazuanna
Try this:

rbind(df, matrix(0, 3, 3, dimnames = list(NULL, names(df


On Fri, Oct 23, 2009 at 11:44 AM, Ashta sewa...@gmail.com wrote:
 Hi all,

 I have the data set  df with three varaibles,

 x1 x2 x3
 1  2   5
 2  4   1
 5  6   0
 1  1   2

 I want to insert more rows ( eg, 3 rows with value  filled with zeros)
 1  2   5
 2  4   1
 5  6   6
 1  1   2
 0  0  0
 0  0  0
 0  0  0

 Can any body help me out?

 Thanks

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.




-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] data frame is killing me! help

2009-10-23 Thread James W. MacDonald



bbslover wrote:



Steve Lianoglou-6 wrote:

Hi,

On Oct 22, 2009, at 2:35 PM, bbslover wrote:


Usage
data(gasoline)
Format
A data frame with 60 observations on the following 2 variables.
octane
a numeric vector. The octane number.
NIR
a matrix with 401 columns. The NIR spectrum

and I see the gasoline data to see below
NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR.1694 nm NIR.1696  
nm

NIR.1698 nm NIR.1700 nm
1 1.242645 1.250789 1.246626 1.250985 1.264189 1.244678 1.245913  
1.221135
2 1.189116 1.223242 1.253306 1.282889 1.215065 1.225211 1.227985  
1.198851
3 1.198287 1.237383 1.260979 1.276677 1.218871 1.223132 1.230321  
1.208742
4 1.201066 1.233299 1.262966 1.272709 1.211068 1.215044 1.232655  
1.206696
5 1.259616 1.273713 1.296524 1.299507 1.226448 1.230718 1.232864  
1.202926
6 1.24109 1.262138 1.288401 1.291118 1.229769 1.227615 1.22763  
1.207576
7 1.245143 1.265648 1.274731 1.292441 1.218317 1.218147 1.73  
1.200446
8 1.222581 1.245782 1.26002 1.290305 1.221264 1.220265 1.227947  
1.188174
9 1.234969 1.251559 1.272416 1.287405 1.211995 1.213263 1.215883  
1.196102


look at this NIR.1686 nm NIR.1688 nm NIR.1690 nm NIR.1692 nm NIR. 
1694 nm

NIR.1696 nm NIR.1698 nm NIR.1700 nm

how can I add letters NIR to my variable, because my 600  
independents never

have NIR as the prefix. however, it is needed to model the plsr.   for
example aa=plsr(y~NIR, data=data ,), the prefix NIR is  
necessary, how

can I do with it?
I'm not really sue that I'm getting you, but if your problem is that  
the column names of your data.frame don't match the variable names  
you'd like to use in your formula, just change the colnames of your  
data.frame to match your formula.


BTW - I have no idea where to get this gasoline data set, so I'm just  
imagining:


eg.
colnames(gasoline) - c('put', 'the', 'variable', 'names', 'that',  
'you', 'want', 'here')


-steve

--
Steve Lianoglou
Graduate Student: Computational Systems Biology
   |  Memorial Sloan-Kettering Cancer Center
   |  Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




thanks for you. but the numbers of indenpendence are so many, it is not easy
to identify them one by one,  is there some better way?


You don't need to identify anything. What you need to do is read the 
help page for the function you want to use, so you (at the very least) 
know how to use the function.


 library(pls)
 data(gasoline)
 fit - plsr(octane~NIR, data=gasoline, validation = CV)
 summary(fit)
Data:   X dimension: 60 401
Y dimension: 60 1
Fit method: kernelpls
Number of components considered: 53

VALIDATION: RMSEP
Cross-validated using 10 random segments.
   (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
CV   1.5431.372   0.3827   0.2522   0.2347   0.2455   0.2281
adjCV1.5431.367   0.3740   0.2497   0.2360   0.2407   0.2243
   7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
CV  0.2311   0.2352   0.24550.25340.27370.28140.2832
adjCV   0.2257   0.2303   0.23950.24730.26460.27050.2726
   14 comps  15 comps  16 comps  17 comps  18 comps  19 comps  20 comps
CV   0.29130.29320.29850.31370.32890.33230.3391
adjCV0.28080.28210.28630.30080.31410.31720.3228
   21 comps  22 comps  23 comps  24 comps  25 comps  26 comps  27 comps
CV   0.34760.33840.33160.32130.31550.31180.3062
adjCV0.33070.32170.31540.30570.30020.29640.2908
   28 comps  29 comps  30 comps  31 comps  32 comps  33 comps  34 comps
CV   0.30330.30340.30740.30830.30940.30870.3105
adjCV0.28810.28810.29170.29260.29360.29290.2946
   35 comps  36 comps  37 comps  38 comps  39 comps  40 comps  41 comps
CV   0.31080.31060.31050.31040.31040.31050.3105
adjCV0.29490.29470.29460.29450.29450.29450.2946
   42 comps  43 comps  44 comps  45 comps  46 comps  47 comps  48 comps
CV   0.31050.31050.31050.31050.31050.31050.3105
adjCV0.29460.29460.29460.29460.29460.29460.2946
   49 comps  50 comps  51 comps  52 comps  53 comps
CV   0.31050.31050.31050.31050.3105
adjCV0.29460.29460.29460.29460.2946

TRAINING: % variance explained
1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps 
8 comps
X 70.9778.5686.15 95.496.1296.9797.32 
  98.1
octane31.9094.6697.71 98.098.6898.9399.06 
  99.1

9 

Re: [R] reset par() within plot layout

2009-10-23 Thread Tom Gottfried
Janke,

Janke ten Holt schrieb:
 Dear list,
 
 I would like to produce a matrix of plots, where par() is reset after
 each plot (see below [simplified] example). When I use layout() to do
 so, I seem to also reset the layout. I have not been able to figure out
 how to prevent this from happening.
 
 Any help is greatly appreciated!
 Janke
 
 Example code:
 #Desired result is a layout of 2 plots: one red and one black
 layout(matrix(1:2, nr=2))
 par.ini - par(no.readonly=TRUE)

look at par.ini: it's a list with all the argument-value pairs for par(). You 
might be able to solve
your problem by removing the appropriate elements from par.ini before calling 
par(par.ini). Do the
following to look which ones need to be kept for the layout:

par()
layout(matrix(1:2, nr=2))
par()

Tom

 par(col=red)
 plot(1:100)
 
 par(par.ini)
 
 plot(1:10)
 
 --
 Janke ten Holt
 Dept. of Psychology/Sociology
 University of Groningen, the Netherlands
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reset par() within plot layout

2009-10-23 Thread Janke ten Holt


Tom Gottfried wrote:
 Janke,
 
 Janke ten Holt schrieb:
 Dear list,

 I would like to produce a matrix of plots, where par() is reset after
 each plot (see below [simplified] example). When I use layout() to do
 so, I seem to also reset the layout. I have not been able to figure out
 how to prevent this from happening.

 Any help is greatly appreciated!
 Janke

 Example code:
 #Desired result is a layout of 2 plots: one red and one black
 layout(matrix(1:2, nr=2))
 par.ini - par(no.readonly=TRUE)
 
 look at par.ini: it's a list with all the argument-value pairs for par(). You 
 might be able to solve
 your problem by removing the appropriate elements from par.ini before calling 
 par(par.ini). Do the
 following to look which ones need to be kept for the layout:
 
 par()
 layout(matrix(1:2, nr=2))
 par()
 
 Tom
 
 par(col=red)
 plot(1:100)

 par(par.ini)

 plot(1:10)

 --
 Janke ten Holt
 Dept. of Psychology/Sociology
 University of Groningen, the Netherlands

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

Yes, that had occured to me too. So I tried:

layout(matrix(1:2, nr=2))
par(no.readonly=TRUE)
plot(1:10)
par(no.readonly=TRUE)

This has differences in
fig
mfg
usr
xaxp
yaxp

But even keeping these back does not solve my problem. So I figured
there must be something else going on that I am unaware of...

btw, your exact suggestion,
 par()
 layout(matrix(1:2, nr=2))
 par()
does not result in any differences.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Functional data analysis - problems with smoothing

2009-10-23 Thread Peter Ehlers

Ben,

 I lose my bet! Good thing my opinion is never worth more than $.02.

Looking at the code for smooth.spline() I see that, although 'argvals'
is the first argument to the function, the error message is generated
by the second argument. [** Spencer, is this intentional? **]

So, Ben, try str(rlz) or class(rlz) and see whether it conforms to
one of the structures described in ?smooth.basis.

 -Peter Ehlers

Benjamin Cheah wrote:

Hi Peter and Spencer (and everyone else out there),

Thank you for your prompt reply to my post re. FDA - FYI - i'm very new to R 
and FDA. Currently half way through a biostats masters degree in Sydney, 
Australia! Absolutely loving it and I hope to include R computing on my future 
CV!

Just to clarify things:

When I entered what you suggested, this is what I obtained with my own 
dataset

str(isi)

 num [1:14] 1 1.5 2 2.5 3 3.5 4 5 7 10 ...

Similarly, for Jim Ramsay's dataset (chapter five of Use R book):

str(age)

 num [1:31] 1 1.25 1.5 1.75 2 3 4 5 6 7 ...

Doesn't this mean that my data object, 'isi' is numeric?

Also, I was looking through Jim Ramsay's datasets - I'm fairly sure that my 
vector, 'isi' was similarly organised to his vector, 'age' in chapter 5 of his 
new Use R FDA book, so I'm fairly certain that the data is numeric but it 
obviously isn't

Hope this email clarified things - apologies if it hasn't - still very new to 
the terminology and overall feel of R.

The R user network is terrific! Thank you for caring!

Ben





From: Peter Ehlers ehl...@ucalgary.ca

Cc: R-help@r-project.org
Sent: Fri, 23 October, 2009 5:07:15 PM
Subject: Re: [R] Functional data analysis - problems with smoothing

The error message is pretty clear: regardless of what
*you* think, R says that 'isi' is not numeric.

Are you sure that 'isi' is not a *factor* object?
I'm willing to bet that it is.

Use str() to check your data.

-Peter Ehlers

Benjamin Cheah wrote:

Hi all,

I'm having major issues with smoothing my functional data

I'm referring to Jim Ramsay's
examples in his books. The following error message keeps appearing,
despite all my data being numeric can anyone kindly offer any suggestions?

isi - vector of argument values - i.e. the independent variable of the curves
rlz - data array
TMSfdPar - functional data parameter.


TMSfdPar = fdPar(TMSbasis, 4, 0.01)
TMSfdsmooth = smooth.basis(isi, rlz, TMSfdPar)

Error in smooth.basis(isi, rlz, TMSfdPar) : 'argvals' is not numeric.

I don't understand why the error message keeps popping up. I've tried playing 
around with Jim Ramsay's datasets and I think my data is organised in a similar 
manner to his, so can't understand what's going on oh, the frustration!

Thanks in advance,

Ben
(the amateur R data analyst and statistician.)



  
__


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.






  
__


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Inserting rows

2009-10-23 Thread John Kane
?rbind

df1 - data.frame(matrix(rep(0,9),nrow=3))
names(df1) - c(x1,x2,x3)

rbind(df,df1)

--- On Fri, 10/23/09, Ashta sewa...@gmail.com wrote:

 From: Ashta sewa...@gmail.com
 Subject: [R] Inserting rows
 To: R help r-help@r-project.org
 Received: Friday, October 23, 2009, 9:44 AM
 Hi all,
 
 I have the data set  df with three varaibles,
 
 x1 x2 x3
 1  2   5
 2  4   1
 5  6   0
 1  1   2
 
 I want to insert more rows ( eg, 3 rows with value 
 filled with zeros)
 1  2   5
 2  4   1
 5  6   6
 1  1   2
 0  0  0
 0  0  0
 0  0  0
 
 Can any body help me out?
 
 Thanks
 
 __
 R-help@r-project.org
 mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained,
 reproducible code.
 


  __
Looking for the perfect gift? Give the gift of Flickr! 

http://www.flickr.com/gift/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Change positions of columns in data frame

2009-10-23 Thread Joel Fürstenberg-Hägg

Hi all,

Probably a simple question, but I just can't find a simple answear in the older 
threads or anywhere else.

I've added some new vectors as columns in a data frame using cbind(). As 
they're all put as the last columns inte the data frame, I would like to move 
them to specific positions. How do you do to change the position of a column in 
a data frame?

I know I can use 
fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), 
Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame 
with the given columns in the specified order, but there must be an easier 
way..?

All the best,

Joel
  
_
Nya Windows 7 - Hitta en dator som passar dig! Mer information. 
http://windows.microsoft.com/shop
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bayesian regression stepwise function?

2009-10-23 Thread Ravi Varadhan
If 'fools rush in where angels fear to tread', then Bayesians 'jump' in
where frequentists fear to 'step'...

Very nice, Chuck!  Definitely one for my list of fortunes.

Ravi.


---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvarad...@jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 




-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Charles C. Berry
Sent: Friday, October 23, 2009 1:18 AM
To: Ben Bolker
Cc: r-help@r-project.org
Subject: Re: [R] Bayesian regression stepwise function?

On Thu, 22 Oct 2009, Ben Bolker wrote:




 Allan.Y wrote:

 Hi everyone,

 I am wondering if there exists a stepwise regression function for the
 Bayesian regression model.  I tried googling, but I couldn't find
 anything.  I know step function exists for regular stepwise regression,
 but nothing for Bayes.


 Why?  That seems so ... un-Bayesian ...

If 'fools rush in where angels fear to tread', then Bayesians 'jump' in 
where frequentists fear to 'step'...

Seriously, there are Bayesian regression approaches that priorize the 
model size (sometimes only implicitly by assigning a prior for the 
inclusion of each candidate regressors). Then they 'jump' between models 
of different sizes.

On CRAN, Package qtlbim (which is specialized to a particular genetics 
problem) implements one such, I think.

Package bqtl does not implement the jumping approach, but does explore a 
model space with differing numbers of regressors for the same (qtl) 
problem.

Perhaps the closest to a general purpose 'stepwise flavored' Bayesian 
regression is implemented in Package BMA, which IIRC borrows step() for 
some of its work.

But CRAN now has more packages than my cortex has neurons, so there are
probably more packages that do something like this. Try

RSiteSearch(jump regression, restric='functions')

and start reading.

HTH,

Chuck





 -- 
 View this message in context:
http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p2601
5081.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


Charles C. Berry(858) 534-2098
 Dept of Family/Preventive
Medicine
E mailto:cbe...@tajo.ucsd.edu   UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] PDF too large, PNG bad quality

2009-10-23 Thread Vijaya Parthiban
Hello Lasse,

Why not try this?

(1) Create 20MB PDF from R
(2) use convert command in linux, examples below:

convert -resize 50% 20mbfile.pdf smallerfile.pdf
convert -resize 75% 20mbfile.pdf image.png

Ghostscript can help you as well for conversion! Using vector formats
(pdf,ps,eps) are good for this purpose, as opposed to scalar (png, jpg,
etc.).

Best,
Parthiban.


2009/10/23 Jim Lemon j...@bitwrit.com.au

 On 10/23/2009 06:07 AM, Lasse Kliemann wrote:

 I wish to save a scatter plot comprising approx. 2 million points
 in order to include it in a LaTeX document.

 Using 'pdf(...)' produces a file of size about 20 MB, which is
 useless.

 Using 'cairo_pdf(...)' produces a smaller file, around 3 MB. This
 is still too large. Not only that the document will be too large,
 but also PDF viewers choke on this. Moreover, Cairo has problems
 with text: by default text looks ugly, like scaled bitmaps. After
 hours of trying different settings, I discovered that choosing a
 different font family can help, e.g.: 'par(family=Mono)'. This
 gives good-looking text. Yet, the problem with the file size
 remains.

 There exists the hint to produdc EPS instead and then convert to
 PDF using 'epstopdf'. The resulting PDF files are slightly
 smaller, but still too large, and PDF viewers still don't like
 it.

 So I gave PNG a try. PNG files are much smaller and PDF viewers
 have no trouble with them. However, fonts look ugly. The same
 trick that worked for Cairo PDF has no effect for PNG. When I
 view the PNGs with a dedicated viewer like 'qiv', even the fonts
 look good. But not when included in LaTeX; I simply use
 '\includegraphics{...}' and run the document through 'pdflatex'.

 I tried both, creating PNG with 'png(...)' and converting from
 PDF to PNG using 'convert' from ImageMagick.

 So my questions are:

 - Is there a way to produce sufficiently lean PDFs directly in R,
   even when the plot comprises several million points?

 - How to produce a PNG that still looks nice when included in a
   LaTeX PDF document?

 Any hints will be greatly appreciated.


 Hi Lasse,
 I may be right off the track, but I can't make much sense of 2 million
 points in a scatterplot. If you are interested in the density of points
 within the plot, you could compute this using something like the bkde2
 function in the KernSmooth package and then plot that using something like
 image.

 Jim

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Change positions of columns in data frame

2009-10-23 Thread John Kane
dataframe xx

x1 x2 x3
1  2   5
2  4   1
5  6   0
1  1   2 

data.frame(xx$x2,xx$x1,xx$x3)

# or


Awkward but works

--- On Fri, 10/23/09, Joel Fürstenberg-Hägg joel_furstenberg_h...@hotmail.com 
wrote:

 From: Joel Fürstenberg-Hägg joel_furstenberg_h...@hotmail.com
 Subject: [R] Change positions of columns in data frame
 To: r-help@r-project.org
 Received: Friday, October 23, 2009, 10:19 AM
 
 Hi all,
 
 Probably a simple question, but I just can't find a simple
 answear in the older threads or anywhere else.
 
 I've added some new vectors as columns in a data frame
 using cbind(). As they're all put as the last columns inte
 the data frame, I would like to move them to specific
 positions. How do you do to change the position of a column
 in a data frame?
 
 I know I can use
 fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID),
 Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create
 a new data frame with the given columns in the specified
 order, but there must be an easier way..?
 
 All the best,
 
 Joel
     
 
       
   
 _
 Nya Windows 7 - Hitta en dator som passar dig! Mer
 information. 
 http://windows.microsoft.com/shop
     [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org
 mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained,
 reproducible code.
 


  __
Be smarter than spam. See how smart SpamGuard is at giving ju
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Nested logit in GEV family

2009-10-23 Thread Min Chen
Hi All,

I know someone has posted similar messages before, but there is
no reply. So I wonder whether there is a way to run a nested logit model in
R. It is in the GEV family; however, the commands fitting GEV don't seem to
work for nested logit. Or maybe I have to write down the log-liklihood
function and let R maximize it?

Thank you!

Best wishes.


Min

-- 
Min Chen
Graduate Student
Department of Agricultural, Food, and Resource Economics
208 Cook Hall
Michigan State University
chenm...@msu.edu / chenmin0...@gmail.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to make R packages?

2009-10-23 Thread Peng Yu
These two are for windows. I have a linux machine. Are there any other
reference for linux?

On Fri, Oct 23, 2009 at 12:08 AM, Jorge Ivan Velez
jorgeivanve...@gmail.com wrote:
 Hi Peng,
 Also, take a look at the following links:
 http://robjhyndman.com/research/Rpackages_notes.pdf
 http://www.biostat.wisc.edu/~kbroman/Rintro/Rwinpack.html
 HTH,
 Jorge

 On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu  wrote:

 I found the following document on making R packages. But it is old.
 I'm wondering if there is more current ones and hopefully more
 complete ones.

 http://biosun1.harvard.edu/courses/individual/bio271/lectures/L6/Rpkg.pdf

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Documenting related data sets

2009-10-23 Thread Michael Friendly
Say I want to include the data from Minard's march on Moscow graphic in 
a package.  They consist of

three data frames, Minard.troops, Minard.cities, Minard.temp.

I've determined that I can document them all in one .Rd file, Minard.Rd, 
that starts:


\name{Minard}
\Rdversion{1.1}
\alias{Minard}
\alias{Minard.cities}
\alias{Minard.troops}
\alias{Minard.temp}
\docType{data}
\title{
Data from Minard's famous graphic map of Napoleon's march on Moscow
}
...
\usage{
data(Minard.troops)
data(Minard.cities)
data(Minard.temp)
}
...

But, AFAICS, it seems I have to save each of these in a separate .RData 
/ .rda file for them
to be found by R CMD check.  Is this correct, or is there some way to 
include them all

in a Minard.RData, via

save(Minard.troops, Minard.cities, Minard.temp, file=Minard.RData)

-Michael


--
Michael Friendly Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.

York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Change positions of columns in data frame

2009-10-23 Thread Peter Ehlers

DF[ ,names(DF)[c(your_preferred_order)]]

 -Peter Ehlers

Joel Fürstenberg-Hägg wrote:

Hi all,

Probably a simple question, but I just can't find a simple answear in the older 
threads or anywhere else.

I've added some new vectors as columns in a data frame using cbind(). As 
they're all put as the last columns inte the data frame, I would like to move 
them to specific positions. How do you do to change the position of a column in 
a data frame?

I know I can use 
fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), 
Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame 
with the given columns in the specified order, but there must be an easier 
way..?

All the best,

Joel
 		 	   		  
_
Nya Windows 7 - Hitta en dator som passar dig! Mer information. 
http://windows.microsoft.com/shop

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] standard errors in bbmle

2009-10-23 Thread Ben Bolker



Peter Dalgaard wrote:
 
 alexander russell wrote:
 Hello,
 Mle2 is a little unforthcoming in the matter of standard errors? Is there
 a
 way to ask the program to supply standard errors along with estimates in
 cases when it doesn't print them 'voluntarily'?
 regards,
 s
 
 What did you do? Which cases? Did you look at the examples? (as in,
 example(mle2) )
 
 

As the package author:  +1.
Perhaps you're looking for summary()?

  Ben Bolker

-- 
View this message in context: 
http://www.nabble.com/standard-errors-in-bbmle-tp26023856p26028127.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] removing random effect from nlme or using varPower() in nls

2009-10-23 Thread Michael A. Gilchrist

Hi Kingsford,

That's exactly what I was looking for!

Thanks.

Mike


On Thu, 22 Oct 2009, Kingsford Jones wrote:


help(gnls, pack=nlme)


hth,
Kingsford Jones

On Thu, Oct 22, 2009 at 4:36 PM, Michael A. Gilchrist mi...@utk.edu wrote:

Hello,

I've been fitting a random effects model using nlme to some data, but I am
discovering that the variation in my random effect is very small.  As a
result, I would like to replace it as a fixed effect (i.e. essentially fit
the same model but with no random effect).

As I understand it I could do this using nls(), but I'm using a number of
options such as weights = varPower() which I am at a loss on how to
implement in that framework.

Is there a way to use nlme but with out a random effect? (a bit absurd, I
realize, but I have the syntax working...)

Alternatively, is there a way to use weights = varPower() with nls?

Any help would be appreciated.

Sincerely,

Mike

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] cut2 once, bin twice...

2009-10-23 Thread Gabor Grothendieck
On Fri, Oct 23, 2009 at 3:58 AM, Dieter Menne
dieter.me...@menne-biomed.de wrote:



 sdanzige wrote:


 I'm using the Hmisc cut2 function to bin a set of data.  It produces bins
 that I like with results like this:

 [96,270]:171
 [69, 96): 54
 [49, 69): 40
 [35, 49): 28
 [28, 35): 14
 [24, 28):  8
 (Other) : 48

 I would like to take a second set of data, and assign it to bins based on
 factors defined by my call to cut 2.


 It used to be quite tricky, but on popular request Brian Ripley has added an
 example how to extract the intervals using regular expression on the bottom
 of the examples for cut (note:cut in base, not cut2 in Hmisc).

 If someone knows of an easier way, please correct me. How about adding this
 information as attribute to the standard cut?


The strapply function in gsubfn can do it with a simpler regular
expression since it extracts based on content rather than delimiters,
which is what you want here:

 # create sample data
 library(gsubfn)
 set.seed(1)
 dat - seq(4, 7, by = 0.05)
 x - sample(dat, 30)
.
 # use cut
 groups - cut(x, breaks = 10)

 # extract interval boundaries using strapply
 strapply(levels(groups), [[:digit:].]+, as.numeric, simplify = TRUE)
 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]  4.0  4.3  4.6  4.9  5.2  5.5  5.8  6.1  6.4   6.7
[2,]  4.3  4.6  4.9  5.2  5.5  5.8  6.1  6.4  6.7   7.0

The above is from

   demo(gsubfn-cut)

For more see the gsubfn home page at http://gsubfn.googlecode.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to find moving averages within each subgroup of a data frame

2009-10-23 Thread Gabor Grothendieck
Try this:

df1$SMA - ave(df1$col4, df1$col1, FUN = function(x) c(NA, (head(x,
-1) + tail(x, -1))/2))

It would also be possible to convert it from long form to wide form
using reshape (or read.zoo in the devel version of zoo), convert that
to a zoo series and the use rollapply in the zoo package.


On Thu, Oct 22, 2009 at 12:28 PM, clair.crossup...@googlemail.com
clair.crossup...@googlemail.com wrote:
 Dear all,

 If I have the following data frame:

 set.seed(21)
 df1 - data.frame(col1=c(rep('a',5), rep('b',5), rep('c',5)), 
 col4=rnorm(1:15))

   col1         col4
 1     a  0.793013171
 2     a  0.522251264
 3     a  1.74641
 4     a -1.271336123
 5     a  2.197389533
 6     b  0.433130777
 7     b -1.570199630
 8     b -0.934905667
 9     b  0.063493345
 10    b -0.002393336
 11    c -2.276781240
 12    c  0.757412225
 13    c -0.548405554
 14    c  0.172549478
 15    c  0.562853068

 How can i make a 2 point moving average within each group? i.e.

 col1    col4    SMA
 a       0.793013171     NA
 a       0.522251264     0.657632218
 a       1.74641     1.134236753
 a       -1.271336123    0.237443059
 a       2.197389533     0.463026705
 b       0.433130777     NA
 b       -1.57019963     -0.568534427
 b       -0.934905667    -1.252552649
 b       0.063493345     -0.435706161
 b       -0.002393336    0.030550005
 c       -2.27678124     NA
 c       0.757412225     -0.759684508
 c       -0.548405554    0.104503336
 c       0.172549478     -0.187928038
 c       0.562853068     0.367701273

 From what i've read, i was thinking it would be something along the
 lines of:

 aggregate(df1$col4, by=list(df1$col1), function(x) {filter(x, rep(1/2,2), 
 sides=1 )} )
 Error in aggregate.data.frame(as.data.frame(x), ...) :
  'FUN' must always return a scalar

 But this tells me (i think) that aggregate should only return a single
 value per group. So what i need, i guess, is something that takes all
 the values in a given group, and returns a vector of the same length.
 Not sure which function to use for that.

 Thanks in advance,

 C xx

 P.S. on a side note, is it possible to extract the values which
 aggregate passes to the function(x) in my example above?

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reset par() within plot layout

2009-10-23 Thread Peter Ehlers



Janke ten Holt wrote:


Tom Gottfried wrote:

Janke,

Janke ten Holt schrieb:

Dear list,

I would like to produce a matrix of plots, where par() is reset after
each plot (see below [simplified] example). When I use layout() to do
so, I seem to also reset the layout. I have not been able to figure out
how to prevent this from happening.

Any help is greatly appreciated!
Janke

Example code:
#Desired result is a layout of 2 plots: one red and one black
layout(matrix(1:2, nr=2))
par.ini - par(no.readonly=TRUE)

look at par.ini: it's a list with all the argument-value pairs for par(). You 
might be able to solve
your problem by removing the appropriate elements from par.ini before calling 
par(par.ini). Do the
following to look which ones need to be kept for the layout:

par()
layout(matrix(1:2, nr=2))
par()

Tom


par(col=red)
plot(1:100)

par(par.ini)

plot(1:10)

--
Janke ten Holt
Dept. of Psychology/Sociology
University of Groningen, the Netherlands

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Yes, that had occured to me too. So I tried:

layout(matrix(1:2, nr=2))
par(no.readonly=TRUE)
plot(1:10)
par(no.readonly=TRUE)

This has differences in
fig
mfg
usr
xaxp
yaxp

But even keeping these back does not solve my problem. So I figured
there must be something else going on that I am unaware of...



I think that what you want is

 layout(matrix(1:2, nr=2))
 opar - par(col=red)
 plot(1:100)
 par(opar)
 plot(1:10)

-Peter Ehlers


btw, your exact suggestion,

par()
layout(matrix(1:2, nr=2))
par()

does not result in any differences.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to make R packages?

2009-10-23 Thread Steve Lianoglou

Hi Pen,

On Oct 23, 2009, at 10:31 AM, Peng Yu wrote:


These two are for windows. I have a linux machine. Are there any other
reference for linux?


The steps are pretty much the same. For instance, in robjhyndman.com  
link, you can simply start from step 2 and go on.


Also, the official Writing R Extensions covers this topic pretty  
thoroughly:


http://cran.r-project.org/doc/manuals/R-exts.html#Creating-R-packages

-steve



On Fri, Oct 23, 2009 at 12:08 AM, Jorge Ivan Velez
jorgeivanve...@gmail.com wrote:

Hi Peng,
Also, take a look at the following links:
http://robjhyndman.com/research/Rpackages_notes.pdf
http://www.biostat.wisc.edu/~kbroman/Rintro/Rwinpack.html
HTH,
Jorge

On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu  wrote:


I found the following document on making R packages. But it is old.
I'm wondering if there is more current ones and hopefully more
complete ones.

http://biosun1.harvard.edu/courses/individual/bio271/lectures/L6/Rpkg.pdf

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.





__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


--
Steve Lianoglou
Graduate Student: Computational Systems Biology
  |  Memorial Sloan-Kettering Cancer Center
  |  Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] reference for mahalanobis {stats}

2009-10-23 Thread Peng Yu
The help on mahalanobis {stats} does not include any reference. I'm
interested in understand why Mahalanobis is defined in its current way
and how to use it. Could somebody point me a good book on this? I have
looked through a few books, but they all give very light explanation
on it.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] cut2 once, bin twice...

2009-10-23 Thread sdanzige


Dieter Menne wrote:
 
 
 It used to be quite tricky, but on popular request Brian Ripley has added
 an example how to extract the intervals using regular expression on the
 bottom of the examples for cut (note:cut in base, not cut2 in Hmisc).
 
 


Thank you, but the regular expression example doesn't seem to work
correctly.

 labs-levels(df$p_bin)

 labs
 [1]  012345  
 [7]  6789   10   11  
[13] 12   13   14   15   16   17  
[19] 18   19   20   [21, 24) [24, 28) [28, 35)
[25] [35, 49) [49, 69) [69, 96) [96,270]

 cbind(lower = as.numeric( sub(\\((.+),.*, \\1, labs) ), upper =
 as.numeric( sub([^,]*,([^]]*)\\], \\1, labs) ))
Warning in cbind(lower = as.numeric(sub(\\((.+),.*, \\1, labs)), upper =
as.numeric(sub([^,]*,([^]]*)\\],  :
  NAs introduced by coercion
Warning in cbind(lower = as.numeric(sub(\\((.+),.*, \\1, labs)), upper =
as.numeric(sub([^,]*,([^]]*)\\],  :
  NAs introduced by coercion
  lower upper
 [1,] 0 0
 [2,] 1 1
 [3,] 2 2
 [4,] 3 3
 [5,] 4 4
 [6,] 5 5
 [7,] 6 6
 [8,] 7 7
 [9,] 8 8
[10,] 9 9
[11,]1010
[12,]1111
[13,]1212
[14,]1313
[15,]1414
[16,]1515
[17,]1616
[18,]1717
[19,]1818
[20,]1919
[21,]2020
[22,]NANA
[23,]NANA
[24,]NANA
[25,]NANA
[26,]NANA
[27,]NANA
[28,]NA   270

--

Any ideas?

Thank you,
-S
-- 
View this message in context: 
http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26027643.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Change positions of columns in data frame

2009-10-23 Thread Phil Spector

Joel -
   Suppose the columns are named x1, x2, x3, x4, and x5.

You can use subscripting:

 x[c('x2','x4','x1','x3','x5')]

or, to save typing the quotes

  subset(x,select=c(x2,x4,x1,x3,x5))


- Phil Spector
 Statistical Computing Facility
 Department of Statistics
 UC Berkeley
 spec...@stat.berkeley.edu



On Fri, 23 Oct 2009, Joel Fürstenberg-Hägg wrote:



Hi all,

Probably a simple question, but I just can't find a simple answear in the older 
threads or anywhere else.

I've added some new vectors as columns in a data frame using cbind(). As 
they're all put as the last columns inte the data frame, I would like to move 
them to specific positions. How do you do to change the position of a column in 
a data frame?

I know I can use 
fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), 
Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame 
with the given columns in the specified order, but there must be an easier 
way..?

All the best,

Joel

_
Nya Windows 7 - Hitta en dator som passar dig! Mer information.
http://windows.microsoft.com/shop
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] cut2 once, bin twice...

2009-10-23 Thread sdanzige


sdanzige wrote:
 
 
 Thank you, but the regular expression example doesn't seem to work
 correctly.
 
 

I wrote a regular expression that does seem to work, so I'll post it here
for anyone else that needs it.

labs-levels(df$p_bin)
cbind(lower=as.numeric(sub([[(],,sub(,.*,,labs))),
upper=as.numeric(sub([])],,sub([[(].*, *,,labs))) )


I fear my inelegance will peg me as a Windows programmer, but so be it... 
-S
-- 
View this message in context: 
http://www.nabble.com/cut2-once%2C-bin-twice...-tp26020736p26028296.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reference for mahalanobis {stats}

2009-10-23 Thread Steve Lianoglou

Hi,

On Oct 23, 2009, at 11:36 AM, Peng Yu wrote:


The help on mahalanobis {stats} does not include any reference. I'm
interested in understand why Mahalanobis is defined in its current way
and how to use it. Could somebody point me a good book on this? I have
looked through a few books, but they all give very light explanation
on it.


The wikipedia page looks pretty good ... the Intuitive Explanation  
section should help:


http://en.wikipedia.org/wiki/Mahalanobis_distance

Googling also provides many more explanations:

http://matlabdatamining.blogspot.com/2006/11/mahalanobis-distance.html

HTH,

-steve

--
Steve Lianoglou
Graduate Student: Computational Systems Biology
  |  Memorial Sloan-Kettering Cancer Center
  |  Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Determining which file is newer

2009-10-23 Thread Dennis Fisher

Colleagues,

I wish to execute a task only if a particular file is newer than a  
second file.  I can access the file modification dates using  
file.info()$mtime.

This yields a time object (? POSIX).
sizeisdir   modemtime   
ctime
Testfile4421FALSE   755 2009-10-23 08:59:09 
2009-10-23 08:59:09
What is the simplest means to compare these two objects?

R2.9.1 in OSX, Windows, Linux

Dennis

Dennis Fisher MD
P  (The P Less Than Company)
Phone: 1-866-PLessThan (1-866-753-7784)
Fax: 1-866-PLessThan (1-866-753-7784)
www.PLessThan.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] opposite estimates from zeroinfl() and hurdle()

2009-10-23 Thread Ben Bolker



Tord Snäll-4 wrote:
 
 Dear all,
 A question related to the following has been asked on R-help before, but 
 I could not find any answer to it. Input will be much appreciated.
 
 I got an unexpected sign of the slope parameter associated with a 
 covariate (diam) using zeroinfl(). It led me to compare the estimates 
 given by zeroinfl() and hurdle():
 [snip]
 

The right thing to do in this case is to poke through the code of hurdle()
and zeroinfl(), but a simple (?) demonstration shows that hurdle()
and zeroinfl() are indeed reporting opposite values :

 hurdle reports  -log(p/(1-p)) = -qlogis(p), where p is the probability
of a zero count:

z = rpois(500,lambda=3)
z = (z[z0])[1:90]
z = c(z,rep(0,10))
hurdle(z~1)  ##
-qlogis(0.1)
## zero coefficient always == -qlogis(0.1)

  zeroinfl reports  log(p/(1-p)), where p is the
zero-inflation:

z = rpois(90,lambda=3)
z = c(z,rep(0,10))
zeroinfl(z~1)  ##
qlogis(0.1)

tmpf = function() {
  z = rpois(90,lambda=3)
  z = c(z,rep(0,10))
  coef(zeroinfl(z~1))[2]
}

rr = replicate(1000,tmpf())

hist(rr,breaks=1000)
summary(rr)
qlogis(0.1)

  Perhaps it would be worth sending an e-mail to the
package maintainers to request a note to this effect in
the documentation, particularly if this a FAQ ...
-- 
View this message in context: 
http://www.nabble.com/opposite-estimates-from-zeroinfl%28%29-and-hurdle%28%29-tp26024735p26029131.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima crashes too

2009-10-23 Thread Alberto Monteiro

Barry Rowlingson wrote:
 
  If you're doing anything in a loop that has the potential to fail
 because of singularities or other conditions when your model can't be
 fitted, you need to stick what you are doing in a 'try' clause. This
 lets you trap errors and do something with them.
 
  Plenty of examples in help(try) or this from me:
 
  for(i in 1:10){
  print(solve(matrix(c(3,3,3,i),2,2)))
  }
 
  This stops the loop at i=3. Now stick it in a try() clause:
 
  for(i in 1:10){
  print(try(solve(matrix(c(3,3,3,i),2,2
  }
 
  and it gives a warning and carries on. If you want your code to do
 something with the failure cases then the help for try() tells you
 what to look for.
 
I will try try.  # irreristible pun

But...

  I'm not sure why your arima produces an error, but I'm assuming the
 numbers are such that the model can't be fitted. I don't really know
 what arima is doing.
 
That's a problem. In this case, I'm fitting an AR(1) time series 
using arima and arma. The documentation does not mention that we 
should  be careful about what numbers are passed to it; on contrary, 
it's perfectly possible (in theory) that an AR(1) time series was
created using an evil phi in 

x[i+1] - phi * x[i] + sigma * n01[i]

and the regression should return the evil phi.

As in this reproducible example:

phi - 1.5
sigma - 0.5
n01 - c(-1, 0.3, -0.3, 0.2, -0.7)
x - rep(0, 6)
for (i in 1:5) x[i+1] - phi * x[i] + sigma * n01[i]
arma(x, order=c(1,0))
arima(x, order=c(1,0,0))

It seems like arma does not force phi to be in the -1:1 range, because
it correctly returns evil phi's, but arima kind of bayesianises
the output, forcing the time series to be stationary.

Probably a longer example may scrash/s stop arima but not arma.

Non-reproducible example (almost surely scrashes/s stops arima but
plays safely in arma):

phi - 1.05
sigma - 0.5
n - 100
n01 - rnorm(n)
x - rep(0, n)
for (i in 1:n) x[i+1] - phi * x[i] + sigma * n01[i]
arma(x, order=c(1,0))
arima(x, order=c(1,0,0))

Increasing phi (or n) will scrash/s stop both arma and arima:


phi - 1.1
sigma - 0.5
n - 100
n01 - rnorm(n)
x - rep(0, n)
for (i in 1:n) x[i+1] - phi * x[i] + sigma * n01[i]
arma(x, order=c(1,0))
arima(x, order=c(1,0,0))

Alberto 'Darth Albmont' Monteiro

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Determining which file is newer

2009-10-23 Thread Henrik Bengtsson
Have a look at your company name and your phone number.  Nudge nudge,
wink wink, say no more... /H

On Fri, Oct 23, 2009 at 9:02 AM, Dennis Fisher fis...@plessthan.com wrote:
 Colleagues,

 I wish to execute a task only if a particular file is newer than a second
 file.  I can access the file modification dates using file.info()$mtime.
 This yields a time object (? POSIX).
                        size    isdir   mode    mtime
           ctime
        Testfile        4421    FALSE   755     2009-10-23 08:59:09
   2009-10-23 08:59:09
 What is the simplest means to compare these two objects?

 R2.9.1 in OSX, Windows, Linux

 Dennis

 Dennis Fisher MD
 P  (The P Less Than Company)
 Phone: 1-866-PLessThan (1-866-753-7784)
 Fax: 1-866-PLessThan (1-866-753-7784)
 www.PLessThan.com

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Determining which file is newer

2009-10-23 Thread Barry Rowlingson
On Fri, Oct 23, 2009 at 5:02 PM, Dennis Fisher fis...@plessthan.com wrote:
 Colleagues,

 I wish to execute a task only if a particular file is newer than a second
 file.  I can access the file modification dates using file.info()$mtime.
 This yields a time object (? POSIX).
                        size    isdir   mode    mtime
           ctime
        Testfile        4421    FALSE   755     2009-10-23 08:59:09
   2009-10-23 08:59:09
 What is the simplest means to compare these two objects?

 Split the string representation of the object into date and time
parts, then split those parts into day, month , year and hour minute
second, taking care of internationalisation and localisation, then
recombine all those parts into a number of seconds since Jan 1 1970,
taking care of leap seconds and the slowing down of the earth due to
the tides of the moon.

 Or just use the usual operators:

  file.info(/etc/passwd)$mtime = file.info(/etc/motd)$mtime
[1] TRUE

Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Change positions of columns in data frame

2009-10-23 Thread Tom Gottfried
You can change the order of nearly everything just by applying an appropriate 
index and assigning to
the original object. Something like

data - data[,order_of_columns]

Tom

Joel Fürstenberg-Hägg schrieb:
 Hi all,
 
 Probably a simple question, but I just can't find a simple answear in the 
 older threads or anywhere else.
 
 I've added some new vectors as columns in a data frame using cbind(). As 
 they're all put as the last columns inte the data frame, I would like to move 
 them to specific positions. How do you do to change the position of a column 
 in a data frame?
 
 I know I can use 
 fieldTrial0809=data.frame(Sample_ID=as.factor(fieldTrial0809$Sample_ID), 
 Plant_ID=as.factor(fieldTrial0809$Plant_ID), ...) to create a new data frame 
 with the given columns in the specified order, but there must be an easier 
 way..?
 
 All the best,
 
 Joel
 
 _
 Nya Windows 7 - Hitta en dator som passar dig! Mer information. 
 http://windows.microsoft.com/shop
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

-- 
Technische Universität München
Department für Pflanzenwissenschaften
Lehrstuhl für Grünlandlehre
Am Hochanger 1
85350 Freising / Germany
Phone: ++49 (0)8161 715324
Fax:   ++49 (0)8161 713243
email: tom.gottfr...@wzw.tum.de
http://www.wzw.tum.de/gruenland

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to make R packages?

2009-10-23 Thread David M Smith
On Fri, Oct 23, 2009 at 12:10 AM, Peng Yu pengyu...@gmail.com wrote:
 I found the following document on making R packages. But it is old.
 I'm wondering if there is more current ones and hopefully more
 complete ones.

Friedrich Leisch made an excellent tutorial for making R packages,
with a tutorial that works in Linux. You can find it here;

http://epub.ub.uni-muenchen.de/6175/

with some additional commentary here:

http://blog.revolution-computing.com/2009/08/creating-r-packages-a-tutorial-draft.html

# David Smith

-- 
David M Smith da...@revolution-computing.com
VP of Community, REvolution Computing  http://blog.revolution-computing.com
Tel: +1 (206) 577-4778 x3203 (Palo Alto, CA, USA)

Download REvolution R free:
http://revolution-computing.com/downloads/revolution-r.php

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Capturing lagged values?

2009-10-23 Thread Arthur Burke
Consider a data structure where header rows supply ID values (school
and grade in the example below) that are not captured in subsequent rows

 
df - read.table(textConnection(
school grade studentid score
101 5 NA NA
NA NA 123 21
NA NA 124 25 
201 6 NA NA
NA NA 231 33), header=TRUE)
closeAllConnections()
 
How do I capture the missing values for school and grade?
 
Art

*

Art Burke

Associate, Evaluation Program

Education Northwest

101 SW Main St, Ste 500

Portland, OR 97204

Phone: 503.275.9592

art.bu...@educationnorthwest.org
mailto:art.bu...@educationnorthwest.org 

http://educationnorthwest.org http://educationnorthwest.org/ 

We have recently changed our name to Education Northwest from
Northwest Regional Educational Laboratory. Please note the new e-mail
and Web addresses in the signature above. You may continue to find us on
the Web at http://www.nwrel.org http://www.nwrel.org/  for the
immediate future as well.




 
 

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] coxph() and survfit()

2009-10-23 Thread tlumley

On Thu, 22 Oct 2009, Mura Tamakou wrote:


Dear All,

I have a question regarding the output of survfit() when I supply a Cox model. 
Lets say for example:


snipped: code that doesn't quite run


we observe that the 95% CIs overlap!! How is this possible since the HR for 
spiders is significant.



It's perfectly natural.

To a good approximation, the p-value for the comparison will be significant 
when the point estimate for each group is outside the confidence interval for 
the other group.

Suppose you had two point estimates with standard error equal to 1.0.  The 
standard error of the difference would be 1.414, so the p-value would be less 
than 0.05 if the two point estimates differ by more  than 1.96*1.414.  The 
confidence intervals will overlap if the point estimates differ by less than 
1.96*2.

-thomas

Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] ggplot2: stat_bin ..count.. with geom_text when NA is present

2009-10-23 Thread Bryan Hanson
One for the ggplot2 gurus...

I have a function which makes a plot just fine if the response vector (res
in the example; fac1 is a factor) has no NA in it.  It plots the data, then
makes a little annotation at the bottom with the data counts using:

p - p + geom_text(aes(x = fac1, y = min(res) - 0.1 * diff(range(res)),
label = paste(n = , ..count.. , sep = )),
color = black, size = 4.0, stat = bin)

If there are NA in the res vector, I get warnings from stat_summary and
geom_point about removing rows; these arise from an earlier part of the
function and the points and error bars all plot.  However, the count
annotation does not appear on the plot when there are NA in res.

Looking at the ggplot2 web site, there is a drop parameter for stat_bin. I
inserted drop = TRUE several places in the snippet above and the function
did not complain but still did not plot the counts.  I looked at the
function bin{ggplot2} which apparently does the work.  There are some
programming tricks there I'm not really familiar with, but generally it
looks like it na.rm or na.omit's in several places, while the drop = TRUE is
carried out as the last step.

So, any suggestions about why the counts don't appear on my plot?  I suppose
I can always clean the data first, but it would be much more practical to do
that in the background during the preparation of the plot.

Thanks as always, Bryan
*
Bryan Hanson
Acting Chair
Professor of Chemistry  Biochemistry
DePauw University, Greencastle IN USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Memory Problems with CSV and Survey Objects

2009-10-23 Thread tlumley



Yes, a 350Mb data frame is a bit big for 32-bit R to handle conveniently.

As you note, the survey package doesn't yet do database-backed replicate-weight 
designs. You can get the same effect yourself without too much work.

First, put the data into a database, such as SQLite.  If you have the data 
frame read in then dbWriteTable will do it.

Now, drop most of the variables, keeping the sampling weights, replicate 
weights, and a couple of other variables.

Create a svrepdesign() with the reduced data set.

When you want to do an analysis, use dbGetQuery() to load the variables you 
need for the analysis, and put them in the $variables component of the 
svrepdesign.

That's exactly what the database-backed functions do for svydesign objects.

[If you only ever want to use a small subset of the variables, it's even 
easier: drop all the extraneous variables and create a svrepdesign with the 
variables you want]

   -thomas

On Fri, 23 Oct 2009, Anthony Damico wrote:


I'm working with a 350MB CSV file on a server that has 3GB of RAM, yet I'm
hitting a memory error when I try to store the data frame into a survey
design object, the R object that stores data for complex sample survey data.

When I launch R, I execute the following line from Windows:
C:\Program Files\R\R-2.9.1\bin\Rgui.exe --max-mem-size=2047M
Anything higher, and I get an error message saying the maximum has been set
to 2047M.

Here are the commands:

library(survey)


#this step takes more than five minutes

data08-read.csv(data08.csv,header=TRUE,nrows=210437)



object.size(data08)

#329877112 bytes

#Looking at Windows Task Manager, Mem Usage for Rgui.exe is already 659,632K


brr.dsgn -svrepdesign( data = data08 , repweights = data08[, grep(

^repwgt , colnames( data08)) ], type = BRR , combined.weights = TRUE ,
weights = data08$mainwgt )
#Error: cannot allocate vector of size 254.5 Mb

#The survey design object does not get created.

#This also causes Windows Task Manager, Mem Usage to spike to 1,748,136K

#And here are some memory diagnostics

memory.limit()

[1] 2047

memory.size()

[1] 1449.06

gc()

  used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells   131148   3.6 593642   15.9  15680924  418.8
Vcells 45479988 347.0  173526492 1324.0 220358611 1681.3

A description of the survey package can be found here:
http://faculty.washington.edu/tlumley/survey/

I tried creating a work-around by using the database-backed survey objects
(DB SO), included in the survey package to conserve memory on larger
datasets like this one.  Unfortunately, I don't think the survey package
supports database connections for replicate weight designs yet, since I've
only been able to get a database connection working after creating a
svydesign object and not a svrepdesign object - and also because neither the
DB SO website nor the svrepdesign help page make any mention of those
parameters.

The DB SOs are described in detail here:
http://faculty.washington.edu/tlumley/survey/svy-dbi.html

Any advice would be truly appreciated.

Thanks,
Anthony Damico

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] how do you know which functions are being debugged?

2009-10-23 Thread Andrew Yee
This is kind of a dumb question:  I know you can use isdebugged() to find
out if a specific function is flagged for debugging, but is there a way to
list all the functions that are flagged for debugging?
Thanks,
Andrew

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Bonferroni with unequal sample sizes

2009-10-23 Thread slreiner
Hello-
   I have run an ANOVA on 4 treatments with unequal sample sizes (n=9,7,10 and 
10).  I want to determine where my sig. differences are between treatments 
using a Bonferroni test, and have run the code:

pairwise.t.test(Wk16, Treatment, p.adf=bonf)

I receive an error message stating that my arguments are of unequal length:

Error in tapply(x, g, mean, na.rm = TRUE) : 
  arguments must have same length

Is there a way to run this test even with unequal sample sizes?

Thanks.
Stephanie Reiner
Andrews Hall 410
Shipping:  Rt. 1208 Greate Road 
Gloucester Point, VA 23062
work: 804-684-7869
cell: 443-286-4795

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bonferroni with unequal sample sizes

2009-10-23 Thread Chuck Cleland
On 10/23/2009 1:30 PM, slrei...@vims.edu wrote:
 Hello-
I have run an ANOVA on 4 treatments with unequal sample sizes (n=9,7,10 
 and 10).  I want to determine where my sig. differences are between 
 treatments using a Bonferroni test, and have run the code:
 
 pairwise.t.test(Wk16, Treatment, p.adf=bonf)
 
 I receive an error message stating that my arguments are of unequal length:
 
 Error in tapply(x, g, mean, na.rm = TRUE) : 
   arguments must have same length
 
 Is there a way to run this test even with unequal sample sizes?

  Unequal sample sizes are not causing the reported error.  The problem
is that 'Wk16' and 'Treatment' do not have the same number of
observations.  If the group sizes are 9, 7, 10, and 10, then 'Wk16' and
'Treatment' should each contain 36 observations.

 Thanks.
 Stephanie Reiner
 Andrews Hall 410
 Shipping:  Rt. 1208 Greate Road 
 Gloucester Point, VA 23062
 work: 804-684-7869
 cell: 443-286-4795
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code. 

-- 
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Yihui Xie
list all objects first; use a loop (explicitly or not) to check
whether (1) your objects are functions (2) functions are debugged

 f = function(x) x
 g = 1
 x = ls()
 debug(f)
 sapply(x[sapply(x, function(i) is.function(get(i)))], isdebugged)
   f
TRUE
 undebug(f)
 sapply(x[sapply(x, function(i) is.function(get(i)))], isdebugged)
f
FALSE

Regards,
Yihui
--
Yihui Xie xieyi...@gmail.com
Phone: 515-294-6609 Web: http://yihui.name
Department of Statistics, Iowa State University
3211 Snedecor Hall, Ames, IA



On Fri, Oct 23, 2009 at 12:28 PM, Andrew Yee y...@post.harvard.edu wrote:
 This is kind of a dumb question:  I know you can use isdebugged() to find
 out if a specific function is flagged for debugging, but is there a way to
 list all the functions that are flagged for debugging?
 Thanks,
 Andrew

        [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Duncan Murdoch

On 10/23/2009 1:28 PM, Andrew Yee wrote:

This is kind of a dumb question:  I know you can use isdebugged() to find
out if a specific function is flagged for debugging, but is there a way to
list all the functions that are flagged for debugging?


No, R doesn't keep any master list, it sets a flag in each one.  So you 
could iterate over all visible objects to find the ones that have the 
flag set (and this is quite slow, there are a lot of places to look), 
but there would always be the chance that one was hiding somewhere you 
didn't look.


Duncan Murdoch

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Yihui Xie
Oops... I forgot to mention that 'envir' (or 'pos') should be
specified in ls()/get() in my last reply if you are looking for
debugged functions in environments other than .GlobalEnv.

Regards,
Yihui
--
Yihui Xie xieyi...@gmail.com
Phone: 515-294-6609 Web: http://yihui.name
Department of Statistics, Iowa State University
3211 Snedecor Hall, Ames, IA



On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote:
 On 10/23/2009 1:28 PM, Andrew Yee wrote:

 This is kind of a dumb question:  I know you can use isdebugged() to find
 out if a specific function is flagged for debugging, but is there a way to
 list all the functions that are flagged for debugging?

 No, R doesn't keep any master list, it sets a flag in each one.  So you
 could iterate over all visible objects to find the ones that have the flag
 set (and this is quite slow, there are a lot of places to look), but there
 would always be the chance that one was hiding somewhere you didn't look.

 Duncan Murdoch

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Gabor Grothendieck
I have often wanted to get such a list as well without having to write
code myself.  If  a function which automatically iterated over all
possible objects and listed out the ones that being debugged were
available in the core it would be useful.  Even if its slow this would
be used at debug time and not a production time so it would not be too
bad.  Even a heuristic that checked likely locations, e.g. just the
global environment, but not all locations might be useful.

On Fri, Oct 23, 2009 at 2:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote:
 On 10/23/2009 1:28 PM, Andrew Yee wrote:

 This is kind of a dumb question:  I know you can use isdebugged() to find
 out if a specific function is flagged for debugging, but is there a way to
 list all the functions that are flagged for debugging?

 No, R doesn't keep any master list, it sets a flag in each one.  So you
 could iterate over all visible objects to find the ones that have the flag
 set (and this is quite slow, there are a lot of places to look), but there
 would always be the chance that one was hiding somewhere you didn't look.

 Duncan Murdoch

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bayesian regression stepwise function?

2009-10-23 Thread JLucke
The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's 
not stepwise.   Ntzoufras shows how to use WinBUGS to conduct Bayesian 
model selection, but again it's not stepwise


Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of 
Statistical Software 7(7), 1--19.
Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ.
Raftery, A. E. (1995), 'Bayesian model selection in social research', 
Sociological Methodology 25, 111-163.







Allan.Y all...@cmu.edu 
Sent by: r-help-boun...@r-project.org
10/22/2009 01:09 PM

To
r-help@r-project.org
cc

Subject
[R]  Bayesian regression stepwise function?







Hi everyone,

I am wondering if there exists a stepwise regression function for the
Bayesian regression model.  I tried googling, but I couldn't find 
anything. 
I know step function exists for regular stepwise regression, but nothing
for Bayes.


Thanks
-- 
View this message in context: 
http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p26013725.html

Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Prediction Error Calculation

2009-10-23 Thread quaildoc

Hello List,

I am fitting a logistic regression model for some presence/absence type
data.  I have numerous covariates I am fitting to explain variation, and I
am using AIC to rank models.  However, I would like to report how well my
best model (s) do at prediction.  I have looked over the archives and the
web and have come up with something that gives me what I think is the mean
prediction error, BUT I am not sure of that. I am sort of unfamiliar with
these types of statistics.  Here is my code:


metrics.global-glm(Type~MPI+IJI+ED+PRD+class2+class3+class5,
family=binomial, data=metrics)## ##Type is my binary response 0 or 1

muhat-metrics.global$fitted.values
##assigns the fitted values a name muhat
global.diag-glm.diag(metrics.global)
##creates a the diagnostic values
cv.err-mean((metrics.global$y-muhat)^2/(1-global.diag$h)^2)
###calculates cv.err
cv.err


My main problem is I am unsure how to interpret what cv.err means for my
model.  I know that h is a leverage statistic for each observation.  I would
appreciate some interpretation clarification.

Thank you.




-- 
View this message in context: 
http://www.nabble.com/Prediction-Error-Calculation-tp26031236p26031236.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Duncan Murdoch

On 10/23/2009 2:27 PM, Gabor Grothendieck wrote:

I have often wanted to get such a list as well without having to write
code myself.  If  a function which automatically iterated over all
possible objects and listed out the ones that being debugged were
available in the core it would be useful.  Even if its slow this would
be used at debug time and not a production time so it would not be too
bad.  Even a heuristic that checked likely locations, e.g. just the
global environment, but not all locations might be useful.


If you look at the code in setBreakpoint, it does a search for functions 
in lots of places (but not everywhere).  You could start from that.


Generally speaking, it's quite hard to think of a way to express what 
the answer would look like to a perfectly general function like you're 
asking for.  Suppose a function lives in the parent environment of the 
environment of an expression that was returned in the result of a call 
to lm that is a local variable in the function that called the function 
where you are asking the question:  how would you return that as an 
answer???  In a language with pointers, you'd just return a pointer to 
it, but R doesn't have those.


Duncan Murdoch



On Fri, Oct 23, 2009 at 2:02 PM, Duncan Murdoch murd...@stats.uwo.ca wrote:

On 10/23/2009 1:28 PM, Andrew Yee wrote:


This is kind of a dumb question:  I know you can use isdebugged() to find
out if a specific function is flagged for debugging, but is there a way to
list all the functions that are flagged for debugging?


No, R doesn't keep any master list, it sets a flag in each one.  So you
could iterate over all visible objects to find the ones that have the flag
set (and this is quite slow, there are a lot of places to look), but there
would always be the chance that one was hiding somewhere you didn't look.

Duncan Murdoch

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] help using R's linprog for LP

2009-10-23 Thread Medha Atre
Hi,

I am using R in one of my courses. I am trying to use R's linprog
package to solve to formulate 2-class classification problem as Linear
programming problem.

For my formulation, I need to set to cvec to all 0s.
I know the points are linearly separable so an optimal solution x
does exist, which satisfies all the constraints.

But given the constraints and setting cvec to all 0s it simply gives
me an x of all 0s.

How can I fix this?

Medha

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Gabor Grothendieck
On Fri, Oct 23, 2009 at 2:34 PM, Duncan Murdoch murd...@stats.uwo.ca wrote:
 On 10/23/2009 2:27 PM, Gabor Grothendieck wrote:

 I have often wanted to get such a list as well without having to write
 code myself.  If  a function which automatically iterated over all
 possible objects and listed out the ones that being debugged were
 available in the core it would be useful.  Even if its slow this would
 be used at debug time and not a production time so it would not be too
 bad.  Even a heuristic that checked likely locations, e.g. just the
 global environment, but not all locations might be useful.

 If you look at the code in setBreakpoint, it does a search for functions in
 lots of places (but not everywhere).  You could start from that.

 Generally speaking, it's quite hard to think of a way to express what the
 answer would look like to a perfectly general function like you're asking
 for.  Suppose a function lives in the parent environment of the environment
 of an expression that was returned in the result of a call to lm that is a
 local variable in the function that called the function where you are asking
 the question:  how would you return that as an answer???  In a language with
 pointers, you'd just return a pointer to it, but R doesn't have those.

I agree its problematic but it would still be useful since after you
have debugged a number of functions its easy to forget which ones are
being debugged.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Andrew Yee
Thanks everyone for the responses.  As a suggestion, wouldn't it be useful,
that when you run isdebugged() without any parameters, it would just report
back the functions that are being flagged for debugging?
Andrew

On Fri, Oct 23, 2009 at 2:12 PM, Yihui Xie xieyi...@gmail.com wrote:

 Oops... I forgot to mention that 'envir' (or 'pos') should be
 specified in ls()/get() in my last reply if you are looking for
 debugged functions in environments other than .GlobalEnv.

 Regards,
 Yihui
 --
 Yihui Xie xieyi...@gmail.com
 Phone: 515-294-6609 Web: http://yihui.name
 Department of Statistics, Iowa State University
 3211 Snedecor Hall, Ames, IA



 On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca
 wrote:
  On 10/23/2009 1:28 PM, Andrew Yee wrote:
 
  This is kind of a dumb question:  I know you can use isdebugged() to
 find
  out if a specific function is flagged for debugging, but is there a way
 to
  list all the functions that are flagged for debugging?
 
  No, R doesn't keep any master list, it sets a flag in each one.  So you
  could iterate over all visible objects to find the ones that have the
 flag
  set (and this is quite slow, there are a lot of places to look), but
 there
  would always be the chance that one was hiding somewhere you didn't look.
 
  Duncan Murdoch
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Duncan Murdoch

On 10/23/2009 2:49 PM, Andrew Yee wrote:

Thanks everyone for the responses.  As a suggestion, wouldn't it be useful,
that when you run isdebugged() without any parameters, it would just report
back the functions that are being flagged for debugging?


Could you give an example of what the answer would look like?  Not all 
functions have names, not all names are unique, ...


Duncan Murdoch


Andrew

On Fri, Oct 23, 2009 at 2:12 PM, Yihui Xie xieyi...@gmail.com wrote:


Oops... I forgot to mention that 'envir' (or 'pos') should be
specified in ls()/get() in my last reply if you are looking for
debugged functions in environments other than .GlobalEnv.

Regards,
Yihui
--
Yihui Xie xieyi...@gmail.com
Phone: 515-294-6609 Web: http://yihui.name
Department of Statistics, Iowa State University
3211 Snedecor Hall, Ames, IA



On Fri, Oct 23, 2009 at 1:02 PM, Duncan Murdoch murd...@stats.uwo.ca
wrote:
 On 10/23/2009 1:28 PM, Andrew Yee wrote:

 This is kind of a dumb question:  I know you can use isdebugged() to
find
 out if a specific function is flagged for debugging, but is there a way
to
 list all the functions that are flagged for debugging?

 No, R doesn't keep any master list, it sets a flag in each one.  So you
 could iterate over all visible objects to find the ones that have the
flag
 set (and this is quite slow, there are a lot of places to look), but
there
 would always be the chance that one was hiding somewhere you didn't look.

 Duncan Murdoch

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.






__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Data format for KSVM

2009-10-23 Thread Noah Silverman

Hi,

I have a process using svm from the e1071 library.  it works.

I want to try using the KSVM library instead.  The same data used wiht 
e1071 gives me an error with KSVM.


My data is a data.frame.

sample code:

svm_formula - formula(y ~ a + B + C)

svm_model - ksvm(formula, data=train_data, type=C-svc, 
kernel=rbfdot, C=1)


I get the following error:

object is not a matrix

So I tried this:

svm_model - ksvm(formula, data=as.matrix(train_data), type=C-svc, 
kernel=rbfdot, C=1, scaled=FALSE)


Now I get this error:
Error in model.fram.definition(data = list(v1 = c(1.1234, -2.3232:
Object is not a matrix

My data was previously scaled with the scale() function so that the mean 
is centered at 0. and the range is {-1,1}


Can anyone provide some suggestions as to why I'm getting an error?

Thanks!

-N

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bayesian regression stepwise function?

2009-10-23 Thread Ravi Varadhan
The stepAIC() function in MASS can be used, with k = log(n), to implement
your suggestion of quasi-Bayesian stepwise selection using the BIC
criterion.

Ravi.


---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvarad...@jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 





-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of jlu...@ria.buffalo.edu
Sent: Friday, October 23, 2009 2:31 PM
To: Allan.Y
Cc: r-help@r-project.org; r-help-boun...@r-project.org
Subject: Re: [R] Bayesian regression stepwise function?

The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's 
not stepwise.   Ntzoufras shows how to use WinBUGS to conduct Bayesian 
model selection, but again it's not stepwise


Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of 
Statistical Software 7(7), 1--19.
Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ.
Raftery, A. E. (1995), 'Bayesian model selection in social research', 
Sociological Methodology 25, 111-163.







Allan.Y all...@cmu.edu 
Sent by: r-help-boun...@r-project.org
10/22/2009 01:09 PM

To
r-help@r-project.org
cc

Subject
[R]  Bayesian regression stepwise function?







Hi everyone,

I am wondering if there exists a stepwise regression function for the
Bayesian regression model.  I tried googling, but I couldn't find 
anything. 
I know step function exists for regular stepwise regression, but nothing
for Bayes.


Thanks
-- 
View this message in context: 
http://www.nabble.com/Bayesian-regression-stepwise-function--tp26013725p2601
3725.html

Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Steve Lianoglou

Hi,

On Oct 23, 2009, at 2:39 PM, Gabor Grothendieck wrote:

On Fri, Oct 23, 2009 at 2:34 PM, Duncan Murdoch  
murd...@stats.uwo.ca wrote:

On 10/23/2009 2:27 PM, Gabor Grothendieck wrote:


I have often wanted to get such a list as well without having to  
write

code myself.  If  a function which automatically iterated over all
possible objects and listed out the ones that being debugged were
available in the core it would be useful.  Even if its slow this  
would
be used at debug time and not a production time so it would not be  
too

bad.  Even a heuristic that checked likely locations, e.g. just the
global environment, but not all locations might be useful.


If you look at the code in setBreakpoint, it does a search for  
functions in

lots of places (but not everywhere).  You could start from that.

Generally speaking, it's quite hard to think of a way to express  
what the
answer would look like to a perfectly general function like you're  
asking
for.  Suppose a function lives in the parent environment of the  
environment
of an expression that was returned in the result of a call to lm  
that is a
local variable in the function that called the function where you  
are asking
the question:  how would you return that as an answer???  In a  
language with
pointers, you'd just return a pointer to it, but R doesn't have  
those.


I agree its problematic but it would still be useful since after you
have debugged a number of functions its easy to forget which ones are
being debugged.


If this is something you just want to use on your own functions that  
are set to `debug`, why not just make write a debug wraps the  
base::debug and does some keeps track of which functions you're  
debugging in some global environment variable (of your creation).


You should probably use something beside substitute here, but:

R .DLIST - character()
R debug - function(fun, text=, condition=NULL) {
  .DLIST - c(.DLIST, substitute(fun))
  base::debug(fun, text=text, condition=condition)
}

R myfun - function(something) {
  a - length(something)
  b - sum(something)
  b
}

R debug(myfun)
R debug(myfun)
R .DLIST
[[1]]
myfun

R myfun(1:10)
debugging in: myfun(1:10)
debug: {
a - length(something)
b - sum(something)
b
}
Browse[2] n
debug: a - length(something)
Browse[2] n
debug: b - sum(something)
Browse[2] n
debug: b
Browse[2] n
exiting from: myfun(1:10)
[1] 55

You can define your own undebug function accordingly.

Would that do what you want?

-steve

--
Steve Lianoglou
Graduate Student: Computational Systems Biology
  |  Memorial Sloan-Kettering Cancer Center
  |  Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] reference for mahalanobis {stats}

2009-10-23 Thread Peng Yu
On Fri, Oct 23, 2009 at 10:59 AM, Steve Lianoglou
mailinglist.honey...@gmail.com wrote:
 Hi,

 On Oct 23, 2009, at 11:36 AM, Peng Yu wrote:

 The help on mahalanobis {stats} does not include any reference. I'm
 interested in understand why Mahalanobis is defined in its current way
 and how to use it. Could somebody point me a good book on this? I have
 looked through a few books, but they all give very light explanation
 on it.

 The wikipedia page looks pretty good ... the Intuitive Explanation section
 should help:

 http://en.wikipedia.org/wiki/Mahalanobis_distance

 Googling also provides many more explanations:

 http://matlabdatamining.blogspot.com/2006/11/mahalanobis-distance.html

I have read the above links. But I still feel that a book might
explain it in a more complete fashion.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how do you know which functions are being debugged?

2009-10-23 Thread Gabor Grothendieck
On Fri, Oct 23, 2009 at 2:52 PM, Duncan Murdoch murd...@stats.uwo.ca wrote:
 On 10/23/2009 2:49 PM, Andrew Yee wrote:

 Thanks everyone for the responses.  As a suggestion, wouldn't it be
 useful,
 that when you run isdebugged() without any parameters, it would just
 report
 back the functions that are being flagged for debugging?

 Could you give an example of what the answer would look like?  Not all
 functions have names, not all names are unique, ...


One could use the formal argument list of the function, i.e. the
output of args, when there is no notion of name since that is always
available.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] help using R's linprog for LP

2009-10-23 Thread Medha Atre
Hi,

I found the reason. By default it puts a condition for x = 0. Is
there a way to get rid of this condition?

Medha


On Fri, Oct 23, 2009 at 2:34 PM, Medha Atre medha.a...@gmail.com wrote:
 Hi,

 I am using R in one of my courses. I am trying to use R's linprog
 package to solve to formulate 2-class classification problem as Linear
 programming problem.

 For my formulation, I need to set to cvec to all 0s.
 I know the points are linearly separable so an optimal solution x
 does exist, which satisfies all the constraints.

 But given the constraints and setting cvec to all 0s it simply gives
 me an x of all 0s.

 How can I fix this?

 Medha


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to apply the Wilcoxon test to a hole table at once?

2009-10-23 Thread Iurie Malai

Hi,

I have a data set:

 Dataset
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17
1 user1  m 22 19 28 24 12 18  9   7   4   5   4   7   5   7   9
2 user2  f 25 19 23 18 18 15  6   8   6   6   7  10   7   7   7
3 user3  f 28 21 24 18 15 12 10   6   7   9   5  10   5   9   5
4 user4  f 26 19 26 21 12 18  6   6   5   1   3   8   6   5   6
5 user5  m 21 22 26 18  9  6  4   6   1   7   2   4   4   6   4
6 user6  m 24  8 25 12 18 12  7   8   4   1   4   6   7   5   6
...
71 user71  m 18  4 10  6  3  6  9   5  10   8   4   5   6   5   5

I can apply the Wilcoxon test on each column one by one, but how to do this
on the hole table at once?

 wilcox.test(X3 ~ X2, alternative=two.sided, data=Dataset)

Wilcoxon rank sum test with continuity correction

data:  X3 by X2 
W = 439, p-value = 0.1291
alternative hypothesis: true location shift is not equal to 0 



I researched on this, but I can't find a solution.
I would really appreciate any help. 

P.S. Excuse my lack of terminology :).

Iurie Malai
Moldova Pedagogical State University
-- 
View this message in context: 
http://www.nabble.com/How-to-apply-the-Wilcoxon-test-to-a-hole-table-at-once--tp26030572p26030572.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Reporting discriminat analysis

2009-10-23 Thread Stats Wolf
Dear all,

I want to report the results of discriminant analysis in a way that
would enable the readers to use my results (functions) for
classification of their own cases. For this it should suffice to
provide discriminant functions along with the cut-off values for
distinguishing the groups. However, I cannot find these cut-off
values. I use the function lda of the MASS package, although I have
also looked at some other procedures for discriminant analysis, and
equally failed.

Thanks,
Stats Wolf

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Bayesian regression stepwise function?

2009-10-23 Thread Frank E Harrell Jr

Ravi Varadhan wrote:

The stepAIC() function in MASS can be used, with k = log(n), to implement
your suggestion of quasi-Bayesian stepwise selection using the BIC
criterion.

Ravi.


Although many statisticians use BIC otherwise, it was only designed to 
compare two pre-specified models.


Frank




---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 


Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvarad...@jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 






-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of jlu...@ria.buffalo.edu
Sent: Friday, October 23, 2009 2:31 PM
To: Allan.Y
Cc: r-help@r-project.org; r-help-boun...@r-project.org
Subject: Re: [R] Bayesian regression stepwise function?

The BIC (Raftery) can be used for quasi-Bayesian model selection, but it's 
not stepwise.   Ntzoufras shows how to use WinBUGS to conduct Bayesian 
model selection, but again it's not stepwise



Ntzoufras, I. (2002), 'Gibbs variable selection using BUGS', Journal of 
Statistical Software 7(7), 1--19.

Ntzoufras, I. (2009), Bayesian modeling using WinBUGS, Wiley, Hoboken, NJ.
Raftery, A. E. (1995), 'Bayesian model selection in social research', 
Sociological Methodology 25, 111-163.








Allan.Y all...@cmu.edu 
Sent by: r-help-boun...@r-project.org

10/22/2009 01:09 PM

To
r-help@r-project.org
cc

Subject
[R]  Bayesian regression stepwise function?







Hi everyone,

I am wondering if there exists a stepwise regression function for the
Bayesian regression model.  I tried googling, but I couldn't find 
anything. 
I know step function exists for regular stepwise regression, but nothing

for Bayes.


Thanks



--
Frank E Harrell Jr   Professor and Chair   School of Medicine
 Department of Biostatistics   Vanderbilt University

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


  1   2   >