Re: [R] write.xls dont find the object in function

2011-12-23 Thread Ronaldo Reis Júnior
Em 20-12-2011 14:09, MacQueen, Don escreveu:
 Or:

 require(xlsx)
 test- function(x){
 +a- data.frame(A=c(1,2),B=c(10,11))
 +write.xlsx(a,file=a.xlsx)
 +  }
   test()
 list.files(patt='xlsx')
 [1] a.xlsx




Thanks,

with the write.xlsx all work.

Inte
Ronaldo

-- 
3ª lei - Na investigação e outros assuntos, o seu orientador está sempre
  certo, na maior parte do tempo.

   --Herman, I. P. 2007. Following the law. NATURE, Vol 445, p. 228.

  Prof. Ronaldo Reis Júnior
|  .''`. UNIMONTES/DBG/Lab. Ecologia Comportamental e Computacional
| : :'  : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia
| `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil
|   `- Fone: (38) 3229-8192 | ronaldo.r...@unimontes.br
| http://www.ppgcb.unimontes.br/lecc | LinuxUser#: 205366


[[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] simple ggplot2 question

2011-12-23 Thread Albert-Jan Roskam
Hello,

I am trying to make a plot using the code below. The plot is divided into two 
parts, using facet_grid. I would like the vertical axis (labelled 'place') to 
be different for each location (=part). So in the upper part, only places 'n' 
through 'z' are shown, while in the lower part, only places 'a' through 'm' are 
shown. I thought 'free_y' would do the trick. I also tried converting variable 
place into class 'factor'.


require(ggplot2)
DF - data.frame(place=letters, value=runif(26), location=c(rep(1, 13), rep(0, 
13)))
qplot(data=DF, x=place, y=value, geom=bar, stat=identity) + 
  coord_flip() + 
  geom_abline(intercept=35, slope=0, colour=red) +
  facet_grid(location ~ ., scales=free_y)
R.version.string # R version 2.10.1 (2009-12-14)


Thank you in advance  merry xmas!
 
Cheers!!
Albert-Jan


~~
All right, but apart from the sanitation, the medicine, education, wine, public 
order, irrigation, roads, a fresh water system, and public health, what have 
the Romans ever done for us?
~~
[[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] GIS operations

2011-12-23 Thread Antonio Rodriges
Hello,

Maybe this will help:

There is rgeos package
with gContains, gCovers, etc.
It is a full implementation of Java Topology Suite
http://cran.r-project.org/web/packages/rgeos/rgeos.pdf

-- 
Kind regards,
Antonio Rodriges

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 Source Code Request Office For National Statistics UK

2011-12-23 Thread Barry Rowlingson
On Wed, Dec 21, 2011 at 1:06 PM, David Samuel
david.sam...@ons.gsi.gov.uk wrote:
 To R Support Team,

 Could you please let me know if ONS would be allowed to view your
 software's source code?

 If proof was ever needed that people don't read software license
agreements, this is it.

 Is there a license that forces you to read the source code before
using the software? I'd license my stuff under that:

THE BAZ PUBLIC LICENSE (hereinafter the BPL)

Under the terms of the BPL you MUST first read all source code and get
a pretty good understanding of what it all does before using THE
SOFTWARE. You MUST report bugs or problems to THE AUTHOR but you MUST
NOT expect them to be fixed. You MUST fix them YOURSELF if you think
you may be technically capable or if YOU think YOU are a better
programmer than THE AUTHOR. Otherwise YOU must SHUT UP.


Barry (hereinbefore Baz)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 transforming a dist

2011-12-23 Thread Taral
Perfect! combn was the trick I needed. Although I'll probably rbind
the stuff instead of building a new frame. :)

On Thu, Dec 22, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com wrote:
 I'm not at all sure what you mean with your matrix: is that supposed to
 be three columns? What about:

 fakedata - data.frame(X=runif(3), Y = runif(3))
 rownames(fakedata) - c(A, B, C)
 dist(fakedata)
          A         B
 B 0.8617733
 C 0.3813032 0.5124284
 data.frame(t(combn(rownames(fakedata), 2)), dist=as.vector(dist(fakedata)))
  X1 X2      dist
 1  A  B 0.8617733
 2  A  C 0.3813032
 3  B  C 0.5124284

 Sarah

 On Thu, Dec 22, 2011 at 7:03 PM, Taral tar...@gmail.com wrote:
 I'd like to convert a dist into a table/matrix/thingy of the form:

 A B value
 A C value
 B C value

 (assuming the dist was over 3 names).

 Is there a way to do this without using a for loop?

 --
 Taral tar...@gmail.com
 Please let me know if there's any further trouble I can give you.
     -- Unknown



 --
 Sarah Goslee
 http://www.functionaldiversity.org



-- 
Taral tar...@gmail.com
Please let me know if there's any further trouble I can give you.
    -- Unknown

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] GIS operations

2011-12-23 Thread Paul Hiemstra
On 12/20/2011 12:54 PM, Antonio Rodriges wrote:
 Hello,

 Is there a way to find points in SpatialPoints which lie inside a given 
 Polygon?

Hi,

Take a look at the over function from the sp-package.

Paul

-- 
Paul Hiemstra, Ph.D.
Global Climate Division
Royal Netherlands Meteorological Institute (KNMI)
Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39
P.O. Box 201 | 3730 AE | De Bilt
tel: +31 30 2206 494

http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Axis manipulation in Stackpoly (Plotrix)

2011-12-23 Thread Jim Lemon

On 12/23/2011 07:52 AM, Ben Neal wrote:


Have made a stacked area plot, but now want to manipulate the x axis: make in 
even increments (50) in order to suppress the tick marks (forming solid bar 
under plot). However, the plot functions do not seem to work, and I cannot find 
documentation for use of xat in stackpoly.

This is a time series of data covering 579 years, and has been successfully 
converted from zoo to matrix. Data plotting fine, jsut want to change axis -

CaribMatrix-as.matrix(ts4Ex)
stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons, 
xaxlab=(1250:2008),
   main=Historical Florida reef fisheries catch by sector, 
axis4=FALSE)

Any assistance appreciated! Thanks very much, Benjamin P Neal, Scripps Inst. of 
Oceanogrpahy


Hi Ben,
In the stackpoly function, you can specify the positions of the axis 
ticks (xat=...) and the labels for those positions (xaxlab=...). You 
don't seem to have supplied the x values, so the values in CaribMatrix 
will be placed display at 1:dim(CaribMatrix)[1]. I'll have to make up 
CaribMatrix...


CaribMatrix-matrix(order(sample(200:500,579,TRUE),nrow=579,ncol=3)
stackpoly(CaribMatrix,stack=TRUE,ylim=c(0,1500),xat=seq(1,600,by=50),
 xaxlab=seq(1450,2000,by=50),staxx=TRUE)

This seems to give a reasonable plot for me.

Jim

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


Re: [R] simple ggplot2 question

2011-12-23 Thread Paul Hiemstra
On 12/23/2011 08:26 AM, Albert-Jan Roskam wrote:
 Hello,

 I am trying to make a plot using the code below. The plot is divided into two 
 parts, using facet_grid. I would like the vertical axis (labelled 'place') to 
 be different for each location (=part). So in the upper part, only places 'n' 
 through 'z' are shown, while in the lower part, only places 'a' through 'm' 
 are shown. I thought 'free_y' would do the trick. I also tried converting 
 variable place into class 'factor'.


 require(ggplot2)
 DF - data.frame(place=letters, value=runif(26), location=c(rep(1, 13), 
 rep(0, 13)))
 qplot(data=DF, x=place, y=value, geom=bar, stat=identity) + 
   coord_flip() + 
   geom_abline(intercept=35, slope=0, colour=red) +
   facet_grid(location ~ ., scales=free_y)
 R.version.string # R version 2.10.1 (2009-12-14)

Hi,

This code using facet_wrap works:

require(ggplot2)
DF - data.frame(place=letters, value=runif(26), location=c(rep(1, 13),
rep(0, 13)))
qplot(data=DF, x=place, y=value, geom=bar, stat=identity) +
  coord_flip() +
  facet_wrap(~ location, scales = free, ncol = 1)

And using facet_grid:

DF - data.frame(place=letters, value=runif(26), location=c(rep(1, 13),
rep(0, 13)))
qplot(data=DF, x=place, y=value, geom=bar, stat=identity) +
  coord_flip() +
  facet_grid(. ~ location, scales = free)

Note that I switched the formula in facet_grid. However, I would expect
your code to also work.

cheers and hope this helps,
Paul



 Thank you in advance  merry xmas!
  
 Cheers!!
 Albert-Jan


 ~~
 All right, but apart from the sanitation, the medicine, education, wine, 
 public order, irrigation, roads, a fresh water system, and public health, 
 what have the Romans ever done for us?
 ~~
   [[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.


-- 
Paul Hiemstra, Ph.D.
Global Climate Division
Royal Netherlands Meteorological Institute (KNMI)
Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39
P.O. Box 201 | 3730 AE | De Bilt
tel: +31 30 2206 494

http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770


[[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] Axis manipulation in Stackpoly (Plotrix)

2011-12-23 Thread Jim Lemon

On 12/23/2011 09:53 PM, Jim Lemon wrote:


CaribMatrix-matrix(order(sample(200:500,579,TRUE),nrow=579,ncol=3)


Oops, that should be

CaribMatrix-matrix(sort(sample(200:500,579,TRUE),nrow=579,ncol=3)

Jim

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


[R] missing value where TRUE/FALSE needed

2011-12-23 Thread Michael Pearmain
Merry Xmas to all,

I am writing a function and curiously this runs sometimes on one data set
and fails on another and i cannot figure out why.
Any help much appreciated.

If i run the code below with
data - iris[ ,1:4]
The code runs fine, but if i run on a large dataset i get the following
error (showing data structures as matrix is large)

 str(cluster.data)
 num [1:9985, 1:811] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:811] 1073949105 1073930585 1073843224 1073792624 ...
#(This is intended to be chr)
 str(iris)
'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 : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 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 : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1 1 1
1 1 1 ...
 str(as.matrix(iris[,1:4]))
 num [1:150, 1:4] 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:4] Sepal.Length Sepal.Width Petal.Length Petal.Width

n.cols - ncol(data)
  n.rows - nrow(data)
  X - as.matrix(data)
  stepsize - 0.05
  c1 - (2 * pi) ** (n.cols / 2)
  c2 - n.rows * (smoothing ** (n.cols + 2))
  c3 - n.rows * (smoothing ** n.cols)

  Kexp - function(sqs){
return (exp((-1 * sqs) / (2 * smoothing ** 2)))
  }

  FindGradient - function(x){
XmY - t(x - t(X))
sqsum - rowSums(XmY * XmY)
K - sapply(sqsum, Kexp)
dens - ((c1 * c3) ** -1) * sum(K)
grad - -1 * ((c1 * c2) ** -1) * colSums(K * XmY)
return (list(gradient = grad,
 density = dens))
  }

  attractors - matrix(0, n.rows, n.cols)
  densities - matrix(0, n.rows)


 density.attractors -
sapply(rep(1:n.rows), function(i) {
  notconverged - TRUE
  # For each row loop through and find the attractor and density value.
  x - (X[i, ])
  iters - as.integer(1)
  # Run gradient ascent for each point to obtain x*
  while(notconverged == TRUE) {
find.gradient - FindGradient(x)
next.x - x + stepsize * find.gradient$gradient
change - sqrt(sum((next.x - x) * (next.x - x)))
notconverged - ifelse(change  tol, TRUE, FALSE)
x - next.x
iters - iters + 1
  }

  # store the attractor and density value
  return(c(densities[i, ] - find.gradient$density,
   attractors[i, ] - x))
})

Error in while (notconverged == TRUE) { :
  missing value where TRUE/FALSE needed


Any help would be great

Mike

[[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] Problem understanding behaviour of mmap package

2011-12-23 Thread Wolfgang Wu
I am trying to solve the following problem using the mmap package:
 
1.) Store some data as a memory mapped file that resides on the disk. 
2.) Open the data and add some new data (or change it) to that file.


1.) works fine using the following code:

cFile -tempfile()
dfrData - data.frame(matrix(data=1, ncol=2, nrow=10))
mmapData - as.mmap(dfrData, file=cFile)
stcData - mmapData$storage.mode
munmap(mmapData)


2.) Now I want to link mmapAfter to the file that should still reside as cFile, 
or? What is the problem here? The error I get is at the second line of the 
code: 
Error in if (!x$bytes) stop(no data to extract) :   missing value where 
TRUE/FALSE needed. 

mmapAfter - mmap(cFile, mode=stcData, extractFUN=as.data.frame)
mmapAfter[];
mmapAfter [1,2] - 2;
mmapAfter [];
munmap(mmapAfter)


I also noted when I execute the following code that 
the  mmapData$filedesc doesn't show the filename I specified in cFile. Is this 
not how it is supposed to be used?


cFile - 'C:/temp/testmapfile'
dfrData - data.frame(matrix(data=1, ncol=2, nrow=10))
mmapData - as.mmap(dfrData, file=cFile)
mmapData$filedesc

Thanks for your help.

Regards,

 
Wolfgang Wu

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lasso based selection for mixed model

2011-12-23 Thread schell
Dear Dieter

There has been some effort examining Lasso-type estimators for mixed models.
It is more involved than with other type of models.

We have studied Gaussian and generalized linear mixed models including a
Lasso-type penalty for the fixed effects and we are now able to provide two
packages: lmmlasso and glmmlasso. Both are available from R-Forge (the
former even from CRAN).

Andreas Groll  provides the glmmLasso package, which is available from CRAN.

Juerg Schelldorfer

--
View this message in context: 
http://r.789695.n4.nabble.com/lasso-based-selection-for-mixed-model-tp887308p4228495.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] ggplot2: behaviour with empty datasets

2011-12-23 Thread Casper Ti. Vector
For example, prepare like this
 df.0 - data.frame(x = 0, y = 0, note = 1)
 df.1 - subset(df.0, note == 1)
 df.2 - subset(df.0, note == 2)

Then a call to
 ggplot() + aes(x = x, y = y) +
   geom_point(data = df.1) + geom_point(data = df.2)
produces the error
 Error in eval(expr, envir, enclos) : object 'x' not found

And a call to
 ggplot(df.1, aes(x = x, y = y, shape = 4)) +
   geom_point() + geom_point(aes(shape = 1), df.2)
plots both a cross (shape = 4) and circle (shape = 1) on position
(0, 0). However, as expected, only a cross should be plotted because
`dt.2' is empty.

So with the current behaviour, if someone writes a script to plot
datasets automatically and distinguish between subsets, he will need
to write different code path for empty and non-empty data subsets.
I think ggplot2 can just choose to plot nothing for empty data
subsets, and only produce error when the dataset (union of a subsets) to
plot is completely empty. This can make scripting easier based on
ggplot2.

-- 
Using GPG/PGP? Please get my current public key (ID: 0xAEF6A134,
valid from 2010 to 2013) from a key server.



signature.asc
Description: Digital signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] missing value where TRUE/FALSE needed

2011-12-23 Thread jim holtman
Does this look similar to the error you are getting:

 while(NA == TRUE) 1
Error in while (NA == TRUE) 1 : missing value where TRUE/FALSE needed

SO 'notconverged' is probably equal to NA.  BTW, what is the value of
'tol'; I do not see it defined.  So when computing 'notconverged' you
have generated an NA.  You can test it to see if this is true.

You can use the following command:

options(error=utils::recover)

and then learn how to use the 'browser' to examine variables when the
error occurs.

On Fri, Dec 23, 2011 at 5:44 AM, Michael Pearmain
michael.pearm...@gmail.com wrote:
 Merry Xmas to all,

 I am writing a function and curiously this runs sometimes on one data set
 and fails on another and i cannot figure out why.
 Any help much appreciated.

 If i run the code below with
 data - iris[ ,1:4]
 The code runs fine, but if i run on a large dataset i get the following
 error (showing data structures as matrix is large)

 str(cluster.data)
  num [1:9985, 1:811] 0 0 0 0 0 0 0 0 0 0 ...
  - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:811] 1073949105 1073930585 1073843224 1073792624 ...
 #(This is intended to be chr)
 str(iris)
 '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 : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 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     : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1 1 1
 1 1 1 ...
 str(as.matrix(iris[,1:4]))
  num [1:150, 1:4] 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
  - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:4] Sepal.Length Sepal.Width Petal.Length Petal.Width

 n.cols - ncol(data)
  n.rows - nrow(data)
  X - as.matrix(data)
  stepsize - 0.05
  c1 - (2 * pi) ** (n.cols / 2)
  c2 - n.rows * (smoothing ** (n.cols + 2))
  c3 - n.rows * (smoothing ** n.cols)

  Kexp - function(sqs){
    return (exp((-1 * sqs) / (2 * smoothing ** 2)))
  }

  FindGradient - function(x){
    XmY - t(x - t(X))
    sqsum - rowSums(XmY * XmY)
    K - sapply(sqsum, Kexp)
    dens - ((c1 * c3) ** -1) * sum(K)
    grad - -1 * ((c1 * c2) ** -1) * colSums(K * XmY)
    return (list(gradient = grad,
                 density = dens))
  }

  attractors - matrix(0, n.rows, n.cols)
  densities - matrix(0, n.rows)


 density.attractors -
    sapply(rep(1:n.rows), function(i) {
      notconverged - TRUE
      # For each row loop through and find the attractor and density value.
      x - (X[i, ])
      iters - as.integer(1)
      # Run gradient ascent for each point to obtain x*
      while(notconverged == TRUE) {
        find.gradient - FindGradient(x)
        next.x - x + stepsize * find.gradient$gradient
        change - sqrt(sum((next.x - x) * (next.x - x)))
        notconverged - ifelse(change  tol, TRUE, FALSE)
        x - next.x
        iters - iters + 1
      }

      # store the attractor and density value
      return(c(densities[i, ] - find.gradient$density,
               attractors[i, ] - x))
    })

 Error in while (notconverged == TRUE) { :
  missing value where TRUE/FALSE needed


 Any help would be great

 Mike

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



-- 
Jim Holtman
Data Munger Guru

What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.

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


Re: [R] Help transforming a dist

2011-12-23 Thread Sarah Goslee
On Fri, Dec 23, 2011 at 1:06 AM, Taral tar...@gmail.com wrote:
 Perfect! combn was the trick I needed. Although I'll probably rbind
 the stuff instead of building a new frame. :)

You should try rbind() with the small example I gave before you use it
on your data. You'll immediately see why I didn't. :)

 On Thu, Dec 22, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com wrote:
 I'm not at all sure what you mean with your matrix: is that supposed to
 be three columns? What about:

 fakedata - data.frame(X=runif(3), Y = runif(3))
 rownames(fakedata) - c(A, B, C)
 dist(fakedata)
          A         B
 B 0.8617733
 C 0.3813032 0.5124284
 data.frame(t(combn(rownames(fakedata), 2)), dist=as.vector(dist(fakedata)))
  X1 X2      dist
 1  A  B 0.8617733
 2  A  C 0.3813032
 3  B  C 0.5124284

 Sarah

 On Thu, Dec 22, 2011 at 7:03 PM, Taral tar...@gmail.com wrote:
 I'd like to convert a dist into a table/matrix/thingy of the form:

 A B value
 A C value
 B C value

 (assuming the dist was over 3 names).

 Is there a way to do this without using a for loop?


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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: behaviour with empty datasets

2011-12-23 Thread Hadley Wickham
See https://github.com/hadley/ggplot2/issues/31 - I totally agree that
it's annoying.
Hadley

PS.  You are more likely to get helpful responses about ggplot2 on the
ggplot mailing list.

On Fri, Dec 23, 2011 at 7:08 AM, Casper Ti. Vector
caspervec...@gmail.com wrote:
 For example, prepare like this
 df.0 - data.frame(x = 0, y = 0, note = 1)
 df.1 - subset(df.0, note == 1)
 df.2 - subset(df.0, note == 2)

 Then a call to
 ggplot() + aes(x = x, y = y) +
   geom_point(data = df.1) + geom_point(data = df.2)
 produces the error
 Error in eval(expr, envir, enclos) : object 'x' not found

 And a call to
 ggplot(df.1, aes(x = x, y = y, shape = 4)) +
   geom_point() + geom_point(aes(shape = 1), df.2)
 plots both a cross (shape = 4) and circle (shape = 1) on position
 (0, 0). However, as expected, only a cross should be plotted because
 `dt.2' is empty.

 So with the current behaviour, if someone writes a script to plot
 datasets automatically and distinguish between subsets, he will need
 to write different code path for empty and non-empty data subsets.
 I think ggplot2 can just choose to plot nothing for empty data
 subsets, and only produce error when the dataset (union of a subsets) to
 plot is completely empty. This can make scripting easier based on
 ggplot2.

 --
    Using GPG/PGP? Please get my current public key (ID: 0xAEF6A134,
 valid from 2010 to 2013) from a key server.


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




-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


[R] error from character NAs in indexing array

2011-12-23 Thread David Winsemius

The help page for ?Extract says:

When extracting, a numerical, logical or character NA index picks an  
unknown element and so returns NA in the corresponding element of a  
logical, integer, numeric, complex or character result, and NULL for a  
list. (It returns 00 for a raw result.] 


However, I seem to be getting different behavior when converting NA's  
to character indices for an array with character dimnames:


str(frailmed)
 num [1:73, 1:2, 1:2] -0.315 -0.31 -0.252 -0.311 -0.29 ...
 - attr(*, dimnames)=List of 3
  ..$ : chr [1:73] 18 19 20 21 ...
  ..$ : chr [1:2] Female Male
  ..$ : chr [1:2] FALSE TRUE

 frailmed[as.character(c(18,18,30)), Female, TRUE]
18 18 30
-0.3477757 -0.3477757 -0.3697707
 frailmed[as.character(c(18,18,30,NA)), Female, TRUE]
Error: subscript out of bounds

I started thinking I might need to use grep or match to convert the  
numeric value for ageLB to a numeric index for [, but after looking  
at the manual is seems I shouldn't be forced to do that.


 as.character(NA)
[1] NA
 frailmed[as.character(NA), Female, TRUE]
Error: subscript out of bounds

I'm guessing that NA_character_ is not being recognized by the  
extraction engine:


 frailmed[NA_character_, Female, TRUE]
Error: subscript out of bounds

Test subset:
test - structure(c(-0.314905119091182, -0.309721294756184,  
-0.252154582680236,
+ 0.218589733943967, 0.227724425528112, 0.219274402692938,  
-0.347775659359612,
+ -0.351245759870099, -0.361973673972107, 0.15524082878,  
0.149765287705234,

+ 0.150206947093079), .Dim = c(3L, 2L, 2L), .Dimnames = list(c(18,
+ 19, 20), c(Female, Male), c(FALSE, TRUE)))
 test
, , FALSE

   Female  Male
18 -0.3149051 0.2185897
19 -0.3097213 0.2277244
20 -0.2521546 0.2192744

, , TRUE

   Female  Male
18 -0.3477757 0.1552556
19 -0.3512458 0.1497653
20 -0.3619737 0.1502069

 test[NA, 1,1]
NA NA NA
  NA   NA   NA
 test[as.character(NA), 1,1]
Error: subscript out of bounds

Is this intended behavior?

I do see from that last successful execution that I may need to attack  
my problem by just processing complete cases, since the return of a  
vector of NA's with length  1 will also cause my logic to fail, but  
this behavior seems different than described.


sessionInfo()
R version 2.14.0 Patched (2011-11-13 r57650)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats graphics  grDevices utils datasets
[7] methods   base

other attached packages:
[1] rms_3.3-2Hmisc_3.9-0  survival_2.36-10 sos_1.3-1
[5] brew_1.0-6   lattice_0.20-0

loaded via a namespace (and not attached):
[1] cluster_1.14.1 grid_2.14.0tools_2.14.0


--

David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 RCOM package

2011-12-23 Thread Uwe Ligges



On 23.12.2011 11:02, KUMAR wrote:

Hi,

 I am using R version 2.14.0 on windows 7. Trying to use RCOM
package and followed example provided at
(http://cran.r-project.org/web/packages/rcom/rcom.pdf) .


Please ask on the mailing list that was set up for that package.
I wonder: Do you have statconnDCOM installed.


Uwe Ligges







 I am able to load the com object as seen in screenshot
2.png.  however I am unable to invoke/get/set property to the object.





Please guide me on the same,



Regards;

Kumar Anand



1.png



2.png




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ggplot2: behaviour with empty datasets

2011-12-23 Thread Casper Ti. Vector
Thanks, I'll post to that list if the encountered problem is just about
ggplot2.
And from the date of the issue on GitHub it seems that I might need to
manually work around the problem some more times :)

On Fri, Dec 23, 2011 at 08:10:05AM -0600, Hadley Wickham wrote:
 See https://github.com/hadley/ggplot2/issues/31 - I totally agree that
 it's annoying.
 Hadley
 
 PS.  You are more likely to get helpful responses about ggplot2 on the
 ggplot mailing list.

-- 
Using GPG/PGP? Please get my current public key (ID: 0xAEF6A134,
valid from 2010 to 2013) from a key server.



signature.asc
Description: Digital signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] nlrob problem

2011-12-23 Thread Martin Maechler
 Pascal A Niklaus pascal.nikl...@ieu.uzh.ch
 on Mon, 19 Dec 2011 20:46:57 +0100 writes:

 Dear all, I am not sure if this mail is for R-help or
 should be sent to R-devel instead, and therefore post to
 both.

 While using nlrob from package 'robustbase', I ran into
 the following problem:

 For psi-functions that can become zero
 (e.g. psi.bisquare), weights in the internal call to nls
 can become zero. Example:

d - data.frame(x=1:5,y=c(2,3,5,10,9))
d.nlrob - nlrob(y ~ k*x,start=list(k=1),data=d,psi=psi.bisquare)

 I think the problem is as follows: After the call to nls, a weighted 
 residual vector is calculated by dividing by sqrt(w). This generates 
 NaN's for weights of zero, which leads to problems in the subsequent 
 call to nls during the next iteration:

  for (iiter in 1:maxit) {
...
  w - psi(resid/Scale, ...)
...
  data$w - sqrt(w)
...
  out - nls( ., data=data ... )
...
  resid - -residuals(out)/sqrt(w)   # NaN for w=0
...
  }

 I wonder whether this problem shouldn't be dealt with by setting 'w' to 
 0 for the NaN cases.

 I can get a fit by calling nlrob with na.action=na.exclude, but I'd have 
 intuitively assumed that na.action applies to the NAs in the data set 
 passed to nlrob, not to the one internally generated by nlrob and passed 
 to nls.

You are right.
The next version of robustbase (0.8-1) will have this fixed.

Note however, that for me, in your example and in other more
interesting ones, as soon as I use a redescending  psi() function,
the IRLS iterations do *not* converge...
but rather strangely ``cycle'' around the true solution.
As a redescending psi() has a non-convex rho(), and non-linear
problems depend quite a bit on feasible/good starting
values, I currently think I would discourage from using such psi().


As others, some possibly more expert in robust non-linear fitting,
may be interested, and have interesting feedback,
I'm crossposting this reply to the R-SIG-robust  mailing
list.
(Further replies should typically only go there!)

 The second 'issue' is that the weights are passed as 'w', whereas the 
 documentation of 'nls' says weights should be given as 'weights' :

 data: an optional data frame in which to evaluate the variables in
‘formula’ and ‘weights’.  Can also be a list or an
environment, but not a matrix.

 I think it would be good to mention in the documentation of 'nls' that 
 weights can be passed as 'w' as well.

You have overlooked that nlrob uses a different *formula* for nls()
 --- nowadays unnecessarily ---
and that formula has  a 'w'.

So this part of your report is not correct.

Martin Maechler, ETH Zurich

 Pascal Niklaus

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


[R] data vector to corresonding percentile ranks

2011-12-23 Thread Steve Jones
I have a problem where I need to calculate the corresponding cohort 
percentile ranks for each of several variables.

Essentially, what I need is a function that will calculate the
distribution-free percentiles from each variable's data vector, returning a
corresponding vector of percentiles:

e.g.:  

percentile.my.data-/function/(my.data)


I tried to make ecdf() perform this task but was unsuccessful.

I'd be grateful for any help or advice...



--
View this message in context: 
http://r.789695.n4.nabble.com/data-vector-to-corresonding-percentile-ranks-tp4228971p4228971.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]

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


[R] Long jobs completing without output

2011-12-23 Thread Brendan Halpin
I've been running a glmer logit on a very large data set (600k obs). 

Running on a 10% subset works correctly, but for the complete data set,
R completes apparently without error, but does not display the results.
Given these jobs take about 200 hours, it's very hard to make progress
by trial and error.

I append the code and the sample and complete output. As is apparent, I
upgraded R during the complete run, but I recall testing on the
subsample with the earlier version too. I am also assuming that
upgrading R will not affect the running process -- is this true?

I'd be grateful for any leads. In the meantime I'll be running with
larger subsamples!

Regards,

Brendan Halpin


- code ---
library(arm)
library(foreign)
mlm - read.dta(../workingdata.dta)
attach(mlm)

gender - as.factor(stu_gend)

yr - year - 1998
failure - (lmer(fail ~
  1 + cao + subj1 + subj2 + subj3 + gender + yr + ageentry + 
as.factor(yrs5)
+ modsize  + meancao + depfemr + (1|deptno) + (1|modinst)  + 
(1|ulid) , 
  na.action = na.exclude, family = binomial (link=logit)))

display(failure, digits=5, detail=TRUE)
--

- output with 10% sample data 
R version 2.14.0 (2011-10-31)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

 library(arm)

arm (Version 1.4-13, built: 2011-6-19)
Working directory is /home/brendan/work/mlmmarks/genderECSR 
 library(foreign)
 mlm - read.dta(../worksample-random1.dta)
 attach(mlm)
 
 gender - as.factor(stu_gend)
 
 yr - year - 1998
 failure - (lmer(fail ~
+   1 + cao + subj1 + subj2 + subj3 + gender + yr + ageentry + 
as.factor(yrs5)
+ + modsize  + meancao + depfemr + (1|deptno) + (1|modinst)  + 
(1|ulid) , na.action = na.exclude, family = binomial (link=logit)))
 
 display(failure, digits=5, detail=TRUE)
glmer(formula = fail ~ 1 + cao + subj1 + subj2 + subj3 + gender + 
yr + ageentry + as.factor(yrs5) + modsize + meancao + depfemr + 
(1 | deptno) + (1 | modinst) + (1 | ulid), family = binomial(link = 
logit), 
na.action = na.exclude)
 coef.est  coef.se   z value   Pr(|z|) 
(Intercept)2.63826   0.97870   2.69568   0.00702
cao   -2.08963   0.11987 -17.43314   0.0
subj1  0.02608   0.23573   0.11064   0.91190
subj2 -0.55668   0.32759  -1.69932   0.08926
subj3 -1.57120   0.30664  -5.12400   0.0
genderM0.36368   0.09188   3.95845   0.8
yr 0.06067   0.01658   3.65996   0.00025
ageentry  -0.00720   0.04338  -0.16598   0.86817
as.factor(yrs5)1  -0.25181   0.05712  -4.40806   0.1
as.factor(yrs5)2  -0.54725   0.07601  -7.20005   0.0
as.factor(yrs5)3  -1.07483   0.08660 -12.41184   0.0
as.factor(yrs5)4  -1.22447   0.14373  -8.51932   0.0
as.factor(yrs5)5  -1.55032   0.31342  -4.94653   0.0
modsize0.03387   0.02533   1.33733   0.18112
meancao1.08747   0.10748  10.11780   0.0
depfemr   -1.49097   0.49350  -3.02122   0.00252

Error terms:
 Groups   NameStd.Dev.
 modinst  (Intercept) 1.14308 
 ulid (Intercept) 1.54030 
 deptno   (Intercept) 0.52497 
 Residual 1.0 
---
number of obs: 63254, groups: modinst, 9076; ulid, 2275; deptno, 26
AIC = 30275.2, DIC = 30237.2
deviance = 30237.2 
 
Loading required package: MASS
Loading required package: Matrix
Loading required package: lattice

Attaching package: ‘Matrix’

The following object(s) are masked from ‘package:base’:

det

Loading required package: lme4

Attaching package: ‘lme4’

The following object(s) are masked from ‘package:stats’:

AIC, BIC

Loading required package: R2WinBUGS
Loading required package: coda

Attaching package: ‘coda’

The following object(s) are masked from ‘package:lme4’:

HPDinterval

Loading required package: abind
Loading required package: foreign

Attaching package: ‘arm’

The following object(s) are masked from ‘package:coda’:

traceplot
--

- output with complete data --
R version 2.13.1 (2011-07-08)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

Re: [R] Help with code

2011-12-23 Thread 1Rnwb
this is how the ouput from the code should be 

structure(list(HTN = 1:10, HTN_FDR = structure(c(4L, 2L, 1L, 
2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc, 
T1D_oc, T1d_w), class = factor), Dyslipidemia = structure(c(3L, 
2L, 1L, 2L, 4L, 1L, 1L, 2L, 4L, 2L), .Label = c(Ctrl_noc, T1D_noc, 
T1D_oc, T1D_w), class = factor), CAD = structure(c(3L, 
2L, 1L, 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc, 
T1D_oc), class = factor), CAD_FDR = structure(c(3L, 2L, 1L, 
2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc, 
T1D_oc), class = factor), Prior_MI = structure(c(3L, 2L, 
1L, 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc, 
T1D_oc), class = factor), t1d_ptype = structure(c(3L, 3L, 
2L, 3L, 1L, 1L, 2L, 3L, 3L, 3L), .Label = c(Ctrl, Ctrl_FDR, 
T1D), class = factor)), .Names = c(HTN, HTN_FDR, Dyslipidemia, 
CAD, CAD_FDR, Prior_MI, t1d_ptype), class = data.frame, row.names
= c(NA, 
-10L))


--
View this message in context: 
http://r.789695.n4.nabble.com/Help-with-code-tp4218989p4228759.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] missing value where TRUE/FALSE needed

2011-12-23 Thread Michael Pearmain
Apologies, I was using top = 0.0001

I had looked at browser and did show notconverged = NA.  But I couldn't
understand why it worked for one and not the other?



On Friday, 23 December 2011, jim holtman jholt...@gmail.com wrote:
 Does this look similar to the error you are getting:

 while(NA == TRUE) 1
 Error in while (NA == TRUE) 1 : missing value where TRUE/FALSE needed

 SO 'notconverged' is probably equal to NA.  BTW, what is the value of
 'tol'; I do not see it defined.  So when computing 'notconverged' you
 have generated an NA.  You can test it to see if this is true.

 You can use the following command:

 options(error=utils::recover)

 and then learn how to use the 'browser' to examine variables when the
 error occurs.

 On Fri, Dec 23, 2011 at 5:44 AM, Michael Pearmain
 michael.pearm...@gmail.com wrote:
 Merry Xmas to all,

 I am writing a function and curiously this runs sometimes on one data set
 and fails on another and i cannot figure out why.
 Any help much appreciated.

 If i run the code below with
 data - iris[ ,1:4]
 The code runs fine, but if i run on a large dataset i get the following
 error (showing data structures as matrix is large)

 str(cluster.data)
  num [1:9985, 1:811] 0 0 0 0 0 0 0 0 0 0 ...
  - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:811] 1073949105 1073930585 1073843224 1073792624
...
 #(This is intended to be chr)
 str(iris)
 '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 : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 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 : Factor w/ 3 levels setosa,versicolor,..: 1 1 1 1 1
1 1
 1 1 1 ...
 str(as.matrix(iris[,1:4]))
  num [1:150, 1:4] 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
  - attr(*, dimnames)=List of 2
  ..$ : NULL
  ..$ : chr [1:4] Sepal.Length Sepal.Width Petal.Length
Petal.Width

 n.cols - ncol(data)
  n.rows - nrow(data)
  X - as.matrix(data)
  stepsize - 0.05
  c1 - (2 * pi) ** (n.cols / 2)
  c2 - n.rows * (smoothing ** (n.cols + 2))
  c3 - n.rows * (smoothing ** n.cols)

  Kexp - function(sqs){
return (exp((-1 * sqs) / (2 * smoothing ** 2)))
  }

  FindGradient - function(x){
XmY - t(x - t(X))
sqsum - rowSums(XmY * XmY)
K - sapply(sqsum, Kexp)
dens - ((c1 * c3) ** -1) * sum(K)
grad - -1 * ((c1 * c2) ** -1) * colSums(K * XmY)
return (list(gradient = grad,
 density = dens))
  }

  attractors - matrix(0, n.rows, n.cols)
  densities - matrix(0, n.rows)


 density.attractors -
sapply(rep(1:n.rows), function(i) {
  notconverged - TRUE
  # For each row loop through and find the attractor and density
value.
  x - (X[i, ])
  iters - as.integer(1)
  # Run gradient ascent for each point to obtain x*
  while(notconverged == TRUE) {
find.gradient - FindGradient(x)
next.x - x + stepsize * find.gradient$gradient
change - sqrt(sum((next.x - x) * (next.x - x)))
notconverged - ifelse(change  tol, TRUE, FALSE)
x - next.x
iters - iters + 1
  }

  # store the attractor and density value
  return(c(densities[i, ] - find.gradient$density,
   attractors[i, ] - x))
})

 Error in while (notconverged == TRUE) { :
  missing value where TRUE/FALSE needed


 Any help would be great

 Mike

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



 --
 Jim Holtman
 Data Munger Guru

 What is the problem that you are trying to solve?
 Tell me what you want to do, not how you want to do it.


[[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] Applying a function

2011-12-23 Thread Joanie Van De Walle

   Hi,

   I need help writing a function

   I capture seal pups mutliple times during the lactation season in order to
   monitor their growth rate. When I release them, the recovery (mother-pup)
   time is not the same for all individuals. I want to know if individuals that
   recover their mother the fastest are the ones with the highest growth rates.

   So, I noted at every release if the pup reovered his mother before we leave
   (yes or no). My dataframe looks like this

   Capture nb individual individual capture motherrecovery growth rate
   1  1  1  n  0.5
   2  1  2  y  0.5
   3  1  3  y  0.5
   4  1  4  y  0.5
   5  1  5  n  0.5
   6  2  1  y  0.3
   7  2  2  y  0.3
   8  3  1  y  0.4
   9  3  2  n  0.4
   10 3  3  y  0.4
   ...

   I want to calculate a rate of mother recovery by individual, i.e. nb of
   recoveries (y)/nb of captures. So for indivial 1 it would be 3/5 = 0.6. I
   want  to write a function that does this for all the individuals in my
   dataframe,  i.e. around 400 individuals (this is why I want to write a
   function, it would be too long by hand)

   Thank you,



   Joanie Van de Walle

   Étudiante à la Maîtrise en Biologie
   Département de Biologie
   1045, avenue de la Médecine Pavillon Vachon,
   bureau 2044 Université Laval
   Québec, Canada, G1V 0A6
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] nlrob problem

2011-12-23 Thread Pascal.Niklaus
Please ignore my previous posting, in the meanwhile I have realised that I
did not look carefully enough how nlrob works (irls), although this is
pretty obvious from the code.

The only issue that remains is the case of zero weight for redescending M
estimators, in which case NaN's are produced. I think it would be better if
one could avoid the NaN's produced by dividing by zero when psi becomes
zero. 

Something along the line of replacing 

   resid - -residuals(out)/sqrt(w)

by 

   resid - ifelse(w0,-residuals(out)/sqrt(w),0);

There probably is a more elegant way to do this. 

This would avoid that one has to use na.action=na.exclude for cases where
there are no NAs in the data set originally passed to nlrob (and the NAs
only occur 'internally' in the code of nlrob).

Hope I haven't missed something else this time...

Pascal Niklaus



--
View this message in context: 
http://r.789695.n4.nabble.com/nlrob-problem-tp4215473p4229016.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] missing value where TRUE/FALSE needed

2011-12-23 Thread jim holtman
Given that the maximum floating point value is:

$double.xmax
[1] 1.797693e+308


and the number you are trying to calculate is 5E323 you are exceeding
the size of numbers you can process.

Have a happy holiday and glad I could help.

On Fri, Dec 23, 2011 at 11:24 AM, Michael Pearmain
michael.pearm...@gmail.com wrote:
 Thanks for al your help with Jim,
 I realise this is really a case of debugging (and so falls to the code
 writer) using your help i have traced the NaN to be from the c1 calculation
 in the FindGradient function.
 there is a place where we use the calculation of c1...

 which at the start was calculated as
 c1 - (2*pi) ** (n.cols/2)

 when i use
 c1 - (2*pi) ** (811/2)
 [1] Inf

 So i think i have found a ceiling on the number of cols i can use to make a
 calcuation.

 Many thanks again for your help, especially on 23rd Dec!

 Mike



 On 23 December 2011 16:01, jim holtman jholt...@gmail.com wrote:

 The next step is to now evaluate the objects used in the calculation;
 e.g., look at the objects in the statement:

  change - sqrt(sum((next.x - x) * (next.x - x)))

 What do each of them show?  Most likely you have something in the data
 that causing this statement to evaluate to NaN (not a number) and this
 is caused by the values of one of these as not being what you think it
 is.  From the value of 'i', it looks like the first iteration (or at
 least when it has the value of 1).  Notice that X is defined as a
 matrix with NAs as the initial values and you are picking up the first
 row (X[i,]) which on the initial pass is an NA.  So start by examining
 the values of all the objects used up to that point and what the
 values of each expression is.

 On Fri, Dec 23, 2011 at 10:38 AM, Michael Pearmain
 michael.pearm...@gmail.com wrote:
  Thanks for guidance;
 
  the return from the dump is below:
 
  So i can see that the NaN is causing NA, as its evaluating to NA in the
  ifelse.
  So i guess my next question is why this happens for my large data set
  but
  not he iris dataset?
 
  Thank you again for your help on this, i really appreciate you taking
  the
  time to help
 
  1: DensitiesAndAttractors(data = cluster.data, smoothing = 0.34, tol =
  1e-05,
  2: denclue_par.R#136: sfSapply(rep(1:n.rows), function(i) {
  3: sapply(x, fun, ..., simplify = simplify, USE.NAMES = USE.NAMES)
  4: lapply(X, FUN, ...)
  5: FUN(1:9985[[1]], ...)
 
  Selection: 5
  Called from: function ()
  {
      if (.isMethodsDispatchOn()) {
          tState - tracingState(FALSE)
          on.exit(tracingState(tState))
      }
      calls - sys.calls()
      from - 0L
      n - length(calls)
      if (identical(sys.function(n), recover))
          n - n - 1L
      for (i in rev(seq_len(n))) {
          calli - calls[[i]]
          fname - calli[[1L]]
          if (!is.na(match(deparse(fname)[1L], c(methods::.doTrace,
              .doTrace {
              from - i - 1L
              break
          }
      }
      if (from == 0L)
          for (i in rev(seq_len(n))) {
              calli - calls[[i]]
              fname - calli[[1L]]
              if (!is.name(fname) || is.na(match(as.character(fname),
                  c(recover, stop, Stop {
                  from - i
                  break
              }
          }
      if (from  0L) {
          if (!interactive()) {
              try(dump.frames())
              cat(gettext(recover called non-interactively; frames
  dumped,
  use debugger() to view\n))
              return(NULL)
          }
          else if (identical(getOption(show.error.messages),
              FALSE))
              return(NULL)
          calls - limitedLabels(calls[1L:from])
          repeat {
              which - menu(calls, title = \nEnter a frame number, or 0
  to
  exit  )
              if (which)
                  eval(substitute(browser(skipCalls = skip), list(skip = 7
  -
                    which)), envir = sys.frame(which))
              else break
          }
      }
      else cat(gettext(No suitable frames for recover()\n))
  }()
  Browse[1] ls()
  [1] change        find.gradient i             iters
  [5] next.x        notconverged  x
  Browse[1] change
  [1] NaN
  Browse[1] i
  [1] 1
  Browse[1]
 
 
  On 23 December 2011 15:22, jim holtman jholt...@gmail.com wrote:
 
  According to your dump, if you look in section 5 (this is the function
  being call by 'sapply') you should do an 'ls()' and see the object
  'notconverged' and if you type its name, I am guessing that is that
  the value NA.  Also look at value of 'change' and 'tol' to determine
  how the 'ifelse' was evaluated.
 
  Let me know how it goes.
 
  On Fri, Dec 23, 2011 at 9:35 AM, Michael Pearmain
  michael.pearm...@gmail.com wrote:
   Hi Jim,
  
   Sorry to email directly, i wasn't sure what i was looking for, can
   you
   give
   me a little guidance then i can post back to the list?  Many thanks
  
   density.data - DensitiesAndAttractors(data = cluster.data,
                                          

Re: [R] data vector to corresonding percentile ranks

2011-12-23 Thread Sarah Goslee
I'm not sure I understand the question, but does quantile() do what you want?

On Fri, Dec 23, 2011 at 10:28 AM, Steve Jones sjone...@jhmi.edu wrote:
 I have a problem where I need to calculate the corresponding cohort
 percentile ranks for each of several variables.

 Essentially, what I need is a function that will calculate the
 distribution-free percentiles from each variable's data vector, returning a
 corresponding vector of percentiles:

 e.g.:

 percentile.my.data-/function/(my.data)


 I tried to make ecdf() perform this task but was unsuccessful.

 I'd be grateful for any help or advice...




-- 
Sarah Goslee
http://www.functionaldiversity.org

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 code

2011-12-23 Thread Sarah Goslee
On Fri, Dec 23, 2011 at 9:25 AM, 1Rnwb sbpuro...@gmail.com wrote:
 this is how the ouput from the code should be

What code might that be?

Readers of the R-help email list have no idea whatsoever what you're asking.

Please include context, as requested in the posting guide.

Sarah


 structure(list(HTN = 1:10, HTN_FDR = structure(c(4L, 2L, 1L,
 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc,
 T1D_oc, T1d_w), class = factor), Dyslipidemia = structure(c(3L,
 2L, 1L, 2L, 4L, 1L, 1L, 2L, 4L, 2L), .Label = c(Ctrl_noc, T1D_noc,
 T1D_oc, T1D_w), class = factor), CAD = structure(c(3L,
 2L, 1L, 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc,
 T1D_oc), class = factor), CAD_FDR = structure(c(3L, 2L, 1L,
 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc,
 T1D_oc), class = factor), Prior_MI = structure(c(3L, 2L,
 1L, 2L, 3L, 1L, 1L, 2L, 3L, 2L), .Label = c(Ctrl_noc, T1D_noc,
 T1D_oc), class = factor), t1d_ptype = structure(c(3L, 3L,
 2L, 3L, 1L, 1L, 2L, 3L, 3L, 3L), .Label = c(Ctrl, Ctrl_FDR,
 T1D), class = factor)), .Names = c(HTN, HTN_FDR, Dyslipidemia,
 CAD, CAD_FDR, Prior_MI, t1d_ptype), class = data.frame, row.names
 = c(NA,
 -10L))



-- 
Sarah Goslee
http://www.functionaldiversity.org

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


Re: [R] data vector to corresonding percentile ranks

2011-12-23 Thread David Winsemius


On Dec 23, 2011, at 10:28 AM, Steve Jones wrote:


I have a problem where I need to calculate the corresponding cohort
percentile ranks for each of several variables.

Essentially, what I need is a function that will calculate the
distribution-free percentiles from each variable's data vector,  
returning a

corresponding vector of percentiles:

e.g.:

percentile.my.data-/function/(my.data)


I tried to make ecdf() perform this task but was unsuccessful.


Unsuccessful? How? Seems like a reasonable strategy:

set.seed(123)
 x - rnorm(1000)
xCdist - ecdf(x)

Seems to give sensible results.

 x[1]
[1] -0.7104066
 100*xCdist(x[1])
[1] 23.4

 x[2]
[1] 0.2568837
 100*xCdist(x[2])
[1] 60


I'd be grateful for any help or advice...


My advice would be to post what code you were trying so that you can  
get help understand what difficulties you need to overcome.


--

David Winsemius, MD
West Hartford, CT

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


Re: [R] data vector to corresonding percentile ranks

2011-12-23 Thread Sarah Goslee
It's far more useful to send your explanation to the list than it is to send it
just to me. I've taken the liberty of doing so.

But this does sound like a job for ecdf() - what did you do, and what went
wrong?

On Fri, Dec 23, 2011 at 12:23 PM, Steven Jones sjone...@jhmi.edu wrote:
 Actually, what I am trying to do is very simple.  I want to take a large 
 vector of data, about 1.5 million observations and return a corresponding 
 vector of percentiles, each corresponding to the elements of the original 
 data vector in a 1:1 correspondence.

 Original data vector:      x--(x1,x2,x3,...xn)
                                ^  ^  ^
                                |  |  |
 Derived percentile vector: p--(p1,p2,p3,...pn)

 Where pn--is the percentile value of xn determined from the rank in the data 
 vector x.

 Ideally, there is a function which could return a vector p from my original 
 data vector x.

 My understanding, which may be incorrect, is that quantile() function returns 
 the value of x for a specified percentile rank in a vector.

 Steve

 -Original Message-
 From: Sarah Goslee [mailto:sarah.gos...@gmail.com]
 Sent: Friday, December 23, 2011 12:11 PM
 To: Steven Jones
 Cc: r-help@r-project.org
 Subject: Re: [R] data vector to corresonding percentile ranks

 I'm not sure I understand the question, but does quantile() do what you want?

 On Fri, Dec 23, 2011 at 10:28 AM, Steve Jones sjone...@jhmi.edu wrote:
 I have a problem where I need to calculate the corresponding cohort
 percentile ranks for each of several variables.

 Essentially, what I need is a function that will calculate the
 distribution-free percentiles from each variable's data vector, returning a
 corresponding vector of percentiles:

 e.g.:

 percentile.my.data-/function/(my.data)


 I tried to make ecdf() perform this task but was unsuccessful.

 I'd be grateful for any help or advice...




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] black and white in qplot? layout 4 graphs in one screen

2011-12-23 Thread Hadley Wickham
You might find the ggplot mailing list a friendlier place to ask
questions about ggplot2.
Hadley

On Wed, Dec 21, 2011 at 2:16 PM, rachaelohde cox.rach...@gmail.com wrote:
 Hello,

 I am trying to plot means and standard errors conditioned by a factor, using
 qplot.  I am successful at getting the bar graph I want with a error bar,
 however I have tried many things and cannot get the bars to change colors.
 Currently showing as red and blue, but need it to be black and white for
 publication.  Any suggestions please?

 Using a data set June, which is str:
 'data.frame':   21 obs. of  6 variables:
  $ BLCK     : Factor w/ 4 levels ,B,I,W: 3 4 2 3 2 2 2 4 3 4 ...
  $ PLOT     : int  3 3 6 1 2 5 1 1 2 2 ...
  $ TRT      : Factor w/ 5 levels ,crop,ten,..: 2 2 2 3 3 3 4 4 4 5 ...
  $ Date     : Factor w/ 8 levels , 6/16/11,..: 2 2 2 2 2 2 2 2 2 2 ...
  $ habitat  : Factor w/ 3 levels ,C,P: 2 2 2 2 2 2 2 2 2 2 ...
  $ Abundance: num  0.333 0 1.333 0 1.667 ...

 Current code is:

 Ab.avg-ddply(June, c(TRT, habitat), function(df)
  return(c(Ab.avg=mean(df$Abundance), Ab.sd=sd(df$Abundance

 avg.plot-qplot(TRT, Ab.avg, fill=factor(habitat),
 data=Ab.avg, geom=bar, position=dodge)

 dodge - position_dodge(width=0.9)
 avg.plot++geom_linerange(aes(ymax=Ab.avg+Ab.sd, ymin=Ab.avg-Ab.sd),
 position=dodge)+theme_bw()

 http://r.789695.n4.nabble.com/file/n4223035/june_bar_graph.png

 Also, would like to plot 4 of these bar graphs (for four dates) on the same
 screen, I cannot get the par() or layout() function to work with qplot.  Is
 there another way?

 Thank you!

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/black-and-white-in-qplot-layout-4-graphs-in-one-screen-tp4223035p4223035.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.



-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] Stacked area plot for time series

2011-12-23 Thread Hadley Wickham
You are more likely to get a helpful response if you provide a
reproducible example - without that I can only guess that you need to
use approx so you get y values at same x values.

Hadley

On Wed, Dec 21, 2011 at 8:13 AM, UncleFish bpn...@ucsd.edu wrote:
 I wish to make a stacked area chart of a time series with three variables.
 The observations are very irregular, covering 500 years, with a few
 historical observations in the range 1500-1850, and then more regular
 observations since 1880 or so.  If I plot this in Excell it truncates the
 time series, but it can make a nice stacked area plot. My question is how
 can I have a correct x axis with a stacked plot in R? I thave tried the data
 in a zoo as well as a ts, and have tried ggplot2 and plotrix, but cannot get
 beyond a lattice plot.

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Stacked-area-plot-for-time-series-tp4221849p4221849.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.



-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] vif function using lm object

2011-12-23 Thread Liviu Andronic
On Wed, Dec 21, 2011 at 9:28 AM, arunkumar akpbond...@gmail.com wrote:
 Hi,

   can anyone please explain why the vif should have more than 2 terms.

 *vif.lm(lmobj) : model contains fewer than 2 terms*

 why it is throwng error if it is one variable.

The _VIF_ is a measure of _multicollinearity_, which occurs between
two or more predictors in a regression model. In other words, you
cannot have collinearity in a single predictor.

Also, to compute the VIF you need to regress one regressor on all the
remaining regressors. If you have a single regressor, then you cannot
estimate a VIF.

See the related articles on Wikipedia.

Regards
Liviu






 -
 Thanks in Advance
        Arun
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/vif-function-using-lm-object-tp4220904p4220904.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.



-- 
Do you know how to read?
http://www.alienetworks.com/srtest.cfm
http://goodies.xfce.org/projects/applications/xfce4-dict#speed-reader
Do you know how to write?
http://garbl.home.comcast.net/~garbl/stylemanual/e.htm#e-mail

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


Re: [R] data vector to corresonding percentile ranks

2011-12-23 Thread R. Michael Weylandt
This seems like what you are interested in:

x - rnorm(50)
ecdf(x)(x)

Michael

On Fri, Dec 23, 2011 at 11:26 AM, Sarah Goslee sarah.gos...@gmail.com wrote:
 It's far more useful to send your explanation to the list than it is to send 
 it
 just to me. I've taken the liberty of doing so.

 But this does sound like a job for ecdf() - what did you do, and what went
 wrong?

 On Fri, Dec 23, 2011 at 12:23 PM, Steven Jones sjone...@jhmi.edu wrote:
 Actually, what I am trying to do is very simple.  I want to take a large 
 vector of data, about 1.5 million observations and return a corresponding 
 vector of percentiles, each corresponding to the elements of the original 
 data vector in a 1:1 correspondence.

 Original data vector:      x--(x1,x2,x3,...xn)
                                ^  ^  ^
                                |  |  |
 Derived percentile vector: p--(p1,p2,p3,...pn)

 Where pn--is the percentile value of xn determined from the rank in the 
 data vector x.

 Ideally, there is a function which could return a vector p from my original 
 data vector x.

 My understanding, which may be incorrect, is that quantile() function 
 returns the value of x for a specified percentile rank in a vector.

 Steve

 -Original Message-
 From: Sarah Goslee [mailto:sarah.gos...@gmail.com]
 Sent: Friday, December 23, 2011 12:11 PM
 To: Steven Jones
 Cc: r-help@r-project.org
 Subject: Re: [R] data vector to corresonding percentile ranks

 I'm not sure I understand the question, but does quantile() do what you want?

 On Fri, Dec 23, 2011 at 10:28 AM, Steve Jones sjone...@jhmi.edu wrote:
 I have a problem where I need to calculate the corresponding cohort
 percentile ranks for each of several variables.

 Essentially, what I need is a function that will calculate the
 distribution-free percentiles from each variable's data vector, returning a
 corresponding vector of percentiles:

 e.g.:

 percentile.my.data-/function/(my.data)


 I tried to make ecdf() perform this task but was unsuccessful.

 I'd be grateful for any help or advice...




 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Problem with RCOM package

2011-12-23 Thread Erich Neuwirth
All questions related to rcom, statconnDCOM, and RExcel
should be posted on the mailing list of the statconnDCOM project.
You can subscribe at
http://rcom.univie.ac.at


On 12/23/2011 11:02 AM, KUMAR wrote:
 Hi,
 
 I am using R version 2.14.0 on windows 7. Trying to use RCOM
 package and followed example provided at
 (http://cran.r-project.org/web/packages/rcom/rcom.pdf) .
 
  
 
 I am able to load the com object as seen in screenshot
 2.png.  however I am unable to invoke/get/set property to the object.
 
  
 
  
 
 Please guide me on the same,
 
  
 
 Regards;
 
 Kumar Anand
 
 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] cast in reshape and reshape2

2011-12-23 Thread Kaiyin Zhong
 library(reshape2)
 x = melt(airquality, id=c('month', 'day'))

With reshape I can cast with multiple functions:

 library(reshape)
 cast(x, month+variable~., c(mean,sd))
   month variable   mean sd
1  5ozone  23.615385  22.224449
2  5  solar.r 181.296296 115.075499
3  5 wind  11.622581   3.531450
4  5 temp  65.548387   6.854870
5  6ozone  29.44  18.207904
6  6  solar.r 190.17  92.882975
7  6 wind  10.27   3.769234
8  6 temp  79.10   6.598589
9  7ozone  59.115385  31.635837
10 7  solar.r 216.483871  80.568344
11 7 wind   8.941935   3.035981
12 7 temp  83.903226   4.315513
13 8ozone  59.961538  39.681210
14 8  solar.r 171.857143  76.834943
15 8 wind   8.793548   3.225930
16 8 temp  83.967742   6.585256
17 9ozone  31.448276  24.141822
18 9  solar.r 167.43  79.118280
19 9 wind  10.18   3.461254
20 9 temp  76.90   8.355671

Is there a way to do the same job with reshape2?
-- 
Kaiyin Zhong
--
Neuroscience Research Institute, Peking University, Beijing, P.R.China 100038

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Applying a function

2011-12-23 Thread Steinar Veka
First, it may be a good idea to use 1 for “yes” and 0 for “no” in the
motherrecovery column.

Then if you name your table e.g tab1,  you may try something like this..

sum(tab1[,2]==1*tab1[,4])/sum(tab1[,2]==1)

2011/12/23 Joanie Van De Walle joanie.van-de-wall...@ulaval.ca


   Hi,

   I need help writing a function

   I capture seal pups mutliple times during the lactation season in order
 to
   monitor their growth rate. When I release them, the recovery (mother-pup)
   time is not the same for all individuals. I want to know if individuals
 that
   recover their mother the fastest are the ones with the highest growth
 rates.

   So, I noted at every release if the pup reovered his mother before we
 leave
   (yes or no). My dataframe looks like this

   Capture nb individual individual capture motherrecovery growth rate
   1  1  1  n  0.5
   2  1  2  y  0.5
   3  1  3  y  0.5
   4  1  4  y  0.5
   5  1  5  n  0.5
   6  2  1  y  0.3
   7  2  2  y  0.3
   8  3  1  y  0.4
   9  3  2  n  0.4
   10 3  3  y  0.4
   ...

   I want to calculate a rate of mother recovery by individual, i.e. nb of
   recoveries (y)/nb of captures. So for indivial 1 it would be 3/5 = 0.6. I
   want  to write a function that does this for all the individuals in my
   dataframe,  i.e. around 400 individuals (this is why I want to write a
   function, it would be too long by hand)

   Thank you,



   Joanie Van de Walle

   Étudiante à la Maîtrise en Biologie
   Département de Biologie
   1045, avenue de la Médecine Pavillon Vachon,
   bureau 2044 Université Laval
   Québec, Canada, G1V 0A6
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Custom XML Readers

2011-12-23 Thread pl.r...@gmail.com
I need to construct a custom XML reader, the files I'm working with are in
funky XML format:

str name=authorPaul H/str
  str name=countryUSA/str
  date name=created_date2010-02-16/date
 
I want to read the file so it looks like:

author = Paul H
country = USA
created_date=2010-02-16

Does any one know how to go about this problem, or know of good references i
could access?

Thanks,
Andy


--
View this message in context: 
http://r.789695.n4.nabble.com/Custom-XML-Readers-tp4229614p4229614.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] Stacked area plot for time series

2011-12-23 Thread Ben Neal
Yes, thanks for the direction. I am new to the list and new to R, and I can see 
that was not a very helpful description. 

I have since resolved my issue. My main issue was that stackpoly does not (I 
think) interact with XOO objects or ts object, and my secondary issue was that 
my time series was highly irregular (with inferred indigenous historical 
consumption rates of fish, separated by large gaps). I resolved this by taking 
my data, making it a zoo to order it, transforming to a ts to interpolate 
missing values, and transforming to a matrix to do the stackpoly plot.

I then resolved the issue of too many tick marks on the x axis by making a new 
axis. I did seem to have difficulty with some plot functions apparently do not 
to function in plotrix (specifically xaxt='n'). 

I am aware that stacked area charts might not an acceptable graphical 
representation of the data in some opinions, and I tend to agree. However, this 
output had to match the format of other graphs in the work, and was thus 
desired in this specific format. 

Thanks for providing this stackpoly tool - very handy! Code follows, and a 
final graphic is attached: 

library(zoo)
library(plotrix)
setwd(~/Documents/Loren/RawData)
CaribData = read.csv(FLA_Loren.csv, header = TRUE)
attach(CaribData)

## Create ZOO to order by time
z1=zoo(Commercial,Year)
z2=zoo(Recreational,Year)
z3=zoo(Subsistence,Year)
z4-merge(z1,z2,z3)

##Convert to ts to interpolate missing NA values
ts1-as.ts(z1)
ts1Ex-na.approx(ts1)
ts2-as.ts(z2)
ts2Ex-na.approx(ts2)
ts3-as.ts(z3)
ts3Ex-na.approx(ts3)

## Merge into one object
ts4Ex-cbind(ts1Ex, ts2Ex, ts3Ex, deparse.level = 1)
##Test plots to examine series
plot(ts1Ex)
plot(ts2Ex)
plot(ts3Ex)
plot(ts4Ex)

## Convert for stacked area plotting and plot
CaribMatrix-as.matrix(ts4Ex)
stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons,
  main=Historical Florida reef fisheries catch by sector, 
axis4=FALSE)

stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons, 
  col.main=red, font.main=4,
  axis4=FALSE,  xat=seq(50,760,by=100), 
  xaxlab=c(1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000))








-Original Message-
From: h.wick...@gmail.com on behalf of Hadley Wickham
Sent: Fri 12/23/2011 9:41 AM
To: UncleFish
Cc: r-help@r-project.org
Subject: Re: [R] Stacked area plot for time series
 
You are more likely to get a helpful response if you provide a
reproducible example - without that I can only guess that you need to
use approx so you get y values at same x values.

Hadley

On Wed, Dec 21, 2011 at 8:13 AM, UncleFish bpn...@ucsd.edu wrote:
 I wish to make a stacked area chart of a time series with three variables.
 The observations are very irregular, covering 500 years, with a few
 historical observations in the range 1500-1850, and then more regular
 observations since 1880 or so.  If I plot this in Excell it truncates the
 time series, but it can make a nice stacked area plot. My question is how
 can I have a correct x axis with a stacked plot in R? I thave tried the data
 in a zoo as well as a ts, and have tried ggplot2 and plotrix, but cannot get
 beyond a lattice plot.

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Stacked-area-plot-for-time-series-tp4221849p4221849.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.



-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] Cuzick's test for trend

2011-12-23 Thread andrewejaffe
I just wrote this up for a project, perhaps 5 years is better than never...

code
cuzick = function(x,z,test.type=c(two.sided, upper, lower)) {

N = length(z)
n = unique(z)

ranks=rank(x)

T = sum(ranks*z)

p = (table(z)/N)
E_Z = sum(unique(z)*p)
E_T = 0.5*N*(N+1)*E_Z

Var_Z = sum(unique(z)^2*p) - E_Z^2
Var_T = N^2*(N+1)/12*Var_Z

Zscore = (T-E_T)/sqrt(Var_T)

if(test.type == two.sided) {
pval = 2*pnorm(-abs(Zscore))
} else if(test.type == upper) {
pval = pnorm(Zscore,lower.tail=F)
} else  pval = pnorm(Zscore,lower.tail=T)

out = data.frame(cbind(Zscore,pval,test.type))
colnames(out) = c(Z,p,testType)
return(out)
}
/code

With data from the Cuzick 1985 paper:
code
z = c(rep(1, 8), rep(2,10), rep(3,9), rep(4,9),rep(5,9))
x = c(0, 0, 1, 1, 2, 2, 4, 9,
0, 0, 5, 7, 8, 11, 13, 23, 25, 97,
2, 3, 6, 9, 10, 11, 11, 12, 21,
0, 3, 5, 6, 10, 19, 56, 100, 132,
2, 4, 6, 6, 6, 7, 18, 39, 60)
 cuzick(x,z,two.sided)
 Z  p  testType
1 2.11 0.0348 two.sided
/code



arin basu-2 wrote
 
 Hi All:
 
 I was looking for, but could not locate in the packages, or in the R
 archive searches if there exists an R implementation of Cuzick's test of
 trend. The test is described as follows:
 
 An extension of the Wilcoxon rank-sum test is developed to handle the
 situation in which a variable is measured for individuals in three or more
 (ordered) groups and a non-parametric test for trend across these groups
 is desired.
 
 Reference:
 
 Cuzick J. A Wilcoxon-type test for trend. Stat Med. 1985
 Jan-Mar;4(1):87-90
 
 Would greatly appreciate your insights. The R version I use is R-2.3.1 
 
 Best,
 Arin Basu
 
 __
 R-help@.ethz mailing list
 https://stat.ethz.ch/mailman/listinfo/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://r.789695.n4.nabble.com/Cuzick-s-test-for-trend-tp807101p4229374.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] Stacked area plot for time series

2011-12-23 Thread Ben Neal
Meant to say ZOO objects . . . 


-Original Message-
From: Ben Neal [mailto:bn...@spg.ucsd.edu]
Sent: Fri 12/23/2011 10:15 AM
To: Hadley Wickham; UncleFish
Cc: r-help@r-project.org
Subject: RE: [R] Stacked area plot for time series
 
Yes, thanks for the direction. I am new to the list and new to R, and I can see 
that was not a very helpful description. 

I have since resolved my issue. My main issue was that stackpoly does not (I 
think) interact with XOO objects or ts object, and my secondary issue was that 
my time series was highly irregular (with inferred indigenous historical 
consumption rates of fish, separated by large gaps). I resolved this by taking 
my data, making it a zoo to order it, transforming to a ts to interpolate 
missing values, and transforming to a matrix to do the stackpoly plot.

I then resolved the issue of too many tick marks on the x axis by making a new 
axis. I did seem to have difficulty with some plot functions apparently do not 
to function in plotrix (specifically xaxt='n'). 

I am aware that stacked area charts might not an acceptable graphical 
representation of the data in some opinions, and I tend to agree. However, this 
output had to match the format of other graphs in the work, and was thus 
desired in this specific format. 

Thanks for providing this stackpoly tool - very handy! Code follows, and a 
final graphic is attached: 

library(zoo)
library(plotrix)
setwd(~/Documents/Loren/RawData)
CaribData = read.csv(FLA_Loren.csv, header = TRUE)
attach(CaribData)

## Create ZOO to order by time
z1=zoo(Commercial,Year)
z2=zoo(Recreational,Year)
z3=zoo(Subsistence,Year)
z4-merge(z1,z2,z3)

##Convert to ts to interpolate missing NA values
ts1-as.ts(z1)
ts1Ex-na.approx(ts1)
ts2-as.ts(z2)
ts2Ex-na.approx(ts2)
ts3-as.ts(z3)
ts3Ex-na.approx(ts3)

## Merge into one object
ts4Ex-cbind(ts1Ex, ts2Ex, ts3Ex, deparse.level = 1)
##Test plots to examine series
plot(ts1Ex)
plot(ts2Ex)
plot(ts3Ex)
plot(ts4Ex)

## Convert for stacked area plotting and plot
CaribMatrix-as.matrix(ts4Ex)
stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons,
  main=Historical Florida reef fisheries catch by sector, 
axis4=FALSE)

stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons, 
  col.main=red, font.main=4,
  axis4=FALSE,  xat=seq(50,760,by=100), 
  xaxlab=c(1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000))








-Original Message-
From: h.wick...@gmail.com on behalf of Hadley Wickham
Sent: Fri 12/23/2011 9:41 AM
To: UncleFish
Cc: r-help@r-project.org
Subject: Re: [R] Stacked area plot for time series
 
You are more likely to get a helpful response if you provide a
reproducible example - without that I can only guess that you need to
use approx so you get y values at same x values.

Hadley

On Wed, Dec 21, 2011 at 8:13 AM, UncleFish bpn...@ucsd.edu wrote:
 I wish to make a stacked area chart of a time series with three variables.
 The observations are very irregular, covering 500 years, with a few
 historical observations in the range 1500-1850, and then more regular
 observations since 1880 or so.  If I plot this in Excell it truncates the
 time series, but it can make a nice stacked area plot. My question is how
 can I have a correct x axis with a stacked plot in R? I thave tried the data
 in a zoo as well as a ts, and have tried ggplot2 and plotrix, but cannot get
 beyond a lattice plot.

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Stacked-area-plot-for-time-series-tp4221849p4221849.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.



-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] Axis manipulation in Stackpoly (Plotrix)

2011-12-23 Thread Ben Neal
Thank you. That helps me understand the issue better. I came up with a similar 
solution . . . but yours is more elegant (I just wrote out the labels . . ).

My solution follows: 

stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons, 
  col.main=red, font.main=4,
  axis4=FALSE,  xat=seq(50,760,by=100), 
  xaxlab=c(1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000))


-Original Message-
From: Jim Lemon [mailto:j...@bitwrit.com.au]
Sent: Fri 12/23/2011 2:53 AM
To: Ben Neal
Cc: r-help@r-project.org
Subject: Re: [R] Axis manipulation in Stackpoly (Plotrix)
 
On 12/23/2011 07:52 AM, Ben Neal wrote:

 Have made a stacked area plot, but now want to manipulate the x axis: make in 
 even increments (50) in order to suppress the tick marks (forming solid bar 
 under plot). However, the plot functions do not seem to work, and I cannot 
 find documentation for use of xat in stackpoly.

 This is a time series of data covering 579 years, and has been successfully 
 converted from zoo to matrix. Data plotting fine, jsut want to change axis -

 CaribMatrix-as.matrix(ts4Ex)
 stackpoly(CaribMatrix, stack=TRUE, xlab=Year,ylab=Catch in tons, 
 xaxlab=(1250:2008),
main=Historical Florida reef fisheries catch by sector, 
 axis4=FALSE)

 Any assistance appreciated! Thanks very much, Benjamin P Neal, Scripps Inst. 
 of Oceanogrpahy

Hi Ben,
In the stackpoly function, you can specify the positions of the axis 
ticks (xat=...) and the labels for those positions (xaxlab=...). You 
don't seem to have supplied the x values, so the values in CaribMatrix 
will be placed display at 1:dim(CaribMatrix)[1]. I'll have to make up 
CaribMatrix...

CaribMatrix-matrix(order(sample(200:500,579,TRUE),nrow=579,ncol=3)
stackpoly(CaribMatrix,stack=TRUE,ylim=c(0,1500),xat=seq(1,600,by=50),
  xaxlab=seq(1450,2000,by=50),staxx=TRUE)

This seems to give a reasonable plot for me.

Jim

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


Re: [R] black and white in qplot? layout 4 graphs in one screen

2011-12-23 Thread baptiste auguie
Hi,

On 22 December 2011 09:16, rachaelohde cox.rach...@gmail.com wrote:
 Hello,

 I am trying to plot means and standard errors conditioned by a factor, using
 qplot.  I am successful at getting the bar graph I want with a error bar,
 however I have tried many things and cannot get the bars to change colors.
 Currently showing as red and blue, but need it to be black and white for
 publication.  Any suggestions please?

Have a look at scale_fill_grey() and grid.arrange()

d= data.frame(x=sample(letters[1:4], 10, replace=TRUE), f=gl(2, 10))
( p = qplot(x, data=d, position=dodge, fill=f) + scale_fill_grey() )
library(gridExtra)
grid.arrange(p, p, p , p, ncol=2)

HTH,

baptiste


 Using a data set June, which is str:
 'data.frame':   21 obs. of  6 variables:
  $ BLCK     : Factor w/ 4 levels ,B,I,W: 3 4 2 3 2 2 2 4 3 4 ...
  $ PLOT     : int  3 3 6 1 2 5 1 1 2 2 ...
  $ TRT      : Factor w/ 5 levels ,crop,ten,..: 2 2 2 3 3 3 4 4 4 5 ...
  $ Date     : Factor w/ 8 levels , 6/16/11,..: 2 2 2 2 2 2 2 2 2 2 ...
  $ habitat  : Factor w/ 3 levels ,C,P: 2 2 2 2 2 2 2 2 2 2 ...
  $ Abundance: num  0.333 0 1.333 0 1.667 ...

 Current code is:

 Ab.avg-ddply(June, c(TRT, habitat), function(df)
  return(c(Ab.avg=mean(df$Abundance), Ab.sd=sd(df$Abundance

 avg.plot-qplot(TRT, Ab.avg, fill=factor(habitat),
 data=Ab.avg, geom=bar, position=dodge)

 dodge - position_dodge(width=0.9)
 avg.plot++geom_linerange(aes(ymax=Ab.avg+Ab.sd, ymin=Ab.avg-Ab.sd),
 position=dodge)+theme_bw()

 http://r.789695.n4.nabble.com/file/n4223035/june_bar_graph.png

 Also, would like to plot 4 of these bar graphs (for four dates) on the same
 screen, I cannot get the par() or layout() function to work with qplot.  Is
 there another way?

 Thank you!

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/black-and-white-in-qplot-layout-4-graphs-in-one-screen-tp4223035p4223035.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] cast in reshape and reshape2

2011-12-23 Thread David Winsemius


On Dec 23, 2011, at 2:58 PM, Kaiyin Zhong wrote:


library(reshape2)
x = melt(airquality, id=c('month', 'day'))


With reshape I can cast with multiple functions:


library(reshape)
cast(x, month+variable~., c(mean,sd))

  month variable   mean sd
1  5ozone  23.615385  22.224449
2  5  solar.r 181.296296 115.075499
3  5 wind  11.622581   3.531450
4  5 temp  65.548387   6.854870
5  6ozone  29.44  18.207904
6  6  solar.r 190.17  92.882975
7  6 wind  10.27   3.769234
8  6 temp  79.10   6.598589
9  7ozone  59.115385  31.635837
10 7  solar.r 216.483871  80.568344
11 7 wind   8.941935   3.035981
12 7 temp  83.903226   4.315513
13 8ozone  59.961538  39.681210
14 8  solar.r 171.857143  76.834943
15 8 wind   8.793548   3.225930
16 8 temp  83.967742   6.585256
17 9ozone  31.448276  24.141822
18 9  solar.r 167.43  79.118280
19 9 wind  10.18   3.461254
20 9 temp  76.90   8.355671

Is there a way to do the same job with reshape2?


Have you looked at the .summarise argument to dcast? That seems to  
deliver the same sort of results one gets with base::aggregate.


--
David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] cast in reshape and reshape2

2011-12-23 Thread David Winsemius


On Dec 23, 2011, at 5:58 PM, David Winsemius wrote:



On Dec 23, 2011, at 2:58 PM, Kaiyin Zhong wrote:


library(reshape2)
x = melt(airquality, id=c('month', 'day'))


With reshape I can cast with multiple functions:


library(reshape)
cast(x, month+variable~., c(mean,sd))

 month variable   mean sd
1  5ozone  23.615385  22.224449
2  5  solar.r 181.296296 115.075499
3  5 wind  11.622581   3.531450
4  5 temp  65.548387   6.854870
5  6ozone  29.44  18.207904
6  6  solar.r 190.17  92.882975
7  6 wind  10.27   3.769234
8  6 temp  79.10   6.598589
9  7ozone  59.115385  31.635837
10 7  solar.r 216.483871  80.568344
11 7 wind   8.941935   3.035981
12 7 temp  83.903226   4.315513
13 8ozone  59.961538  39.681210
14 8  solar.r 171.857143  76.834943
15 8 wind   8.793548   3.225930
16 8 temp  83.967742   6.585256
17 9ozone  31.448276  24.141822
18 9  solar.r 167.43  79.118280
19 9 wind  10.18   3.461254
20 9 temp  76.90   8.355671

Is there a way to do the same job with reshape2?


Have you looked at the .summarise argument to dcast? That seems to  
deliver the same sort of results one gets with base::aggregate.


Actually I see after looking at examples on the plyr-reshape- 
googlegroups group that it is not '.summarise' but rather 'summarise'.  
Unfortunately there are no links in the help pages that seem to  
describe its proper use ... a not uncommon failing for that package in  
my experience.




--

David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Latent class multinomial (or conditional) logit using R?

2011-12-23 Thread adanmartinez
Hi everyone?

Does anybody know how can I estimate a
Latent class multinomial (or conditional) logit using R?

I have tried flexmix, poLCA, and
they do not seem to support this model.
thanks in advance

adan

--
View this message in context: 
http://r.789695.n4.nabble.com/Latent-class-multinomial-or-conditional-logit-using-R-tp4230083p4230083.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] Latent class multinomial (or conditional) logit using R?

2011-12-23 Thread Christian Hennig
You may have a look at lcmixed/flexmixedruns in package fpc. This provides 
flexmix methods that can be used to fit latent class models with locally 
independent multinomials (no logits, though, and I'm not sure that for 
categorical data only this is much different from what poLCA offers).


Best regards,
Christian

On Fri, 23 Dec 2011, adanmartinez wrote:


Hi everyone?

Does anybody know how can I estimate a
Latent class multinomial (or conditional) logit using R?

I have tried flexmix, poLCA, and
they do not seem to support this model.
thanks in advance

adan

--
View this message in context: 
http://r.789695.n4.nabble.com/Latent-class-multinomial-or-conditional-logit-using-R-tp4230083p4230083.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.



*** --- ***
Christian Hennig
University College London, Department of Statistical Science
Gower St., London WC1E 6BT, phone +44 207 679 1698
chr...@stats.ucl.ac.uk, www.homepages.ucl.ac.uk/~ucakche

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 code

2011-12-23 Thread Rui Barradas
Hello,

There are so many people posting answers that I'm curious and decided to try
one.

I don't know if this is it but it doesn't give an error and it reformats
your data
according to the rules in your original code.

#
nr - dim(c1)[1]
nc - dim(c1)[2]
c2 - NULL
c2_row - rep(, nc-1)
for(i in 1:nr){
ptype- as.character(c1$t1d_ptype[i])
stype- substr(ptype,1,4)
num_comp - sum(c1[i,] == Y)
for(j in 1:(nc-1)){
c2_row[j] - 
if(num_comp  0){
c1_tmp - as.integer(c1[i ,j])
if(ptype == T1D){
if(c1_tmp == 1) c2_row[j] - T1D_oc
if(c1_tmp == 2) c2_row[j] - T1D_w
}
if(stype == Ctrl){
if(c1_tmp == 1) c2_row[j] - Ctrl_oc
if(c1_tmp == 2) c2_row[j] - Ctrl_w
}
}
else{
if(ptype == T1D)  c2_row[j] - T1D_noc
if(stype == Ctrl) c2_row[j] - Ctrl_noc
}
}
c2 - rbind(c2, c(c2_row, ptype))
}
c2 - data.frame(c2)
colnames(c2) - colnames(c1)
c2

I bet there's a way to work on entire objects. This is C-like.

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/Help-with-code-tp4218989p4229987.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] Applying a function

2011-12-23 Thread Rui Barradas

Joanie Van De Walle wrote
 
 Hi,
 
I need help writing a function
 
I capture seal pups mutliple times during the lactation season in order
 to
monitor their growth rate. When I release them, the recovery
 (mother-pup)
time is not the same for all individuals. I want to know if individuals
 that
recover their mother the fastest are the ones with the highest growth
 rates.
 
So, I noted at every release if the pup reovered his mother before we
 leave
(yes or no). My dataframe looks like this
 
Capture nb individual individual capture motherrecovery growth rate
1  1  1  n  0.5
2  1  2  y  0.5
3  1  3  y  0.5
4  1  4  y  0.5
5  1  5  n  0.5
6  2  1  y  0.3
7  2  2  y  0.3
8  3  1  y  0.4
9  3  2  n  0.4
10 3  3  y  0.4
...
 
I want to calculate a rate of mother recovery by individual, i.e. nb of
recoveries (y)/nb of captures. So for indivial 1 it would be 3/5 = 0.6.
 I
want  to write a function that does this for all the individuals in my
dataframe,  i.e. around 400 individuals (this is why I want to write a
function, it would be too long by hand)
 
Thank you,
 
 
 __
 R-help@ mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

Hello,

Here is a one line function

recovery.rate - function(x) unlist(lapply(split(tab1, tab1[,2]),
function(x) mean(x[[4]]==y)))

It works with that table. I hope it helps

Merry Christmas

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/Applying-a-function-tp4229277p4230100.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] if statement problem

2011-12-23 Thread reena
Hello,

I want to do fisher test for the rows in data file which has value less than
5 otherwise chi square test .The p values from both test should be stored in
one resulted file. but there is some problem with bold if statement. I don't
know how 
implement this line properly.


x = cbind(obs1,obs2,exp1,exp2)
a = matrix(c(0,0,0,0), ncol=2, byrow =TRUE)#matrix with initialized
values

for (i in 1: length(x[,1]))
{
  *if((x[i,1] || x[i,2] || x[i,3] || x[i,4])  5)*
 {
 a[1,1]  - x[i,1];
 a[1,2] - x[i,2];
 a[2,1] - x[i,3];
 a[2,2] - x[i,4];
 result - fisher.test(a)
 write.table(result[[p.value]],file=results.txt,
sep=\n, append=TRUE, col.names=FALSE, row.names=FALSE);
   }

 else
{
a[1,1] - x[i,1];
a[1,2] - x[i,2];
a[2,1] - x[i,3];
a[2,2] - x[i,4];
result - chisq.test(a)
write.table(result[[p.value]],file=results.txt,
sep=\n, append=TRUE, col.names=FALSE, row.names=FALSE);}
}

Regards
R

--
View this message in context: 
http://r.789695.n4.nabble.com/if-statement-problem-tp4230026p4230026.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] Trouble converting hourly data into daily data

2011-12-23 Thread Sean Baumgarten
Hi Jean,

Thanks for the help. I couldn't quite get the results I needed with the
merge command, but I ended up using the following work-around:

Weather - read.csv(Weather.csv)
Weather$diff.time - abs(.5 - Weather$TimeNumeric)
agg - aggregate(diff.time ~ Date, data = Weather, FUN = which.min)
n.obs - cumsum(rle(as.double(Weather$Date))$lengths)
n.obs - c(0, n.obs[1:(length(n.obs) - 1)])
noon.ind - agg$diff.time + n.obs
subset - Weather[noon.ind,]

Cheers,
Sean

On Mon, Dec 19, 2011 at 6:03 AM, Jean V Adams jvad...@usgs.gov wrote:


 Sean Baumgarten wrote on 12/14/2011 06:38:08 PM:

  Hello,
 
  I have a data frame with hourly or sub-hourly weather records that span
  several years, and from that data frame I'm trying to select only the
  records taken closest to noon for each day. Here's what I've done so far:
 
  #Add a column to the data frame showing the difference between noon and
 the
  observation time (I converted time to a 0-1 scale so 0.5 represents
 noon):
  data$Diff_from_noon - abs(0.5-data$Time)
 
  #Find the minimum value of Diff_from_noon for each Date:
  aggregated - aggregate(Diff_from_noon ~ Date, data, FUN=min)
 
 
  The problem is that the aggregated data frame only has two columns:
 Date
  and Diff_from_noon. I can't figure out how to get the columns with the
  actual weather variables to carry over from the original data frame.
 
  Any suggestions you have would be much appreciated.
 
  Thanks,
  Sean


 You don't provide any example data, so I will use data from R datasets,
 airquality.  After using the aggregate() function to find the minimum Day
 for each Month, merge the resulting data frame with the original data frame
 to see all the columns corresponding to the selected minimums.

  aggregated - aggregate(Day ~ Month, airquality, FUN=min)
  aggregated
   Month Day
 1 5   1
 2 6   1
 3 7   1
 4 8   1
 5 9   1
  merge(aggregated, airquality)
   Month Day Ozone Solar.R Wind Temp
 1 5   141 190  7.4   67
 2 6   1NA 286  8.6   78
 3 7   1   135 269  4.1   84
 4 8   139  83  6.9   81
 5 9   196 167  6.9   91

 For your data, the code would look like this:
 aggregated - aggregate(Diff_from_noon ~ Date, data, FUN=min)
 merge(aggregated, data)

 I recommend that you use a name other than data for your data frame,
 since data() is a built in R function.

 Jean

[[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] Applying a function

2011-12-23 Thread Rui Barradas
Sorry, in the function body, NO 'tab1', use 'x' only:

recovery.rate - function(x) unlist(lapply(split(x, x[,2]), function(x)
mean(x[[4]]==y)))

The error is because 'tab1' existed in the environment and the function
would find it.

This time, tested after removing 'tab1'.

Rui Barradas



--
View this message in context: 
http://r.789695.n4.nabble.com/Applying-a-function-tp4229277p4230106.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] if statement problem

2011-12-23 Thread Rui Barradas

reena wrote
 
 Hello,
 
 I want to do fisher test for the rows in data file which has value less
 than 5 otherwise chi square test .The p values from both test should be
 stored in one resulted file. but there is some problem with bold if
 statement. I don't know how 
 implement this line properly.
 
 
 x = cbind(obs1,obs2,exp1,exp2)
 a = matrix(c(0,0,0,0), ncol=2, byrow =TRUE)#matrix with
 initialized values
 
 for (i in 1: length(x[,1]))
 {
   *if((x[i,1] || x[i,2] || x[i,3] || x[i,4])  5)*

 

Hello,
Try

*if(any(x[i,] 5))*

Merry Christmas
Rui Barradas



--
View this message in context: 
http://r.789695.n4.nabble.com/if-statement-problem-tp4230026p4230135.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] Latent class multinomial (or conditional) logit using R?

2011-12-23 Thread adanmartinez
Hi Christian,
Thanks a lot for the suggestion.
As you expected, lcmixed does not provide more functionality than poLCA for
the case of
only categorical data...
Do you know whether there is an explanation (some technical difficulty I 
may be missing, I am Economist so it is very likely) of why I can fit a
Poisson, a logit latent class
model but not a conditional logit?

In fact, I started using R when I realize that my (very rudimentary) EM
script written in Matlab
had tons of problems to converge.
Now I start wondering if there is some specific situation I should be aware
of.

Thanks a lot again

--
View this message in context: 
http://r.789695.n4.nabble.com/Latent-class-multinomial-or-conditional-logit-using-R-tp4230083p4230194.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] Custom XML Readers

2011-12-23 Thread Ben Tupper
Hi Andy,

On Dec 23, 2011, at 2:51 PM, pl.r...@gmail.com wrote:

 I need to construct a custom XML reader, the files I'm working with are in
 funky XML format:
 
 str name=authorPaul H/str
  str name=countryUSA/str
  date name=created_date2010-02-16/date
 
 I want to read the file so it looks like:
 
 author = Paul H
 country = USA
 created_date=2010-02-16
 
 Does any one know how to go about this problem, or know of good references i
 could access?
 


Have you tried Duncan Temple Lang's XML package for R?  It works very well for 
parsing and building XML formatted data.

http://www.omegahat.org/RSXML/

Cheers,
Ben



 Thanks,
 Andy
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Custom-XML-Readers-tp4229614p4229614.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.

Ben Tupper
Bigelow Laboratory for Ocean Sciences
180 McKown Point Rd. P.O. Box 475
West Boothbay Harbor, Maine   04575-0475 
http://www.bigelow.org

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] cast in reshape and reshape2

2011-12-23 Thread Hadley Wickham
On Fri, Dec 23, 2011 at 1:58 PM, Kaiyin Zhong kindlych...@gmail.com wrote:
 library(reshape2)
 x = melt(airquality, id=c('month', 'day'))

 With reshape I can cast with multiple functions:

 library(reshape)
 cast(x, month+variable~., c(mean,sd))
   month variable       mean         sd
 1      5    ozone  23.615385  22.224449
 2      5  solar.r 181.296296 115.075499
 3      5     wind  11.622581   3.531450
 4      5     temp  65.548387   6.854870
 5      6    ozone  29.44  18.207904
 6      6  solar.r 190.17  92.882975
 7      6     wind  10.27   3.769234
 8      6     temp  79.10   6.598589
 9      7    ozone  59.115385  31.635837
 10     7  solar.r 216.483871  80.568344
 11     7     wind   8.941935   3.035981
 12     7     temp  83.903226   4.315513
 13     8    ozone  59.961538  39.681210
 14     8  solar.r 171.857143  76.834943
 15     8     wind   8.793548   3.225930
 16     8     temp  83.967742   6.585256
 17     9    ozone  31.448276  24.141822
 18     9  solar.r 167.43  79.118280
 19     9     wind  10.18   3.461254
 20     9     temp  76.90   8.355671

 Is there a way to do the same job with reshape2?

No, this is something that reshape2 does not support.

Hadley

-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] cast in reshape and reshape2

2011-12-23 Thread Hadley Wickham
 Have you looked at the .summarise argument to dcast? That seems to deliver
 the same sort of results one gets with base::aggregate.


 Actually I see after looking at examples on the plyr-reshape-googlegroups
 group that it is not '.summarise' but rather 'summarise'. Unfortunately
 there are no links in the help pages that seem to describe its proper use
 ... a not uncommon failing for that package in my experience.

?plyr::summarise seems pretty helpful to me.  If you can do better,
please submit a patch - they are very much appreciated.

Hadley

-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
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.


Re: [R] Help creating a symmetric matrix?

2011-12-23 Thread Rui Barradas

Matt Considine wrote
 
 Hi,
 I am trying to work with the output of the MINE analysis routine found at
http://www.exploredata.net
 
 Specifically, I am trying to read the results into a matrix (ideally an 
 n x n x 6 matrix, but I'll settle right now for getting one column into 
 a matrix.)
 
 The problem I have is not knowing how to take what amounts to being one 
 half of a symmetric matrix - excluding the diagonal - and getting it 
 into a matrix.  I have tried using lower.tri as found here
https://stat.ethz.ch/pipermail/r-help/2008-September/174516.html
 but it appears to only partially fill in the matrix.  My code and an 
 example of the output is below.  Can anyone point me to an example that 
 shows how to create a matrix with this sort of input?
 
 Thank you in advance,
 Matt
 
 #v-newx[,3]
 #or, for the sake of this example
 v-c(0.33740, 0.26657, 0.23388, 0.23122, 0.21476, 0.20829, 0.20486, 
 0.19439, 0.19237,
 0.18633, 0.17298, 0.17174, 0.16822, 0.16480, 0.15027)
 z-diag(6)
 ind - lower.tri(z)
 z[ind] - t(v)[ind]
 
 z
  [,1][,2] [,3] [,4] [,5] [,6]
 [1,] 1.0 0.00000
 [2,] 0.26657 1.00000
 [3,] 0.23388 0.192371000
 [4,] 0.23122 0.18633   NA100
 [5,] 0.21476 0.17298   NA   NA10
 [6,] 0.20829 0.17174   NA   NA   NA1
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@ mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

Hello,

Aren't you complicating?

In the last line of your code, why use 'v[ind]' if 'ind' indexes the matrix,
not the vector?

z-diag(6)
ind - lower.tri(z)
z[ind] - v#This works
z

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/Help-creating-a-symmetric-matrix-tp4227301p4230335.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] cast in reshape and reshape2

2011-12-23 Thread David Winsemius


On Dec 23, 2011, at 9:23 PM, Hadley Wickham wrote:

Have you looked at the .summarise argument to dcast? That seems to  
deliver

the same sort of results one gets with base::aggregate.



Actually I see after looking at examples on the plyr-reshape- 
googlegroups
group that it is not '.summarise' but rather 'summarise'.  
Unfortunately
there are no links in the help pages that seem to describe its  
proper use

... a not uncommon failing for that package in my experience.


?plyr::summarise seems pretty helpful to me.  If you can do better,
please submit a patch - they are very much appreciated.


My failure to find it stemmed from it not being mentioned in any way  
in package reshape2's help files, but maybe I was mistaken that it was  
meant to be used in that context?


--

David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Optimising timeboxing in xts

2011-12-23 Thread Hasan Diwan
I don't know if timeboxing is the correct term to use to accomplish
what I'm attempting, so allow me to explain. I have a set n of tagged
observations in time series t. What I'm interested in is taking i
seconds before and after every n. My code is below:
# observations.xts is an xts time series and arg is the number of
seconds to for the timebox
timeboxes - sapply(c(1:as.numeric(last(index(observations.xts))-arg)),
function(x) { return(observations.xts[index(x) - arg, index(x)+arg])},
simplify=TRUE)
What I'd like to see in timeboxes is observations.xts elements with
timestamps within the range of abs(arg).
-- 
Sent from my mobile device
Envoyait de mon portable

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