date:20090909

Try this:

rowsum(Income, format(Date, '%m-%Y'))

or

tapply(Income, format(Date, '%m-%Y'), sum)

On Wed, Sep 9, 2009 at 7:51 AM, clair.crossup...@googlemail.com 
clair.crossup...@googlemail.com wrote:

 Dear all,

 Lets say I have a data frame as follows:


  Date - as.Date(c('2006-08-23', '2006-08-30', '2006-09-06', '2006-09-13',
 '2006-09-20'))
  Income - c(73.79, 72.46, 76.32, 72.43, 72.62)
  data.frame(Date, Income)
Date Income
 1 2006-08-23  73.79
 2 2006-08-30  72.46
 3 2006-09-06  76.32
 4 2006-09-13  72.43
 5 2006-09-20  72.62
 

 is there a way to group the data by month (summing the values in each
 month), i.e.

 Date   Income
 2006-08  146.25
 2006-09  221.37

 Thanks in advance,
 C.C.

 __
 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

[[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] Help with data containing date

Just one other addition.  If by monthly summary you mean a summary
that has an entry for each year/month combo rather than just one
entry for all Januaries, one entry for all Februaries etc. then try this
which works in both the Date and chron cases:

aggregate(z, as.yearmon, mean)


On Tue, Sep 8, 2009 at 11:36 PM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 Try this:

 Lines - Date             A    B        C       D           E
  1978-10-22   18  20.64   0.0   0.176     -1.76
  1978-10-23   15  17.06   0.4   0.147      2.52
  1978-10-24    3   7.588   0.0   0.068     -6.86
  1978-10-25    9  11.491   0.0   0.102    -1.01
  1978-10-26   13  14.98   1.4    0.130     1.26 

 library(zoo); library(chron)
 z - read.zoo(textConnection(Lines), header = TRUE, FUN = as.chron)
 aggregate(z, months, mean)
 aggregate(z, years, mean)

 chron was used above to take advantage of the convenient months()
 and years() functions that are available for it although it could be
 done using Date class with only slightly more complexity:

 library(zoo)
 z - read.zoo(textConnection(Lines), header = TRUE)
 aggregate(z, format(time(z), %m), mean)
 aggregate(z, format(time(z), %Y), mean)

 Read the three zoo vignettes and the article on dates in R News 4/1 as well
 as the references to that article.

 On Tue, Sep 8, 2009 at 11:23 PM, Subodh Acharyashoeb...@gmail.com wrote:
 Hello Everyone,I think this is a very simple problem, I have been struggling
 with it for a few days now.
 I have a 10-year daily data  in the following format.


  Date             A    B        C       D           E
  1978-10-22   18  20.64   0.0   0.176     -1.76
  1978-10-23   15  17.06   0.4   0.147      2.52
  1978-10-24    3   7.588   0.0   0.068     -6.86
  1978-10-25    9  11.491   0.0   0.102    -1.01
  1978-10-26   13  14.98   1.4    0.130     1.26

 I want to calculate the monthly and Annual average averages of A, B, C, D,
 and E, for the 10 years.

 I tried to use the xts package to convert the data into a time series object
 but was not able to even change it into the time series object.
 Any help would be highly appreciated

 Thank you in advance.

 --
 Subodh Acharya
 University of Florida
 Gainesville, FL.

        [[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] Problem with print()

2009-09-09 Thread Amparo Albalate


Dear R-users,
I´m having for the first time a problem while printing out values in the 
screen,
I have a function wich takes quite a long time to execute, and I thought 
it´d be useful to insert a print statement inside the main loop to keep 
control of the current iteration,
however, for some reason, the prints() are not shown in real time, and 
only when I stop the execution, all prints  up to the current iteration 
are shown. This never happened before,

but we observed this problem today (in a machine with vista).
Maybe any of you has experienced the same problem and/or has any idea?

Thanks in advance!

__
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] gridlines in contour plot

2009-09-09 Thread FMH

Dear All,

Someone suggesting me to use the filled.contour function to plot the image 
together with the color index, and an example from the help menu is show below. 
 
#
require(grDevices) # for colours
filled.contour(volcano, color = terrain.colors, asp = 1)# simple
x - 10*1:nrow(volcano)
y - 10*1:ncol(volcano)
filled.contour(x, y, volcano, color = terrain.colors, 
    plot.title = title(main = The Topography of Maunga Whau,
    xlab = Meters North, ylab = Meters West),
    plot.axes = { axis(1, seq(100, 800, by = 100))
  axis(2, seq(100, 600, by = 100)) },
    key.title = title(main=Height\n(meters)),
    key.axes = axis(4, seq(90, 190, by = 10)))# maybe also asp=1
#
 
I am still trying to draw several gridlines with 10m gaps, which might 
represent the height of the mountain on the image but still can't manage it. 
 
Could someone please give some hints?
 
Thank you
Fir





__
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] Very basic question regarding plot.Design...

2009-09-09 Thread Petar Milin


Thank you very much for the help!
Now, I have an additional question regarding transcending from Design to 
rms: Before, it was possible to plot interaction, and get lines with

   plot(ols2, gender=NA, marital.status=NA, xlab='gender',
   ylab='anxiety', conf.int=FALSE)
Now, I am lost how to do that. If model is:
   ols2 - ols(anxiety ~ marital.status * gender)
Then, I need:
   p2 - Predict(ols2, marital.status=., gender=.)
And plot(p2) does not give me lines, but points with conf.int. Also, 
axes are transposed as to what old plot(ols2) would give.
In addition, if I try old function (plot(ols2), above), error message 
returns:

   Error in predictDesign(fit, adj, type = x, non.slopes = non.slopes)
   could not find function Varcov

Best,
PM

Frank E Harrell Jr wrote:

Petar Milin wrote:

I would like to have a line on this plot, instead of two points:

x1 = rnorm(100, 10, 2.5)
x2 = rnorm(100, 26, 3.2)
x1 = as.data.frame(x1)
x2 = as.data.frame(x2)
colnames(x1) = 'anxiety'
colnames(x2) = 'anxiety'
x1$gender = 'male'
x2$gender = 'female'
dat = rbind(x1, x2)

require(Design)

attach(dat)
d=datadist(gender)
options(datadist=d)

ols1 - ols(anxiety ~ gender, data=dat, x=T, y=T)

plot(ols1, gender=NA, xlab=gender, ylab=anxiety,
ylim=c(5,30), conf.int=FALSE)

detach(dat)


Your code has many problems and inefficiencies in it.  Here are some 
suggested fixes and the commands needed using the new rms package:


require(rms)
n - 100
anxiety - c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2))
gender  - c(rep('male', n), rep('female',n))
d - datadist(gender); options(datadist='d')
f - ols(anxiety ~ gender)
p - Predict(f, gender=.)
# For this case could also do p - Predict(f); plot(p) would give a
# vertical dot chart
p   # print estimates
plot(p) # horizontal dot chart; preferred for categorical predictors
# To take control using lines:
with(p, plot(1:2, yhat, type='l', xlab='gender numeric'))

Frank



Thanks!
PM

Frank E Harrell Jr wrote:

Petar Milin wrote:

Hello ALL!
I have a problem to plot factor (lets say gender) as a line, or at 
least both line and point, from ols model:

ols1 - ols(Y ~ gender, data=dat, x=T, y=T)
plot(ols1, gender=NA, xlab=gender, ylab=Y,
ylim=c(5,30), conf.int=FALSE)

If I convert gender into discrete numeric predictor, and use 
forceLines=TRUE, plot is not nice and true, since it shows values 
between 1 and 2.


Thanks!
PM



Petar,

forceLines seems to be doing what it was intended to do.  I'm not 
clear on why you need a line, though.  If you provide self-contained 
code and data that replicate your problem I may be able to help more, 
or you can try a new package I'm about to announce.


Frank







__
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] Forecast - How to create variables with summary() results parameters

Try this:

accuracy(xfor)

The model information can be obtained by:

xfor$model

On Wed, Sep 9, 2009 at 8:11 AM, Pedro Souto pmmso...@gmail.com wrote:

 Hi,
 I would like to create variables in R containing parameters of
 summary(*Forecast
 Results*).
 Using the following code:

 library(forecast)
 data - AirPassengers
 xets - ets(data, model=ZZZ, damped=NULL)
 xfor - forecast(xets,h=12, level=c(80,95))
 summary(xfor)

 the output is:

 Forecast method: ETS(M,A,M)

 Model Information:
 ETS(M,A,M)

 Call:
  ets(y = data, model = ZZZ, damped = NULL)

  Smoothing parameters:
alpha = 0.5901
beta  = 0.0058
gamma = 1e-04

  Initial states:
l = 126.9791
b = 1.6483
s = 0.8865 0.7928 0.9226 1.0582 1.2186 1.2371
   1.1069 0.9818 0.985 1.0149 0.8946 0.901

  sigma:  0.0367

 AIC AICc  BIC
 1391.174 1395.457 1438.691

 In-sample error measures:
ME   RMSEMAEMPE   MAPE   MASE
  1.1644124 10.9944227  8.0668541  0.2032631  2.9361762  0.3119416

 Forecasts:
 Point ForecastLo 80Hi 80Lo 95Hi 95
 Jan 1961   444.9979 424.0692 465.9267 412.9901 477.0057
 Feb 1961   444.1187 419.8324 468.4049 406.9760 481.2613
 Mar 1961   506.4125 475.2920 537.5330 458.8178 554.0072
 Apr 1961   493.9890 460.5937 527.3844 442.9153 545.0628

 How can I associate to a variable for example the In sample error
 measures
 parameters?

 Thanks,
 Pedro

[[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.




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

[[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] gstat---2 basic plot questions

2009-09-09 Thread ONKELINX, Thierry

Dear Steve,

Here are a fex examples for empirical variograms using ggplot2

library(gstat)
library(ggplot2)
data(meuse)
coordinates(meuse) = ~x+y
vgIso - variogram(log(zinc)~x+y, meuse)
vgIso$id - iso
vgAniso - variogram(log(zinc)~x+y, meuse, alpha=c(0,45,90,135))
vgAniso$id - aniso

Empirical - rbind(vgIso, vgAniso)
ggplot(Empirical, aes(x = dist, y = gamma, weight = np / (dist ^ 2),
colour = id, shape = factor(dir.hor))) + geom_smooth() + geom_point()
ggplot(Empirical, aes(x = dist, y = gamma, weight = np / (dist ^ 2),
colour = id)) + geom_smooth() + geom_point() + facet_wrap(~dir.hor)

newVgIso - merge(vgIso[, -4], expand.grid(dist = vgIso$dist, dir.hor =
unique(vgAniso$dir.hor)))
Empirical2 - rbind(newVgIso, vgAniso)
ggplot(Empirical2, aes(x = dist, y = gamma, weight = np / (dist ^ 2),
colour = id)) + geom_smooth() + geom_point() + facet_wrap(~dir.hor)

HTH,

Thierry 




ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-Oorspronkelijk bericht-
Van: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
Namens bevenro
Verzonden: dinsdag 8 september 2009 20:23
Aan: r-help@r-project.org
Onderwerp: [R] gstat---2 basic plot questions


Hi all--

I'm new to R, statistics and programming, so sorry if this is a really
basic question!

I have plotted a directional variogram, and I want to

a. overlay the omni-directional line over each directional panel

b. display the directional variograms in a single panel with a legend
that associated each line to each degree measurement.

The line I'm using is

plot(DirectionalVar,multipanel=FALSE,auto.key=TRUE)

this gives me the single panel, but I don't know how to bind the legend
to the lines (it's tricky to tell which is which) or to overlay the
omni-directional line.

Just looking to keep this as simple as I can--

Thanks!

Steve
--
View this message in context:
http://www.nabble.com/gstat---2-basic-plot-questions-tp25351553p25351553
.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.

Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.

__
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 character to factor in data frame

You could use lapply indeed of sapply:

irisf[,index] - lapply(irisf[,index], as.factor)
str(irisf)

On Wed, Sep 9, 2009 at 8:06 AM, Petr PIKAL petr.pi...@precheza.cz wrote:

 Dear all

 I have a simple problem which I thought is easy to solve but what I tried
 did not work. I want to change character variables to factor in data
 frame. It goes easily from factor to character, but I am stuck in how to
 do backwards conversion.

 Here is an example

 irisf-iris
 irisf[,2]-factor(irisf[,2]) # create second factor

 str(irisf)
 'data.frame':   150 obs. of  5 variables:
  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
  $ Sepal.Width : Factor w/ 23 levels 2,2.2,2.3,..: 15 10 12 11 16 19
 14 14 9 11 ...
  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
  $ Species : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1 1
 1 1 1 1 ...
 index-sapply(irisf, is.factor)
 irisf[,index]-sapply(irisf[,index], as.character)

 str(irisf)
 'data.frame':   150 obs. of  5 variables:
  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
  $ Sepal.Width : chr  3.5 3 3.2 3.1 ...
  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
  $ Species : chr  setosa setosa setosa setosa ...

 I hoped that backwards conversion would be strightforward but...

 irisf[,index]-sapply(irisf[,index], as.factor)
 str(irisf)
 'data.frame':   150 obs. of  5 variables:
  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
  $ Sepal.Width : chr  3.5 3 3.2 3.1 ...
  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
  $ Species : chr  setosa setosa setosa setosa ...

 I want to get both ch columns converted to factor but sapply produces
 character variables (which is documented as it produces matrix or vector)

  R.Version()
 $platform
 [1] i386-pc-mingw32

 $arch
 [1] i386

 $os
 [1] mingw32

 $system
 [1] i386, mingw32

 $status
 [1] Under development (unstable)

 $major
 [1] 2

 $minor
 [1] 10.0

 $year
 [1] 2009

 $month
 [1] 07

 $day
 [1] 15

 $`svn rev`
 [1] 48932

 $language
 [1] R

 $version.string
 [1] R version 2.10.0 Under development (unstable) (2009-07-15 r48932)

 Regards
 Petr

 __
 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

[[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] Facing an error during SVM prediction

2009-09-09 Thread Abbas R. Ali

Hi
I am facing an error after using ksvm() and/or svm() when I can call predict() 
it is giving me the error: 
Error in .local(object, ...) : test vector does not match model !
My dimentions of trainingset is 134 x 95 and validationset is 66 x 94

sample code of prediction:
model.ksvm = ksvm(as.factor(Disposition) ~ ., data = trainingset, kernel = 
rbfdot, kpar = list(sigma=0.05), C = 5, cross = 3)
pred_valid = predict(model.ksvm, validationset, type = decision)  
 
thanks
 


  
[[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] Problem with print()

See ?flush.console

for (i in 1:1e2){
print(i)
flush.console()
}

On Wed, Sep 9, 2009 at 7:28 AM, Amparo Albalate
amparo.albal...@uni-ulm.dewrote:

 Dear R-users,
 Iæ´ having for the first time a problem while printing out values in the
 screen,
 I have a function wich takes quite a long time to execute, and I thought
 itæ² be useful to insert a print statement inside the main loop to keep
 control of the current iteration,
 however, for some reason, the prints() are not shown in real time, and only
 when I stop the execution, all prints  up to the current iteration are
 shown. This never happened before,
 but we observed this problem today (in a machine with vista).
 Maybe any of you has experienced the same problem and/or has any idea?

 Thanks in advance!


 __
 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

[[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] Problem with print()

If you are running with RGUI, change the mode to nonbuffered output
under the Misc tab or use flush.console()

On Wed, Sep 9, 2009 at 6:28 AM, Amparo
Albalateamparo.albal...@uni-ulm.de wrote:
 Dear R-users,
 I惴 having for the first time a problem while printing out values in the
 screen,
 I have a function wich takes quite a long time to execute, and I thought it悲
 be useful to insert a print statement inside the main loop to keep control
 of the current iteration,
 however, for some reason, the prints() are not shown in real time, and only
 when I stop the execution, all prints  up to the current iteration are
 shown. This never happened before,
 but we observed this problem today (in a machine with vista).
 Maybe any of you has experienced the same problem and/or has any idea?

 Thanks in advance!


 __
 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.





-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

__
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] predict-fuction for metaMDS (vegan)

2009-09-09 Thread Kim Vanselow

Dear r-Community,
Step1: I would like to calculate a NMDS (package vegan, function metaMDS) with 
species data.
Step2: Then I want to plot environmental variables over it, using function 
envfit.
The Problem: One of these environmental variables is cos(EXPOSURE). But for 
flat releves there is no exposure. The value is missing and I can't call it 0 
as 0 stands for east and west. Therefore I kicked all releves with missing 
environmental variables. Both, metaMDS and envfit then work without problems.
Now I want to bring the releves with missing environmetal variables back to my 
ordination-plot.

Gavin Simpson gave me the advice to use the predict-function for the same 
missing value problem when I was calculating a cca. This worked without problem.

As my species data was recorded in Braun-Blanquet-numbers (ordinal scale) I 
would prefer to calculate a NMDS. Does anybody know a similar function to the 
predict function which works with NMDS or does anybody know how to modify the 
predict function so that it will work also for NMDS?

Thank you very much!

Kim
 
 Original-Nachricht 
 Datum: Fri, 04 Sep 2009 18:11:09 +0100
 Von: Gavin Simpson gavin.simp...@ucl.ac.uk
 An: Kim Vanselow vanse...@gmx.de
 CC: r-help@r-project.org
 Betreff: Re: [R] NA in cca (vegan)

 On Fri, 2009-09-04 at 17:15 +0200, Kim Vanselow wrote:
  Dear all,
  I would like to calculate a cca (package vegan) with species and
  environmental data. One of these environmental variables is
  cos(EXPOSURE).
  The problem: for flat releves there is no exposure. The value is
  missing and I can't call it 0 as 0 stands for east and west.
  The cca does not run with missing values. What can I do to make vegan
  cca ignoring these missing values?
  Thanks a lot,
  Kim

 Hi Kim,

 This is timely as Jari Oksanen (lead developer on vegan) has been
 looking into making this happen automatically in vegan ordination
 functions. The solution for something like cca is very simple but it
 gets more complicated when you might like to allow features like
 na.exclude etc and have all the functions that operate on objects of
 class cca work nicely.

 For the moment, you should just process your data before it goes into
 cca. Here I assume that you have two data frames; i) Y is the species
 data, and ii) X the environmental data. Further I assume that only one
 variable in X has missings, lets call this Exposure:

 ## dummy data
 set.seed(1234)
 ## 20 samples of 10 species
 Y - data.frame(matrix(rpois(20*10, 2), ncol = 10))
 ## 20 samples and 5 env variables
 X - data.frame(matrix(rnorm(20*5), ncol = 5))
 names(X) - c(paste(Var, 1:4, sep = ), Exposure)
 ## simulate some NAs in Exposure
 X$Exposure[sample(1:20, 3)] - NA
 ## show X
 X

 ## Now create a new variable indicating which are missing
 miss - with(X, is.na(Exposure))

 ## now create new X and Y omitting these rows
 Y2 - Y[!miss, ]
 X2 - X[!miss, ]

 ## Now submit to CCA
 mod - cca(Y2 ~ ., data = X2)
 mod

 ## plot it
 plot(mod, display = c(sites,bp), scaling = 3)

 ## It'd be nice to get predictions for the 3 samples we missed out
 pred - predict(mod, newdata = Y[miss, ], type = wa, scaling = 3)

 ## add these points to the ordination:
 points(pred[, 1:2], col = red, cex = 1.5)

 HTH

 G
 --
 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
  Dr. Gavin Simpson [t] +44 (0)20 7679 0522
  ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
  Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
  Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
  UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

 __
 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.
-- 
Jetzt kostenlos herunterladen: Internet Explorer 8 und Mozilla Firefox 3 -
sicherer, schneller und einfacher! http://portal.gmx.net/de/go/atbrowser

__
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] R-help Digest, Vol 79, Issue 9

2009-09-09 Thread Ista Zahn

 From: Carlos Alzola calz...@cox.net
 To: r-help@r-project.org
 Date: Tue, 8 Sep 2009 23:19:11 -0400
 Subject: [R] SAS vs. R in web application
 Good evening,

 I have been asked to investigate the pros and cons of using SAS vs. R in a 
 web application. Either SAS or R would be the engine used to make some very 
 simple calculations and to produce graphs, preferably in png format.

 The advantages of R are pretty obvious as there would be no licensing issues. 
 The only drawback I can see is that when calling it in batch (using R CMD 
 BATCH), a DOS window appears. Thus I have some basic questions:

 a) Is it possible to have R operate in the background without the DOS window 
 appear? How?
 b) Is it correct that there will be no licensing issues?
 c) What would be an efficient way to run it? I am thinking of having R 
 running in the client's local machine and upload the results to a central 
 server.

I don't think you want to do this. Much better to have the R process
run on the server. My advice is to use the apache2 web server on a
unix-alike system, along with the RApache module
(http://biostat.mc.vanderbilt.edu/rapache/). There are other ways to
use  R in a web page, but in my experience RApache is easier to set
up/use and is more  robust. My website at http://yourpsyche.org uses
RAapche, and I've been very pleased with it.

 If using SAS, would the model described in c) above be the best way to design 
 it, or would it be better to upload the raw data to the server and have SAS 
 perform the calculations there. Would this option require a multi-user SAS 
 license? (I know, I should check with SAS Institute, but I thought I'd ask 
 anyway. Someone in the list may have done something similar).

No experience with SAS, but again I think you're better off  uploading
the data and performing calculations server-side, regardless of the
program you use for the calculations.

 Thanks in advance for any suggestions.

 Carlos Alzola

__
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] Very basic question regarding plot.Design...

2009-09-09 Thread Frank E Harrell Jr


Petar Milin wrote:

Thank you very much for the help!
Now, I have an additional question regarding transcending from Design to 
rms: Before, it was possible to plot interaction, and get lines with

   plot(ols2, gender=NA, marital.status=NA, xlab='gender',
   ylab='anxiety', conf.int=FALSE)
Now, I am lost how to do that. If model is:
   ols2 - ols(anxiety ~ marital.status * gender)
Then, I need:
   p2 - Predict(ols2, marital.status=., gender=.)
And plot(p2) does not give me lines, but points with conf.int. Also, 
axes are transposed as to what old plot(ols2) would give.


You're right, the plot is not optimal for the case of two factor 
predictors.  Here is a way to take control:


m - sample(c('a','b'), 2*n, TRUE)
d -  datadist(d, m)
f - ols(anxiety ~ gender*m)
p - Predict(f, m=., gender=.)
plot(p)  # default plot
Key(.5, .5)

im - as.integer(p$m)
xYplot(Cbind(yhat,lower,upper) ~ im, groups=gender,
   type='l', data=p, ylim=c(5,30), xlab='m',
   scales=list(x=list(at=1:2, labels=levels(p$m))),
   label.curve=list(offset=unit(.3,in)))

For the future I'll see if I can incorporate this automatically in 
plot.Predict with there are only 2 predictors being varied and they are 
both factors.



In addition, if I try old function (plot(ols2), above), error message 
returns:

   Error in predictDesign(fit, adj, type = x, non.slopes = non.slopes)
   could not find function Varcov


That is because we need to update Design to be compatible with the 
latest Hmisc.


Frank



Best,
PM

Frank E Harrell Jr wrote:

Petar Milin wrote:

I would like to have a line on this plot, instead of two points:

x1 = rnorm(100, 10, 2.5)
x2 = rnorm(100, 26, 3.2)
x1 = as.data.frame(x1)
x2 = as.data.frame(x2)
colnames(x1) = 'anxiety'
colnames(x2) = 'anxiety'
x1$gender = 'male'
x2$gender = 'female'
dat = rbind(x1, x2)

require(Design)

attach(dat)
d=datadist(gender)
options(datadist=d)

ols1 - ols(anxiety ~ gender, data=dat, x=T, y=T)

plot(ols1, gender=NA, xlab=gender, ylab=anxiety,
ylim=c(5,30), conf.int=FALSE)

detach(dat)


Your code has many problems and inefficiencies in it.  Here are some 
suggested fixes and the commands needed using the new rms package:


require(rms)
n - 100
anxiety - c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2))
gender  - c(rep('male', n), rep('female',n))
d - datadist(gender); options(datadist='d')
f - ols(anxiety ~ gender)
p - Predict(f, gender=.)
# For this case could also do p - Predict(f); plot(p) would give a
# vertical dot chart
p   # print estimates
plot(p) # horizontal dot chart; preferred for categorical predictors
# To take control using lines:
with(p, plot(1:2, yhat, type='l', xlab='gender numeric'))

Frank



Thanks!
PM

Frank E Harrell Jr wrote:

Petar Milin wrote:

Hello ALL!
I have a problem to plot factor (lets say gender) as a line, or at 
least both line and point, from ols model:

ols1 - ols(Y ~ gender, data=dat, x=T, y=T)
plot(ols1, gender=NA, xlab=gender, ylab=Y,
ylim=c(5,30), conf.int=FALSE)

If I convert gender into discrete numeric predictor, and use 
forceLines=TRUE, plot is not nice and true, since it shows values 
between 1 and 2.


Thanks!
PM



Petar,

forceLines seems to be doing what it was intended to do.  I'm not 
clear on why you need a line, though.  If you provide self-contained 
code and data that replicate your problem I may be able to help 
more, or you can try a new package I'm about to announce.


Frank










--
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.

Re: [R] SRS Required sample size for simulated data

2009-09-09 Thread Robert Baer

My guess is that by table you are really looking for a 'dataframe'. 
Further, I guess that an indexed format is most useful to you.

Try:

# Create the requested dataframe
dat=data.frame(type=c(rep(Hypermarket,10),rep(Supermarket,15),rep(Minimarket,20),rep(Cornershop,20),rep(Spazashop,35)),
value=c(Hypermarket,Supermarket,Minimarket,Cornershop,Spazashop))#  Print out first 15 rows. dat[1:15,]R helpersPlease help me 
combine the simulated data to a form of table where:Hypermarket have 10 rows, supermarket have 15 rows,..., spazashops with35 
rows.Hypermarket - rnorm(10, mean=2, sd=7000)Supermarket - rnorm(15, mean=12000, sd=4000)Minimarket - rnorm(20, mean=1, 
sd=4000)Cornershop - rnorm(20, mean= 8000, sd=3000)Spazashop - rnorm(35, mean= 7000, sd=3000)i am trying to write a code to simulate 
data such that i have 10% ofhypermarkets, 15% of supermarkets,etc.[[alternative HTML version deleted]]-Inline Attachment 
follows-__r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do 
read the posting guide http://www.R-project.org/po!
sting-guide.htmland 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 
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] predict-fuction for metaMDS (vegan)

2009-09-09 Thread Gavin Simpson

On Wed, 2009-09-09 at 14:43 +0200, Kim Vanselow wrote:
 Dear r-Community,
 Step1: I would like to calculate a NMDS (package vegan, function
 metaMDS) with species data.
 Step2: Then I want to plot environmental variables over it, using
 function envfit.
 The Problem: One of these environmental variables is cos(EXPOSURE).
 But for flat releves there is no exposure. The value is missing and I
 can't call it 0 as 0 stands for east and west. Therefore I kicked all
 releves with missing environmental variables. Both, metaMDS and envfit
 then work without problems.
 Now I want to bring the releves with missing environmetal variables
 back to my ordination-plot.
 
 Gavin Simpson gave me the advice to use the predict-function for the
 same missing value problem when I was calculating a cca. This worked
 without problem.

Hi Kim,

Apologies for not replying (yet) to your private email to me (asking
essentially the same as here) --- it was a busy, late, working day for
me yesterday.

Also note that Jari has commented on how you are coding your Exposure
variable; I glossed over that bit when providing my response and you
should probably rethink your coding along the lines Jari suggested.

There isn't a predict function for metaMDS() because there aren't rules
or relationships that would allow you to do what predict.cca does but
for a nMDS. In the CCA case we were estimating the locations in
ordination space for the left-out samples on the basis of their species
composition and computing their site score as a weighted average of the
species scores determined from the data that went into building the CCA.

In nMDS all you have is the species information (at first) which is
converted to dissimilarities and thence to nMDS configurations. In this
sense, none of the information is missing (in your case), so you don't
need to leave out the samples with missing exposures. So far so good.
(Note the point about rethinking your coding, above).

What you now want to do is superimpose upon that plot the env data where
you might have missingness. envfit doesn't allow missings and it is not
immediately clear to me how it might be modified to do so, certainly not
without a lot of changes.

Instead, ordisurf() may be more useful, but you will loose the nice
plot all vectors on the ordination at once feature of envfit (you gain
a lot with ordisurf though as there is no reason to assume anything is
linear across an nMDS configuration).

A cursory look at the guts of ordisurf indicates that it can handle
missing values in the (env) variable you wish to overlay onto the nMDS
plot; the data is passed to mgcv::gam usig the formula interface which
deals with the NA.

If obj contains your nMDS model and Y is a matrix of env variables
containing Exposure, then something like this (untested)

ordisurf(obj, Y$Exposure)

will plot the ordination and surface in a single step. You might want to
control things a bit more so take a look at ?ordisurf.

 
 As my species data was recorded in Braun-Blanquet-numbers (ordinal
 scale) I would prefer to calculate a NMDS. Does anybody know a similar
 function to the predict function which works with NMDS or does anybody
 know how to modify the predict function so that it will work also for
 NMDS?

You could also try capscale() also in vegan. This is like CCA and RDA
but allows you to use any dissimilarity coefficient. It is a bit like a
constrained Principal Coordinates Analysis. This can use the rda method
for predict to do what you did with the CCA earlier.

HTH

G

 
 Thank you very much!
 
 Kim
  
  Original-Nachricht 
  Datum: Fri, 04 Sep 2009 18:11:09 +0100
  Von: Gavin Simpson gavin.simp...@ucl.ac.uk
  An: Kim Vanselow vanse...@gmx.de
  CC: r-help@r-project.org
  Betreff: Re: [R] NA in cca (vegan)
 
  On Fri, 2009-09-04 at 17:15 +0200, Kim Vanselow wrote:
   Dear all,
   I would like to calculate a cca (package vegan) with species and
   environmental data. One of these environmental variables is
   cos(EXPOSURE).
   The problem: for flat releves there is no exposure. The value is
   missing and I can't call it 0 as 0 stands for east and west.
   The cca does not run with missing values. What can I do to make vegan
   cca ignoring these missing values?
   Thanks a lot,
   Kim
 
  Hi Kim,
 
  This is timely as Jari Oksanen (lead developer on vegan) has been
  looking into making this happen automatically in vegan ordination
  functions. The solution for something like cca is very simple but it
  gets more complicated when you might like to allow features like
  na.exclude etc and have all the functions that operate on objects of
  class cca work nicely.
 
  For the moment, you should just process your data before it goes into
  cca. Here I assume that you have two data frames; i) Y is the species
  data, and ii) X the environmental data. Further I assume that only one
  variable in X has missings, lets call this Exposure:
 
  ## dummy data
  set.seed(1234)
  ## 20 samples of 10 species

Re: [R] How to reduce memory demands in a function?

2009-09-09 Thread Richard Gunton


Dear Jim,

Thanks for telling me about gc() - that may help.

Here's some more detail on my function.  It's to implement a methodology for 
classifying plant communities into functional groups according to the values of 
species traits, and comparing the usefulness of alternative classifications by 
using them to summarise species abundance data and then correlating 
site-differences based on these abundance data with site-differences based on 
environmental variables.  The method is described in Pillar et al (2009), J. 
Vegetation Science 20: 334-348.

First the function produces a set of classifications of species by applying the 
function agnes() from
the library cluster to all possible combinations of the variables of a 
species-by-traits dataframe Q.  It also prepares a distance matrix dR based on 
environmental variables for the same sites.  Then there is a loop that takes 
each ith classification in turn and summarises a raw-data dataframe of species 
abundances into a shorter dataframe Xi, by
grouping clusters of its rows according to the classification. It then 
calculates a distance matrix (dXi) based on this summary of abundances, and 
another distance matrix (dQi) based on corresponding variables of the matrix Q 
directly.  Finally in the loop, mantel.partial() from the library vegan is 
used to run a partial mantel
test between dXi and dR,
conditioned on dQi.  The argument
permutations is set to zero, and only the mantel statistic is
stored.  

The loop also contains a forward stepwise selection procedure so that not all 
classifications are actually used. After all classifications using (e.g.) a 
single variable have been used, the variable(s) involved in the best 
classification(s) are specified for inclusion in the next cycle of the loop.  I 
wonder how lucid that all was...

I began putting together the main parts of the code, but I fear it takes so 
much explanation (not to mention editing for transparency) that it may not be 
worth the effort unless someone is really committed to following it through - 
it's about 130 lines in total.  I could still do this if it's likely to be 
worthwhile...

However, the stepwise procedure was only fully implemented after I sent my 
first email.  Now that none of the iterative output is stored except the final 
mantel statistic and essential records of which classifications were used, the 
memory demand has decreased.  The problem that now presents itself is simply 
that the function still takes a very long time to run (e.g. 10 hours to go 
through 7 variables stepwise, with distance matrices of dimension 180 or so).

Two parts of the code that feel clumsy to me already are:
unstack(stack(by(X,f,colSums)))# to reduce a dataframe X to a dataframe 
with fewer rows by summing within sets of rows defined by the factor f
and
V - list();  for(i in 1:n) V[[i]] - grep(pattern[i], x);  names(V) - 1:q;  V 
- stack(V);  V[,1] # to get the indices of multiple matches from x (which is a 
vector of variable names some of which may be repeated) for each of several 
character strings in pattern

The code also involves relating columns of dataframes to each
other using character-matching - e.g. naming columns with paste() -ed strings 
of all the variables used to create them, and then subsequently strsplit() -ing 
these names so that columns can be selected that contain any of the
specified variable names. 

I'm grateful for any advice!

Thanks,   Richard.



From: jim holtman jholt...@gmail.com
To: Richard Gunton r.gun...@talk21.com
Sent: Tuesday, 8 September, 2009 2:08:51 PM
Subject: Re: [R] How to reduce memory demands in a function?

Can you at least post what the function is doing and better yet,
provide commented, minimal, self-contained, reproducible code.  You
can put call to memory.size() to see how large things are growing,
delete temporary objects when not needed, make calls to gc(), ... ,
but it is hard to tell what without an example.


On Mon, Sep 7, 2009 at 4:16 AM, Richard Guntonr..gun...@talk21.com wrote:

I've written a function that regularly throws the cannot allocate vector of 
size X Kb error, since it contains a loop that creates large numbers of big 
distance matrices. I'd be very grateful for any simple advice on how to reduce 
the memory demands of my function.  Besides increasing memory.size to the 
maximum available, I've tried reducing my dist objects to 3 sig. fig.s (not 
sure if that made any difference), I've tried the distance function daisy() 
from package cluster instead of dist(), and I've avoided storing unnecessary 
intermediary objects as far as possible by nesting functions in the same 
command.  I've even tried writing each of my dist() objects to a text file, one 
line for each, and reading them in again one at a time as and when required, 
using scan() - and although this seemed to avoid the memory problem, it ran so 
slowly that it wasn't much use for someone with deadlines to meet...

I

[R] Matrix multiplication and random numbers

2009-09-09 Thread RFish


Dear All

I new to using R and am struggling with some matrix multiplication. 

I have two matrices, one containing random numbers, these are multiplied
together to get another matrix which is different each time. When I put in
another for loop to repeat this process a multiple times the matrices are
all the same. I’m sure there is a way to keep the randomness of the
different matrices but I think I’m putting the for loop in the wrong place.
Also when I try and save the matrices using write.table only the first one
is saved

The code so far is below, any help would be greatly appreciated 

Cheers

Tom 

InitialPop-matrix(c(500,0,0,0,0,0,0))

matmult-function(InitialPop,N){
mat3-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)

for (i in 1:N){
PVAmatrix-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)
mat3-mat3%*%PVAmatrix
}

ans-mat3 %*% InitialPop

return(ans)

}

z-matmult(InitialPop,1)
for (i in 1:4) print(z[,1])
write.table((z),c:\\Documents and Settings\\...) 

-- 
View this message in context: 
http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25365184.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] Derivative of nonparametric curve

2009-09-09 Thread spencerg

 This may be overkill for your application, but you might be 
interested in the fda package, for which a new book appeared a couple 
of months ago:  Functional Data Analysis with R and Matlab (Springer 
Use R! series, by Ramsay, Hooker and Graves;  I'm the third author).  
The package includes a scripts subdirectory with R code to recreate 
all but one of the 76 figures in the book.  [To find this scripts 
directory, use system.file('scripts', package='fda').] 



 Functional data analysis generalizes spline smoothing in two 
important ways: 



  (1) It supports the use of an arbitrary finite basis set to 
approximate elements of a function space;  spline smoothing uses splines 
only, usually cubic splines.  The first derivative of a cubic spline is 
piecewise quadratic, and the second derivative is piecewise linear.  If 
you want something smoother than linear, you need at least a quartic 
spline, and Ramsay has recommended quintics -- degree 5 polynomials = 
order 6 spline. 



  (2) It allows the curve to be smoothed using an arbitrary 
linear differential operator, not just the second derivative.  This can 
be important if you have theory saying that the truth should follow a 
particular differential equation.  Otherwise, if you want to estimate 
the second derivative, Ramsay has recommended smoothing with the fourth 
derivative rather than the second.  (In any event, smoothing is achieved 
by penalized least squares with the penalty being proportional to the 
integral of the square of the chosen linear differential operator.) 



 To reinforce this second point, chapter 11 of Functional Data 
Analysis with R and Matlab describes functional differential 
analysis, which will estimate non-constant coefficients in a 
differential equation model. 



 Hope this helps. 
 Spencer Graves



Liaw, Andy wrote:

From: Rolf Turner
  

On 8/09/2009, at 9:07 PM, FMH wrote:



Dear All,

I'm looking for a way on computing the derivative of first and  
second order of a smoothing curve produced by a nonprametric  
regression. For instance, if we run the R script below, a smooth  
nonparametric regression curve is produced.


provide.data(trawl)
Zone92   - (Year == 0  Zone == 1)
Position - cbind(Longitude - 143, Latitude)
dimnames(Position)[[2]][1] - Longitude - 143
sm.regression(Longitude, Score1, method = aicc, col = red,   
model = linear)


Could someone please give some hints on the way to find the  
derivative on the curve at some points ?
  

See

?smooth.spline
and
?predict.smooth.spline



Since sm.regression() (from the sm package, I presume) uses kernel
methods, a kernel-based estimator of derivatives is available in the
KernSmooth package.

Andy
 
  

cheers,

Rolf Turner

##
Attention:\ This e-mail message is privileged and 
confid...{{dropped:9}}


__
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.



Notice:  This e-mail message, together with any attachme...{{dropped:12}}

__
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] How to reduce memory demands in a function?

If it is taking 10 hours to run, you might want to use Rprof on a
small subset to see where time is being spent.  Also can you use
matrices instead of dataframe because there is a cost in accessing
them if you are just using them in a matrix-like way.  Take a look at
the output from Rprof and that may help identify places to optimize
the code.  Also how large are the distance matrices that you are
creating?  Depending on you input, these can be large.

On Wed, Sep 9, 2009 at 9:24 AM, Richard Guntonr.gun...@talk21.com wrote:

 Dear Jim,

 Thanks for telling me about gc() - that may help..

 Here's some more detail on my function.  It's to implement a methodology for
 classifying plant communities into functional groups according to the values
 of species traits, and comparing the usefulness of alternative
 classifications by using them to summarise species abundance data and then
 correlating site-differences based on these abundance data with
 site-differences based on environmental variables.  The method is described
 in Pillar et al (2009), J. Vegetation Science 20: 334-348.

 First the function produces a set of classifications of species by applying
 the function agnes() from the library cluster to all possible combinations
 of the variables of a species-by-traits dataframe Q.  It also prepares a
 distance matrix dR based on environmental variables for the same sites.
 Then there is a loop that takes each ith classification in turn and
 summarises a raw-data dataframe of species abundances into a shorter
 dataframe Xi, by grouping clusters of its rows according to the
 classification. It then calculates a distance matrix (dXi) based on this
 summary of abundances, and another distance matrix (dQi) based on
 corresponding variables of the matrix Q directly.  Finally in the loop,
 mantel.partial() from the library vegan is used to run a partial mantel
 test between dXi and dR, conditioned on dQi.  The argument permutations is
 set to zero, and only the mantel statistic is stored.

 The loop also contains a forward stepwise selection procedure so that not
 all classifications are actually used. After all classifications using
 (e.g.) a single variable have been used, the variable(s) involved in the
 best classification(s) are specified for inclusion in the next cycle of the
 loop. I wonder how lucid that all was...

 I began putting together the main parts of the code, but I fear it takes so
 much explanation (not to mention editing for transparency) that it may not
 be worth the effort unless someone is really committed to following it
 through - it's about 130 lines in total.  I could still do this if it's
 likely to be worthwhile...

 However, the stepwise procedure was only fully implemented after I sent my
 first email.  Now that none of the iterative output is stored except the
 final mantel statistic and essential records of which classifications were
 used, the memory demand has decreased.  The problem that now presents itself
 is simply that the function still takes a very long time to run (e.g. 10
 hours to go through 7 variables stepwise, with distance matrices of
 dimension 180 or so).

 Two parts of the code that feel clumsy to me already are:
 unstack(stack(by(X,f,colSums)))    # to reduce a dataframe X to a dataframe
 with fewer rows by summing within sets of rows defined by the factor f
 and
 V - list();  for(i in 1:n) V[[i]] - grep(pattern[i], x);  names(V) -
 1:q;  V - stack(V);  V[,1] # to get the indices of multiple matches from x
 (which is a vector of variable names some of which may be repeated) for each
 of several character strings in pattern

 The code also involves relating columns of dataframes to each other using
 character-matching - e.g. naming columns with paste() -ed strings of all the
 variables used to create them, and then subsequently strsplit() -ing these
 names so that columns can be selected that contain any of the specified
 variable names.

 I'm grateful for any advice!

 Thanks,   Richard.

 
 From: jim holtman jholt...@gmail.com
 To: Richard Gunton r.gun...@talk21.com
 Sent: Tuesday, 8 September, 2009 2:08:51 PM
 Subject: Re: [R] How to reduce memory demands in a function?

 Can you at least post what the function is doing and better yet,
 provide commented, minimal, self-contained, reproducible code.  You
 can put call to memory.size() to see how large things are growing,
 delete temporary objects when not needed, make calls to gc(), ... ,
 but it is hard to tell what without an example.


 On Mon, Sep 7, 2009 at 4:16 AM, Richard Guntonr.gun...@talk21.com wrote:

 I've written a function that regularly throws the cannot allocate vector of
 size X Kb error, since it contains a loop that creates large numbers of big
 distance matrices. I'd be very grateful for any simple advice on how to
 reduce the memory demands of my function.  Besides increasing memory.size to
 the maximum available, I've tried reducing my dist objects to 3

Re: [R] Matrix multiplication and random numbers

I am not sure what you mean by being the same each time.  If I make
successive calls to the function, I get different results:

 z
 [,1]
[1,]   0.
[2,]   0.
[3,] 201.6382
[4,]   0.
[5,]   0.
[6,]   0.
[7,]   0.
 matmult(InitialPop,1)
[,1]
[1,]   0.000
[2,]   0.000
[3,] 130.998
[4,]   0.000
[5,]   0.000
[6,]   0.000
[7,]   0.000
 matmult(InitialPop,1)
 [,1]
[1,]   0.
[2,]   0.
[3,] 219.5943
[4,]   0.
[5,]   0.
[6,]   0.
[7,]   0.
 matmult(InitialPop,1)
 [,1]
[1,]   0.
[2,]   0.
[3,] 200.0077
[4,]   0.
[5,]   0.
[6,]   0.
[7,]   0.
 matmult(InitialPop,1)
 [,1]
[1,]   0.
[2,]   0.
[3,] 256.1736
[4,]   0.
[5,]   0.
[6,]   0.
[7,]   0.


On Wed, Sep 9, 2009 at 9:36 AM, RFishtomworthing...@talk21.com wrote:

 Dear All

 I new to using R and am struggling with some matrix multiplication.

 I have two matrices, one containing random numbers, these are multiplied
 together to get another matrix which is different each time. When I put in
 another for loop to repeat this process a multiple times the matrices are
 all the same. I’m sure there is a way to keep the randomness of the
 different matrices but I think I’m putting the for loop in the wrong place.
 Also when I try and save the matrices using write.table only the first one
 is saved

 The code so far is below, any help would be greatly appreciated

 Cheers

 Tom

 InitialPop-matrix(c(500,0,0,0,0,0,0))

 matmult-function(InitialPop,N){
 mat3-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)

 for (i in 1:N){
 PVAmatrix-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)
 mat3-mat3%*%PVAmatrix
 }

 ans-mat3 %*% InitialPop

 return(ans)

 }

 z-matmult(InitialPop,1)
 for (i in 1:4) print(z[,1])
 write.table((z),c:\\Documents and Settings\\...)

 --
 View this message in context: 
 http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25365184.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.




-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

__
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] barplot with lines instead of bars

2009-09-09 Thread John Kane

A clumsy way but it seems to work

data - data.frame(cbind(k = 0:3, fk = c(11, 20,7,2), f0k = c(13.72, 
  17.64, 7.56, 1.08), fkest = c(11.85, 17.78, 8.89, 1.48)))
d - t(data[,2:4])
# barplot(d, beside=TRUE)

xps1 - xps2 - c(.95,1,1.05, 1.95, 2, 2.05, 2.95, 3, 3.05, 3.95, 4, 4.05)
yps1  - rep(0, 12)
yps2 - d

plot(1, 1, xlim=c(0,5),ylim=c(min(d), max(d)), type=n, xaxt=n,
  xlab=Hi There, ylab=Skinny bars)
  
arrows(xps1,yps1, xps2,yps2, angle=0, col=c(red,blue,green), lwd=2)

==

--- On Tue, 9/8/09, rafamoral rafa_moral2...@yahoo.com.br wrote:

 From: rafamoral rafa_moral2...@yahoo.com.br
 Subject: Re: [R] barplot with lines instead of bars
 To: r-help@r-project.org
 Received: Tuesday, September 8, 2009, 2:12 PM
 
 How can I draw thin bars in a barplot?
 Rafael
 
 
 hadley wrote:
  
  What's the difference between a line and a thin bar?
  Hadley
  
  On Tue, Sep 8, 2009 at 12:17 PM, rafamoralrafa_moral2...@yahoo.com.br
  wrote:
 
  I'm sorry, but I think I was misunderstood. What I
 need is something like
  this:
 
  http://img525.imageshack.us/img525/2818/imagemyu.jpg
 
  Lines instead of bars
 
  Thanks!
 
  Rafael.
 
 
  ONKELINX, Thierry wrote:
 
  Here is a solutions using ggplot2 and reshape
 
  library(reshape)
  library(ggplot2)
  data - data.frame(k = 0:3, fk = c(11,
 20,7,2), f0k = c(13.72, 17.64,
  7.56, 1.08), fkest = c(11.85, 17.78, 8.89,
 1.48))
  Molten - melt(data, id.vars = k)
  ggplot(Molten, aes(x = k, y = value, colour =
 variable)) + geom_line()
 
  HTH,
 
  Thierry
 
 
 
  ir. Thierry Onkelinx
  Instituut voor natuur- en bosonderzoek /
 Research Institute for Nature
  and
  Forest
  Cel biometrie, methodologie en kwaliteitszorg
 / Section biometrics,
  methodology and quality assurance
  Gaverstraat 4
  9500 Geraardsbergen
  Belgium
  tel. + 32 54/436 185
  thierry.onkel...@inbo.be
  www.inbo.be
 
  To call in the statistician after the
 experiment is done may be no more
  than asking him to perform a post-mortem
 examination: he may be able to
  say what the experiment died of.
  ~ Sir Ronald Aylmer Fisher
 
  The plural of anecdote is not data.
  ~ Roger Brinner
 
  The combination of some data and an aching
 desire for an answer does not
  ensure that a reasonable answer can be
 extracted from a given body of
  data.
  ~ John Tukey
 
  -Oorspronkelijk bericht-
  Van: r-help-boun...@r-project.org
 [mailto:r-help-boun...@r-project.org]
  Namens Rafael Moral
  Verzonden: dinsdag 8 september 2009 16:45
  Aan: r-help
  Onderwerp: [R] barplot with lines instead of
 bars
 
  Dear useRs,
 
  I want to plot the following barplot with
 lines instead of bars. Is
  there
  a way?
 
  data - data.frame(cbind(k = 0:3, fk =
 c(11, 20,7,2), f0k = c(13.72,
  17.64, 7.56, 1.08), fkest = c(11.85, 17.78,
 8.89, 1.48)))
  d - t(data[,2:4])
  barplot(d, beside=TRUE)
 
  Regards,
  Rafael.
 
 
 
 
 
  [[elided Yahoo spam]]
 
        [[alternative HTML version deleted]]
 
 
  Druk dit bericht a.u.b. niet onnodig af.
  Please do not print this message
 unnecessarily.
 
  Dit bericht en eventuele bijlagen geven enkel
 de visie van de schrijver
  weer
  en binden het INBO onder geen enkel beding,
 zolang dit bericht niet
  bevestigd is
  door een geldig ondertekend document. The
 views expressed in  this
  message
  and any annex are purely those of the writer
 and may not be regarded as
  stating
  an official position of INBO, as long as the
 message is not confirmed by
  a
  duly
  signed document.
 
 
 __
  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.
 
 
 
  --
  View this message in context:
  http://www.nabble.com/barplot-with-lines-instead-of-bars-tp25347695p25350500.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.
 
  
  
  
  -- 
  http://had.co.nz/
  
  __
  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.
  
  
 
 -- 
 View this message in context:

Re: [R] predict-fuction for metaMDS (vegan)

2009-09-09 Thread Dylan Beaudette

On Wed, Sep 9, 2009 at 5:43 AM, Kim Vanselowvanse...@gmx.de wrote:
 Dear r-Community,
 Step1: I would like to calculate a NMDS (package vegan, function metaMDS) 
 with species data.
 Step2: Then I want to plot environmental variables over it, using function 
 envfit.
 The Problem: One of these environmental variables is cos(EXPOSURE). But for 
 flat releves there is no exposure. The value is missing and I can't call it 0 
 as 0 stands for east and west. Therefore I kicked all releves with missing 
 environmental variables. Both, metaMDS and envfit then work without problems.


Hi,

How about using something like modeled solar radiation instead of a
transformed aspect angle? It is fairly simple to do in a GIS (like
GRASS), and can give you much more information that does not come from
a circular distribution.

Cheers,
Dylan



 Now I want to bring the releves with missing environmetal variables back to 
 my ordination-plot.

 Gavin Simpson gave me the advice to use the predict-function for the same 
 missing value problem when I was calculating a cca. This worked without 
 problem.

 As my species data was recorded in Braun-Blanquet-numbers (ordinal scale) I 
 would prefer to calculate a NMDS. Does anybody know a similar function to the 
 predict function which works with NMDS or does anybody know how to modify the 
 predict function so that it will work also for NMDS?

 Thank you very much!

 Kim

  Original-Nachricht 
 Datum: Fri, 04 Sep 2009 18:11:09 +0100
 Von: Gavin Simpson gavin.simp...@ucl.ac.uk
 An: Kim Vanselow vanse...@gmx.de
 CC: r-help@r-project.org
 Betreff: Re: [R] NA in cca (vegan)

 On Fri, 2009-09-04 at 17:15 +0200, Kim Vanselow wrote:
  Dear all,
  I would like to calculate a cca (package vegan) with species and
  environmental data. One of these environmental variables is
  cos(EXPOSURE).
  The problem: for flat releves there is no exposure. The value is
  missing and I can't call it 0 as 0 stands for east and west.
  The cca does not run with missing values. What can I do to make vegan
  cca ignoring these missing values?
  Thanks a lot,
  Kim

 Hi Kim,

 This is timely as Jari Oksanen (lead developer on vegan) has been
 looking into making this happen automatically in vegan ordination
 functions. The solution for something like cca is very simple but it
 gets more complicated when you might like to allow features like
 na.exclude etc and have all the functions that operate on objects of
 class cca work nicely.

 For the moment, you should just process your data before it goes into
 cca. Here I assume that you have two data frames; i) Y is the species
 data, and ii) X the environmental data. Further I assume that only one
 variable in X has missings, lets call this Exposure:

 ## dummy data
 set.seed(1234)
 ## 20 samples of 10 species
 Y - data.frame(matrix(rpois(20*10, 2), ncol = 10))
 ## 20 samples and 5 env variables
 X - data.frame(matrix(rnorm(20*5), ncol = 5))
 names(X) - c(paste(Var, 1:4, sep = ), Exposure)
 ## simulate some NAs in Exposure
 X$Exposure[sample(1:20, 3)] - NA
 ## show X
 X

 ## Now create a new variable indicating which are missing
 miss - with(X, is.na(Exposure))

 ## now create new X and Y omitting these rows
 Y2 - Y[!miss, ]
 X2 - X[!miss, ]

 ## Now submit to CCA
 mod - cca(Y2 ~ ., data = X2)
 mod

 ## plot it
 plot(mod, display = c(sites,bp), scaling = 3)

 ## It'd be nice to get predictions for the 3 samples we missed out
 pred - predict(mod, newdata = Y[miss, ], type = wa, scaling = 3)

 ## add these points to the ordination:
 points(pred[, 1:2], col = red, cex = 1.5)

 HTH

 G
 --
 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
  Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
  ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
  Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
  Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
  UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

 __
 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.
 --
 Jetzt kostenlos herunterladen: Internet Explorer 8 und Mozilla Firefox 3 -
 sicherer, schneller und einfacher! http://portal.gmx.net/de/go/atbrowser

 __
 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

Re: [R] Scan and read.table

2009-09-09 Thread Juliet Hannah

Do you run into problems if you use something like:

cc - rep(numeric,9)
mydata - read.table(yourdata,header=TRUE,colClasses=cc,skip=1,nrows=numRows)

__
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] strange results in summary and IQR functions

2009-09-09 Thread John Kane

Have a look at ?IQR
  Note that this function computes the quartiles using the quantile function 
rather than following Tukey's recommendations, i.e., IQR(x) = quantile(x,3/4) - 
quantile(x,1/4).

It looks like boxplot() gives the results you expect.  

tt - boxplot(x)
tt







--- On Tue, 9/8/09, Chunhao Tu tu_chun...@yahoo.com wrote:

 From: Chunhao Tu tu_chun...@yahoo.com
 Subject: [R]  strange results in summary and IQR functions
 To: r-help@r-project.org
 Received: Tuesday, September 8, 2009, 11:09 AM
 
 Dear R users,
 Something is strange in summary and IQR. Suppose, I have a
 data set and I
 would like to find the Q1, Q2, Q3 and IQR.  
  
 x-c(2,4,11,12,13,15,31,31,37,47)
  summary(x)
    Min. 1st Qu.  Median   
 Mean 3rd Qu.    Max. 
    2.00   11.25   14.00   20.30   31.00   47.00
 
  IQR(x)
 [1] 19.75
 However, I test the same data set in SAS proc univariate,
 and SAS shows
 that Q1=11, Q2=14 and Q3=31. I think most of us agree that
 Q1 is 11 not
 11.25. 
 
 Could someone please explain to me why R shows Q1=11.25 not
 11?
 
 Many Thanks
 Tu
 
 
 -- 
 View this message in context: 
 http://www.nabble.com/strange-results-in-summary-and-IQR-functions-tp25348079p25348079.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.
 


  __
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.

Re: [R] Help-R-graphic

2009-09-09 Thread Arif Chandra


I am Arif. I have made program code for Vector Auto Regressive in terms
of completing my undergraduate program using R. I have an important
question related to my project.
If I have:
data(Canada)
var.2c - VAR(Canada, p = 2, type = const)
var.2c.stabil - stability(var.2c, type = OLS-CUSUM)
I want to get the value of plot(var.2c.stabil). Can you help me what should 
I do or write so the result can occur?
 
It means if I have source code:
data(Canada)
x=acf(Canada)
 
I will get the value of acf if I write x in R.
 



Thank you for answering my question
_
See all the ways you can stay connected to friends and family

[[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] weird vector size cannot be NA/NaN problem

2009-09-09 Thread Werner Wernersen

Hi,

I have a weird problem with my data but I cannot really locate it and cannot 
make a small example data set do reproduce the problem.
I basically divide one numerical column of a data frame with another. When I 
run describe() on that column, I get
Error in vector(integer, length) : vector size cannot be NA/NaN

The two original columns comprise many zero values and I think the particular 
row which causes the problem ends up with an Inf. 

I fear that other functions might also not be reliable for this data and that I 
have to repair it somehow. It actually happens with many of the variables from 
the various tables which originally are in Stata .dta format and which I import 
with the foreign library.

I am working on windows XP and R 2.9.1

I assume that this description is too vague but if anybody has an idea, I would 
appreciate it.

Thanks so much,
  Werner




__
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] Very basic question regarding plot.Design...

2009-09-09 Thread Petar Milin


Sorry, I hope this will be the last (for now, at least).

Following your advices, I did:
 n - 100
 anxiety - c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2),
  rnorm(100, 25, 3.1), rnorm(100, 10, 2.6))
 m.status - c(rep('married', n*2), rep('not married', n*2))
 gender  - c(rep('male', n), rep('female', n),
  rep('male', n), rep('female',n))

 im - as.integer(p2$m.status)
 m.status - as.factor(m.status)
 gender - as.factor(gender)

 require(rms)

 d - datadist(gender, m.status); options(datadist='d')

 ols1 - ols(anxiety ~ m.status)
 p1 - Predict(ols1, m.status=.)
 ols2 - ols(anxiety ~ m.status * gender)
 p2 - Predict(ols2, m.status=., gender=.)

 pdf('figs/anova.pdf', height=6, width=8)
 par(mfrow=c(1,2), mar=c(4,4,1,1), cex=1.2)
 with(p1, plot(1:2, yhat, type='l', xlab='marital status',
  ylab='anxiety', ylim=c(5,30), lwd=1.2))
 xyplot(yhat ~ im, groups=gender,
  type='l', data=p2, ylim=c(5,30), xlab='marital status',
  ylab='anxiety', scales=list(x=list(at=1:2,
  labels=levels(p2$m.status))), col=c('black','black'),
  lwd=1.2, label.curve=list(offset=unit(.3,in)))
 par(mfrow=c(1,1))
 dev.off()

Graphs are great! However, now, they does not appear in one panel.

Sorry, again, for asking so many questions. I am trying to wrap this.

Best,
PM

__
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: mixing colour and linetype in geom_line

2009-09-09 Thread Benoit Boulinguiez

Hi all,
 
I try to represent a multiple curve graphic where the x-axis is the
temperature and the different y-axes are the different X (X22,X43,X44...)
some X corresponds to the same molecule (22 and 44 are for CO2 for instance)
so I use the same colour for them.
 
I wanna mix the linetype with the colour to be able to visually see the
difference between X43 and X45
The best I have done up to now is this code but it crashes with :Error in
col2rgb(colour, TRUE) : invalid color name 'AA'
 
if I skip the linetype in geom it works perfectly of course
 
 
THT_N2_ATGMS_plot-ggplot(THT_N2_ATGMS,aes(x=Temp)) +
 
# geom_line(aes(y=X22,colour=CO2)) +
# geom_line(aes(y=X44,colour=CO2)) +
 
 geom_line(aes(y=X43,colour=AA,linetype=43)) +
 geom_line(aes(y=X45,colour=AA,linetype=45)) 
 
 
 
data set looks like:
 Temp  X22  X43  X44  X45  X48  X58  X60
125.97650 4.62e-12 6.14e-11 3.93e-10 1.29e-11 2.05e-10 6.78e-12 9.31e-12
226.57025 4.73e-12 6.11e-11 3.91e-10 1.28e-11 2.05e-10 6.80e-12 9.43e-12
327.16400 4.62e-12 6.04e-11 3.91e-10 1.27e-11 2.05e-10 6.87e-12 9.28e-12
427.75650 4.57e-12 6.06e-11 3.90e-10 1.27e-11 2.03e-10 6.79e-12 9.26e-12
528.34892 4.64e-12 6.01e-11 3.89e-10 1.28e-11 2.03e-10 6.75e-12 9.18e-12
628.94142 4.71e-12 6.01e-11 3.88e-10 1.27e-11 2.02e-10 6.67e-12 9.24e-12
729.53125 4.70e-12 5.93e-11 3.87e-10 1.27e-11 2.02e-10 6.65e-12 9.20e-12
830.12233 4.64e-12 5.91e-11 3.86e-10 1.26e-11 2.01e-10 6.63e-12 9.11e-12
930.71483 4.71e-12 5.91e-11 3.85e-10 1.25e-11 2.01e-10 6.59e-12 9.17e-12
10   31.30983 4.51e-12 5.88e-11 3.85e-10 1.23e-11 2.00e-10 6.54e-12 9.00e-12
11   31.90233 4.52e-12 5.80e-11 3.85e-10 1.23e-11 1.99e-10 6.60e-12 9.09e-12
12   32.49475 4.47e-12 5.80e-11 3.83e-10 1.24e-11 1.99e-10 6.57e-12 8.95e-12
13   33.08458 4.58e-12 5.79e-11 3.83e-10 1.23e-11 1.98e-10 6.57e-12 9.02e-12
14   33.67575 4.56e-12 5.75e-11 3.82e-10 1.22e-11 1.97e-10 6.43e-12 8.89e-12
15   34.26558 4.53e-12 5.74e-11 3.80e-10 1.23e-11 1.97e-10 6.42e-12 8.97e-12
16   34.85933 4.54e-12 5.70e-11 3.80e-10 1.22e-11 1.96e-10 6.48e-12 8.93e-12
17   35.45183 4.55e-12 5.67e-11 3.79e-10 1.21e-11 1.96e-10 6.38e-12 8.86e-12
18   36.04683 4.54e-12 5.66e-11 3.78e-10 1.20e-11 1.95e-10 6.37e-12 8.90e-12
19   36.63933 4.53e-12 5.61e-11 3.77e-10 1.19e-11 1.94e-10 6.45e-12 8.85e-12
20   37.23042 4.50e-12 5.57e-11 3.77e-10 1.21e-11 1.94e-10 6.35e-12 8.82e-12
21   37.82417 4.55e-12 5.58e-11 3.76e-10 1.19e-11 1.94e-10 6.29e-12 8.86e-12
22   38.41400 4.57e-12 5.57e-11 3.75e-10 1.20e-11 1.92e-10 6.30e-12 8.80e-12

 

Regards/Cordialement

-
Benoit Boulinguiez
Ph.D student
Ecole de Chimie de Rennes (ENSCR) Bureau 1.20 
Equipe CIP UMR CNRS 6226 Sciences Chimiques de Rennes
Avenue du Général Leclerc 
CS 50837 
35708 Rennes CEDEX 7 
Tel 33 (0)2 23 23 80 83
Fax 33 (0)2 23 23 81 20
 http://www.ensc-rennes.fr/ http://www.ensc-rennes.fr/ 

 

[[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] more efficient vectorization of a function ?

2009-09-09 Thread Carlos Hernandez

dear All,
i'm using the following two functions:

share.vector - function (vec1)
{
  vec1 - vec1 - max(vec1,na.rm=TRUE) -0.1  ## this line avoids overflow
  vec1 - exp(vec1)
  vec2 - vec1/(1+sum(vec1,na.rm=TRUE))
  vec2
}

share.matrix - function (mat1)
{
  out1 - apply(mat1,2,share.vector)
  return(out1)
}

vec1 is a vector (of numeric data, usually small numbers), mat1 is a matrix
with many vec1's

is there another way to program them such that they are more efficient (in
terms of time)?

i appreciate any hints or advice.

best regards,
Carlos

[[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] barplot with lines instead of bars

2009-09-09 Thread S Ellison

What is wrong with using plot(..., type=h)?

data - data.frame(cbind(k = 0:3, fk = c(11, 20,7,2), f0k = c(13.72, 
  17.64, 7.56, 1.08), fkest = c(11.85, 17.78, 8.89, 1.48)))
d - t(data[,2:4])

plot(rep(1:4,3)+rep(seq(-0.1,0.1,0.1), 4), as.vector(d), col=rep(1:3,
each=4),type=h, lwd=3, axes=F)
box()
axis(2)
axis(1, at=1:4, labels=1:4)
legend(2.7,20, legend=row.names(d), col=1:3, lwd=3)


 John Kane jrkrid...@yahoo.ca 09/09/2009 15:30:39 
A clumsy way but it seems to work

data - data.frame(cbind(k = 0:3, fk = c(11, 20,7,2), f0k = c(13.72, 
  17.64, 7.56, 1.08), fkest = c(11.85, 17.78, 8.89, 1.48)))
d - t(data[,2:4])
# barplot(d, beside=TRUE)

xps1 - xps2 - c(.95,1,1.05, 1.95, 2, 2.05, 2.95, 3, 3.05, 3.95, 4,
4.05)
yps1  - rep(0, 12)
yps2 - d

plot(1, 1, xlim=c(0,5),ylim=c(min(d), max(d)), type=n, xaxt=n,
  xlab=Hi There, ylab=Skinny bars)
  
arrows(xps1,yps1, xps2,yps2, angle=0, col=c(red,blue,green),
lwd=2)

==

--- On Tue, 9/8/09, rafamoral rafa_moral2...@yahoo.com.br wrote:

 From: rafamoral rafa_moral2...@yahoo.com.br
 Subject: Re: [R] barplot with lines instead of bars
 To: r-help@r-project.org 
 Received: Tuesday, September 8, 2009, 2:12 PM
 
 How can I draw thin bars in a barplot?
 Rafael
 
 
 hadley wrote:
  
  What's the difference between a line and a thin bar?
  Hadley
  
  On Tue, Sep 8, 2009 at 12:17 PM,
rafamoralrafa_moral2...@yahoo.com.br
  wrote:
 
  I'm sorry, but I think I was misunderstood. What I
 need is something like
  this:
 
  http://img525.imageshack.us/img525/2818/imagemyu.jpg 
 
  Lines instead of bars
 
  Thanks!
 
  Rafael.
 
 
  ONKELINX, Thierry wrote:
 
  Here is a solutions using ggplot2 and reshape
 
  library(reshape)
  library(ggplot2)
  data - data.frame(k = 0:3, fk = c(11,
 20,7,2), f0k = c(13.72, 17.64,
  7.56, 1.08), fkest = c(11.85, 17.78, 8.89,
 1.48))
  Molten - melt(data, id.vars = k)
  ggplot(Molten, aes(x = k, y = value, colour =
 variable)) + geom_line()
 
  HTH,
 
  Thierry
 
 


  ir. Thierry Onkelinx
  Instituut voor natuur- en bosonderzoek /
 Research Institute for Nature
  and
  Forest
  Cel biometrie, methodologie en kwaliteitszorg
 / Section biometrics,
  methodology and quality assurance
  Gaverstraat 4
  9500 Geraardsbergen
  Belgium
  tel. + 32 54/436 185
  thierry.onkel...@inbo.be 
  www.inbo.be 
 
  To call in the statistician after the
 experiment is done may be no more
  than asking him to perform a post-mortem
 examination: he may be able to
  say what the experiment died of.
  ~ Sir Ronald Aylmer Fisher
 
  The plural of anecdote is not data.
  ~ Roger Brinner
 
  The combination of some data and an aching
 desire for an answer does not
  ensure that a reasonable answer can be
 extracted from a given body of
  data.
  ~ John Tukey
 
  -Oorspronkelijk bericht-
  Van: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] 
  Namens Rafael Moral
  Verzonden: dinsdag 8 september 2009 16:45
  Aan: r-help
  Onderwerp: [R] barplot with lines instead of
 bars
 
  Dear useRs,
 
  I want to plot the following barplot with
 lines instead of bars. Is
  there
  a way?
 
  data - data.frame(cbind(k = 0:3, fk =
 c(11, 20,7,2), f0k = c(13.72,
  17.64, 7.56, 1.08), fkest = c(11.85, 17.78,
 8.89, 1.48)))
  d - t(data[,2:4])
  barplot(d, beside=TRUE)
 
  Regards,
  Rafael.
 
 
 
 


  [[elided Yahoo spam]]
 
[[alternative HTML version deleted]]
 
 
  Druk dit bericht a.u.b. niet onnodig af.
  Please do not print this message
 unnecessarily.
 
  Dit bericht en eventuele bijlagen geven enkel
 de visie van de schrijver
  weer
  en binden het INBO onder geen enkel beding,
 zolang dit bericht niet
  bevestigd is
  door een geldig ondertekend document. The
 views expressed in  this
  message
  and any annex are purely those of the writer
 and may not be regarded as
  stating
  an official position of INBO, as long as the
 message is not confirmed by
  a
  duly
  signed document.
 
 
 __
  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.
 
 
 
  --
  View this message in context:
 
http://www.nabble.com/barplot-with-lines-instead-of-bars-tp25347695p25350500.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,

Re: [R] Very basic question regarding plot.Design...

2009-09-09 Thread Frank E Harrell Jr


Petar Milin wrote:

Sorry, I hope this will be the last (for now, at least).

Following your advices, I did:
 n - 100
 anxiety - c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2),
  rnorm(100, 25, 3.1), rnorm(100, 10, 2.6))
 m.status - c(rep('married', n*2), rep('not married', n*2))
 gender  - c(rep('male', n), rep('female', n),
  rep('male', n), rep('female',n))

 im - as.integer(p2$m.status)
 m.status - as.factor(m.status)
 gender - as.factor(gender)

 require(rms)

 d - datadist(gender, m.status); options(datadist='d')

 ols1 - ols(anxiety ~ m.status)
 p1 - Predict(ols1, m.status=.)
 ols2 - ols(anxiety ~ m.status * gender)
 p2 - Predict(ols2, m.status=., gender=.)

 pdf('figs/anova.pdf', height=6, width=8)
 par(mfrow=c(1,2), mar=c(4,4,1,1), cex=1.2)
 with(p1, plot(1:2, yhat, type='l', xlab='marital status',
  ylab='anxiety', ylim=c(5,30), lwd=1.2))
 xyplot(yhat ~ im, groups=gender,
  type='l', data=p2, ylim=c(5,30), xlab='marital status',
  ylab='anxiety', scales=list(x=list(at=1:2,
  labels=levels(p2$m.status))), col=c('black','black'),
  lwd=1.2, label.curve=list(offset=unit(.3,in)))
 par(mfrow=c(1,1))
 dev.off()

Graphs are great! However, now, they does not appear in one panel.


Lattice graphics do not respect par(mfrow=).  You have to use multiple 
panels within the graph, or use the lattice print method to compose a 
page with multiple lattice graphs.


Frank



Sorry, again, for asking so many questions. I am trying to wrap this.

Best,
PM




--
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.

Re: [R] optim() argument scoping: passing parameter values into user's subfunction

2009-09-09 Thread Bert Gunter

Inline below.

Bert Gunter
Genentech Nonclinical Biostatistics

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Ben Bolker
Sent: Tuesday, September 08, 2009 7:29 PM
To: r-help@r-project.org
Subject: Re: [R] optim() argument scoping: passing parameter values into
user's subfunction

Szumiloski, John wrote:

 Dear useRs,

 I have a complicated function to be optimized with optim(), and whose
 parameters are passed to another function within its evaluation.  This
 function allows for the parameters to enter as arguments to various
 probability distribution functions.

 However, I am violating some scoping convention, as somewhere within the
 hierarchy of calls a variable is not visible.  I will give a genericized
 example here.

 myFxn - function(parms, Y, phi, other args) {body of function}
 ### I want to optimize this over its first argument

  optim(par=numeric(2), fn=myFxn, ### end of named args, next are
 all in ...
  Y=data,  other args, 
  phi=function(r) pnorm(r, mean=parms[2], sd=parms[1])
  )

This is complicated, but I believe this is what's happening (corrections
from wiser folks would be appreciated if I misstate the situation below):

parms is visible only within the scope of myFxn; by lexical scoping, when
phi is called within myFxn, it goes looking for parms within the environment
it was _defined_, which is optim's, and doesn't find it. Hence the error
message.  Ben's solution below works because it explicitly makes myFxn phi's
closure, thereby giving it access to parms. 

An alternative is to pass parms as a named argument to myFxn. Then it will
look for parms in the _calling_ environment, which will be that of myFxn,
where it will find parms. That is, modify Ben's original version to:

myFxn2 - function(parms, Y,phi,X) {
 ## phi(1) ## not necessary
  sum((Y-(parms[1]+parms[2]*X))^2)
}

set.seed(1001)
X = 1:10
Y = 1+2*X+rnorm(10)

optim(par=c(1,2), fn=myFxn2,
  phi = function(r,parms) pnorm(r, mean=parms[2], sd=parms[1]),
  Y=Y, X=X)

This works just fine.

Given the complexity of the original query, I'm not sure this helps, but
maybe ...

--

 Error in pnorm(r, mean = parms[2], sd = parms[1]) : 
   object 'parms' not found
 debugger()## options(error=expression(dump.frames())) in
 .RProfile
 Message:  Error in pnorm(r, mean = parms[2], sd = parms[1]) : 
   object 'parms' not found
 Available environments had calls:
 1: optim(par = numeric(2), fn = myFxn, Y=data, other args, ..
 2: function (par) 
 3: fn(par, ...)
 4: ifelse(logical vector, phi(Y), other stuff within myFxn
 5: phi(Y)
 6: pnorm(r, mean = parms[2], sd = parms[1])

 Now, when using the debugger in environments 1 and 2, I can see the
 value of 'par' correctly.  I cannot access environment 4 as it just
 returns the original error message.  Trying to access environment 3
 gives the (to me) cryptic  'Error in get(.obj, envir =
 dump[[.selection]]) :  argument ... is missing, with no default' and
 returns to the top level without debugging.  

 I will try to explain to the best of my ability what I think is
 happening here.  Environments 2 and 3 are from the first lines of
 optim(), where it is building an internal function to evaluate the
 candidate parameter values.  When accessing environment 3, it seems like
 when it fills out the ... argument of fn(), it is passing
 phi=function(r) pnorm(r, mean=parms[2], sd=parms[1]) but upon trying to
 evaluate the variable 'parms', it cannot see it in the search path.
 When actually running the original call, 'parms' is apparently not
 evaluated yet, but is once the pnorm call is hit.  It appears the
 'parms' variable is being evaluated before the fn(par) is evaluated into
 myFxn(parms=par).   

 A point which is probably, but not certainly, irrelevant:  myFxn() has
 ... as its final argument, so as to pass tuning arguments to
 integrate().  The function being integrated contains phi(), as well as
 other stuff, of a dummy variable.  The calls in the debugging tree,
 however, are *not* those involving integrate().  

 It would probably be possible to include some
 substitute/eval/as.name/etc. constructions within the function myFxn, in
 order to avoid this problem, but as myFxn already involves numeric
 integrations whose integrand involves the optimized parameters
 themselves, computing on the language at each step of the optimization
 seems like a bad idea.  

 My question: is there a straightforward and efficient way of
 constructing this function and optimizing it with optim(), allowing for
 the argument phi to pass an arbitrary distribution function whose
 parameters are the global ones being optimized?  In particular, in the
 third environment of the debugger tree, is there a way to force the
 fn(par, ...)  myFxn(parms=par, ...) evaluation before the ... get
 evaluated?

 Thanks, John

 John  Szumiloski,  Ph.D.

 Senior Biometrician

Re: [R] Matrix multiplication and random numbers

2009-09-09 Thread Chris Stubben



RFish wrote:
 
 I new to using R and am struggling with some matrix multiplication. 
 

I'm not sure what you're trying to print, but you could place this vector in
an expression

mat3-expression(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0))

# and then evaluate to get a new matrix each time
matrix(eval(mat3), nrow=7)

#I think this may be easier to follow. First create a matrix of zeros, stick
in fertilities and then add random survival probabilities each time


mat3-diag(0,7)
#fertilities
mat3[1,3:7]-c(1.9, 4.8, 9.7, 18, 32.6)
# random survival on sub-diagonal
mat3[row(mat3)==col(mat3)+1]-rnorm(6,0.6021,0.0987)


# and if you want to project the population over 10 time steps in a loop ?

n-matrix(c(500,0,0,0,0,0,0))

popsize - matrix(numeric(7 * 10), nrow = 7)
 for (i in 1:10) {
popsize[, i] - n
 mat3[row(mat3)==col(mat3)+1]-rnorm(6,0.6021,0.0987)
 n - mat3 %*% n
}


 [,1] [,2] [,3] [,4]  [,5]  [,6]   [,7] 
[,8]
[1,]  500   0.   0. 531.6256 709.89940 940.19337 1697.52862
3403.6610
[2,]0 352.5116   0.   0. 298.97874 424.71160  561.32525
1027.1605
[3,]0   0. 279.8029   0.   0.0 231.45988  316.83352 
424.8883
[4,]0   0.   0. 147.8957   0.0   0.0  136.36804 
220.7370
[5,]0   0.   0.   0.  96.92715   0.00.0 
108.6551
[6,]0   0.   0.   0.   0.0  69.875270.0   
0.
[7,]0   0.   0.   0.   0.0   0.0   65.86229   
0.


Chris Stubben

-- 
View this message in context: 
http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25369495.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] barplot with lines instead of bars

2009-09-09 Thread John Kane

My bad memory? I forgot that option existed.

--- On Wed, 9/9/09, S Ellison s.elli...@lgc.co.uk wrote:

 From: S Ellison s.elli...@lgc.co.uk
 Subject: Re: [R] barplot with lines instead of bars
 To: r-help@r-project.org, John Kane jrkrid...@yahoo.ca
 Received: Wednesday, September 9, 2009, 12:02 PM
 What is wrong with using plot(...,
 type=h)?

 data - data.frame(cbind(k = 0:3, fk = c(11, 20,7,2),
 f0k = c(13.72, 
           17.64, 7.56, 1.08),
 fkest = c(11.85, 17.78, 8.89, 1.48)))
 d - t(data[,2:4])

 plot(rep(1:4,3)+rep(seq(-0.1,0.1,0.1), 4), as.vector(d),
 col=rep(1:3,
 each=4),type=h, lwd=3, axes=F)
 box()
 axis(2)
 axis(1, at=1:4, labels=1:4)
 legend(2.7,20, legend=row.names(d), col=1:3, lwd=3)

  John Kane jrkrid...@yahoo.ca
 09/09/2009 15:30:39 
 A clumsy way but it seems to work

 data - data.frame(cbind(k = 0:3, fk = c(11, 20,7,2),
 f0k = c(13.72, 
           17.64, 7.56, 1.08),
 fkest = c(11.85, 17.78, 8.89, 1.48)))
 d - t(data[,2:4])
 # barplot(d, beside=TRUE)

 xps1 - xps2 - c(.95,1,1.05, 1.95, 2, 2.05, 2.95, 3,
 3.05, 3.95, 4,
 4.05)
 yps1  - rep(0, 12)
 yps2 - d

 plot(1, 1, xlim=c(0,5),ylim=c(min(d), max(d)), type=n,
 xaxt=n,
           xlab=Hi There,
 ylab=Skinny bars)

 arrows(xps1,yps1, xps2,yps2, angle=0,
 col=c(red,blue,green),
 lwd=2)

 ==

 --- On Tue, 9/8/09, rafamoral rafa_moral2...@yahoo.com.br
 wrote:

  From: rafamoral rafa_moral2...@yahoo.com.br
  Subject: Re: [R] barplot with lines instead of bars
  To: r-help@r-project.org

  Received: Tuesday, September 8, 2009, 2:12 PM

  How can I draw thin bars in a barplot?
  Rafael

  hadley wrote:

   What's the difference between a line and a thin
 bar?
   Hadley

   On Tue, Sep 8, 2009 at 12:17 PM,
 rafamoralrafa_moral2...@yahoo.com.br
   wrote:

   I'm sorry, but I think I was misunderstood.
 What I
  need is something like
   this:

   http://img525.imageshack.us/img525/2818/imagemyu.jpg 

   Lines instead of bars

   Thanks!

   Rafael.

   ONKELINX, Thierry wrote:

   Here is a solutions using ggplot2 and
 reshape

   library(reshape)
   library(ggplot2)
   data - data.frame(k = 0:3, fk =
 c(11,
  20,7,2), f0k = c(13.72, 17.64,
   7.56, 1.08), fkest = c(11.85, 17.78,
 8.89,
  1.48))
   Molten - melt(data, id.vars = k)
   ggplot(Molten, aes(x = k, y = value,
 colour =
  variable)) + geom_line()

   HTH,

   Thierry

   ir. Thierry Onkelinx
   Instituut voor natuur- en bosonderzoek /
  Research Institute for Nature
   and
   Forest
   Cel biometrie, methodologie en
 kwaliteitszorg
  / Section biometrics,
   methodology and quality assurance
   Gaverstraat 4
   9500 Geraardsbergen
   Belgium
   tel. + 32 54/436 185
   thierry.onkel...@inbo.be

   www.inbo.be 

   To call in the statistician after the
  experiment is done may be no more
   than asking him to perform a post-mortem
  examination: he may be able to
   say what the experiment died of.
   ~ Sir Ronald Aylmer Fisher

   The plural of anecdote is not data.
   ~ Roger Brinner

   The combination of some data and an
 aching
  desire for an answer does not
   ensure that a reasonable answer can be
  extracted from a given body of
   data.
   ~ John Tukey

   -Oorspronkelijk bericht-
   Van: r-help-boun...@r-project.org

  [mailto:r-help-boun...@r-project.org]

   Namens Rafael Moral
   Verzonden: dinsdag 8 september 2009
 16:45
   Aan: r-help
   Onderwerp: [R] barplot with lines instead
 of
  bars

   Dear useRs,

   I want to plot the following barplot
 with
  lines instead of bars. Is
   there
   a way?

   data - data.frame(cbind(k = 0:3, fk
 =
  c(11, 20,7,2), f0k = c(13.72,
   17.64, 7.56, 1.08), fkest = c(11.85,
 17.78,
  8.89, 1.48)))
   d - t(data[,2:4])
   barplot(d, beside=TRUE)

   Regards,
   Rafael.

   [[elided Yahoo spam]]

    [[alternative HTML version deleted]]

   Druk dit bericht a.u.b. niet onnodig af.
   Please do not print this message
  unnecessarily.

   Dit bericht en eventuele bijlagen geven
 enkel
  de visie van de schrijver
   weer
   en binden het INBO onder geen enkel
 beding,
  zolang dit bericht niet
   bevestigd is
   door een geldig ondertekend document.
 The
  views expressed in  this
   message
   and any annex are purely those of the
 writer
  and may not be regarded as
   stating
   an official position of INBO, as long as
 the
  message is not confirmed by
   a
   duly
   signed document.

  __
   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

Re: [R] vector of Vectors

2009-09-09 Thread Giovanni Petris


The code below seems to contradict your claim that vectors cannot
contain vectors.

 zz - vector(mode = list, length = 3)
 zz
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

Best,
Giovanni

 Date: Tue, 08 Sep 2009 23:00:54 -0500 (CDT)
 From: sclar...@illinois.edu
 Sender: r-help-boun...@r-project.org
 Precedence: list
 
 I just learned that vectors can't contain vectors, which frankly simply 
 confuses me (can you imagine arrays or lists in any other language not being 
 able to contain arrays or lists?). At any rate, I need to create a data 
 structure (size to be determined at runtime) which I will instantiate and 
 then fill with vectors using a for loop. Clearly I'm new to the R game. Data 
 Frame doesnt seem to be the right tool clearly.
 
 __
 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.
 
 

-- 

Giovanni Petris  gpet...@uark.edu
Associate Professor
Department of Mathematical Sciences
University of Arkansas - Fayetteville, AR 72701
Ph: (479) 575-6324, 575-8630 (fax)
http://definetti.uark.edu/~gpetris/

__
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] (no subject)

2009-09-09 Thread Jari Oksanen

Hello Kim  Gav,

Gavin Simpson gavin.simpson at ucl.ac.uk writes:

 
 On Wed, 2009-09-09 at 14:43 +0200, Kim Vanselow wrote:
  Dear r-Community,
  Step1: I would like to calculate a NMDS (package vegan, function
  metaMDS) with species data.
  Step2: Then I want to plot environmental variables over it, using
  function envfit.
  The Problem: One of these environmental variables is cos(EXPOSURE).
  But for flat releves there is no exposure. The value is missing and I
  can't call it 0 as 0 stands for east and west. Therefore I kicked all
  releves with missing environmental variables. Both, metaMDS and envfit
  then work without problems.
  Now I want to bring the releves with missing environmetal variables
  back to my ordination-plot.
  
  Gavin Simpson gave me the advice to use the predict-function for the
  same missing value problem when I was calculating a cca. This worked
  without problem.


 
 Also note that Jari has commented on how you are coding your Exposure
 variable; I glossed over that bit when providing my response and you
 should probably rethink your coding along the lines Jari suggested.

Yes, and Dylan Beaudette in gave some more concrete ideas in this thread 
(provided it is the Exposure to sun, and not to, say, wind or pollution source).
 
 There isn't a predict function for metaMDS() because there aren't rules
 or relationships that would allow you to do what predict.cca does but
 for a nMDS. In the CCA case we were estimating the locations in
 ordination space for the left-out samples on the basis of their species
 composition and computing their site score as a weighted average of the
 species scores determined from the data that went into building the CCA.

Actually, you can add new points to NMDS. I once wrote a function for one
user who asked this, but I did not have it in vegan, because its use needed
great technical skill, and was too tricky for a general package. For instance, 
you need to have a rectangular dissimilarity matrix between new points and 
old points (which you can find with Gav's distance() function in analogue). 

 
 What you now want to do is superimpose upon that plot the env data where
 you might have missingness. envfit doesn't allow missings and it is not
 immediately clear to me how it might be modified to do so, certainly not
 without a lot of changes.

True. Seems to be bad design. You should change four functions and the
user interface to have this. Even in that case it would be na.rm (remove rows
with any missing data), and if  you want to do that, you can do it by hand:

scor - scores(mymetaMDSresult)
keep - complete.cases(myenvdata)
envfit(scor[keep,] ~., myenvdata[keep,])
 
 Instead, ordisurf() may be more useful, but you will loose the nice
 plot all vectors on the ordination at once feature of envfit (you gain
 a lot with ordisurf though as there is no reason to assume anything is
 linear across an nMDS configuration).
 
 A cursory look at the guts of ordisurf indicates that it can handle
 missing values in the (env) variable you wish to overlay onto the nMDS
 plot; the data is passed to mgcv::gam usig the formula interface which
 deals with the NA.

Yes, this is true (in most cases).
 

 
 You could also try capscale() also in vegan. This is like CCA and RDA
 but allows you to use any dissimilarity coefficient. It is a bit like a
 constrained Principal Coordinates Analysis. This can use the rda method
 for predict to do what you did with the CCA earlier.
 
However, it does not handle missing values directly (not even in the R-Forge
version, and there are no immediate plans to change this). So you must
remove missing data rows manually.

Cheers, Jari Oksanen

__
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] Matrix multiplication and random numbers

2009-09-09 Thread RFish


Sorry I probably wasn't clear with my description. The reason i put for loop
in was that I want to do the matrix multiplication about 1000 times to get a
1000 different matrices. Therefore I was hoping the for loop would be able
to automate this then use write.table to write to an external document

Cheers

Tom 

   

jholtman wrote:
 
 I am not sure what you mean by being the same each time.  If I make
 successive calls to the function, I get different results:
 
 z
  [,1]
 [1,]   0.
 [2,]   0.
 [3,] 201.6382
 [4,]   0.
 [5,]   0.
 [6,]   0.
 [7,]   0.
 matmult(InitialPop,1)
 [,1]
 [1,]   0.000
 [2,]   0.000
 [3,] 130.998
 [4,]   0.000
 [5,]   0.000
 [6,]   0.000
 [7,]   0.000
 matmult(InitialPop,1)
  [,1]
 [1,]   0.
 [2,]   0.
 [3,] 219.5943
 [4,]   0.
 [5,]   0.
 [6,]   0.
 [7,]   0.
 matmult(InitialPop,1)
  [,1]
 [1,]   0.
 [2,]   0.
 [3,] 200.0077
 [4,]   0.
 [5,]   0.
 [6,]   0.
 [7,]   0.
 matmult(InitialPop,1)
  [,1]
 [1,]   0.
 [2,]   0.
 [3,] 256.1736
 [4,]   0.
 [5,]   0.
 [6,]   0.
 [7,]   0.
 
 
 On Wed, Sep 9, 2009 at 9:36 AM, RFishtomworthing...@talk21.com wrote:

 Dear All

 I new to using R and am struggling with some matrix multiplication.

 I have two matrices, one containing random numbers, these are multiplied
 together to get another matrix which is different each time. When I put
 in
 another for loop to repeat this process a multiple times the matrices are
 all the same. I’m sure there is a way to keep the randomness of the
 different matrices but I think I’m putting the for loop in the wrong
 place.
 Also when I try and save the matrices using write.table only the first
 one
 is saved

 The code so far is below, any help would be greatly appreciated

 Cheers

 Tom

 InitialPop-matrix(c(500,0,0,0,0,0,0))

 matmult-function(InitialPop,N){
 mat3-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)

 for (i in 1:N){
 PVAmatrix-matrix(c(0,rnorm(1,0.6021,0.0987),0,0,0,0,0,0,0,rnorm(1,0.6021,0.0987),0,0,0,0,1.9,0,0,rnorm(1,0.6021,0.0987),0,0,0,4.8,0,0,0,rnorm(1,0.6021,0.0987),0,0,9.7,0,0,0,0,rnorm(1,0.6021,0.0987),0,18,0,0,0,0,0,rnorm(1,0.6021,0.0987),32.6,0,0,0,0,0,0),nrow=7)
 mat3-mat3%*%PVAmatrix
 }

 ans-mat3 %*% InitialPop

 return(ans)

 }

 z-matmult(InitialPop,1)
 for (i in 1:4) print(z[,1])
 write.table((z),c:\\Documents and Settings\\...)

 --
 View this message in context:
 http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25365184.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.

 
 
 
 -- 
 Jim Holtman
 Cincinnati, OH
 +1 513 646 9390
 
 What is the problem that you are trying to solve?
 
 __
 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/Matrix-multiplication-and-random-numbers-tp25365184p25366609.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] R code for creating and appending to frequency table

2009-09-09 Thread drlucyasher


Apologies for what might seem like an simple question.

I have written a model which gives me a frequency distribution for a
particular score within a set. What I now want to do is loop this so that I
get many different frequency distributions and append them to a table with a
collum which specifies which loop the frequency distribution is from. What I
wantto end up with would look something like this:

 Score
 1  2  3  4  5
1   5  6  7  4  1
2   2  5  8  1  1
3   3  8  9  2  1
4   etc 
5
.

Without the loop I can get the model to report the frequency distributions
using the c(table()) command.

SO then I added the loop. I added a command to specify that this object fdt
was controlled by time (t). My original object which calculates the freq
distribution for a single loop is fd 

t- 50 #number of simulations
fdt- numeric (t)#define terms affected by time
for (l in 1:t) #loop start stop

and then at the end of the loop have tried various commands:

e.g.
fdt[l]- c(table(fd))
fdt[l] - fd   c(table(fdt))
c(table(fdt[l]))- fd 
 c(table(fdt[l]))- c(table(fd))
and some other commands using append and assign.

None of these produces what I am looking for and some mess up the rest of my
caculations. When I tried the second of these I get the warning Warning:
number of items to replace is not a multiple of replacement length.

Any help would be much appreciated.

-- 
View this message in context: 
http://www.nabble.com/R-code-for-creating-and-appending-to-frequency-table-tp25368860p25368860.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] lag a data.frame column?

2009-09-09 Thread Mark Knecht

Sometimes it's the simple things...

Why doesn't this lag X$x by 3 and place it in X$x1? (i.e. - Na's in
the first 3 rows and then values showing up...)

The help page does talk about time series. If lag doesn't work on
data.frame columns then what would be the right function to use to lag
by a variable amount?

Thanks,
Mark


X=data.frame(x=seq(1:10))
X$x1=lag(X$x, 3)
X

__
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] R code for creating and appending to frequency table

If the result of the 'table' is a different length each time, then you
will have to consider the use of a list.  If it is the same length,
then something like this should work

results - matrix(ncol=50, nr=10)
for (i in 1:50) results[,i] - your computation for table

On Wed, Sep 9, 2009 at 12:40 PM, drlucyasherlas...@rvc.ac.uk wrote:

 Apologies for what might seem like an simple question.

 I have written a model which gives me a frequency distribution for a
 particular score within a set. What I now want to do is loop this so that I
 get many different frequency distributions and append them to a table with a
 collum which specifies which loop the frequency distribution is from. What I
 wantto end up with would look something like this:

     Score
     1  2  3  4  5
 1   5  6  7  4  1
 2   2  5  8  1  1
 3   3  8  9  2  1
 4   etc
 5
 .

 Without the loop I can get the model to report the frequency distributions
 using the c(table()) command.

 SO then I added the loop. I added a command to specify that this object fdt
 was controlled by time (t). My original object which calculates the freq
 distribution for a single loop is fd

 t- 50 #number of simulations
 fdt- numeric (t)#define terms affected by time
 for (l in 1:t) #loop start stop

 and then at the end of the loop have tried various commands:

 e.g.
 fdt[l]- c(table(fd))
 fdt[l] - fd       c(table(fdt))
 c(table(fdt[l]))- fd
  c(table(fdt[l]))- c(table(fd))
 and some other commands using append and assign.

 None of these produces what I am looking for and some mess up the rest of my
 caculations. When I tried the second of these I get the warning Warning:
 number of items to replace is not a multiple of replacement length.

 Any help would be much appreciated.

 --
 View this message in context: 
 http://www.nabble.com/R-code-for-creating-and-appending-to-frequency-table-tp25368860p25368860.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.




-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

__
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] lag a data.frame column?

2009-09-09 Thread Achim Zeileis


On Wed, 9 Sep 2009, Mark Knecht wrote:


Sometimes it's the simple things...

Why doesn't this lag X$x by 3 and place it in X$x1?


It does.


(i.e. - Na's in the first 3 rows and then values showing up...)


Because this is not how the ts class handles lags.

What happens is that X$x is transformed to ts
  as.ts(X$x)
which is now a regular series with frequency 1 starting at 1 and ending at 
10. If you apply lag(), the data is not modified at all, just the time 
index is shifted

  lag(as.ts(X$x), 3)
Thus it does not create any NAs or - even worse - throws away observations 
(which is not necessary because the frequency time series is known and the 
time index can be extended).


BTW: You almost surely wanted lag(..., -3). Personally, I also don't find 
this intuitive but it's how things are (as documented on the man page).



The help page does talk about time series. If lag doesn't work on
data.frame columns then what would be the right function to use to lag
by a variable amount?


That depends what you want to do. If your data really is a time series, 
then using a time series class (such as ts, or zoo etc.) would 
probably be preferable. This would probably also get you further benefits 
for data processing.


If for some reason you can't do that, it shouldn't be too difficult to 
write a function that does what you want for your personal use

  mylag - function(x, k) c(rep(NA, k), x[1:(length(x)-k)])
which assumes that k is a positive integer and length(x)  k.

Best,
Z

__
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] lag a data.frame column?

Try this:

 x - 1:10
 head(c(rep(NA, 3), x), -3)
 [1] NA NA NA  1  2  3  4  5  6  7
 tail(c(x, rep(NA, 3)), -3)
 [1]  4  5  6  7  8  9 10 NA NA NA

depending on which direction you want.

On Wed, Sep 9, 2009 at 1:43 PM, Mark Knecht markkne...@gmail.com wrote:
 Sometimes it's the simple things...

 Why doesn't this lag X$x by 3 and place it in X$x1? (i.e. - Na's in
 the first 3 rows and then values showing up...)

 The help page does talk about time series. If lag doesn't work on
 data.frame columns then what would be the right function to use to lag
 by a variable amount?

 Thanks,
 Mark


 X=data.frame(x=seq(1:10))
 X$x1=lag(X$x, 3)
 X

 __
 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] fitting nonlinear model

2009-09-09 Thread Bill Hyman

I only have 4 data points and want to fit a curve. It does not work in modreg 
due to too few data. Do you have any idea? Many 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] fitting nonlinear model

2009-09-09 Thread Bill Hyman

My data look like:

Np year
962
915
897
85   10



- Original Message 
From: Bill Hyman billhym...@yahoo.com
To: r-help@r-project.org
Sent: Wednesday, September 9, 2009 11:23:26 AM
Subject: [R] fitting nonlinear model

I only have 4 data points and want to fit a curve. It does not work in modreg 
due to too few data. Do you have any idea? Many 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.

__
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] Writing R Scripts and passing command line arguments

2009-09-09 Thread Abhishek Pratap

Hi All
Thanks for the pointers. reading them I see that what I intend to can
certainly be done without much pain. However with my test scripts I am not
able to fully understand Rscript. The ?Rscript option doesnt print out a lot
to help me get the minute details.

Any good example/s or some info on Rscript will help.

Thanks,
-Abhi

On Mon, Sep 7, 2009 at 2:37 PM, cls59 ch...@sharpsteen.net wrote:

Abhishek Pratap wrote:

1. What's the best way to pass command line arguments to R scripts ?

As Gabor mentioned, the commandArgs function and the getopt package provide
some excellent starting points for this.

Abhishek Pratap wrote:

2. How to execute R scripts from command line ? When I use R CMD BATCH I
see
no output on the screen

I believe R CMD BATCH dumps all of it's output to a file ending in .Rout.
If
you want more control over input and output to your script then Rscript is
the utility to use.

Abhishek Pratap wrote:

3. What does R --slave --vanilla do ?

If I recall correctly, --vanilla makes it so that R does not waste time
trying to load a previously saved session and also makes it so that R does
not try to save history and environment variables when it exits. Vanilla
also disables the loading of options from profile files such as
~/.Rprofile.

I think --slave makes R shut up about it's self and run more quietly than
it
normally does.

Hope that helps!

-Charlie

-
Charlie Sharpsteen
Undergraduate
Environmental Resources Engineering
Humboldt State University
--
View this message in context:
http://www.nabble.com/Writing-R-Scripts-and-passing-command-line-arguments-tp25334067p25334648.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.

Re: [R] ggplot2: mixing colour and linetype in geom_line

2009-09-09 Thread Felipe Carrillo

Is this what you want?

x - structure(list(Temp = c(25.9765, 26.57025, 27.164, 27.7565, 28.34892, 
28.94142, 29.53125, 30.12233, 30.71483, 31.30983, 31.90233, 32.49475, 
33.08458, 33.67575, 34.26558, 34.85933, 35.45183, 36.04683, 36.63933, 
37.23042, 37.82417, 38.414), X22 = c(4.62e-12, 4.73e-12, 4.62e-12, 
4.57e-12, 4.64e-12, 4.71e-12, 4.7e-12, 4.64e-12, 4.71e-12, 4.51e-12, 
4.52e-12, 4.47e-12, 4.58e-12, 4.56e-12, 4.53e-12, 4.54e-12, 4.55e-12, 
4.54e-12, 4.53e-12, 4.5e-12, 4.55e-12, 4.57e-12), X43 = c(6.14e-11, 
6.11e-11, 6.04e-11, 6.06e-11, 6.01e-11, 6.01e-11, 5.93e-11, 5.91e-11, 
5.91e-11, 5.88e-11, 5.8e-11, 5.8e-11, 5.79e-11, 5.75e-11, 5.74e-11, 
5.7e-11, 5.67e-11, 5.66e-11, 5.61e-11, 5.57e-11, 5.58e-11, 5.57e-11
), X44 = c(3.93e-10, 3.91e-10, 3.91e-10, 3.9e-10, 3.89e-10, 3.88e-10, 
3.87e-10, 3.86e-10, 3.85e-10, 3.85e-10, 3.85e-10, 3.83e-10, 3.83e-10, 
3.82e-10, 3.8e-10, 3.8e-10, 3.79e-10, 3.78e-10, 3.77e-10, 3.77e-10, 
3.76e-10, 3.75e-10), X45 = c(1.29e-11, 1.28e-11, 1.27e-11, 1.27e-11, 
1.28e-11, 1.27e-11, 1.27e-11, 1.26e-11, 1.25e-11, 1.23e-11, 1.23e-11, 
1.24e-11, 1.23e-11, 1.22e-11, 1.23e-11, 1.22e-11, 1.21e-11, 1.2e-11, 
1.19e-11, 1.21e-11, 1.19e-11, 1.2e-11), X48 = c(2.05e-10, 2.05e-10, 
2.05e-10, 2.03e-10, 2.03e-10, 2.02e-10, 2.02e-10, 2.01e-10, 2.01e-10, 
2e-10, 1.99e-10, 1.99e-10, 1.98e-10, 1.97e-10, 1.97e-10, 1.96e-10, 
1.96e-10, 1.95e-10, 1.94e-10, 1.94e-10, 1.94e-10, 1.92e-10), 
X58 = c(6.78e-12, 6.8e-12, 6.87e-12, 6.79e-12, 6.75e-12, 
6.67e-12, 6.65e-12, 6.63e-12, 6.59e-12, 6.54e-12, 6.6e-12, 
6.57e-12, 6.57e-12, 6.43e-12, 6.42e-12, 6.48e-12, 6.38e-12, 
6.37e-12, 6.45e-12, 6.35e-12, 6.29e-12, 6.3e-12), X60 = c(9.31e-12, 
9.43e-12, 9.28e-12, 9.26e-12, 9.18e-12, 9.24e-12, 9.2e-12, 
9.11e-12, 9.17e-12, 9e-12, 9.09e-12, 8.95e-12, 9.02e-12, 
8.89e-12, 8.97e-12, 8.93e-12, 8.86e-12, 8.9e-12, 8.85e-12, 
8.82e-12, 8.86e-12, 8.8e-12)), .Names = c(Temp, X22, 
X43, X44, X45, X48, X58, X60), class = data.frame, row.names = 
c(1, 
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20, 21, 22))

attach(x)
 xmelt - melt(x,id.vars=Temp);head(xmelt)
 qplot(Temp,value,colour=variable,linetype=variable,data=xmelt,geom=line)

Felipe D. Carrillo  
Supervisory Fishery Biologist  
Department of the Interior  
US Fish  Wildlife Service  
California, USA


--- On Wed, 9/9/09, Benoit Boulinguiez benoit.boulingu...@ensc-rennes.fr 
wrote:

 From: Benoit Boulinguiez benoit.boulingu...@ensc-rennes.fr
 Subject: [R]  ggplot2: mixing colour and linetype in geom_line
 To: r-help@r-project.org
 Date: Wednesday, September 9, 2009, 8:24 AM
 Hi all,
  
 I try to represent a multiple curve graphic where the
 x-axis is the
 temperature and the different y-axes are the different X
 (X22,X43,X44...)
 some X corresponds to the same molecule (22 and 44 are for
 CO2 for instance)
 so I use the same colour for them.
  
 I wanna mix the linetype with the colour to be able to
 visually see the
 difference between X43 and X45
 The best I have done up to now is this code but it crashes
 with :Error in
 col2rgb(colour, TRUE) : invalid color name 'AA'
  
 if I skip the linetype in geom it works perfectly of
 course
  
  
 THT_N2_ATGMS_plot-ggplot(THT_N2_ATGMS,aes(x=Temp)) +
  
 # geom_line(aes(y=X22,colour=CO2)) +
 # geom_line(aes(y=X44,colour=CO2)) +
  
  geom_line(aes(y=X43,colour=AA,linetype=43)) +
  geom_line(aes(y=X45,colour=AA,linetype=45)) 
  
  
  
 data set looks like:
          Temp   
   X22      X43     
 X44      X45      X48 
     X58      X60
 1    25.97650 4.62e-12 6.14e-11 3.93e-10 1.29e-11
 2.05e-10 6.78e-12 9.31e-12
 2    26.57025 4.73e-12 6.11e-11 3.91e-10 1.28e-11
 2.05e-10 6.80e-12 9.43e-12
 3    27.16400 4.62e-12 6.04e-11 3.91e-10 1.27e-11
 2.05e-10 6.87e-12 9.28e-12
 4    27.75650 4.57e-12 6.06e-11 3.90e-10 1.27e-11
 2.03e-10 6.79e-12 9.26e-12
 5    28.34892 4.64e-12 6.01e-11 3.89e-10 1.28e-11
 2.03e-10 6.75e-12 9.18e-12
 6    28.94142 4.71e-12 6.01e-11 3.88e-10 1.27e-11
 2.02e-10 6.67e-12 9.24e-12
 7    29.53125 4.70e-12 5.93e-11 3.87e-10 1.27e-11
 2.02e-10 6.65e-12 9.20e-12
 8    30.12233 4.64e-12 5.91e-11 3.86e-10 1.26e-11
 2.01e-10 6.63e-12 9.11e-12
 9    30.71483 4.71e-12 5.91e-11 3.85e-10 1.25e-11
 2.01e-10 6.59e-12 9.17e-12
 10   31.30983 4.51e-12 5.88e-11 3.85e-10
 1.23e-11 2.00e-10 6.54e-12 9.00e-12
 11   31.90233 4.52e-12 5.80e-11 3.85e-10
 1.23e-11 1.99e-10 6.60e-12 9.09e-12
 12   32.49475 4.47e-12 5.80e-11 3.83e-10
 1.24e-11 1.99e-10 6.57e-12 8.95e-12
 13   33.08458 4.58e-12 5.79e-11 3.83e-10
 1.23e-11 1.98e-10 6.57e-12 9.02e-12
 14   33.67575 4.56e-12 5.75e-11 3.82e-10
 1.22e-11 1.97e-10 6.43e-12 8.89e-12
 15   34.26558 4.53e-12 5.74e-11 3.80e-10
 1.23e-11 1.97e-10 6.42e-12 8.97e-12
 16   34.85933 4.54e-12 5.70e-11 3.80e-10
 1.22e-11 1.96e-10 6.48e-12 8.93e-12
 17   35.45183 4.55e-12 5.67e-11 3.79e-10
 1.21e-11 1.96e-10 6.38e-12 8.86e-12
 18   36.04683 4.54e-12 5.66e-11 3.78e-10
 1.20e-11 1.95e-10 6.37e-12 8.90e-12
 19   36.63933 4.53e-12 5.61e-11 3.77e-10

Re: [R] eps file with embedded font

2009-09-09 Thread Ted Harding

Thanks, Paul.
In fact I had been hoping to lure you to the surface, from the
12740-km ocean depths which you inhabit, during our hours of
darkness!

Your suggested modification of the options in the embedFonts
command indeed produces an EPS file which displays without clipping.

However, when these files are imported into a PS document using
'groff', the blanking problem described at the end of my included
orifginal posting (below) persisted.

So I posted a summary of the problem to the 'groff' mailing-list.
This can be found in the thread
  [Groff] EPS importing problem, Ted Harding, 2009/09/09
at
  http://lists.gnu.org/archive/html/groff/2009-09/threads.html

There is a reply by Tadziu Hoffmann:

   Attached are two EPS files:
   
 test1.eps
 test1-EMB2.eps
   
  [snip]

  File test1-EMB2.eps[*] contains a call to setpagedevice
  (through setpagesize), which is a bin no-no for EPS files.
  Does the problem still occur If you delete the line that says
  720 360 /letter setpagesize?

[*] which was generated by your modified embedFonts()

and, as my reply confirms, commenting-out that line resolved the
problem (though there a minor issue that the BoundingBox is not
what is really wanted, but that can be resolved by editing the
EPS file). I'm not sure where the insertion of the call to
setpagedevice arose from, but I think it is down to 'gs'.

There was a much longer private reply from a Robert Herrmann,
which I shall send you privately since it may be of interest,
along with the PS filoe I got after implementing Tadziu Hoffmann.

Meanwhile, those who are interested by the issues arising in this
thread, raised by Simone Gabbriellini, may find the above useful.

Ted.

On 08-Sep-09 23:51:27, Paul Murrell wrote:
 Hi
 Thanks for the further analysis on this Ted.  I think the problem is 
 that, with such a wide plot, you are running into the default paper 
 size.  If you look at the EPS produced by ghostscript, you will see a 
 line like this ...
 
 612 792 /letter setpagesize
 
 ... and notice that the value 612 corresponds to the unexpected right 
 hand clipping margin.  So a possible solution is to specify a paper
 size 
 or, more generally, a device size, that is larger (especially wider) 
 than the original plot.  Here's an example for this particular case ...
 
 embedFonts(file=test1.eps,
 outfile=test1-EMB.eps,
 options=-dDEVICEWIDTHPOINTS=720 -dDEVICEHEIGHTPOINTS=360)
 
 Now, rather than clipping the output to the default paper size, the 
 result is clipped to the edges of the plot.
 
 Hope that helps.
 
 Paul
 
 p.s.  Another useful tip that helps in these situations is to get 
 ghostscript to just calculate a bounding box for your plot.  For
 example ...
 
   gs -dNOPAUSE -dBATCH -q -sDEVICE=bbox test1.eps
 %%BoundingBox: 4 18 691 336
 %%HiResBoundingBox: 4.968000 18.71 690.134885 335.214341
 
 ... which shows that ghostscript can produce the right bounding box, if
 it ignores the default paper size for output.
 
 
 Ted Harding wrote:
 I am going back to Simone's original query (though this will
 split the thread) because subsequent replies did not include
 his original. Some comments interspersed below; the main
 response at the end.
 
 I have had some private correspondence with Simone, who sent me
 two of his files that exhibit the problem, and this has enabled
 me to form an idea of where the trouble may lie. It would seem
 that either there is something seriously wrong with the function
 embedFonts(), or with ghostscript when executing the command
 generated by embedFonts(), or with both. I shal descibe the results
 of my esperminets, in the hope that some expert in embedFonts(),
 or in ghostscript, or in both, can make a useful suggestion.
 
 On 04-Sep-09 14:01:44, Simone Gabbriellini wrote:
 Dear list,
 I am trying to make eps file with embedded font.
 I use:

 postscript(ranking-exp-all.eps, horizontal=TRUE, onefile=FALSE,
paper=special, height=8, width=12, family=Helvetica)
 # plot stuff
 dev.off()

 since R does not embed font, I then use:

 embedFonts(file=indegdistr.eps, outfile=indegdistrEMB.eps,
fontpaths=System/Library/Fonts)
 
 I think Simone intended to use a different filename here, probably
 
   embedFonts(file=ranking-exp-all.eps,
  outfile=ranking-exp-all-EMB.eps,
  fontpaths=System/Library/Fonts)
 
 in line with his previous command above. This is not important here.
 
 the problem is that the second file, with font embedded, is cutted
 near the end, and the very last part of the plots and the border are
 off the page...
 
 In fact the bottom of the graphic is slightly clipped when viewed
 in ghostview, and the righthand side is severely clipped.
 
 I use R 2.8.1 on a Mac OSX
 any help more than welcome,
 regards,
 Simone
 
 Ths problem Simone encountered can be reproduced in a simple way
 as follows.
 
   postscript(file=test1.eps,family=Times,horizontal=FALSE,

Re: [R] Writing R Scripts and passing command line arguments

Put this line in a file called test.R:
cat(command args are:\n); print(commandArgs())

and then call it like this from the shell or Windows cmd line assuming
that Rscript is in your path:

Rscript test.R abc def

On Wed, Sep 9, 2009 at 2:37 PM, Abhishek Pratap abhishek@gmail.com wrote:
Hi All
Thanks for the pointers. reading them I see that what I intend to can
certainly be done without much pain. However with my test scripts I am not
able to fully understand Rscript. The ?Rscript option doesnt print out a lot
to help me get the minute details.

Any good example/s or some info on Rscript will help.

Thanks,
-Abhi

On Mon, Sep 7, 2009 at 2:37 PM, cls59 ch...@sharpsteen.net wrote:

Abhishek Pratap wrote:

1. What's the best way to pass command line arguments to R scripts ?

As Gabor mentioned, the commandArgs function and the getopt package provide
some excellent starting points for this.

Abhishek Pratap wrote:

2. How to execute R scripts from command line ? When I use R CMD BATCH I
see
no output on the screen

I believe R CMD BATCH dumps all of it's output to a file ending in .Rout.
If
you want more control over input and output to your script then Rscript is
the utility to use.

Abhishek Pratap wrote:

3. What does R --slave --vanilla do ?

I think --slave makes R shut up about it's self and run more quietly than
it
normally does.

Hope that helps!

-Charlie

[[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] fitting nonlinear model

2009-09-09 Thread cls59




Bill Hyman wrote:
 
 My data look like:
 
 Np year
 962
 915
 897
 85   10
 
 

And which equation are you trying to fit to this data?

-Charlie

-
Charlie Sharpsteen
Undergraduate
Environmental Resources Engineering
Humboldt State University
-- 
View this message in context: 
http://www.nabble.com/fitting-nonlinear-model-tp25370697p25371220.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] fitting nonlinear model

2009-09-09 Thread milton ruser

Hi Bill,

I am not sure what you want, but...

mydf-read.table(stdin(), head=T, sep=,)
Np,year
96,2
91,5
89,7
85,10

plot(Np~year, data=mydf)
mymodel-lm(Np~year, data=mydf)
abline(mymodel, col=red)

bests
milton

On Wed, Sep 9, 2009 at 2:27 PM, Bill Hyman billhym...@yahoo.com wrote:

 My data look like:

 Np year
 962
 915
 897
 85   10



 - Original Message 
 From: Bill Hyman billhym...@yahoo.com
 To: r-help@r-project.org
 Sent: Wednesday, September 9, 2009 11:23:26 AM
 Subject: [R] fitting nonlinear model

 I only have 4 data points and want to fit a curve. It does not work in
 modreg due to too few data. Do you have any idea? Many 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.htmlhttp://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.htmlhttp://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] fitting nonlinear model

2009-09-09 Thread Bill Hyman

Hi Milton,

Thanks for your help. Actually, I would like to fit a non-linear fashion. For 
some data like below, 'lm' may not work very well. Do you have idea? Thanks 
again!

Bill



0 1 
2 0.9 
5 0.5 
7 0.1 
10 0.01 





From: milton ruser milton.ru...@gmail.com

Cc: R-help@r-project.org
Sent: Wednesday, September 9, 2009 12:01:52 PM
Subject: Re: [R] fitting nonlinear model


Hi Bill,
 
I am not sure what you want, but...
 
mydf-read.table(stdin(), head=T, sep=,)
Np,year
96,2
91,5
89,7
85,10
 
plot(Np~year, data=mydf)
mymodel-lm(Np~year, data=mydf)
abline(mymodel, col=red)
 
bests 
milton




My data look like:

Np year
962
915
897
85   10




- Original Message 

To: r-help@r-project.org
Sent: Wednesday, September 9, 2009 11:23:26 AM
Subject: [R] fitting nonlinear model

[[elided Yahoo spam]]

__
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.




  
[[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] How to return/show the indexes of unusual points in boxplot?

2009-09-09 Thread guox

I am wondering if you know how to return by function or show in boxplot,
the indexes of unusual points, such as,
points that are outside the box or in [Q3+1.5IQR,Max].
For example,

 boxplot(c(4,rnorm(20),8))

There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8 in
the boxplot
or return them by function?

Thanks,

-james

__
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] fitting nonlinear model

2009-09-09 Thread cls59




Bill Hyman wrote:
 
 Hi Milton,
 
 Thanks for your help. Actually, I would like to fit a non-linear fashion.
 For some data like below, 'lm' may not work very well. Do you have idea?
 Thanks again!
 
 

That's why information equation you are trying to fit is very important. For
example, the BOD data set in R is:

BOD
  Time demand
118.3
22   10.3
33   19.0
44   16.0
55   15.6
67   19.8

BOD demand can be modeled as a function of Time using the following
equation:

demand = BODu * ( 1 - exp( -K * Time ) )

Where BODu and K are the unknown parameters of the model. One way of doing a
non-linear fit in R is to use nls(), the nonlinear least-squares function:

model - nls( demand ~ BODu * ( 1 - exp( -K * Time ) ), 
data = BOD, 
start = list( BODu = max( BOD[['demand']]), k = 0.1 ) 
)

Note that with nls(), it is necessary to provide starting guesses for the
parameters as a list using the start parameters of the nls function.

Hope this helps!

-Charlie

P.S. 
Providing an example of the equation you are trying to fit to your data will
help us provide an answer that is more specific to your situation.


-
Charlie Sharpsteen
Undergraduate
Environmental Resources Engineering
Humboldt State University
-- 
View this message in context: 
http://www.nabble.com/fitting-nonlinear-model-tp25370697p25371862.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] The code behind the function

2009-09-09 Thread Chunhao Tu


Hi R users,
I have a question. How can I see the code behind the function. For example,

 boxplot
function (x, ...) 
UseMethod(boxplot)
environment: namespace:graphics

I really would like to see how people code this. Could someone please show
me how to see the code behind the function?
Many Thanks
Tu


-- 
View this message in context: 
http://www.nabble.com/The-code-behind-the-function-tp25370743p25370743.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] fitting nonlinear model

2009-09-09 Thread milton ruser

May be this:

mydf-read.table(stdin(), head=T, sep=,)
Xvar,Yvar
0,1
2,0.9
5,0.5
7,0.1
10,0.01
mymodel-loess(Yvar~Xvar, data=mydf)

plot(mymodel)
mydf.new-data.frame(cbind(Xvar=seq(from=0, to=10, by=0.1)))
mydf.new$Yvar-predict(mymodel, newdata=mydf.new)
lines(Yvar~Xvar, data=mydf.new, lty=2)

cheers

milton

On Wed, Sep 9, 2009 at 3:01 PM, milton ruser milton.ru...@gmail.com wrote:

 Hi Bill,

 I am not sure what you want, but...

 mydf-read.table(stdin(), head=T, sep=,)
 Np,year
 96,2
 91,5
 89,7
 85,10

 plot(Np~year, data=mydf)
 mymodel-lm(Np~year, data=mydf)
 abline(mymodel, col=red)

 bests
 milton

   On Wed, Sep 9, 2009 at 2:27 PM, Bill Hyman billhym...@yahoo.com wrote:

 My data look like:

 Np year
 962
 915
 897
 85   10



 - Original Message 
 From: Bill Hyman billhym...@yahoo.com
 To: r-help@r-project.org
 Sent: Wednesday, September 9, 2009 11:23:26 AM
 Subject: [R] fitting nonlinear model

 I only have 4 data points and want to fit a curve. It does not work in
 modreg due to too few data. Do you have any idea? Many 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.htmlhttp://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.htmlhttp://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] Combining simulated data

2009-09-09 Thread Schalk Heunis

Kabeli

I think this is what you want:

Hypermarket - matrix(rnorm(10, mean=2, sd=7000))
Supermarket - matrix(rnorm(15, mean=12000, sd=4000))
Minimarket   - matrix(rnorm(20, mean=1, sd=4000))
Cornershop  - matrix(rnorm(20, mean= 8000, sd=3000))
Spazashop   - matrix(rnorm(35, mean= 7000, sd=3000))
rbind(Hypermarket,Supermarket,Minimarket,Cornershop,Spazashop)

This creates a matrix with one column and 100 rows.  You can add names
also by doing something like this for e.g. Spazashop
Spazashop   - matrix(rnorm(35, mean= 7000, sd=3000), dimnames =
list(paste(Spazashop,1:35)))

Hope this helps.

Schalk Heunis



On Wed, Sep 9, 2009 at 2:26 AM, KABELI MEFANEkabelimef...@yahoo.co.uk wrote:
 R helpers

 Please help me combine the simulated data to a form of table where: 
 Hypermarket have 10 rows, supermarket have 15 rows,..., spazashops with 
 35 rows.

 Hypermarket - rnorm(10, mean=2, sd=7000)
 Supermarket - rnorm(15, mean=12000, sd=4000)
 Minimarket   - rnorm(20, mean=1, sd=4000)
 Cornershop  - rnorm(20, mean= 8000, sd=3000)
 Spazashop   - rnorm(35, mean= 7000, sd=3000)

 i am trying to write a code to simulate data such that i have 10% of 
 hypermarkets, 15% of supermarkets,etc.




        [[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] The code behind the function

2009-09-09 Thread Erik Iverson

--
How can I see the code behind the function. For example,

 boxplot
function (x, ...) 
UseMethod(boxplot)
environment: namespace:graphics

I really would like to see how people code this.
--

That *is* the code.  boxplot is a generic function, which calls another 
function based on what class of object you pass to it. You can find out which 
classes the boxplot method works on by using 

 methods(boxplot) 

[1] boxplot.default  boxplot.formula*

   Non-visible functions are asterisked

Next, try typing 

 boxplot.default

to see the implementation for the default function, and then since 
boxplot.formula is hidden, usually

getAnywhere(boxplot.formula) 

will retrieve it.

Section 10.9 of An Introduction to R explains this more:

http://cran.r-project.org/doc/manuals/R-intro.html#Object-orientation

Best Regards,
Erik Iverson   

__
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 return/show the indexes of unusual points in boxplot?

Try this:

bp - boxplot(c(4,rnorm(20),8))
bp$out


On Wed, Sep 9, 2009 at 4:25 PM, g...@ucalgary.ca wrote:

 I am wondering if you know how to return by function or show in boxplot,
 the indexes of unusual points, such as,
 points that are outside the box or in [Q3+1.5IQR,Max].
 For example,

  boxplot(c(4,rnorm(20),8))

 There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8 in
 the boxplot
 or return them by function?

 Thanks,

 -james

 __
 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

[[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] The code behind the function

2009-09-09 Thread cls59




Chunhao Tu wrote:
 
 Hi R users,
 I have a question. How can I see the code behind the function. For
 example,
 
 boxplot
 function (x, ...) 
 UseMethod(boxplot)
 environment: namespace:graphics
 
 I really would like to see how people code this. Could someone please show
 me how to see the code behind the function?
 Many Thanks
 Tu
 
 


When you type a function name and see something like:

UseMethod(boxplot)

This means the the function you are inquiring about is a generic function-
that is it examines the class of the object passed to it and then calls a
class function. These class functions, also known as methods, all have the
form:

boxplot.class-name

So, a quick way to see which methods may be available is to use the methods
function()

 methods( boxplot )
[1] boxplot.default  boxplot.formula* boxplot.matrix

This shows that there are visible methods for objects of class matrix and
stats and a default method that is used in all other cases. Sometimes there
may be other methods that have been hidden in away in package namespaces
where it is harder to find them as in the case of boxplot.formula.

For visible methods, simply type their name to see their contents. I suspect
you will be most interested in boxplot.default:

 boxplot.default
[output omitted]

For non-visible methods, like boxplot.formula, use the getAnywhere function:

 getAnywhere( 'boxplot.formula' )
[output ommitted]

Hope this helps!

-Charlie

-
Charlie Sharpsteen
Undergraduate
Environmental Resources Engineering
Humboldt State University
-- 
View this message in context: 
http://www.nabble.com/The-code-behind-the-function-tp25370743p25372134.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 to return/show the indexes of unusual points in boxplot?

boxplot returns a dataframe that has the values in it at $out:

 x - boxplot(c(4,rnorm(20),8))
 x
$stats
   [,1]
[1,] -1.5364498
[2,] -0.5282799
[3,] -0.1398736
[4,]  0.3065579
[5,]  1.3430388

$n
[1] 22

$conf
   [,1]
[1,] -0.4210947
[2,]  0.1413474

$out
[1] 4 8

$group
[1] 1 1

$names
[1] 1



On Wed, Sep 9, 2009 at 3:25 PM, g...@ucalgary.ca wrote:
 I am wondering if you know how to return by function or show in boxplot,
 the indexes of unusual points, such as,
 points that are outside the box or in [Q3+1.5IQR,Max].
 For example,

 boxplot(c(4,rnorm(20),8))

 There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8 in
 the boxplot
 or return them by function?

 Thanks,

 -james

 __
 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.




-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

__
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] lag a data.frame column?

2009-09-09 Thread Mark Knecht

On Wed, Sep 9, 2009 at 11:09 AM, Achim Zeileis
achim.zeil...@wu-wien.ac.at wrote:
 On Wed, 9 Sep 2009, Mark Knecht wrote:

 Sometimes it's the simple things...

 Why doesn't this lag X$x by 3 and place it in X$x1?

 It does.

 (i.e. - Na's in the first 3 rows and then values showing up...)

 Because this is not how the ts class handles lags.

 What happens is that X$x is transformed to ts
  as.ts(X$x)
 which is now a regular series with frequency 1 starting at 1 and ending at
 10. If you apply lag(), the data is not modified at all, just the time index
 is shifted
  lag(as.ts(X$x), 3)
 Thus it does not create any NAs or - even worse - throws away observations
 (which is not necessary because the frequency time series is known and the
 time index can be extended).

 BTW: You almost surely wanted lag(..., -3). Personally, I also don't find
 this intuitive but it's how things are (as documented on the man page).

 The help page does talk about time series. If lag doesn't work on
 data.frame columns then what would be the right function to use to lag
 by a variable amount?

 That depends what you want to do. If your data really is a time series, then
 using a time series class (such as ts, or zoo etc.) would probably be
 preferable. This would probably also get you further benefits for data
 processing.

 If for some reason you can't do that, it shouldn't be too difficult to write
 a function that does what you want for your personal use
  mylag - function(x, k) c(rep(NA, k), x[1:(length(x)-k)])
 which assumes that k is a positive integer and length(x)  k.

 Best,
 Z



Thank you very much for the explanation. It is helpful.

I think the function is probably the best answer for me short-term as
I don't know much about time series - I need to learn - and I have a
lot of data.frames where the function can help. Over time if I learn
about ts and zoo maybe that will make my code a bit better.

Thanks,
Mark

__
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 return/show the indexes of unusual points in boxplot?

2009-09-09 Thread guox

Thank you very much.
But it seems that
x$out returns the values not the indexes of the  values (1,22).
-james

 boxplot returns a dataframe that has the values in it at $out:

 x - boxplot(c(4,rnorm(20),8))
 x
 $stats
[,1]
 [1,] -1.5364498
 [2,] -0.5282799
 [3,] -0.1398736
 [4,]  0.3065579
 [5,]  1.3430388

 $n
 [1] 22

 $conf
[,1]
 [1,] -0.4210947
 [2,]  0.1413474

 $out
 [1] 4 8

 $group
 [1] 1 1

 $names
 [1] 1



 On Wed, Sep 9, 2009 at 3:25 PM, g...@ucalgary.ca wrote:
 I am wondering if you know how to return by function or show in
 boxplot,
 the indexes of unusual points, such as,
 points that are outside the box or in [Q3+1.5IQR,Max].
 For example,

 boxplot(c(4,rnorm(20),8))

 There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8
 in
 the boxplot
 or return them by function?

 Thanks,

 -james

 __
 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.




 --
 Jim Holtman
 Cincinnati, OH
 +1 513 646 9390

 What is the problem that you are trying to solve?




__
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 return/show the indexes of unusual points in boxplot?

use 'match' to get the indices;

match(x$out, yourData)


On Wed, Sep 9, 2009 at 4:07 PM, g...@ucalgary.ca wrote:
 Thank you very much.
 But it seems that
 x$out returns the values not the indexes of the  values (1,22).
 -james

 boxplot returns a dataframe that has the values in it at $out:

 x - boxplot(c(4,rnorm(20),8))
 x
 $stats
            [,1]
 [1,] -1.5364498
 [2,] -0.5282799
 [3,] -0.1398736
 [4,]  0.3065579
 [5,]  1.3430388

 $n
 [1] 22

 $conf
            [,1]
 [1,] -0.4210947
 [2,]  0.1413474

 $out
 [1] 4 8

 $group
 [1] 1 1

 $names
 [1] 1



 On Wed, Sep 9, 2009 at 3:25 PM, g...@ucalgary.ca wrote:
 I am wondering if you know how to return by function or show in
 boxplot,
 the indexes of unusual points, such as,
 points that are outside the box or in [Q3+1.5IQR,Max].
 For example,

 boxplot(c(4,rnorm(20),8))

 There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8
 in
 the boxplot
 or return them by function?

 Thanks,

 -james

 __
 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.




 --
 Jim Holtman
 Cincinnati, OH
 +1 513 646 9390

 What is the problem that you are trying to solve?









-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem that you are trying to solve?

__
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 return/show the indexes of unusual points in boxplot?

2009-09-09 Thread Jorge Ivan Velez

Hi James,
If that's the case, then

set.seed(123)
x - c(4,rnorm(20),8)
bp - boxplot(x)
bp$out
# [1] 4 8
which(x %in% bp$out)
# [1]  1 22

is what you are looking for.  :-)

HTH,
Jorge


On Wed, Sep 9, 2009 at 4:07 PM, g...@ucalgary.ca wrote:

 Thank you very much.
 But it seems that
 x$out returns the values not the indexes of the  values (1,22).
 -james

  boxplot returns a dataframe that has the values in it at $out:
 
  x - boxplot(c(4,rnorm(20),8))
  x
  $stats
 [,1]
  [1,] -1.5364498
  [2,] -0.5282799
  [3,] -0.1398736
  [4,]  0.3065579
  [5,]  1.3430388
 
  $n
  [1] 22
 
  $conf
 [,1]
  [1,] -0.4210947
  [2,]  0.1413474
 
  $out
  [1] 4 8
 
  $group
  [1] 1 1
 
  $names
  [1] 1
 
 
 
  On Wed, Sep 9, 2009 at 3:25 PM, g...@ucalgary.ca wrote:
  I am wondering if you know how to return by function or show in
  boxplot,
  the indexes of unusual points, such as,
  points that are outside the box or in [Q3+1.5IQR,Max].
  For example,
 
  boxplot(c(4,rnorm(20),8))
 
  There are 2 unusual points 4 and 8. How to show the indexes of 4 and 8
  in
  the boxplot
  or return them by function?
 
  Thanks,
 
  -james
 
  __
  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.
 
 
 
 
  --
  Jim Holtman
  Cincinnati, OH
  +1 513 646 9390
 
  What is the problem that you are trying to solve?
 
 
 

 __
 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] Stats help with calculating between and within subject variance and confidence intervals

2009-09-09 Thread Paul


Hello.
I'm trying to find a way in R to calculate between and within subject 
variances and confidence intervals for some analytical method 
development data.


I've found a reference to a method in Burdick, R. K.  Graybill, F. A. 
1992, Confidence Intervals on variance components, CRC Press.  This 
example is for Balanced Data confidence interval calculation from Pg 
62.  The data are fill weights from bottles sampled randomly from a 
sample of four filling machines.  There are 12 values, and the 
confidence intervals are for 1-2a = 95%.  I  have got the same results 
as the book but using slightly different fomulae (see variables for H1, 
G1 and H12 and G12).


I'd appreciate any help, and any comments on whether their is a better 
way to do this.


Thanks

Paul.

 BGBottles
  Machine weight
11  14.23
21  14.96
31  14.85
42  16.46
52  16.74
62  15.94
73  14.98
83  14.88
93  14.87
10   4  15.94
11   4  16.07
12   4  14.91

 BGaov-aov(weight~Machine,data=BGBottles)

 summary(BGaov)
   Df Sum Sq Mean Sq F value   Pr(F)   
Machine  3 5.3294  1.7765  9.7702 0.004733 **
Residuals8 1.4546  0.1818
---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 BGaov
Call:
  aov(formula = weight ~ Machine, data = BGBottles)

Terms:
Machine Residuals
Sum of Squares  5.329425  1.454600
Deg. of Freedom3 8

Residual standard error: 0.4264094
Estimated effects may be unbalanced

 BGaov-aov(weight~Machine+Error(Machine),data=BGBottles)

 summary(BGaov)

Error: Machine
   Df Sum Sq Mean Sq
Machine  3 5.3294  1.7765

Error: Within
 Df  Sum Sq Mean Sq F value Pr(F)
Residuals  8 1.45460 0.18182   


 BGaov

Call:
aov(formula = weight ~ Machine + Error(Machine), data = BGBottles)

Grand Mean: 15.4025

Stratum 1: Machine

Terms:
Machine
Sum of Squares  5.329425
Deg. of Freedom3

Estimated effects may be unbalanced

Stratum 2: Within

Terms:
   Residuals
Sum of Squares 1.4546
Deg. of Freedom 8

Residual standard error: 0.4264094

 confint-95

 a-(1-(confint/100))/2

 #str(BGaov) # get structure of BGAOV object

 #str(summary(BGaov)) # get structure of BGaov summary

 #summary(aov.1.e)[[2]][[1]]$Df # Could also use this syntax

 grandmean-as.vector(BGaov$(Intercept)[[1]][1]) # Grand Mean (I think)

 within-summary(BGaov)$Error: Within[[1]]$Mean Sq  # S2^2Mean 
Square Value for Within Machine = 0.1819


 dfMachine-summary(BGaov)$Error: Machine[[1]]$Df  # DF for within = 3

 dfWithin-summary(BGaov)$Error: Within[[1]]$Df  # DF for within = 8

 machine-summary(BGaov)$Error: Machine[[1]]$Mean Sq # S1^2Mean 
Square for Machine


 between-(machine-within)/((dfWithin/(dfMachine+1))+1) # (S1^2-S2^2)/J

 total-between+within

 between # Between Run Variance
[1] 0.53155

 within # Within Run Variance
[1] 0.181825

 total # Total Variance
[1] 0.713375

 betweenCV-sqrt(between)/grandmean * 100 # Between Machine CV%

 withinCV-sqrt(within)/grandmean * 100 # Within Machine CV%

 totalCV-sqrt(total)/grandmean * 100 # Total CV%

 #within confidence intervals (Chi Squared Method)

 withinLCB-within/qf(1-a,8,Inf) # Within LCB

 withinUCB-within/qf(a,8,Inf) # Within UCB

 #Between Confidence Intervals (Modified Large Sample Method)

 n1-dfMachine

 n2-dfWithin

 G1-1-(1/qf(1-a,n1,Inf)) # According to Burdick and Graybill this 
should be a


 G2-1-(1/qf(1-a,n2,Inf))

 H1-(1/qf(a,n1,Inf))-1  # and this should be 1-a, but my results 
don't agree


 H2-(1/qf(a,n2,Inf))-1

 
G12-((qf(1-a,n1,n2)-1)^2-(G1^2*qf(1-a,n1,n2)^2)-(H2^2))/qf(1-a,n1,n2) # 
again, should be a, not 1-a


 H12-((1-qf(a,n1,n2))^2-H1^2*qf(a,n1,n2)^2-G2^2)/qf(a,n1,n2) # again, 
should be 1-a, not a


 Vu-H1^2*machine^2+G2^2*within^2+H12*machine*within

 Vl-G1^2*machine^2+H2^2*within^2+G12*within*machine

 betweenLCB-(machine-within-sqrt(Vl))/J # Betwen LCB

 betweenUCB-(machine-within+sqrt(Vu))/J # Between UCB

 #Total Confidence Intervals (Graybill-Wang Method)

 y-(machine+(J-1)*within)/J

 totalLCB-y-(sqrt(G1^2*machine^2+G2^2*(J-1)^2*within^2)/J) # Total LCB

 totalUCB-y+(sqrt(H1^2*machine^2+H2^2*(J-1)^2*within^2)/J) # Total UCB

 result-data.frame(Name=c(within, between, 
total),CV=c(withinCV,betweenCV,totalCV),LCB=c(sqrt(withinLCB)/grandmean*100,sqrt(betweenLCB)/grandmean*100,sqrt(totalLCB)/grandmean*100),UCB=c(sqrt(withinUCB)/grandmean*100,sqrt(betweenUCB)/grandmean*100,sqrt(totalUCB)/grandmean*100))


 result
Name   CV  LCB   UCB
1  within 2.768443 1.869964  5.303702
2 between 4.733483 2.123816 18.551051
3   total 5.483625 3.590745 18.772389

__
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 return/show the indexes of unusual points in boxplot?

2009-09-09 Thread guox

Thanks for your help.
-jmaes
 Hi James,
 If that's the case, then

 set.seed(123)
 x - c(4,rnorm(20),8)
 bp - boxplot(x)
 bp$out
 # [1] 4 8
 which(x %in% bp$out)
 # [1]  1 22

 is what you are looking for.  :-)

 HTH,
 Jorge


 On Wed, Sep 9, 2009 at 4:07 PM, g...@ucalgary.ca wrote:

 Thank you very much.
 But it seems that
 x$out returns the values not the indexes of the  values (1,22).
 -james

  boxplot returns a dataframe that has the values in it at $out:
 
  x - boxplot(c(4,rnorm(20),8))
  x
  $stats
 [,1]
  [1,] -1.5364498
  [2,] -0.5282799
  [3,] -0.1398736
  [4,]  0.3065579
  [5,]  1.3430388
 
  $n
  [1] 22
 
  $conf
 [,1]
  [1,] -0.4210947
  [2,]  0.1413474
 
  $out
  [1] 4 8
 
  $group
  [1] 1 1
 
  $names
  [1] 1
 
 
 
  On Wed, Sep 9, 2009 at 3:25 PM, g...@ucalgary.ca wrote:
  I am wondering if you know how to return by function or show in
  boxplot,
  the indexes of unusual points, such as,
  points that are outside the box or in [Q3+1.5IQR,Max].
  For example,
 
  boxplot(c(4,rnorm(20),8))
 
  There are 2 unusual points 4 and 8. How to show the indexes of 4 and
 8
  in
  the boxplot
  or return them by function?
 
  Thanks,
 
  -james
 
  __
  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.
 
 
 
 
  --
  Jim Holtman
  Cincinnati, OH
  +1 513 646 9390
 
  What is the problem that you are trying to solve?
 
 
 

 __
 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] Combining simulated data

2009-09-09 Thread KABELI MEFANE

Thank you very much, now i can proceed to my next task. THANK YOU 

--- On Wed, 9/9/09, Schalk Heunis schalk.heu...@enerweb.co.za wrote:

From: Schalk Heunis schalk.heu...@enerweb.co.za
Subject: Re: [R] Combining simulated data
To: KABELI MEFANE kabelimef...@yahoo.co.uk
Cc: R-help@r-project.org
Date: Wednesday, 9 September, 2009, 8:45 PM

Kabeli

I think this is what you want:

Hypermarket - matrix(rnorm(10, mean=2, sd=7000))
Supermarket - matrix(rnorm(15, mean=12000, sd=4000))
Minimarket   - matrix(rnorm(20, mean=1, sd=4000))
Cornershop  - matrix(rnorm(20, mean= 8000, sd=3000))
Spazashop   - matrix(rnorm(35, mean= 7000, sd=3000))
rbind(Hypermarket,Supermarket,Minimarket,Cornershop,Spazashop)

This creates a matrix with one column and 100 rows.  You can add names
also by doing something like this for e.g. Spazashop
Spazashop   - matrix(rnorm(35, mean= 7000, sd=3000), dimnames =
list(paste(Spazashop,1:35)))

Hope this helps.

Schalk Heunis

On Wed, Sep 9, 2009 at 2:26 AM, KABELI MEFANEkabelimef...@yahoo.co.uk wrote:
 R helpers

 Please help me combine the simulated data to a form of table where: 
 Hypermarket have 10 rows, supermarket have 15 rows,..., spazashops with 
 35 rows.

 Hypermarket - rnorm(10, mean=2, sd=7000)
 Supermarket - rnorm(15, mean=12000, sd=4000)
 Minimarket   - rnorm(20, mean=1, sd=4000)
 Cornershop  - rnorm(20, mean= 8000, sd=3000)
 Spazashop   - rnorm(35, mean= 7000, sd=3000)

 i am trying to write a code to simulate data such that i have 10% of 
 hypermarkets, 15% of supermarkets,etc.

        [[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] Joining Characters in R {issue with paste}

2009-09-09 Thread milton ruser

You not tryed ?paste  :-)

paste(a,b,sep=)

bests

milton

On Wed, Sep 9, 2009 at 5:08 PM, Abhishek Pratap abhishek@gmail.comwrote:

 Hi Guys
 I am want to join to strings in R. I am using paste but not getting
 desirable result.

 For the sake of clarity, a quick example:

  a=Bio
  b=iology
  paste(a,b)
 [1] Bio iology


 *There is a SPACE in the word biology which is what I dont want *


 Thanks,
 -Abhi

[[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.htmlhttp://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] Joining Characters in R {issue with paste}

2009-09-09 Thread Erik Iverson

And did you read the help file, ?paste , paying attention to the arguments and 
their descriptions, specifically the sep argument?  Presumably, you want, 

paste(a, b, sep = )

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Abhishek Pratap
Sent: Wednesday, September 09, 2009 4:09 PM
To: r-help@r-project.org
Subject: [R] Joining Characters in R {issue with paste}

Hi Guys
I am want to join to strings in R. I am using paste but not getting
desirable result.

For the sake of clarity, a quick example:

 a=Bio
 b=iology
 paste(a,b)
[1] Bio iology

*There is a SPACE in the word biology which is what I dont want *

Thanks,
-Abhi

[[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] Joining Characters in R {issue with paste}

2009-09-09 Thread Abhishek Pratap

I did try ?paste and paste(a,b,separator=). same result
Thanks,
-Abhi

On Wed, Sep 9, 2009 at 5:11 PM, milton ruser milton.ru...@gmail.com wrote:

 You not tryed ?paste  :-)

 paste(a,b,sep=)

 bests

 milton

 On Wed, Sep 9, 2009 at 5:08 PM, Abhishek Pratap abhishek@gmail.comwrote:

 Hi Guys
 I am want to join to strings in R. I am using paste but not getting
 desirable result.

 For the sake of clarity, a quick example:

  a=Bio
  b=iology
  paste(a,b)
 [1] Bio iology


 *There is a SPACE in the word biology which is what I dont want *


 Thanks,
 -Abhi

[[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.htmlhttp://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] Joining Characters in R {issue with paste}

2009-09-09 Thread Abhishek Pratap

guys I am sorry may be it is that part of the day where my eyes are not able
to read.
using separator instead of sep.

Thanks a lot,
-Abhi

On Wed, Sep 9, 2009 at 5:12 PM, Erik Iverson eiver...@nmdp.org wrote:

 And did you read the help file, ?paste , paying attention to the arguments
 and their descriptions, specifically the sep argument?  Presumably, you
 want,

 paste(a, b, sep = )

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
 On Behalf Of Abhishek Pratap
 Sent: Wednesday, September 09, 2009 4:09 PM
 To: r-help@r-project.org
 Subject: [R] Joining Characters in R {issue with paste}

 Hi Guys
 I am want to join to strings in R. I am using paste but not getting
 desirable result.

 For the sake of clarity, a quick example:

  a=Bio
  b=iology
  paste(a,b)
 [1] Bio iology


 *There is a SPACE in the word biology which is what I dont want *


 Thanks,
 -Abhi

 [[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] Joining Characters in R {issue with paste}

2009-09-09 Thread Clint Bowman


?paste and look at sep.

On Wed, 9 Sep 2009, Abhishek Pratap wrote:


Hi Guys
I am want to join to strings in R. I am using paste but not getting
desirable result.

For the sake of clarity, a quick example:


a=Bio
b=iology
paste(a,b)

[1] Bio iology


*There is a SPACE in the word biology which is what I dont want *


Thanks,
-Abhi

[[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.



--
Clint BowmanINTERNET:   cl...@ecy.wa.gov
Air Quality Modeler INTERNET:   cl...@math.utah.edu
Department of Ecology   VOICE:  (360) 407-6815
PO Box 47600FAX:(360) 407-7534
Olympia, WA 98504-7600

__
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] Joining Characters in R {issue with paste}

2009-09-09 Thread Rolf Turner



On 10/09/2009, at 9:08 AM, Abhishek Pratap wrote:


Hi Guys
I am want to join to strings in R. I am using paste but not getting
desirable result.

For the sake of clarity, a quick example:


a=Bio
b=iology
paste(a,b)

[1] Bio iology


*There is a SPACE in the word biology which is what I dont want *


(1) RTFM

(2) paste(a,b,sep=)

cheers,

Rolf Turner

##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

__
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] Stats help with calculating between and within subject variance and confidence intervals

2009-09-09 Thread Bert Gunter

Paul:

If these data are real -- or at least a reasonable facsimile -- then even
though the machines might be considered random -- i.e. a sample from a
potential population of machines -- there are too few of them to get a
reasonable estimate of their variance. Better to treat them as fixed and
just do a standard oneway anova with a single within = residual error
component.

Incidentally, this situation is quite common in gauge RR = analytical
method development studies. While the problem is unavoidable -- you only
have so many machines or operators or whatever -- it is unfortunate that
standard statistical references do not point out that it's just basically
silly to try to estimate variance components with so few df, the resulting
confidence intervals, as here, being so wide as to be useless (reflecting
the inadequate information).

Incidentally, a better way to get at this when there are sufficient df is
via the lme() function in the nlme package -- it will work with unbalanced
data and not just in the balanced data situation. But there would be a
considerable learning curve required, I realize.

Cheers,

Bert Gunter
Genentech Nonclinical Biostatistics

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Paul
Sent: Wednesday, September 09, 2009 1:38 PM
To: R-help@r-project.org
Subject: [R] Stats help with calculating between and within subject variance
and confidence intervals

Hello.
I'm trying to find a way in R to calculate between and within subject 
variances and confidence intervals for some analytical method 
development data.

I've found a reference to a method in Burdick, R. K.  Graybill, F. A. 
1992, Confidence Intervals on variance components, CRC Press.  This 
example is for Balanced Data confidence interval calculation from Pg 
62.  The data are fill weights from bottles sampled randomly from a 
sample of four filling machines.  There are 12 values, and the 
confidence intervals are for 1-2a = 95%.  I  have got the same results 
as the book but using slightly different fomulae (see variables for H1, 
G1 and H12 and G12).

I'd appreciate any help, and any comments on whether their is a better 
way to do this.

Thanks

Paul.

  BGBottles
   Machine weight
11  14.23
21  14.96
31  14.85
42  16.46
52  16.74
62  15.94
73  14.98
83  14.88
93  14.87
10   4  15.94
11   4  16.07
12   4  14.91

  BGaov-aov(weight~Machine,data=BGBottles)

  summary(BGaov)
Df Sum Sq Mean Sq F value   Pr(F)   
Machine  3 5.3294  1.7765  9.7702 0.004733 **
Residuals8 1.4546  0.1818
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  BGaov
Call:
   aov(formula = weight ~ Machine, data = BGBottles)

Terms:
 Machine Residuals
Sum of Squares  5.329425  1.454600
Deg. of Freedom3 8

Residual standard error: 0.4264094
Estimated effects may be unbalanced

  BGaov-aov(weight~Machine+Error(Machine),data=BGBottles)

  summary(BGaov)

Error: Machine
Df Sum Sq Mean Sq
Machine  3 5.3294  1.7765

Error: Within
  Df  Sum Sq Mean Sq F value Pr(F)
Residuals  8 1.45460 0.18182   

  BGaov

Call:
aov(formula = weight ~ Machine + Error(Machine), data = BGBottles)

Grand Mean: 15.4025

Stratum 1: Machine

Terms:
 Machine
Sum of Squares  5.329425
Deg. of Freedom3

Estimated effects may be unbalanced

Stratum 2: Within

Terms:
Residuals
Sum of Squares 1.4546
Deg. of Freedom 8

Residual standard error: 0.4264094

  confint-95

  a-(1-(confint/100))/2

  #str(BGaov) # get structure of BGAOV object

  #str(summary(BGaov)) # get structure of BGaov summary

  #summary(aov.1.e)[[2]][[1]]$Df # Could also use this syntax

  grandmean-as.vector(BGaov$(Intercept)[[1]][1]) # Grand Mean (I think)

  within-summary(BGaov)$Error: Within[[1]]$Mean Sq  # S2^2Mean 
Square Value for Within Machine = 0.1819

  dfMachine-summary(BGaov)$Error: Machine[[1]]$Df  # DF for within = 3

  dfWithin-summary(BGaov)$Error: Within[[1]]$Df  # DF for within = 8

  machine-summary(BGaov)$Error: Machine[[1]]$Mean Sq # S1^2Mean 
Square for Machine

  between-(machine-within)/((dfWithin/(dfMachine+1))+1) # (S1^2-S2^2)/J

  total-between+within

  between # Between Run Variance
[1] 0.53155

  within # Within Run Variance
[1] 0.181825

  total # Total Variance
[1] 0.713375

  betweenCV-sqrt(between)/grandmean * 100 # Between Machine CV%

  withinCV-sqrt(within)/grandmean * 100 # Within Machine CV%

  totalCV-sqrt(total)/grandmean * 100 # Total CV%

  #within confidence intervals (Chi Squared Method)

  withinLCB-within/qf(1-a,8,Inf) # Within LCB

  withinUCB-within/qf(a,8,Inf) # Within UCB

  #Between Confidence Intervals (Modified Large Sample Method)

  n1-dfMachine

  n2-dfWithin

  G1-1-(1/qf(1-a,n1,Inf)) # According to Burdick and Graybill this 
should be a

Re: [R] Recursion is slow

2009-09-09 Thread Martin Morgan

Hi Bryan --

Bryan Keller wrote:
 Bill,
 
 Thanks for the great tips for speeding up functions in your response.  Those 
 are really useful for me.  Even with the improvements the recursion is still 
 many times slower than I need it to be in order to make it useful in R.  I 
 may just have to suck it up and call a compiled language.

Probably you've moved on, but it struck me that identical values of A()
are being calculated repeatedly, so one can perhaps 'memoize' them
(record that they've already been calculated). I'm not sure that I have
this exactly right, or that I've memoized in the right place, but this

Cmemo - function(T, m) {
C - function(lt, m) {
if (lt == 1L) {
R - as.numeric(0 = m  m = T[1])
} else {
R - 0
for (u in 0:T[lt]) {
if (is.na(memo[lt-1L, offset + m-u]))
memo[lt-1L, offset + m-u] - C(lt-1L, m-u)
R - R + memo[lt-1L, offset + m-u]
}
}
R
}
offset - sum(T[-1]) - m + 1L
nrow - length(T) - 1L
memo - matrix(rep(NA_real_, nrow * (offset + m)), nrow=nrow)
C(length(T), m)
}

seems to give consistent answers quickly (A1 is a tidied version of your
original recursion, B is Bill Dunlap's)

 system.time(ansA1 - A1( c(20,40,35,21,31), 100))
   user  system elapsed
 12.144   0.008  12.202
 system.time(ansB - B( c(20,40,35,21,31), 100))
   user  system elapsed
  0.272   0.000   0.270
 system.time(ansC - Cmemo( c(20,40,35,21,31), 100))
   user  system elapsed
  0.064   0.000   0.066
 all.equal(ansA1, ansB)
[1] TRUE
 all.equal(ansA1, ansC)
[1] TRUE

It's a little memory inefficient, but seems to scale fairly well.

Martin

 
 Bryan
 
 -
 Bryan Keller, Doctoral Student/Project Assistant
 Educational Psychology - Quantitative Methods
 The University of Wisconsin - Madison
 
 - Original Message -
 From: William Dunlap wdun...@tibco.com
 Date: Thursday, September 3, 2009 6:41 pm
 Subject: RE: [R] Recursion is slow
 To: Bryan Keller bskel...@wisc.edu, r-help@r-project.org
 
 Bill Dunlap
 TIBCO Software Inc - Spotfire Division
 wdunlap tibco.com  

 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] On Behalf Of Bryan Keller
 Sent: Thursday, September 03, 2009 2:11 PM
 To: r-help@r-project.org
 Subject: [R] Recursion is slow

 The following recursion is about 120 times faster in C#.  I 
 know R is not known for its speed with recursions but I'm 
 wondering if anyone has a tip about how to speed things up in R.

 #T is a vector and m is a number between 1 and sum(T)

 A - function(T,m) {
 lt - length(T)

 if (lt == 1) {
 if (0 = m  m = T[1]) {
 return(1)
 } else {
 return(0)
 }
 }
 R - 0
 for (u in 0:T[lt]) {
 R - (R+(A(T[1:(lt-1)],(m-u
 }
 return(R)
 }

 For starters, remove all repeated calculations from
 the for loop.  Then vectorize the length(T)==2 case
 to avoid a lot of calls to the scalar case.  Finally,
 I noticed that the answer was independent of the
 order of T and that putting it in reverse sorted order
 was the faster order, so I did that.  I use default
 values of arguments so you don't have to supply
 derived values when you start but the recursions
 don't have to recompute some things.

 B - function(T,m,lt=length(T), Tsorted = rev(sort(T))) {
if (lt == 1L) {
   R - as.integer(m = Tsorted  0L = m)
} else if (lt == 2L) {
   mu - m - (0:Tsorted[2L])
   R - sum(mu = Tsorted[1L]  0L = mu)
} else {
   R - 0L
   lt1 - lt-1L
   T1 - Tsorted[1:lt1]
   for (mu in m-(0:Tsorted[lt])) {
  R - R + B(unused, mu, lt1, T1)
   }
}
R
 }

 E.g.,

 system.time(A( c(20,40,35,21,31), 100))
user  system elapsed 
   13.230.08   12.93 
 system.time(B( c(20,40,35,21,31), 100))
user  system elapsed 
0.350.000.33 
 A( c(20,40,35,21,31), 100)
 [1] 193363
 B( c(20,40,35,21,31), 100)
 [1] 193363

 Bill Dunlap
 TIBCO Software Inc - Spotfire Division
 wdunlap tibco.com 

 -
 Bryan Keller, Doctoral Student/Project Assistant
 Educational Psychology - Quantitative Methods
 The University of Wisconsin - Madison

 __
 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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide

[R] Xyplot, multi line title via main, all lines left justified

2009-09-09 Thread Afshartous, David


All,

Below is an xyplot plot with multiple panels and a title produced via main:

library(lattic)
data.ex = data.frame(y = rnorm(10), t = rep(1:5, 2), group = rep(c(0,1),
each = 5))

xyplot(y ~ t | as.factor(group), data = data.ex,
main = list(Put figure caption here xx
want this line left justified ))

I must be mis-interpreting the help description for main under xyplot, and
the descriptions of just, hjust, and vjust in textGrob, since manipulation
of these arguments do not seem to produce the desired result of having the
second line of text left-justified.

Any help much appreciated.

David

__
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] DOE in R?

2009-09-09 Thread Carlos J. Gil Bellosta

Hello,

I think I misread the link you sent me yesterday. In any case, the
reason why SAS generates a model with 30 trials whereas R produces
another one with 25 is that, if the parameter for the number of trials
is not specified in the proc/function call, both systems adhere to
different defaults: number of cols in the design matrix plus 10 for SAS
and plus 5 for R. You can check proc optex and optFederov documentation
for details.

Beyond that, SAS implements several optimization methods, among which,
Federov's. R implements just (one version of) Federov's. Finally, the
non-deterministic nature of the search for optimality (and local optima)
may lead to discrepancies between the outputs. You will have to work
through the different proc/function inputs to make sure you are running
an analogous algorithm on both systems.

Best regards,

Carlos J. Gil Bellosta
http://www.datanalytics.com




On Sun, 2009-09-06 at 20:33 -0400, b miner wrote:
 I attempted to use the package algdesign. I used the following code.
 However, the results were very much not matching the reference I noted
 (which is located at
 http://www2.sas.com/proceedings/sugi31/196-31.pdf). Instead of 30
 design points, I received 25 and those that were returned, only half
 or so matched the reference. I am inexperienced with optimal designs
 so I dont know if I am doing something wrong, the package is not the
 correct one for the task or a combination. Here is the code in case
 anyone has any input:
 
 #CODE:
 
 runif(1) #for a bug in program (assumes random seed object exists)
 
 dat-gen.factorial(c(3,3,3,3,2),varNames=c(Intro,Duration,GOTO,Fee,Color))
 dat #show design plan
 
 desD-optFederov(~Intro+Duration+GOTO+Fee+Color
 +quad(Intro)+quad(Duration)+quad(GOTO)+quad(Fee)+Intro*Duration
 +Intro*GOTO+Intro*Fee+Intro*Color+Duration*GOTO+Duration*Fee
 +Duration*Color+GOTO*Fee+GOTO*Color
 +Fee*Color,dat,crit=D,maxIteration=1000,eval=TRUE)
 
 #D
 desD$D
 
 #design
 desD$design
 design-desD$design
 
 
 
 
  Subject: Re: [R] DOE in R?
  From: c...@datanalytics.com
  To: b_mi...@live.com
  CC: r-help@r-project.org
  Date: Sun, 6 Sep 2009 14:57:36 +0200
  
  Hello,
  
  This is your starting point:
  
  http://cran.r-project.org/web/views/ExperimentalDesign.html
  
  Best regards,
  
  Carlos J. Gil Bellosta
  http://www.datanalytics.com
  
  
  On Thu, 2009-09-03 at 17:38 -0700, B_miner wrote:
   Hello!
   
   
   This is not a topic I am well versed in but required to become
 well versed
   in...I welcome any assistance!
   
   Using R, I want to create an optimal design for an experiment.
 I'll be
   analyzing the results with logistic regression or some generalized
 linear
   model. I am thinking that the algdesign package can help (but no
 idea where
   to start?).
   
   I'm presenting an example here that I have seen the answer to (in
 SAS) in
   order to make sure I would have gotten it *right*.
   
   There are 5 factors: 4 are quantitative with three levels each and
 1 is
   qualitative with two levels.
   
   Factor and levels:
   
   Intro: 0, 1.99, 2.99
   Duration: 6, 9 ,12
   GOTO: 3.99, 4.99, 5.99
   Fee: 0, 15, 45
   Color: Red, White
   
   
   In order to screen these factors, I would want to get a design
 where I could
   evaluate all main effects, all first order interactions and the
 squared
   terms of Intro, Duration, GOTO and FEE (for example Intro*Intro).
   
   Looking for the D-optimal design.
   
   Is this something that R can provide?
   
   
   These are, according to the SAS paper I read the following:
   Obs intro duration goto fee color
   1 0.00 6 3.99 0 WHITE
   2 0.00 6 3.99 45 RED
   3 0.00 6 5.99 0 RED
   4 0.00 6 5.99 45 WHITE
   5 0.00 9 3.99 45 RED
   6 0.00 9 4.99 15 WHITE
   7 0.00 9 5.99 0 RED
   8 0.00 12 3.99 0 RED
   9 0.00 12 3.99 45 WHITE
   10 0.00 12 5.99 0 WHITE
   11 0.00 12 5.99 45 RED
   12 0.00 12 5.99 45 WHITE
   13 1.99 6 3.99 15 RED
   14 1.99 6 4.99 45 WHITE
   15 1.99 6 5.99 0 WHITE
   16 1.99 9 5.99 45 RED
   17 1.99 12 3.99 0 WHITE
   18 1.99 12 5.99 15 RED
   19 2.99 6 3.99 0 WHITE
   20 2.99 6 3.99 45 WHITE
   21 2.99 6 4.99 0 RED
   22 2.99 6 5.99 15 WHITE
   23 2.99 6 5.99 45 RED
   24 2.99 9 3.99 0 RED
   25 2.99 12 3.99 15 WHITE
   26 2.99 12 3.99 45 RED
   27 2.99 12 4.99 0 WHITE
   28 2.99 12 4.99 45 RED
   29 2.99 12 5.99 0 RED
   30 2.99 12 5.99 45 WHITE
   
   
  
  
  
  
 
 
 __
 Windows Live: Make it easier for your friends to see what you’re up to
 on Facebook. Find out more.

__
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] Recursion is slow

2009-09-09 Thread William Dunlap

Note that the memoization optimization is orthogonal
to the simpler optimizations I suggested (don't repeat
calculations in loops, sort the input vector, do a little
math to avoid so many recursions to the leaves).  You
can combine them to get

Cmemo1 - function(T, m) {
C - function(lt, m) {
if (lt == 1L) {
R - as.numeric(0 = m  m = T[1])
} else if (lt == 2L) {
mu - m - (0:T[2L])
R - sum(mu = T[1L]  0L = mu)
} else {
R - 0
lt1 - lt-1L
for (mu in m-(0:T[lt])) {
if (is.na(memo[lt1, offset + mu]))
memo[lt-1L, offset + mu] - C(lt1, mu)
R - R + memo[lt1, offset + mu]
}
}
R
}
T - rev(sort(T))
m - as.integer(m)
offset - sum(T[-1]) - m + 1L
nrow - length(T) - 1L
memo - matrix(rep(NA_real_, nrow * (offset + m)), nrow=nrow)
C(length(T), m)
}

 system.time(z1-Cmemo1(c(10,30,20,45,3,100,21,37,50,20,37,46,90),409))
   user  system elapsed 
   0.440.000.37 
 system.time(z-Cmemo(c(10,30,20,45,3,100,21,37,50,20,37,46,90),409))
   user  system elapsed 
   1.240.011.19 
 identical(z,z1)
[1] TRUE

Bill Dunlap
TIBCO Software Inc - Spotfire Division
wdunlap tibco.com  

 -Original Message-
 From: Bryan Keller [mailto:bskel...@wisc.edu] 
 Sent: Wednesday, September 09, 2009 3:25 PM
 To: Martin Morgan
 Cc: William Dunlap; r-help@r-project.org
 Subject: Re: [R] Recursion is slow
 
 Thanks Martin!  It seems to be right on and it is blazing 
 fast!  I greatly appreciate the responses from you and Bill 
 both.  As a beginning user of R it is really helpful to be 
 able to compare my code with yours and try to figure out your tricks.
 
 Bryan
 
 -
 Bryan Keller, Doctoral Student/Project Assistant
 Educational Psychology - Quantitative Methods
 The University of Wisconsin - Madison
 
 - Original Message -
 From: Martin Morgan mtmor...@fhcrc.org
 Date: Wednesday, September 9, 2009 4:40 pm
 Subject: Re: [R] Recursion is slow
 To: Bryan Keller bskel...@wisc.edu
 Cc: William Dunlap wdun...@tibco.com, r-help@r-project.org
 
 
  Hi Bryan --
  
  Bryan Keller wrote:
   Bill,
   
   Thanks for the great tips for speeding up functions in your 
  response.  Those are really useful for me.  Even with the 
 improvements 
  the recursion is still many times slower than I need it to 
 be in order 
  to make it useful in R.  I may just have to suck it up and call a 
  compiled language.
  
  Probably you've moved on, but it struck me that identical 
 values of A()
  are being calculated repeatedly, so one can perhaps 'memoize' them
  (record that they've already been calculated). I'm not sure 
 that I have
  this exactly right, or that I've memoized in the right 
 place, but this
  
  Cmemo - function(T, m) {
  C - function(lt, m) {
  if (lt == 1L) {
  R - as.numeric(0 = m  m = T[1])
  } else {
  R - 0
  for (u in 0:T[lt]) {
  if (is.na(memo[lt-1L, offset + m-u]))
  memo[lt-1L, offset + m-u] - C(lt-1L, m-u)
  R - R + memo[lt-1L, offset + m-u]
  }
  }
  R
  }
  offset - sum(T[-1]) - m + 1L
  nrow - length(T) - 1L
  memo - matrix(rep(NA_real_, nrow * (offset + m)), nrow=nrow)
  C(length(T), m)
  }
  
  seems to give consistent answers quickly (A1 is a tidied 
 version of your
  original recursion, B is Bill Dunlap's)
  
   system.time(ansA1 - A1( c(20,40,35,21,31), 100))
 user  system elapsed
   12.144   0.008  12.202
   system.time(ansB - B( c(20,40,35,21,31), 100))
 user  system elapsed
0.272   0.000   0.270
   system.time(ansC - Cmemo( c(20,40,35,21,31), 100))
 user  system elapsed
0.064   0.000   0.066
   all.equal(ansA1, ansB)
  [1] TRUE
   all.equal(ansA1, ansC)
  [1] TRUE
  
  It's a little memory inefficient, but seems to scale fairly well.
  
  Martin
  
   
   Bryan
   
   -
   Bryan Keller, Doctoral Student/Project Assistant
   Educational Psychology - Quantitative Methods
   The University of Wisconsin - Madison
   
   - Original Message -
   From: William Dunlap wdun...@tibco.com
   Date: Thursday, September 3, 2009 6:41 pm
   Subject: RE: [R] Recursion is slow
   To: Bryan Keller bskel...@wisc.edu, r-help@r-project.org
   
   Bill Dunlap
   TIBCO Software Inc - Spotfire Division
   wdunlap tibco.com  
  
   -Original Message-
   From: r-help-boun...@r-project.org 
   [mailto:r-help-boun...@r-project.org] On Behalf Of Bryan Keller
   Sent: Thursday, September 03, 2009 2:11 PM
   To: r-help@r-project.org
   Subject: [R] Recursion is slow
  
   The following recursion is about 120 times faster in C#.  I 
   know R is not known for its speed with recursions but I'm 
   wondering if anyone has a tip about how to speed things up in R.
  
   #T is a vector and m is a number between 1 and sum(T)

Re: [R] Help on percentage of random numbers for different classes

2009-09-09 Thread KABELI MEFANE

R-list
 
I am sorry for asking this stupid question, but i have been running in circles. 
I want to randomly generate a scaling point of between 1 and 10, for say 
hundred entries, where the first 10% percent is has rates between 2 and 7, the 
next 15% 3 and 7, 20% between 3 and 9, 20% between 3 and  10, 35% between 5 and 
10. The problem is that i can only generate the usual 100 using runif function
 
  y-c(ceiling(10*runif(100)))
 y
  [1] 10  8  5  2  4  1  6  7  1  6  8  8  8  9  7  7  8  8  2  7  3 10  1  7  1
 [26] 10  4  8  8  8  9  3  7  8  4  6  7  2  3  1  9  8  2  6  7  4  8  8  9  7
 [51]  6  5  4  1  8  7  9  8 10  5  3  7  5  5  4  4  7  4 10  4  9  1  5 10 10
 [76]  5  5 10  7  3  4  4  9 10  6  2  6  6  6  3  8  2  2  4  4 10  6  9  4  3

I just want to try to avoid small numbers as much as possible. I am open to 
suggestions, please please please.
 
Kabeli


  
[[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] Very basic question regarding plot.Design...

2009-09-09 Thread Frank E Harrell Jr

I just added a new option nlines to plot.Predict to make it easy to get 
interaction line plots for categorical predictors.  To get the new 
version type the following after running require(rms):


source('http://biostat.mc.vanderbilt.edu/tmp/plot.Predict.s')

Here are some examples:

require(rms)
n - 100
set.seed(1)
gender - c(rep('male', n), rep('female',n))
m - sample(c('a','b'), 2*n, TRUE)
d -  datadist(gender, m); options(datadist='d')
anxiety - runif(2*n) + .2*(gender=='female') + .4*(gender=='female'  
m=='b')

tapply(anxiety, llist(gender,m), mean)
f - ols(anxiety ~ gender*m)
p - Predict(f, gender=., m=.)
plot(p) # horizontal dot chart; usually preferred for categorical 
predictors

Key(.5, .5)
plot(p, ~gender, groups='m', nlines=TRUE)
plot(p, ~m, groups='gender', nlines=TRUE)
plot(p, ~gender|m, nlines=TRUE)


Frank

Petar Milin wrote:

Sorry, I hope this will be the last (for now, at least).

Following your advices, I did:
 n - 100
 anxiety - c(rnorm(n, 10, 2.5), rnorm(n, 26, 3.2),
  rnorm(100, 25, 3.1), rnorm(100, 10, 2.6))
 m.status - c(rep('married', n*2), rep('not married', n*2))
 gender  - c(rep('male', n), rep('female', n),
  rep('male', n), rep('female',n))

 im - as.integer(p2$m.status)
 m.status - as.factor(m.status)
 gender - as.factor(gender)

 require(rms)

 d - datadist(gender, m.status); options(datadist='d')

 ols1 - ols(anxiety ~ m.status)
 p1 - Predict(ols1, m.status=.)
 ols2 - ols(anxiety ~ m.status * gender)
 p2 - Predict(ols2, m.status=., gender=.)

 pdf('figs/anova.pdf', height=6, width=8)
 par(mfrow=c(1,2), mar=c(4,4,1,1), cex=1.2)
 with(p1, plot(1:2, yhat, type='l', xlab='marital status',
  ylab='anxiety', ylim=c(5,30), lwd=1.2))
 xyplot(yhat ~ im, groups=gender,
  type='l', data=p2, ylim=c(5,30), xlab='marital status',
  ylab='anxiety', scales=list(x=list(at=1:2,
  labels=levels(p2$m.status))), col=c('black','black'),
  lwd=1.2, label.curve=list(offset=unit(.3,in)))
 par(mfrow=c(1,1))
 dev.off()

Graphs are great! However, now, they does not appear in one panel.

Sorry, again, for asking so many questions. I am trying to wrap this.

Best,
PM




--
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.

Re: [R] Help on percentage of random numbers for different classes

2009-09-09 Thread Nordlund, Dan (DSHS/RDA)

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
 Behalf Of KABELI MEFANE
 Sent: Wednesday, September 09, 2009 4:06 PM
 To: R-help@r-project.org
 Subject: Re: [R] Help on percentage of random numbers for different classes

 R-list

 I am sorry for asking this stupid question, but i have been running in 
 circles. I want
 to randomly generate a scaling point of between 1 and 10, for say hundred 
 entries,
 where the first 10% percent is has rates between 2 and 7, the next 15% 3 
 and 7,
 20% between 3 and 9, 20% between 3 and  10, 35% between 5 and 10. The
 problem is that i can only generate the usual 100 using runif function

   y-c(ceiling(10*runif(100)))
  y
   [1] 10  8  5  2  4  1  6  7  1  6  8  8  8  9  7  7  8  8  2  7  3 10  1  
 7  1
  [26] 10  4  8  8  8  9  3  7  8  4  6  7  2  3  1  9  8  2  6  7  4  8  8  
 9  7
  [51]  6  5  4  1  8  7  9  8 10  5  3  7  5  5  4  4  7  4 10  4  9  1  5 10 
 10
  [76]  5  5 10  7  3  4  4  9 10  6  2  6  6  6  3  8  2  2  4  4 10  6  9  
 4  3

 I just want to try to avoid small numbers as much as possible. I am open to
 suggestions, please please please.

 Kabeli

If I understand you correctly, this might do what you want.  n is a vector with 
the number samples you want in each range, x is then minimum for each range and 
y is the maximum.  The function, s, samples from a given range, a specified 
number of times.  mapply applies the function using the first, second, ... 
elements in turn returning a list with the samples.  

n - c(10,15,20,20,35)
x - c(2,3,3,3,5)
y - c(7,7,9,10,10)

s - function(mn, mx, n) {sample(mn:mx, n, replace=TRUE)}
unlist(mapply(s, x, y, n))

Hope this is helpful,

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA  98504-5204

__
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] Read.csv in R with dynamic file (1st) argument

2009-09-09 Thread Steven Kang

Dear R users,


I have numerous data sets (csv files) saved in the folder which has the same
name as individual data.
(i.e data x1 saved in x1 folder, data x2 in x2 folder etc)

I would like to read in the desired data set name using 'scan' function and
assign this inputted value to an object so that it can be used in the
'read.csv' function.

For example,

x - scan()
1: 0708
2:
Read 1 item

dat - read.csv(D://R//Data//x//x.csv, head=TRUE, sep = ,)
Error in file(file, r) : cannot open the connection
In addition: Warning message:
In file(file, r) :
  cannot open file ('D://R//Data//x//x.csv': No such file or directory


 In SAS, this can be done by assigning x as a macro variable, is there an
equivalent function in R?

Moreover, is there any way of using the inputted variable which starts with
number 0 (i.e 0102) in the 'read.csv' function?

Appreciate for providing your expertise in resolving this problem.

[[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] Read.csv in R with dynamic file (1st) argument

Try this:

read.csv(sprintf(D://R//Data//%04d//%04d.csv, x, x), header = TRUE)

On Wed, Sep 9, 2009 at 9:32 PM, Steven Kang stochastick...@gmail.comwrote:

 Dear R users,


 I have numerous data sets (csv files) saved in the folder which has the
 same
 name as individual data.
 (i.e data x1 saved in x1 folder, data x2 in x2 folder etc)

 I would like to read in the desired data set name using 'scan' function and
 assign this inputted value to an object so that it can be used in the
 'read.csv' function.

 For example,

 x - scan()
 1: 0708
 2:
 Read 1 item

 dat - read.csv(D://R//Data//x//x.csv, head=TRUE, sep = ,)
 Error in file(file, r) : cannot open the connection
 In addition: Warning message:
 In file(file, r) :
  cannot open file ('D://R//Data//x//x.csv': No such file or directory


  In SAS, this can be done by assigning x as a macro variable, is there an
 equivalent function in R?

 Moreover, is there any way of using the inputted variable which starts with
 number 0 (i.e 0102) in the 'read.csv' function?

 Appreciate for providing your expertise in resolving this problem.

[[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.




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

[[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] inline function in apply

2009-09-09 Thread bwgoudey


I've been trying to filling in the missing variables in a matrix with the
mean from the column they are in. So the following code does this just fine. 

#3X3 matrix where the middle number is missing. 
data=matrix(c(1:4,NA,6:9), 3, 3)

#replace missing values in an vector with the mean of the vector
fill_in_NA-function (x)
{ 
x[is.na(x)]-sum(x[!is.na(x)])/length(x[!is.na(x)])
return(x)
}
#replace the missing value with 5 (mean of 4 and 6)
apply(data, 2, fill_in_NA)


I'm curious as to whether or not I can reduce the function with a single
inline function call (I'm aware that it will be less readable). My initial
thought was something like

apply(data, 2, function(x)
(x[is.na(x)]-sum(x[!is.na(x)])/length(x[!is.na(x)])))

but this returns a single vector. The problem is that the x in my inline
function doesn't seem to refer to what I thought it did. Could anyway one
suggest some appropriate code or possible provide me with a better
understanding of what my current inline function is actually doing?

Thanks in advance


 
-- 
View this message in context: 
http://www.nabble.com/inline-function-in-apply-tp25375733p25375733.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] inline function in apply