Re: [R] ftables package, zero rows

2008-09-02 Thread Marc Flockerzi

thanks, but not quite what i wanted.
to be more precise: i want the whole zero-lines to be deleted, including the 
attribute. in my final table only rows with non-zero rows should remain.
regarding my small toy example:
this is what i have:
>test

>  c a b c d
>  a b
>  1 4 2 0 0 0
>  5 0 1 0 1
>  6 0 0 1 0
>  7 0 0 0 0
>  8 0 0 0 0
>  9 0 0 0 0
>  2 4 0 0 0 0
>  5 0 0 0 0
>  6 0 0 0 0
>  7 0 1 0 1
>  8 0 1 0 1
>  9 1 0 0 0


this is what i want:

>  ca b c d
>  ab
>  1  4 2 0 0 0
>  5 0 1 0 1
>  6 0 0 1 0
> 2   7 0 1 0 1
>  8 0 1 0 1
> 9 1 0 0 0

best regards
marc


> -Ursprüngliche Nachricht-
> Von: "Henrique Dallazuanna" <[EMAIL PROTECTED]>
> Gesendet: 02.09.08 18:33:53
> An: "Marc Flockerzi" <[EMAIL PROTECTED]>
> CC: r-help@r-project.org
> Betreff: Re: [R] ftables package, zero rows


> 
> Try this:
> 
> test[test == 0] <- ''
> test
> 
> 
> On Tue, Sep 2, 2008 at 1:03 PM, Marc Flockerzi <[EMAIL PROTECTED]>
> wrote:
> dear all,
> 
>  i'm just about to do some straightforward contingency tables using 
> ftables (and ctab() for percents).
> 
>  the problem:
>  factor "a" are regions, factor "b" are subregions.
>  every region "a" consists of some subregions "b", but obviously not 
> every subregion "b" is part of every region "a".
>  if i use the ftable() function, the table contains a lot of zero 
> rows which i don't want in my output.
> 
>  minimal example:
> 
>  a <- c(1,1,1,1,1,2,2,2,2,2)
>  > b <- c(4,5,6,5,4,7,8,9,8,7)
>  > c <- c("a","b","c","d","a","b","b","a","d","d")
>  > A <- cbind(a,b,c)
>  > A
>  a b c
>  [1,] "1" "4" "a"
>  [2,] "1" "5" "b"
>  [3,] "1" "6" "c"
>  [4,] "1" "5" "d"
>  [5,] "1" "4" "a"
>  [6,] "2" "7" "b"
>  [7,] "2" "8" "b"
>  [8,] "2" "9" "a"
>  [9,] "2" "8" "d"
>  [10,] "2" "7" "d"
>  > test <- ftable(a,b,c)
>  > test
>  c a b c d
>  a b
>  1 4 2 0 0 0
>  5 0 1 0 1
>  6 0 0 1 0
>  7 0 0 0 0
>  8 0 0 0 0
>  9 0 0 0 0
>  2 4 0 0 0 0
>  5 0 0 0 0
>  6 0 0 0 0
>  7 0 1 0 1
>  8 0 1 0 1
>  9 1 0 0 0
> 
>  my question: how can i "delete" the zero rows and preserve the 
> structure and attributes of the original table?
>  simply doing something like:
>  test2 <- test[test>0]
>  obviously only returns the non-zero values, but not the nice 
> structure and attributes of the original table.
> 
>  to do it by hand is not an option as the original table has like 
> 2000 rows, 1500 of which are zero...
> 
>  thanks in advance
>  marc
>  __
>  "Hostage" mit Bruce Willis kostenlos anschauen!
> 
>  __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
>  and provide commented, minimal, self-contained, reproducible code.
> 
> 
> 
> -- 
> Henrique Dallazuanna
> Curitiba-Paraná-Brasil
> 25° 25' 40" S 49° 16' 22" O
> 
> 


__
"Hostage" mit Bruce Willis kostenlos anschauen!

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

2008-09-02 Thread Aval Sarri
Does  ?I() help?
See
? I()

Regards
Aval



On Wed, Sep 3, 2008 at 9:30 AM,  <[EMAIL PROTECTED]> wrote:
> R Users:
>
> I'm wondering:
>
>  Why does my logical vector becomes a numeric vector when stuffed into a data 
> frame?  How do I change this so that it's retained as a logical data type?  
> I've tried a couple of things but to no avail.
>
> Here's my example code:
>
> # Exercise 4-1 in Non-Detects and Data Analysis. Dennis Helsel.  2005.
>
> # "Create two new variables in the Interval Endpoints format, StartCu and 
> EndCu,
> # that will contain the same information given by the current variables
>
> library(NADA)
> data(CuZn)
> names(CuZn)
>
> StartCu <- ifelse(CuZn$CuCen == "TRUE", 0, CuZn$Cu)
> EndCu <- CuZn$Cu
>
> CuIEP = data.frame(cbind(Cu = CuZn$Cu, CuCen = CuZn$CuCen, StartCu, EndCu))
>
> class(CuZn$CuCen)
> #returns "logical"
> class(CuIEP$CuCen)
> #returns "numeric"
>
> CuIEP2 = data.frame(cbind(Cu = CuZn$Cu, CuCen = as.logical(CuZn$CuCen), 
> StartCu, EndCu))
> class(CuIEP2$CuCen)
> #returns "numeric"
>
> CuIEP3 = data.frame(cbind(Cu = CuZn$Cu, CuCen = I(CuZn$CuCen), StartCu, 
> EndCu))
> class(CuIEP3$CuCen)
> #returns "numeric"
>
> I think that I might be missing something fairly fundamental about data 
> coercion in R.   ... would love to figure this out.
>
> Any assistance would be appreciated.
>
> Thanks much,
>
> Matt Findley
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] nls.control()

2008-09-02 Thread Dieter Menne
Ranney, Steven  montana.edu> writes:

> 
> and I'm trying to fit a simple Von Bertalanffy growth curve with program:
> 
> VonB=nls(TL~Linf*(1-exp(-k*(age-t0))), data=box5.4,
> start=list(Linf=1000, k=0.1, t0=0.1), trace=TRUE)
...
> Everything works as it should right up until the confint(VonB) statement.  
When I ask
> for the confidence intervals, R goes through the process then stops and gives 
> 
> "Error in prof$getProfile() : 
>   step factor 0.000488281 reduced below 'minFactor' of 0.000976562".
> 
> However, when I use nls.control=(minFactor=0.01) OR 
nls.control(minFactor=0.01), 

In my experience, it never helps to reduce the control factors once this error 
comes up; it is a problem that by profiling you hit a boundary where 
convergence problems are serious.

Suggestions: first, try with some smaller p-value, e.g. 0.8 or 0.7 instead of 
the default of 0.95; profiling will stay in safer region.

It may be possible to use algorithm "port" in combination with lower/upper, but 
I am not sure if this is honoured by profile.

Dieter

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lattice command equivalent to points

2008-09-02 Thread Deepayan Sarkar
On 9/2/08, Steven McKinney <[EMAIL PROTECTED]> wrote:
>
>
>  This is close, but maybe not optimal lattice coding.
>  I haven't yet figured out how to suppress the x axis
>  labeling.
>
>  bwplot(yield ~ 1|year, panel = function(x, y, ...){panel.bwplot(x, y, ..., 
> pch = "|"); panel.points(x, mean(y), ...,  pch=17)}, data = barley, 
> horizontal = FALSE, xlab = "")
>

A direct translation would be

bwplot(yield ~ year, data = barley, horizontal = FALSE,
   panel = function(x, y, ...) {
   panel.bwplot(x, y, ..., pch = "|");
   mean.values <- tapply(y, x, mean)
   panel.points(1:2, mean.values, ..., pch = 17)
   })

and that seems to work.

-Deepayan

>  Steve McKinney


>
>
>
>  -Original Message-
>  From: [EMAIL PROTECTED] on behalf of Surai
>  Sent: Tue 9/2/2008 5:53 PM
>  To: r-help@r-project.org
>  Subject: [R] lattice command equivalent to points
>
>  Hello,
>  This is my first post to R-help so forgive me if I'm committing a fatal 
> error somehow.  Feel free to let me know.
>
>  My question is this:  I would like to add a triangle to represent the mean 
> in a box plot generated with the bwplot command.   How do I do this?  I am 
> able to do this using the boxplot and points command but points does not work 
> with bwplot.If you run the following code in R, you'll see that I'm 
> trying to reproduce graphs from method 2 by modifying code from method 1.
>
>  Thank you,
>  Surai Jones
>
>
>  library(lattice)
>  attach(barley)
>
>  #method 1
>  bwplot(yield~year, pch="|")
>
>  #method 2
>  boxplot(yield~year)
>  mean.values<-tapply(yield,year, mean)
>  points(1:2, mean.values, pch = 17)
>
>
>
> [[alternative HTML version deleted]]
>
>
> __
>  R-help@r-project.org mailing list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>  and provide commented, minimal, self-contained, reproducible code.
>

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


[R] basic dataframe question

2008-09-02 Thread Matthew.Findley
R Users:

I'm wondering:

 Why does my logical vector becomes a numeric vector when stuffed into a data 
frame?  How do I change this so that it's retained as a logical data type?  
I've tried a couple of things but to no avail.

Here's my example code:

# Exercise 4-1 in Non-Detects and Data Analysis. Dennis Helsel.  2005.

# "Create two new variables in the Interval Endpoints format, StartCu and EndCu,
# that will contain the same information given by the current variables

library(NADA)
data(CuZn)
names(CuZn)

StartCu <- ifelse(CuZn$CuCen == "TRUE", 0, CuZn$Cu)
EndCu <- CuZn$Cu

CuIEP = data.frame(cbind(Cu = CuZn$Cu, CuCen = CuZn$CuCen, StartCu, EndCu))

class(CuZn$CuCen)
#returns "logical"
class(CuIEP$CuCen)
#returns "numeric"

CuIEP2 = data.frame(cbind(Cu = CuZn$Cu, CuCen = as.logical(CuZn$CuCen), 
StartCu, EndCu))
class(CuIEP2$CuCen)
#returns "numeric"

CuIEP3 = data.frame(cbind(Cu = CuZn$Cu, CuCen = I(CuZn$CuCen), StartCu, EndCu))
class(CuIEP3$CuCen)
#returns "numeric"

I think that I might be missing something fairly fundamental about data 
coercion in R.   ... would love to figure this out.

Any assistance would be appreciated.

Thanks much,

Matt Findley

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lattice command equivalent to points

2008-09-02 Thread Steven McKinney


This is close, but maybe not optimal lattice coding.
I haven't yet figured out how to suppress the x axis
labeling.

bwplot(yield ~ 1|year, panel = function(x, y, ...){panel.bwplot(x, y, ..., pch 
= "|"); panel.points(x, mean(y), ...,  pch=17)}, data = barley, horizontal = 
FALSE, xlab = "")

Steve McKinney


-Original Message-
From: [EMAIL PROTECTED] on behalf of Surai
Sent: Tue 9/2/2008 5:53 PM
To: r-help@r-project.org
Subject: [R] lattice command equivalent to points
 
Hello,
This is my first post to R-help so forgive me if I'm committing a fatal error 
somehow.  Feel free to let me know.
 
My question is this:  I would like to add a triangle to represent the mean in a 
box plot generated with the bwplot command.   How do I do this?  I am able to 
do this using the boxplot and points command but points does not work with 
bwplot.    If you run the following code in R, you'll see that I'm trying to 
reproduce graphs from method 2 by modifying code from method 1.
 
Thank you,
Surai Jones
 
 
library(lattice)
attach(barley)
 
#method 1
bwplot(yield~year, pch="|") 
 
#method 2
boxplot(yield~year)
mean.values<-tapply(yield,year, mean)
points(1:2, mean.values, pch = 17)


  
[[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] Free SQL Database with R

2008-09-02 Thread hadley wickham
sqlite

But for a sensible answer, you need to define what you mean by best.

Hadley

On Tue, Sep 2, 2008 at 8:16 AM, Chibisi Chima-Okereke <[EMAIL PROTECTED]> wrote:
> Hi all,
>
> could someone tell me which is the best user-friendly free/open source sql
> database system to use with R, particularly as a back-end for a data-hungry
> web page?
>
> Thanks in advance.
>
> Kind Regards
>
> Chibisi
>
>[[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.
>



-- 
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] Polychoric and tetrachoric correlation

2008-09-02 Thread John Fox
Dear Andy,

Yes, the tetrachoric correlation is a special case of the polychoric
correlation when both factors are dichotomous.

The 95-percent confidence interval that you suggest might be adequate if the
sample size is sufficiently large and the correlation isn't too close to 0
or 1, but it is probably not in general terribly trustworthy.

I hope this helps,
 John 



Original message:

Andy Fugard a.fugard at ed.ac.uk 
Mon Sep 1 19:25:54 CEST 2008 

Hi there,

Am I correct to believe that tetrachoric correlation is a special case 
of polychoric correlation when there are only two levels to the ordered 
factor?  Thus it should be okay to use hetcor from the polycor package 
to build a matrix of correlations for binary variables?

If this is true, how can one estimate 95% confidence intervals for the 
correlations?  My guess would be

   mat = hetcor(dataframe)

   mat$correlation - (1.96 * mat$std.errors)
   mat$correlation + (1.96 * mat$std.errors)

Thanks,

Andy

-- 
Andy Fugard, Postgraduate Research Student
Psychology (Room S6), The University of Edinburgh,
   7 George Square, Edinburgh EH8 9JZ, UK
+44 (0)78 123 87190   http://figuraleffect.googlepages.com/

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox

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

2008-09-02 Thread John Fox
Dear Nguyen D Nguyen,

The tcltk package is part of the standard R distribution. I'm not sure why
you got this error, and you don't provide any details about your friend's
system (OS, version of R, version of the Rcmdr package). To see why the
problem lies, you might try loading just the tcltk package via
library(tcltk). If that doesn't work, then something must be broken in the R
installation or perhaps Tcl/Tk isn't installed on the system.

I hope this helps,
 John

--

Original message:

Nguyen Dinh Nguyen n.nguyen at garvan.org.au 
Mon Sep 1 02:09:28 CEST 2008 
Dear all,
A friend of mine, who just installed Rcmdr package. When calling the
package, the error as following comes up. And actually, I've checked in CRAN
packages, the tcltk package doesn't exist anymore, but tcltk2.
Any help is appreciated
Regards
Nguyen D Nguyen
Garvan Institute of Medical Research
Sydney, Australia

library(Rcmdr)
Loading required package: tcltk
Loading Tcl/Tk interface ...Error in structure(.External("dotTclObjv", objv,
PACKAGE = "tcltk"), class = "tclObj") :
  [tcl] invalid command name "se".

Error : .onLoad failed in 'loadNamespace' for 'tcltk'
Error: package 'tcltk' could not be loaded



--
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] converting values of a dataframe to numeric (where possible)

2008-09-02 Thread jim holtman
Try this as a solution:

> df <- data.frame(a=letters[15:17], b=c("21","NA","23"), c=10:12, d=15:17)
> # convert to numeric
> x <- as.matrix(df)
> mode(x) <- "numeric"
Warning message:
In eval(expr, envir, enclos) : NAs introduced by coercion
> cbind(df, RTot=rowSums(x, na.rm=TRUE))
  a  b  c  d RTot
1 o 21 10 15   46
2 p NA 11 16   27
3 q 23 12 17   52
>


On Tue, Sep 2, 2008 at 5:50 PM, Ralikwen <[EMAIL PROTECTED]> wrote:
>
> Hi,
> I am new to R.
> I have a dataframe with many columns some of which contain genuine strings
> some numerical values as strings. The dataframe is created by cast so I have
> no control over the resulting data type.
> I want to attach columns as aggregates of other columns to the dataframe.
> Here is the solution that I came up with (after a lot of struggle):
>
> castNum <- function(n) {
>x<-as.numeric(as.character(n))
>if (is.na(x)){
>  return(n)
>}else{
>  return(x)
>}
> }
> df <- data.frame(a=letters[15:17], b=c("21","NA","23"), c=10:12, d=15:17)
> cbind(df,RTot=rowSums(as.data.frame(lapply(df, function(x)
> castNum(x)))[2:4],na.rm=TRUE))
>
> This works, but is full of warnings and looks extremely ugly.
> Could you direct me how to achieve the same result in a more elegant way?
>
> Thx.
> Balázs
>
>
>
> --
> View this message in context: 
> http://www.nabble.com/converting-values-of-a-dataframe-to-numeric-%28where-possible%29-tp19279139p19279139.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

What is the problem that you are trying to solve?

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


Re: [R] how to reduce stress value in isoMDS?

2008-09-02 Thread 陈武
Sorry, wrong code. The right one here:

library(MASS)
cl<-read.table("e:/data.txt",header=T,sep=",")
row.names(cl)<-colnames(cl)
cm<-as.matrix(cl)
loc<-sammon(cm)
jpeg(filename="e:/plot.gif",width = 480, height = 480, units = "px",
pointsize = 12, quality = 75, bg = "white", res = NA, restoreConsole = TRUE)
plot(loc$points,type="p")
text(loc$points,rownames(cl),cex=1,pos=1,offset=1)
dev.off()

And "e:/data.txt" contains a 40*40 dissimilarity matrix. Thanks for you
advices!


2008/9/3, ³ÂÎä <[EMAIL PROTECTED]>:
>
> I apply isoMDS to my data, but the result turns out to be bad as the stress
> value stays around 31! Yeah, 31 ,not 3.1... I don't know if I ignore
> something before recall isoMDS.
> My code as follow:
>
> m <- read.table("e:/tsdata.txt",header=T,sep=",")
> article_number <- ts(m, start = 2004,end=2008, frequency = 1
> ,names=colnames(m))
> jpeg(filename="e:/tsmap.gif",width = 480, height = 480, units = "px",
> pointsize = 12, quality = 75, bg = "white", res = NA, restoreConsole = TRUE)
> plot(article_number, plot.type="single",
> lty=c(1,1,1,1,1),col=c(1,2,3,4,5),las=1)
> max<-range(m)
> x<-c(2004,2004,2004,2004,2004)
>
> y<-c(max[2]*0.96,max[2]*0.9199,max[2]*0.88,max[2]*0.84,max[2]*0.7999)
> points(x,y,col=c(1,2,3,4,5),pch=15)
> text(x,y,colnames(m),pos=4,offset=0.4)
> dev.off()
>
> A 40*40 matrix in "e:/tsdata.txt". How should I do to improve the effect?
> Thank you!
>

[[alternative HTML version deleted]]

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


[R] how to reduce stress value in isoMDS?

2008-09-02 Thread 陈武
I apply isoMDS to my data, but the result turns out to be bad as the stress
value stays around 31! Yeah, 31 ,not 3.1... I don't know if I ignore
something before recall isoMDS.
My code as follow:

m <- read.table("e:/tsdata.txt",header=T,sep=",")
article_number <- ts(m, start = 2004,end=2008, frequency = 1
,names=colnames(m))
jpeg(filename="e:/tsmap.gif",width = 480, height = 480, units = "px",
pointsize = 12, quality = 75, bg = "white", res = NA, restoreConsole = TRUE)
plot(article_number, plot.type="single",
lty=c(1,1,1,1,1),col=c(1,2,3,4,5),las=1)
max<-range(m)
x<-c(2004,2004,2004,2004,2004)
y<-c(max[2]*0.96,max[2]*0.9199,max[2]*0.88,max[2]*0.84,max[2]*0.7999)
points(x,y,col=c(1,2,3,4,5),pch=15)
text(x,y,colnames(m),pos=4,offset=0.4)
dev.off()

A 40*40 matrix in "e:/tsdata.txt". How should I do to improve the effect?
Thank you!

[[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] lattice command equivalent to points

2008-09-02 Thread Surai
Hello,
This is my first post to R-help so forgive me if I'm committing a fatal error 
somehow.  Feel free to let me know.
 
My question is this:  I would like to add a triangle to represent the mean in a 
box plot generated with the bwplot command.   How do I do this?  I am able to 
do this using the boxplot and points command but points does not work with 
bwplot.    If you run the following code in R, you'll see that I'm trying to 
reproduce graphs from method 2 by modifying code from method 1.
 
Thank you,
Surai Jones
 
 
library(lattice)
attach(barley)
 
#method 1
bwplot(yield~year, pch="|") 
 
#method 2
boxplot(yield~year)
mean.values<-tapply(yield,year, mean)
points(1:2, mean.values, pch = 17)


  
[[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] R Newbie: quantmod and zoo: Warning in rbind.zoo(...) : column names differ

2008-09-02 Thread Achim Zeileis

On Wed, 3 Sep 2008, Aval Sarri wrote:


Hello;

I am trying following but getting a warning  message : Warning in
rbind.zoo(...) : column names differ, no matter whatever I do.

Also I do not want to specify column names manually, since I am just
writing a wrapper function around getSymbols to get chunks of data
from various sources - oanda, dividends etc.

I tried giving col.names = T/F, header = T/F and skip = 1 but no help.

I think problem is that getSymbols returns a zoo objects whose
index/first column is null.  I write these data in a file and read it
again using read.zoo; when I try to append to these data (using one
more call to getSymbols) it throws this warning message.


The data in objects part1 and part3 are both matrices of dimension 10x1 
and 11x1, respectively, with column name USD.EUR. part2, on the other 
hand, contains a vector of length 10 and consequently has no column names.

Hence, column matching does not work and the c()/rbind() method complains.

The easiest thing is to turn part3 into a vector and then combine
  c(part2, part3[,1])

The reason for part2 being different is the following: zoo seems to drop 
the dimension of 1-column series in a few places, e.g., in read.zoo(). 
I just added a drop = TRUE argument to read.zoo() in the devel version.


hth,
Z




library("quantmod")
options(warn = 1)

part1<-getSymbols(Symbols="USD/EUR", src="oanda", from="2008-01-01",
to="2008-01-10", auto.assign=F, return.class="zoo")
print(dimnames(part1))
write.zoo(part1,"USDEUR", col.names=T) # writes as

part2 <- read.zoo("USDEUR", header=T)
print (dimnames(part2)) # dinames or attributes

part3<-getSymbols(Symbols="USD/EUR", src="oanda", from="2008-01-21",
to="2008-01-31", auto.assign=F, return.class="zoo")
print(dimnames(part3))

allpart <- c(part2, part3)
cat ("-allparts\n")
print(dimnames(allpart))

write.zoo (allpart, "USDEUR.all", col.names=T)


I am not sure how to handle this, please kindly provide some pointer.

Thanks and Regards
-Aval.

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

2008-09-02 Thread Greg Snow
OK, there are a couple of approaches that you can take.

A simple approach is to use the 'locator' function to find the coordinates 
where you want to put the numbers by clicking on the plot, then use that result 
with the 'text' function.  This is fine for a single plot with just a few 
outlier to add, but does not work if you want the process to be fully automated.

You can use the boxplot function with argument 'plot=FALSE' to compute the 
stats, then pass that to the 'bxp' function to do the actual plotting.  The bxp 
function has an 'at' argument so that you can specify where you want the 
individual boxes plotted at (x-coordinate), then use that information along 
with the 'text' function to plot the values.  Or use the bxp function again, 
but capture the return value which is the x-coordinates of the boxes.

If you want to plot the numbers along the top or bottom of the plotting area, 
then 'par("usr")' returns a vector of length 4, the last 2 values are the 
y-coordinates for the bottom of the plotting area and the top respectively.

You can use the function 'updateusr' from the TeachingDemos package to change 
the coordinates of a plot after the fact so that the x-coordinates match with 
what you expect (using information from locator or
par('usr')), but this is probably overkill compared to the above methods.

Hope this helps, if it doesn't, give a bit more detail.


From: Sherri Heck [EMAIL PROTECTED]
Sent: Tuesday, September 02, 2008 5:11 PM
To: Greg Snow
Cc: r-help@r-project.org
Subject: Re: [R] boxplot - label outliers

Hi Greg,

I have the values of the outliers from the boxplot stats.  I am just
having a difficult time adding the values to the plot in the appropriate
places.
Hope this is clearer.

Thanks!

sherri

Greg Snow wrote:
> Are you trying to find the values of the outliers? Or just how to add the 
> text to the plot in the appropriate place? Or both?
>
> --
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> [EMAIL PROTECTED]
> (801) 408-8111
>
>
>
>
>> -Original Message-
>> From: [EMAIL PROTECTED]
>> [mailto:[EMAIL PROTECTED] On Behalf Of Sherri Heck
>> Sent: Tuesday, September 02, 2008 3:38 PM
>> To: r-help@r-project.org
>> Subject: [R] boxplot - label outliers
>>
>> Hi All-
>>
>> I have 24 boxplots on one graph.  I do not have the whiskers
>> extending to the outliers, but I would like to label the
>> maximum value of each outlier above the whiskers.  I have the
>> stats but am having trouble figuring out how to label the whiskers.
>>
>> Any suggestions would be great!
>>
>> sherri
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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] R Newbie: quantmod and zoo: Warning in rbind.zoo(...) : column names differ

2008-09-02 Thread Aval Sarri
Hello;

I am trying following but getting a warning  message : Warning in
rbind.zoo(...) : column names differ, no matter whatever I do.

Also I do not want to specify column names manually, since I am just
writing a wrapper function around getSymbols to get chunks of data
from various sources - oanda, dividends etc.

I tried giving col.names = T/F, header = T/F and skip = 1 but no help.

I think problem is that getSymbols returns a zoo objects whose
index/first column is null.  I write these data in a file and read it
again using read.zoo; when I try to append to these data (using one
more call to getSymbols) it throws this warning message.



library("quantmod")
options(warn = 1)

part1<-getSymbols(Symbols="USD/EUR", src="oanda", from="2008-01-01",
to="2008-01-10", auto.assign=F, return.class="zoo")
print(dimnames(part1))
write.zoo(part1,"USDEUR", col.names=T) # writes as

part2 <- read.zoo("USDEUR", header=T)
print (dimnames(part2)) # dinames or attributes

part3<-getSymbols(Symbols="USD/EUR", src="oanda", from="2008-01-21",
to="2008-01-31", auto.assign=F, return.class="zoo")
print(dimnames(part3))

allpart <- c(part2, part3)
cat ("-allparts\n")
print(dimnames(allpart))

write.zoo (allpart, "USDEUR.all", col.names=T)


I am not sure how to handle this, please kindly provide some pointer.

Thanks and Regards
-Aval.

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

2008-09-02 Thread Sherri Heck

Hi Greg,

I have the values of the outliers from the boxplot stats.  I am just 
having a difficult time adding the values to the plot in the appropriate 
places.

Hope this is clearer.

Thanks!

sherri

Greg Snow wrote:

Are you trying to find the values of the outliers? Or just how to add the text 
to the plot in the appropriate place? Or both?

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



  

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Sherri Heck
Sent: Tuesday, September 02, 2008 3:38 PM
To: r-help@r-project.org
Subject: [R] boxplot - label outliers

Hi All-

I have 24 boxplots on one graph.  I do not have the whiskers
extending to the outliers, but I would like to label the
maximum value of each outlier above the whiskers.  I have the
stats but am having trouble figuring out how to label the whiskers.

Any suggestions would be great!

sherri

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] converting values of a dataframe to numeric (where possible)

2008-09-02 Thread Ralikwen

Hi,
I am new to R.
I have a dataframe with many columns some of which contain genuine strings
some numerical values as strings. The dataframe is created by cast so I have
no control over the resulting data type.
I want to attach columns as aggregates of other columns to the dataframe.
Here is the solution that I came up with (after a lot of struggle):

castNum <- function(n) {
x<-as.numeric(as.character(n))
if (is.na(x)){  
  return(n)  
}else{
  return(x)
}
}
df <- data.frame(a=letters[15:17], b=c("21","NA","23"), c=10:12, d=15:17)
cbind(df,RTot=rowSums(as.data.frame(lapply(df, function(x)
castNum(x)))[2:4],na.rm=TRUE))

This works, but is full of warnings and looks extremely ugly.
Could you direct me how to achieve the same result in a more elegant way?

Thx.
Balázs



-- 
View this message in context: 
http://www.nabble.com/converting-values-of-a-dataframe-to-numeric-%28where-possible%29-tp19279139p19279139.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] question about GLM

2008-09-02 Thread Ben Bolker
Aiste Aistike  gmail.com> writes:

> 
> Hello R-users,
> 
> I do not have much knowledge about generalized linear models therefore my
> question maybe quite stupid.
> 
> I have data from 20 towns with their population and number of people with an
> illness from those towns. I would like to use glm function in R so I
> can calculate proportions of ill people (and later on produce confidence
> intervals). I also want to compare those with original proportions of ill
> people.
> 
> If I use:
> 
>   model1 <- glm(ill ~ offset(log(total)), family = poisson)
> # ill - number of people with illness
> #total - total number of people
> 
> with predict.glm I could get number of people (count data), but not the
> proportions. If the obtained number I divide by 'total', I get the same
> proportion for everyone. But what I want is a way of modeling proportions.
> This probably requires to fit a different model but my lack of knowledge
> isn't helping here.
> 

  Not stupid -- but -- wouldn't a binomial model

 glm(cbind(ill,total-ill) ~ 1, family=binomial)

 make more sense?

  Read ?predict.glm carefully to determine whether you
are predicting responses on the linear predictor (=log-odds)
scale or the original scale

  Ben Bolker

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


Re: [R] How to write the log function in R

2008-09-02 Thread Ben Bolker
toh  yahoo.com> writes:

> 
> 
> Hi,
> 
> I need help to write the following log-likelihood function in R:
> log(L) = sum [(y_k - y_k-1).log(m(t_k) - m(t_k-1)) - (m(t_k) - m(t_k-1)) -
> log((y_k -y_k-1)!)]
> 
> I need to write this log-likelihood function in order to find the parameters
> by MLE method. 
> 

   if you have a vector y_k then diff(y_k) should give you a vector
corresponding to (y_k - y_k-1)  [I assume you mean y_{k-1} rather
than {y_k}-1 above!]

  log(x!) is lfactorial(x) in R

  sum() works as expected

  the rest should be pretty easy ...

  Ben Bolker

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

2008-09-02 Thread Ranney, Steven
All - 

I have data:

TL  age
388 4
418 4
438 4
428 5
539 10
432 4
444 7
421 4
438 4
419 4
463 6
423 4
...
[truncated]

and I'm trying to fit a simple Von Bertalanffy growth curve with program:

#Creates a Von Bertalanffy growth model
VonB=nls(TL~Linf*(1-exp(-k*(age-t0))), data=box5.4,
start=list(Linf=1000, k=0.1, t0=0.1), trace=TRUE)
#Scatterplot of the data
plot(TL~age, data=box5.4, pch=19, xlab="Age (yr)", ylab="Total length (mm)")
#Incorporates the model VonB onto the called scatterplot
mod1a=seq(2, 11, by=0.001)
lines(mod1a, predict(VonB, list(age = mod1a)), col="blue", lwd=2)

#summarizes the model VonB with the model formula, parameter estimates, std. 
errors, 
#and a correlation matrix
summary(VonB, corr=T)
#Provides 95% confidence intervals for the parameter estimates
confint(VonB)

Everything works as it should right up until the confint(VonB) statement.  When 
I ask
for the confidence intervals, R goes through the process then stops and gives 
 
"Error in prof$getProfile() : 
  step factor 0.000488281 reduced below 'minFactor' of 0.000976562".

However, when I use nls.control=(minFactor=0.01) OR 
nls.control(minFactor=0.01), 
I get the same error.  I'm at a loss.  I'm not sure what else I can do to have 
R reduce 
the step size if it's not nls.control().

Thanks, 

SR

Steven H. Ranney
 
Graduate Research Assistant (Ph.D) 
USGS Montana Cooperative Fishery Research Unit 
Montana State University 
P.O. Box 173460 
Bozeman, MT 59717-3460 

phone: (406) 994-6643 
fax: (406) 994-7479


[[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] boxplot - label outliers

2008-09-02 Thread Greg Snow
Are you trying to find the values of the outliers? Or just how to add the text 
to the plot in the appropriate place? Or both?

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Sherri Heck
> Sent: Tuesday, September 02, 2008 3:38 PM
> To: r-help@r-project.org
> Subject: [R] boxplot - label outliers
>
> Hi All-
>
> I have 24 boxplots on one graph.  I do not have the whiskers
> extending to the outliers, but I would like to label the
> maximum value of each outlier above the whiskers.  I have the
> stats but am having trouble figuring out how to label the whiskers.
>
> Any suggestions would be great!
>
> sherri
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] boxplot - label outliers

2008-09-02 Thread Sherri Heck

Hi All-

I have 24 boxplots on one graph.  I do not have the whiskers extending 
to the outliers, but I would like to label the maximum value of each 
outlier above the whiskers.  I have the stats but am having trouble 
figuring out how to label the whiskers.


Any suggestions would be great!

sherri

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] annotating individual panels produced by xyplot

2008-09-02 Thread Martin Brown
That helped -- got it working.  Thanks! -mjb

On Tue, Sep 2, 2008 at 1:28 PM, Deepayan Sarkar
<[EMAIL PROTECTED]> wrote:
> See ?packet.number
>
> -Deepayan
>

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

2008-09-02 Thread Daniel Malter
With that you should probably get advice from your local stats department.
Although you describe your procedure, we do not know your data. And in
particular, we do not know what you do in R. 

Just from inspecting your graph, it looks that your estimated function
undershoots/overshoots the fitted values systematically for certain
intervals of the fit. For example, over the entire last part of the fitted
curve, the actual data points lie predominantly above the fitted curve and
for a long interval before that they lie predominantly below the fitted
curve. This should not be so, which indicates that your fitted function,
despite its relative fit, may not reflect your data generating process well.

Regarding fixing the function in the first observation/data point: That's
wrong. This point would then carry an infinitely greater amount of
information than all the other points (because you assume zero error for
this point). Just imagine you would have a second point like this somewhere
else on the timeline. Then you could perfectly fit your nonlinear function
with two data points. You could only do that if your first point is
nonstochastic, i.e. if there is no error and you would get the EXACT same
value at that point in time every time you run your experiment.

Again, I think it's a question the definition of your function.

Best,
Daniel

-
cuncta stricte discussurus
-

-Ursprüngliche Nachricht-
Von: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Im
Auftrag von LuriFax
Gesendet: Tuesday, September 02, 2008 8:06 AM
An: r-help@r-project.org
Betreff: [R] Help with nonlinear regressional


Dear All,

I am doing experiments in live plant tissue using a laser confocal
microscope. The method is called "fluorescence recovery after
photo-bleaching"  (FRAP) and here follows a short summary:

1. Record/ measure fluorescence intensity in a defined, round region of
interest (ROI, in this case a small spot) to determine the initial intensity
value before the bleaching. This pre-bleach value is also used for
normalising the curve (pre-bleach is then set to 1).

2. Bleach this ROI (with high laser intensity).

3. Record/ measure the recovery of fluorescence over time in the ROI until
it reaches a steady state (a plateau).
.
n. Fit the measured intensity for each time point and mesure the half time
(the timepoint which the curve has reached half the plateau), and more...

The recovery of fluorescence in the ROI is used as a measurement of protein
diffusion in the time range of the experiment. A steep curve means that the
molecules has diffused rapidly into the observed ROI and vice versa.



When I do a regressional curve fit without any constraints I get a huge
deviation from the measured value and the fitted curve at the first data
point in the curve (se the bottom picture).

My question is simply: can I constrain the fitting so that the first point
used in fitting is equal to the measured first point? Also, is this method
of fitting statistically justified / a correct way of doing it when it comes
to statistical error? 

Since the first point in the curve is critical for the calculation of the
halftime I get a substantial deviation when I compare the halftime from a
"automatically" fitted curve (let software decide) and a fitting with a
constrained first-point (y0).

I assume that all measured values have the same amount of noise and
therefore it seems strange that the first residual deviates that strongly
(the curve fit is even not in the range of the standard deviation of the
first point). 


I will greatly appreciate some feedback. Thank you.

---
http://www.nabble.com/file/p19268931/CurveFit_SigmaPlot.png
--
View this message in context:
http://www.nabble.com/Help-with-nonlinear-regressional-tp19268931p19268931.h
tml
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] annotating individual panels produced by xyplot

2008-09-02 Thread Deepayan Sarkar
On Tue, Sep 2, 2008 at 1:23 PM, Martin Brown <[EMAIL PROTECTED]> wrote:
> Hi all,
>
> I'm new to R and lattice and panel functions.  I've produced a lattice
> graph using xyplot.  Now I would like to add various, and *different*,
> annotations to each individual panel.  I tried something like this,
> using panel.text  ...
>
> 
> xyplot(dent ~ era | vegzone,
>   groups=seedtype,
>panel=function(...) {
>panel.xyplot(...)
>panel.text(6.5,50,labels="my 1st annotation")
>panel.text(6.5,70,labels="my 2nd annotation")  
> } )
> 
>
> ... but that didn't work, as both annotations appeared in all panels.
> I want to control which panel each annotation appears in.  I feel
> there must be a way to assign the contents of panel text to individual
> panels, or to make a list of coordinates and annotations and have them
> meted out over the various panels.  However, I've searched through the
> R-help archives and just can't figure it out.
>
> Any tips would be appreciated.  Thanks so much!

See ?packet.number

-Deepayan

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] annotating individual panels produced by xyplot

2008-09-02 Thread Martin Brown
Hi all,

I'm new to R and lattice and panel functions.  I've produced a lattice
graph using xyplot.  Now I would like to add various, and *different*,
annotations to each individual panel.  I tried something like this,
using panel.text  ...


xyplot(dent ~ era | vegzone,
   groups=seedtype,
panel=function(...) {
panel.xyplot(...)
panel.text(6.5,50,labels="my 1st annotation")
panel.text(6.5,70,labels="my 2nd annotation")  
} )


... but that didn't work, as both annotations appeared in all panels.
I want to control which panel each annotation appears in.  I feel
there must be a way to assign the contents of panel text to individual
panels, or to make a list of coordinates and annotations and have them
meted out over the various panels.  However, I've searched through the
R-help archives and just can't figure it out.

Any tips would be appreciated.  Thanks so much!

Martin John Brown
Portland, Oregon, USA
http://martinjohnbrown.net

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] two lattice graphs in one object

2008-09-02 Thread Deepayan Sarkar
On Tue, Sep 2, 2008 at 6:24 AM, Andreas Krause <[EMAIL PROTECTED]> wrote:
>
> When I create a lattice/Trellis type graph, I typically write a function that 
> returns the graph, as in
> do.graph <- function(x, y, ...)
>  {
>require(lattice)
>return(xyplot(y~x, ...))
>  }
>
> My question today is this:
> If I want two graphs on one page, one way of achieving it is to print the 
> objects into defined areas, as in
>
> gr1 <- xyplot(rnorm(111) ~ runif(111))
> gr2 <- xyplot(runif(111) ~ runif(111))
> print(gr1, pos=c(0, 0, 1, 0.5), more=T)
> print(gr2, pos=c(0, 0.5, 1, 1), more=F)
>
> Instead of using the print method, can I create a single trellis object that 
> contains those two "sub-graphs"?
> I do not think so, given what I know about the design of these objects.
> I am hoping for a pleasant surprise though.

Well, you cannot return it as a single "trellis" object. However, you
could always return it as a list with multiple "trellis" objects, and
they will just get printed one by one. You can attach print()
arguments to the objects themselves, so that takes care of the layout.
For example, try

plist <- {
gr1 <-
xyplot(rnorm(111) ~ runif(111), main = "Plot A",
   plot.args = list(position = c(0, 0, 1, 0.5), more = TRUE))
gr2 <-
xyplot(runif(111) ~ runif(111), main = "Plot B",
   plot.args = list(position = c(0, 0.5, 1, 1), more = FALSE))
list(gr1, gr2)
}

print(plist)


Actually you will see some output on the console:

> print(plist)
[[1]]

[[2]]

This is from the print method for lists. If you want to avoid that,
you can always set a class on the list you return and write a
corresponding print() method.

-Deepayan

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

2008-09-02 Thread rkevinburton
Thank you. 
I am not saying the data is wrong. I can do somethiing like:
y = tread + seasonal + remainder
and it gives me back the original data almost exactly.
I just don't know how to interpret it. The data is clearly not periodic but I 
was expecting to get more information about the function that was indicated in 
the seasonal component. Something similar to the impulse response to a function 
generates values at basically all frequencies but different amplitudes. There 
is something different in the response to this function than say what would be 
expected from a Fourier analysis of frequencies.

Kevin

 stephen sefick <[EMAIL PROTECTED]> wrote: 
> .15+.52 #seasonal (.01*52) I think because you said it was periodic
> [1] 0.67
> > .8+.67 #seasonal + trend + positive remainder
> [1] 1.47
> now if you look at the little bit that is in the remainder being
> negative then you can probably subtract about .4ish which is close to
> 1  which is the value of the time series in question, I think.
> 
> Is this example periodic?  Is your data periodic?
> 
> On Tue, Sep 2, 2008 at 12:21 PM,  <[EMAIL PROTECTED]> wrote:
> > There was a typo. I wnated to form an array so it should be:
> >
> > y <- numeric(365)
> >
> > Now you should be able to reproduce it.
> >
> > Kevin
> >
> >  stephen sefick <[EMAIL PROTECTED]> wrote:
> >> I can't reproduce this because the data has two points 0 and one at
> >> the ends of the data set, and I get an na.fail error.  There is no
> >> periodic part to this data- it doesn't seem because there are only two
> >> points.
> >>
> >> stephen
> >>
> >> On Tue, Sep 2, 2008 at 11:38 AM,  <[EMAIL PROTECTED]> wrote:
> >> > I don't understand the output of stl. As a simple example:
> >> >
> >> > y <- numeric(1:365)
> >> > y[250] = 1
> >> >
> >> > stl <- stl(ts(y, frequency=7), s.window="periodic")
> >> >
> >> > This returns without error but the results are puzzling to me. If you 
> >> > plot the results it is probably easiest to visualize what I mean.
> >> >
> >> > plot(stl)
> >> >
> >> > This shows the original data (a single spike at 250). A trend (which 
> >> > also shows a bump at 250). It is the rest that I have a question on. For 
> >> > the "seasonal" component it seems to show a sinusoid like wave with a 
> >> > period roughly a week (7 days) long all with the same amplitude. I can't 
> >> > see how a single spike can generate a "seasonal" component that is 
> >> > periodic for every period in the data. Finally the "remainder" portion 
> >> > of the data generated seems to show just what I want, a representation 
> >> > of the input. But if this is ruly the remainder (data - (trend + 
> >> > seasonal)) then shouldn't it have all entries close to zero? Please help 
> >> > me with my misunderstanding if you have any experience with stl.
> >> >
> >> > Finally it has been suggested that in order to find an overall formula 
> >> > to represent the data a model will need to be constructed. I 
> >> > unfortunately don't have any experience in developing a model. Any hints 
> >> > on where to start?
> >> >
> >> > Thank you.
> >> >
> >> > Kevin
> >> >
> >> > __
> >> > R-help@r-project.org mailing list
> >> > https://stat.ethz.ch/mailman/listinfo/r-help
> >> > PLEASE do read the posting guide 
> >> > http://www.R-project.org/posting-guide.html
> >> > and provide commented, minimal, self-contained, reproducible code.
> >> >
> >>
> >>
> >>
> >> --
> >> Stephen Sefick
> >> Research Scientist
> >> Southeastern Natural Sciences Academy
> >>
> >> Let's not spend our time and resources thinking about things that are
> >> so little or so large that all they really do for us is puff us up and
> >> make us feel like gods. We are mammals, and have not exhausted the
> >> annoying little problems of being mammals.
> >>
> >>   -K. Mullis
> >
> >
> 
> 
> 
> -- 
> Stephen Sefick
> Research Scientist
> Southeastern Natural Sciences Academy
> 
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods. We are mammals, and have not exhausted the
> annoying little problems of being mammals.
> 
>   -K. Mullis

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


Re: [R] More help with stl?

2008-09-02 Thread stephen sefick
what kind of data is this?

On Tue, Sep 2, 2008 at 3:10 PM,  <[EMAIL PROTECTED]> wrote:
> The data is real. The fact that there are a bunch of zeros and only one value 
> of 1 is just the way things are. I have over 20,000 data sets and some are 
> like this. Admittedly this is not periodic but ideally you should see all 
> frequencies at various amplitudes, remniscent of the impulse response to a 
> system. I was not expecting long string of sinusoids all at the same 
> amplitude.
>
> Thank you for your imput. This will help in my understanding.
>
> Kevin
>
>  Ryan Hafen <[EMAIL PROTECTED]> wrote:
>>
>>
>> Trying your example:
>>
>> y <- numeric(365) y y[250] = 1 y
>>
>> y.stl <- stl(ts(y, frequency=7), s.window="periodic")
>>
>> First of all, pay attention to the axes on your plot - the scales are
>> different for each panel. Your seasonal component is quite small in magnitude
>> compared to everything else.
>>
>> Also, if you are unsure that data = seasonal + trend + remainder, just try
>>
>> apply(y.stl$time.series, 1, sum)
>>
>> which adds up the three components. This will get you back your original time
>> series.
>>
>> The problem with your example is that you need to be giving sensible data to
>> the stl procedure. How does data with a bunch of zeros and one 1 represent
>> anything with weekly periodicity? For example, try the following plot:
>>
>> library(lattice) xyplot(y ~ 1:365 | factor(rep(1:7, 53)[1:365]))
>>
>> This groups your data into all Mondays, all Tuesdays, etc. Do you see 
>> anything
>> here indicating periodicity?
>>
>> It was your specification of frequency=7 that created the cyclical pattern 
>> you
>> see in the seasonal component. The STL procedure has a step where it smooths,
>> in this case, each day of the week, and then strings each of those fitted
>> values back together. In the case of this data, it gets a positive value for
>> day 5 (refer to lattice plot above), and hence the seasonal pattern you see.
>>
>> If you read the documentation, you will see that s.window="periodic" causes
>> the mean to be taken for each day of the week, which forces the day of the
>> week to be constant. If you would like the seasonal to be able to adapt, try
>> something like:
>>
>> y.stl <- stl(ts(y, frequency=7), s.window=9, s.degree=1) plot(y.stl)
>>
>> This will use local linear fitting to each week day and allow the seasonal to
>> evolve over time. You can play with s.window argument.
>>
>> If you provided this example just as an example, hopefully this explanation
>> will be helpful. However, if this is what your data actually looks like, I
>> don't see how stl will do you any good.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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.
>



-- 
Stephen Sefick
Research Scientist
Southeastern Natural Sciences Academy

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

-K. Mullis

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


[R] question about GLM

2008-09-02 Thread Aiste Aistike
Hello R-users,

I do not have much knowledge about generalized linear models therefore my
question maybe quite stupid.

I have data from 20 towns with their population and number of people with an
illness from those towns. I would like to use glm function in R so I
can calculate proportions of ill people (and later on produce confidence
intervals). I also want to compare those with original proportions of ill
people.

If I use:

  model1 <- glm(ill ~ offset(log(total)), family = poisson)
# ill - number of people with illness
#total - total number of people

with predict.glm I could get number of people (count data), but not the
proportions. If the obtained number I divide by 'total', I get the same
proportion for everyone. But what I want is a way of modeling proportions.
This probably requires to fit a different model but my lack of knowledge
isn't helping here.

I would appreciate any help.


Aiste

[[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] More help with stl?

2008-09-02 Thread rkevinburton
The data is real. The fact that there are a bunch of zeros and only one value 
of 1 is just the way things are. I have over 20,000 data sets and some are like 
this. Admittedly this is not periodic but ideally you should see all 
frequencies at various amplitudes, remniscent of the impulse response to a 
system. I was not expecting long string of sinusoids all at the same amplitude. 

Thank you for your imput. This will help in my understanding.

Kevin

 Ryan Hafen <[EMAIL PROTECTED]> wrote: 
> 
> 
> Trying your example:
> 
> y <- numeric(365) y y[250] = 1 y
> 
> y.stl <- stl(ts(y, frequency=7), s.window="periodic")
> 
> First of all, pay attention to the axes on your plot - the scales are
> different for each panel. Your seasonal component is quite small in magnitude
> compared to everything else.
> 
> Also, if you are unsure that data = seasonal + trend + remainder, just try
> 
> apply(y.stl$time.series, 1, sum)
> 
> which adds up the three components. This will get you back your original time
> series.
> 
> The problem with your example is that you need to be giving sensible data to
> the stl procedure. How does data with a bunch of zeros and one 1 represent
> anything with weekly periodicity? For example, try the following plot:
> 
> library(lattice) xyplot(y ~ 1:365 | factor(rep(1:7, 53)[1:365]))
> 
> This groups your data into all Mondays, all Tuesdays, etc. Do you see anything
> here indicating periodicity?
> 
> It was your specification of frequency=7 that created the cyclical pattern you
> see in the seasonal component. The STL procedure has a step where it smooths,
> in this case, each day of the week, and then strings each of those fitted
> values back together. In the case of this data, it gets a positive value for
> day 5 (refer to lattice plot above), and hence the seasonal pattern you see.
> 
> If you read the documentation, you will see that s.window="periodic" causes
> the mean to be taken for each day of the week, which forces the day of the
> week to be constant. If you would like the seasonal to be able to adapt, try
> something like:
> 
> y.stl <- stl(ts(y, frequency=7), s.window=9, s.degree=1) plot(y.stl)
> 
> This will use local linear fitting to each week day and allow the seasonal to
> evolve over time. You can play with s.window argument.
> 
> If you provided this example just as an example, hopefully this explanation
> will be helpful. However, if this is what your data actually looks like, I
> don't see how stl will do you any good.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] how to increase row number

2008-09-02 Thread ram basnet
Dear Chuck Cleland,
 
Thank you very much. I got the answer what i searching for.
Thanks again.
 
sincerely,
Ram Kumar Basnet

--- On Tue, 9/2/08, Chuck Cleland <[EMAIL PROTECTED]> wrote:

From: Chuck Cleland <[EMAIL PROTECTED]>
Subject: Re: [R] how to increase row number
To: [EMAIL PROTECTED]
Cc: r-help@r-project.org
Date: Tuesday, September 2, 2008, 9:43 AM

On 9/2/2008 12:30 PM, ram basnet wrote:
> Hi
>  
> I think this may be simple task but creating trouble to me. 
> In my calculation, i have got large number of rows (5546) and column 182
but i am not able to see all the rows. May be somebody have idea how i can see
all the rows ? I will be greatful if any body help me.
> Thanks in advance

  See the max.print argument to options().

> options(max.print = 999)

  should allow you to "see" all of the elements.

?options

> Sincerely,
> Ram Kumar Basnet
> Graduate student
> Wageningen University
> Netherlands
>   
>   [[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.

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




  
[[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] More help with stl?

2008-09-02 Thread Ryan Hafen


Trying your example:

y <- numeric(365) y y[250] = 1 y

y.stl <- stl(ts(y, frequency=7), s.window="periodic")

First of all, pay attention to the axes on your plot - the scales are
different for each panel. Your seasonal component is quite small in magnitude
compared to everything else.

Also, if you are unsure that data = seasonal + trend + remainder, just try

apply(y.stl$time.series, 1, sum)

which adds up the three components. This will get you back your original time
series.

The problem with your example is that you need to be giving sensible data to
the stl procedure. How does data with a bunch of zeros and one 1 represent
anything with weekly periodicity? For example, try the following plot:

library(lattice) xyplot(y ~ 1:365 | factor(rep(1:7, 53)[1:365]))

This groups your data into all Mondays, all Tuesdays, etc. Do you see anything
here indicating periodicity?

It was your specification of frequency=7 that created the cyclical pattern you
see in the seasonal component. The STL procedure has a step where it smooths,
in this case, each day of the week, and then strings each of those fitted
values back together. In the case of this data, it gets a positive value for
day 5 (refer to lattice plot above), and hence the seasonal pattern you see.

If you read the documentation, you will see that s.window="periodic" causes
the mean to be taken for each day of the week, which forces the day of the
week to be constant. If you would like the seasonal to be able to adapt, try
something like:

y.stl <- stl(ts(y, frequency=7), s.window=9, s.degree=1) plot(y.stl)

This will use local linear fitting to each week day and allow the seasonal to
evolve over time. You can play with s.window argument.

If you provided this example just as an example, hopefully this explanation
will be helpful. However, if this is what your data actually looks like, I
don't see how stl will do you any good.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] aov or lme effect size calculation

2008-09-02 Thread Greg Trafton
Sorry about that.  My problem is computational, not statistical and  
exactly as you say:  I don't quite know how to get the correct  
variance component from either aov or lme.  the way to compute partial  
eta squared is:


partial-eta-squared = SS(effect) / (SS(effect) + SS(error))

AOV gives Sum Squares for both effects and the interaction, but lme  
doesn't even give that in default format.


thanks,
greg

On Sep 2, 2008, at 11:43 AM, Doran, Harold wrote:


Greg

You haven't really explained what your problem is. If it is conceptual
(i.e., how do I do it) this is not really the right place for in-depth
statistical advice, but it is often given. OTOH, if your problem is
computational, please explain what that is? For example, maybe you  
know

how to compute eta-squared, but you want to extract the variance
component and you can't figure that out.

Without more info, it is hard to help. Now, with that said, with lme  
(or
mixed models) you have multiple variance components, so how would  
you go

about computing eta-squared anyhow?


-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Greg Trafton
Sent: Tuesday, September 02, 2008 10:25 AM
To: r-help@r-project.org
Subject: [R] aov or lme effect size calculation

(A repost of this request with a bit more detail)

Hi, All.  I'd like to calculate effect sizes for aov or lme
and seem to have a bit of a problem.  partial-eta squared
would be my first choice, but I'm open to suggestions.

I have a completely within design with 2 conditions
(condition and palette).

Here is the aov version:


fit.aov <- (aov(correct ~ cond * palette + Error(subject),

data=data))

summary(fit.aov)


Error: subject
 Df  Sum Sq Mean Sq F value Pr(>F) Residuals 15
0.17326 0.01155

Error: Within
Df  Sum Sq Mean Sq F valuePr(>F)
cond  1 0.32890 0.32890  52.047 4.906e-09 ***
palette   1 0.21971 0.21971  34.768 4.447e-07 ***
cond:palette  1 0.50387 0.50387  79.735 1.594e-11 ***
Residuals45 0.28437 0.00632
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

and here is the lme version:


fm1 <- lme(correct ~ cond * palette, random=~1 | subject,

data=data)  > anova(fm1)
numDF denDF  F-value p-value
(Intercept)  145 4031.042  <.0001
cond 145   52.047  <.0001
palette  145   34.768  <.0001
cond:palette 145   79.735  <.0001

Thanks so much!
Greg

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] What questions can be asked given this information?

2008-09-02 Thread JohnDoe365

This is a general statistic question, but maybe ...

Suppose there are 1000 persons.

In 2002, 400 out of those 1000 persons participated in a questionaire and
330 said they own a homepage.
In 2008, 530 out of those 1000 persons participated in a questionaire and
490 said they own a homepage.

What kind of information can I draw from this data?


Given the sample size of 400 resp. 530 out of 1000, can I give an answer to
'significance'?

Kind regards,

Johann
-- 
View this message in context: 
http://www.nabble.com/What-questions-can-be-asked-given-this-information--tp19275609p19275609.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] aov or lme effect size calculation

2008-09-02 Thread Doran, Harold
Greg

Upgrade your packages to the supported versions (lme4 and Matrix), and
use lmer and not lme.

### Example
> example(lmer)
> anova(fm1)
Analysis of Variance Table
 Df Sum Sq Mean Sq F value
Days  1  30032   30032  45.854

Your method for eta-squared with a mixed model is another story,
however.


> -Original Message-
> From: Greg Trafton [mailto:[EMAIL PROTECTED] 
> Sent: Tuesday, September 02, 2008 1:57 PM
> To: Doran, Harold
> Cc: r-help@r-project.org
> Subject: Re: [R] aov or lme effect size calculation
> 
> Sorry about that.  My problem is computational, not 
> statistical and exactly as you say:  I don't quite know how 
> to get the correct variance component from either aov or lme. 
>  the way to compute partial eta squared is:
> 
> partial-eta-squared = SS(effect) / (SS(effect) + SS(error))
> 
> AOV gives Sum Squares for both effects and the interaction, 
> but lme doesn't even give that in default format.
> 
> thanks,
> greg
> 
> On Sep 2, 2008, at 11:43 AM, Doran, Harold wrote:
> 
> > Greg
> >
> > You haven't really explained what your problem is. If it is 
> conceptual 
> > (i.e., how do I do it) this is not really the right place 
> for in-depth 
> > statistical advice, but it is often given. OTOH, if your problem is 
> > computational, please explain what that is? For example, maybe you 
> > know how to compute eta-squared, but you want to extract 
> the variance 
> > component and you can't figure that out.
> >
> > Without more info, it is hard to help. Now, with that said, 
> with lme 
> > (or mixed models) you have multiple variance components, so 
> how would 
> > you go about computing eta-squared anyhow?
> >
> >> -Original Message-
> >> From: [EMAIL PROTECTED]
> >> [mailto:[EMAIL PROTECTED] On Behalf Of Greg Trafton
> >> Sent: Tuesday, September 02, 2008 10:25 AM
> >> To: r-help@r-project.org
> >> Subject: [R] aov or lme effect size calculation
> >>
> >> (A repost of this request with a bit more detail)
> >>
> >> Hi, All.  I'd like to calculate effect sizes for aov or 
> lme and seem 
> >> to have a bit of a problem.  partial-eta squared would be my first 
> >> choice, but I'm open to suggestions.
> >>
> >> I have a completely within design with 2 conditions (condition and 
> >> palette).
> >>
> >> Here is the aov version:
> >>
> >>> fit.aov <- (aov(correct ~ cond * palette + Error(subject),
> >> data=data))
> >>> summary(fit.aov)
> >>
> >> Error: subject
> >>  Df  Sum Sq Mean Sq F value Pr(>F) Residuals 15
> >> 0.17326 0.01155
> >>
> >> Error: Within
> >> Df  Sum Sq Mean Sq F valuePr(>F)
> >> cond  1 0.32890 0.32890  52.047 4.906e-09 ***
> >> palette   1 0.21971 0.21971  34.768 4.447e-07 ***
> >> cond:palette  1 0.50387 0.50387  79.735 1.594e-11 ***
> >> Residuals45 0.28437 0.00632
> >> ---
> >> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >>
> >> and here is the lme version:
> >>
> >>> fm1 <- lme(correct ~ cond * palette, random=~1 | subject,
> >> data=data)  > anova(fm1)
> >> numDF denDF  F-value p-value
> >> (Intercept)  145 4031.042  <.0001
> >> cond 145   52.047  <.0001
> >> palette  145   34.768  <.0001
> >> cond:palette 145   79.735  <.0001
> >>
> >> Thanks so much!
> >> Greg
> >>
> >> __
> >> R-help@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/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] Help with stl

2008-09-02 Thread stephen sefick
#from the ?stl examples
g<-stl(nottem, "per")
g
#in g are the residuals, sesonal trend, and the remainder.  It looks
like you are going to have to #model this decompostion to get at what
you want.

Stephen

On Tue, Sep 2, 2008 at 12:50 AM, Ryan <[EMAIL PROTECTED]> wrote:
>   charter.net> writes:
>
>>
>> I just realized after some tips and a little digging that what I was trying 
>> to
> do "manually" has already been
>> done. I was trying to fit my data using 'lm' then taking the "residual" data
> and trying to do a spectral
>> estimate (for seasonality) usiing fft and then passing the "residual" of all
> of that to arima to get the
>> irregular portion of the time series forecast equation. I found 'stl' that
> advertises that it does all of
>> this (and probably better than my efforts). The only problem that I found was
> that it takes a time-series
>> object. Time-series objects don't handle missing observations. Say I only 
>> have
> two observations for the
>> year (enough to fit a line through). I can easily fit a line through the
> observation points but if I add in
>> zeros for the missing observations least squares will un
>>  doubtedly throw my observations out and fit the "wrong" line. The other 
>> steps
> will more than likely have
>> residual data so I don't have to worry for these steps about missing data. 
>> Has
> anyon!
>>  e used 'stl' with missing observations? If so what sort of tricks did you 
>> use
> to get around the missing data
>> and time series?
>>
>> Thank you.
>>
>> Kevin
>>
>> __
>> R-help  r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
> Someone can correct me if I am wrong, but R's STL implementation
> uses the Fortran version of stl, which in interest of speed does not
> allow for missing observations.  If you use the stl implemented in S,
> you should be able to have missing observations.
>
> If you read the original stl paper, the algorithm is described and is
> pretty straightforward to implement or do your own variation on.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Stephen Sefick
Research Scientist
Southeastern Natural Sciences Academy

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

-K. Mullis

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


Re: [R] Text

2008-09-02 Thread Nair, Murlidharan T
Dani,
You can put text on your graphs using the locator function. Here is an example

hist(rnorm(1000,0,1) #plots a histogram for you
xy.cord<-locator(1) #Once you give this click on the graph to get the 
coordinates
text(xy.cord, "Text you want to write")

Hope this helps.
Cheers../Murli




-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Greg Snow
Sent: Tuesday, September 02, 2008 1:38 PM
To: Dani Valverde; R Help
Subject: Re: [R] Text

In recent versions of R there are the functions grconvertX and grconvertY that 
will convert between the different coordiate systems.  You can use them to 
convert from the device or ndc systems to user coordinates, then place the text 
using the text function.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Dani Valverde
> Sent: Monday, September 01, 2008 11:17 AM
> To: R Help
> Subject: [R] Text
>
> Hello,
> I would like to place some text that should appear in the
> same position of the graphic device, preferably outside the
> plotting area, regardless the x and y axes limits. How can I do this?
> Best,
>
> Dani
>
> --
> Daniel Valverde Saubí
>
> Grup de Biologia Molecular de Llevats
> Facultat de Veterinària de la Universitat Autònoma de
> Barcelona Edifici V, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
>
> Centro de Investigación Biomédica en Red en Bioingeniería,
> Biomateriales y Nanomedicina (CIBER-BBN)
>
> Grup d'Aplicacions Biomèdiques de la RMN Facultat de
> Biociències Universitat Autònoma de Barcelona Edifici Cs, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
> +34 93 5814126
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

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


Re: [R] multinomial estimation output stat question - not R question

2008-09-02 Thread Greg Snow
Mark,

There are a couple of possible things that could be going on here:

In regular ANOVA cases you can have a situation where you have 3 groups, A, B, 
and C where A and C are significantly different from each other, but B lies 
between them in such a way that we cannot say that B is significantly different 
from A or C (variation is large enough that the mean of B could equal that of A 
or C).  Clearly B cannot equal both A and C if C and A are not the same, it is 
just a matter of lack of evidence.  B could be the same as A, or it could be 
the same as C, or it could be something different from either.

So in your case it could be that there is evidence to show a significant 
difference between positive and negative, but not enough to show how neutral 
compares to them.

You could try refitting your model with a different baseline to see if there is 
a significant difference between the new baseline and one of the other levels 
of the factor.

Another possibility is that in some cases of logistic regressions (and that 
could easily carry over to multinomial regressions) you get a large coefficient 
that is very meaningful, but due to the flattness of the likelihood in that 
region, the wald test overestimates the variance by quite a bit and results in 
non-significant conclusions.  Look at the size of the coefficient and the size 
of the standard error estimate, if both are large, then this could be the case 
and you should ignore the wald test and look more at other types of tests.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of
> [EMAIL PROTECTED]
> Sent: Monday, September 01, 2008 6:31 PM
> To: r-help@r-project.org
> Subject: [R] multinomial estimation output stat question -
> not R question
>
> I am estimating a multinomial model with two quantitative
> predictors, X1 and X2, and 3 responses. The responses are
> called neutral, positive and negative with neutral being the
> baseline. There are actually many models being estimated
> because I estimate the model over time and also for various
> parameter sets but that's not important. When I estimate a
> model, since neutral is the baseline and there is no
> interaction term, I get back coefficients
>
> X1 negative
> X2 negative
>
> X1 positive
> X2 positive
>
> Usually the signs of the coefficients are what I would
> expect. Also, I've read about Anova so I think that I kind of
>  understand what that is doing. But, what I'm confused about
> is the following: In some of the models, I can get back wald
> statistics for X1 say, where both the X1 negative Wald stat
> and the X1 positive Wald stat are not significant.
> Yet, the pvalue from the Anova for the X1 variable overall is
> significant ? Is this possible ? I think I'm not
> understanding  the Anova output as well as I thought because,
> to me, that seems inconsistent ?
>
> I understand that the Wald statistics for the particular
> variables are kind of analogous to the t-stats in a regular
> regression in that they are a function of the decrease in
> deviance conditional on all the rest of the variables being
> in the model. The pvalue in the Anova table I thought was
> kind of doing the same thing except not differentiating
> between the factors and just calculating the decrease in
> deviance due to
> X1 overall without regard to the particular factor response ?
>
> If I'm right in my interpretation of the Anova output, then
> can that still happen ?
>
> If I'm wrong about my interpretation, and it can happen, can
> someone tell me where to look for an explanation on why that
> can happen and possibly explain where my interpretation is
> wrong ? I just want to understand my output as best as I can.
>
> If it can't happen, then it's puzzling because it is happening.
>
> Thanks for any insights, comments or references. The output
> is not easily reproducible or else I would reproduce it here.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Text

2008-09-02 Thread Greg Snow
In recent versions of R there are the functions grconvertX and grconvertY that 
will convert between the different coordiate systems.  You can use them to 
convert from the device or ndc systems to user coordinates, then place the text 
using the text function.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
[EMAIL PROTECTED]
(801) 408-8111



> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Dani Valverde
> Sent: Monday, September 01, 2008 11:17 AM
> To: R Help
> Subject: [R] Text
>
> Hello,
> I would like to place some text that should appear in the
> same position of the graphic device, preferably outside the
> plotting area, regardless the x and y axes limits. How can I do this?
> Best,
>
> Dani
>
> --
> Daniel Valverde Saubí
>
> Grup de Biologia Molecular de Llevats
> Facultat de Veterinària de la Universitat Autònoma de
> Barcelona Edifici V, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
>
> Centro de Investigación Biomédica en Red en Bioingeniería,
> Biomateriales y Nanomedicina (CIBER-BBN)
>
> Grup d'Aplicacions Biomèdiques de la RMN Facultat de
> Biociències Universitat Autònoma de Barcelona Edifici Cs, Campus UAB
> 08193 Cerdanyola del Vallès- SPAIN
> +34 93 5814126
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] receiving "Error: cannot allocate vector of size 1.5 Gb"

2008-09-02 Thread Prof Brian Ripley
Please study the rw-FAQ.  With a 2GB address space your chance of getting 
a 1.5GB contiguous block is essentially zero.


On Tue, 2 Sep 2008, Hayes, Daniel wrote:


Dear all,



In my attempt to run the below modelling command in R 2.7.0 under windows XP 
(4GB RAM with /3GB switch set) I receive the following error:

Error: cannot allocate vector of size 1.5 Gb



I have searched a bit and have tried adding: --max-mem-size=3071M  to the 
command line (when set to 3G I get the error that 3072M is too much)

I also run:


memory.size()


[1] 11.26125


memory.size(max=T)


[1] 13.4375



Modelling script:

model.females <- quote(gamlss(WAZ11~cs(sqrtage,df=12)+country, 
sigma.formula=~cs(sqrtage,df=3)+country,

 nu.formula=~cs(sqrtage,df=1), tau.formula=~cs(sqrtage,df=1),

   data=females, family=BCPE, control=con))

fit.females <- eval(model.females)



the females  (1,654KB) that is being modelled by the GAMLSS package contains 
158,533 observations

I have further installed various memory optimization programs under XP but to 
no avail.

I believe that I perhaps need to set the Vcells and Ncells but am not sure 
which nor to what limits.

Any other help in maximizing my RAM usage in R would be great



I am quite a novice so please excuse any obvious mistakes or omissions.

Thank you in advance for your help



Dr. Daniel Hayes


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



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

2008-09-02 Thread Duncan Murdoch

On 9/2/2008 12:42 PM, Raphael Fraser wrote:

Can R plot 3D graphs? If so, how would you plot x^2 + y^2 + z^2 = 1?


There are several ways.  Plotting a function z = f(x,y) is easiest; 
persp or image or contour can do that.  Plotting a general surface is 
harder, but there are examples of how to do this in ?persp3d in the rgl 
package.


Duncan Murdoch

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


Re: [R] how to increase row number

2008-09-02 Thread Chuck Cleland
On 9/2/2008 12:30 PM, ram basnet wrote:
> Hi
>  
> I think this may be simple task but creating trouble to me. 
> In my calculation, i have got large number of rows (5546) and column 182 but 
> i am not able to see all the rows. May be somebody have idea how i can see 
> all the rows ? I will be greatful if any body help me.
> Thanks in advance

  See the max.print argument to options().

> options(max.print = 999)

  should allow you to "see" all of the elements.

?options

> Sincerely,
> Ram Kumar Basnet
> Graduate student
> Wageningen University
> Netherlands
>   
>   [[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.

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

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


Re: [R] More help with stl?

2008-09-02 Thread stephen sefick
.15+.52 #seasonal (.01*52) I think because you said it was periodic
[1] 0.67
> .8+.67 #seasonal + trend + positive remainder
[1] 1.47
now if you look at the little bit that is in the remainder being
negative then you can probably subtract about .4ish which is close to
1  which is the value of the time series in question, I think.

Is this example periodic?  Is your data periodic?

On Tue, Sep 2, 2008 at 12:21 PM,  <[EMAIL PROTECTED]> wrote:
> There was a typo. I wnated to form an array so it should be:
>
> y <- numeric(365)
>
> Now you should be able to reproduce it.
>
> Kevin
>
>  stephen sefick <[EMAIL PROTECTED]> wrote:
>> I can't reproduce this because the data has two points 0 and one at
>> the ends of the data set, and I get an na.fail error.  There is no
>> periodic part to this data- it doesn't seem because there are only two
>> points.
>>
>> stephen
>>
>> On Tue, Sep 2, 2008 at 11:38 AM,  <[EMAIL PROTECTED]> wrote:
>> > I don't understand the output of stl. As a simple example:
>> >
>> > y <- numeric(1:365)
>> > y[250] = 1
>> >
>> > stl <- stl(ts(y, frequency=7), s.window="periodic")
>> >
>> > This returns without error but the results are puzzling to me. If you plot 
>> > the results it is probably easiest to visualize what I mean.
>> >
>> > plot(stl)
>> >
>> > This shows the original data (a single spike at 250). A trend (which also 
>> > shows a bump at 250). It is the rest that I have a question on. For the 
>> > "seasonal" component it seems to show a sinusoid like wave with a period 
>> > roughly a week (7 days) long all with the same amplitude. I can't see how 
>> > a single spike can generate a "seasonal" component that is periodic for 
>> > every period in the data. Finally the "remainder" portion of the data 
>> > generated seems to show just what I want, a representation of the input. 
>> > But if this is ruly the remainder (data - (trend + seasonal)) then 
>> > shouldn't it have all entries close to zero? Please help me with my 
>> > misunderstanding if you have any experience with stl.
>> >
>> > Finally it has been suggested that in order to find an overall formula to 
>> > represent the data a model will need to be constructed. I unfortunately 
>> > don't have any experience in developing a model. Any hints on where to 
>> > start?
>> >
>> > Thank you.
>> >
>> > Kevin
>> >
>> > __
>> > R-help@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide 
>> > http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>> >
>>
>>
>>
>> --
>> Stephen Sefick
>> Research Scientist
>> Southeastern Natural Sciences Academy
>>
>> Let's not spend our time and resources thinking about things that are
>> so little or so large that all they really do for us is puff us up and
>> make us feel like gods. We are mammals, and have not exhausted the
>> annoying little problems of being mammals.
>>
>>   -K. Mullis
>
>



-- 
Stephen Sefick
Research Scientist
Southeastern Natural Sciences Academy

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

-K. Mullis

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


[R] 3D Graphs

2008-09-02 Thread Raphael Fraser
Can R plot 3D graphs? If so, how would you plot x^2 + y^2 + z^2 = 1?

Thanks
Raphael

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] graphics: per mil plus delta notation

2008-09-02 Thread Henrique Dallazuanna
Try this:

plot(1, ylab=expression({delta}^13*C~'\211'~VPDB))

2008/9/2 Nora Hanson <[EMAIL PROTECTED]>

> Hello All,
>
> I am looking for some help with was is probably very simple:
>
> I want the y-axis title '¦Ä13C¡ë VPDB' and similarly for nitrogen on a
> secondary axis. I am able to get the delta notation and per mil notation
> separately but one messes up the other with the inclusion of quotations:
>
> ylab=expression("{delta}^13*C \211 VPDB")
> gives the per mil just fine but not the delta
> ylab=expression({delta}^13*C \211 VPDB)
> gives an error message for the inclusion of the backslash, and
> ylab=expression({delta}^13*C "\211 VPDB")
> gives an error message for the inclustion of the quotation marks.
>
> Apologies for the simple question, any help is very much appreciated!
>
> Thank you in advance,
>
> Nora
>
>
>
> --
> Nora Hanson
> Gatty Marine Institute
> Sea Mammal Research Unit
> University of St. Andrews
> St. Andrews
> Fife KY16 9AL
> Scotland
> Mobile: 07846140350
> [EMAIL PROTECTED]
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Henrique Dallazuanna
Curitiba-Paran¨¢-Brasil
25¡ã 25' 40" S 49¡ã 16' 22" O

[[alternative HTML version deleted]]

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


Re: [R] receiving "Error: cannot allocate vector of size 1.5 Gb"

2008-09-02 Thread Rory.WINSTON
See

http://cran.r-project.org/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-uses_0021


Rory Winston
RBS Global Banking & Markets
Office: +44 20 7085 4476

-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Hayes, Daniel
Sent: 02 September 2008 14:48
To: r-help@r-project.org
Subject: [R] receiving "Error: cannot allocate vector of size 1.5 Gb"

Dear all,



In my attempt to run the below modelling command in R 2.7.0 under windows XP 
(4GB RAM with /3GB switch set) I receive the following error:

Error: cannot allocate vector of size 1.5 Gb



I have searched a bit and have tried adding: --max-mem-size=3071M  to the 
command line (when set to 3G I get the error that 3072M is too much)

I also run:

> memory.size()

[1] 11.26125

> memory.size(max=T)

[1] 13.4375



Modelling script:

model.females <- quote(gamlss(WAZ11~cs(sqrtage,df=12)+country, 
sigma.formula=~cs(sqrtage,df=3)+country,

  nu.formula=~cs(sqrtage,df=1), tau.formula=~cs(sqrtage,df=1),

data=females, family=BCPE, control=con))

fit.females <- eval(model.females)



the females  (1,654KB) that is being modelled by the GAMLSS package contains 
158,533 observations

I have further installed various memory optimization programs under XP but to 
no avail.

I believe that I perhaps need to set the Vcells and Ncells but am not sure 
which nor to what limits.

Any other help in maximizing my RAM usage in R would be great



I am quite a novice so please excuse any obvious mistakes or omissions.

Thank you in advance for your help



Dr. Daniel Hayes


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

***
The Royal Bank of Scotland plc. Registered in Scotland No 90312. Registered 
Office: 36 St Andrew Square, Edinburgh EH2 2YB. 
Authorised and regulated by the Financial Services Authority 

This e-mail message is confidential and for use by the=2...{{dropped:22}}

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

2008-09-02 Thread Henrique Dallazuanna
Try this:

test[test == 0] <- ''
test

On Tue, Sep 2, 2008 at 1:03 PM, Marc Flockerzi <[EMAIL PROTECTED]>wrote:

> dear all,
>
> i'm just about to do some straightforward contingency tables using ftables
> (and ctab() for percents).
>
> the problem:
> factor "a" are regions, factor "b" are subregions.
> every region "a" consists of some subregions "b", but obviously not every
> subregion "b" is part of every region "a".
> if i use the ftable() function, the table contains a lot of zero rows which
> i don't want in my output.
>
> minimal example:
>
> a <- c(1,1,1,1,1,2,2,2,2,2)
> > b <- c(4,5,6,5,4,7,8,9,8,7)
> > c <- c("a","b","c","d","a","b","b","a","d","d")
> > A <- cbind(a,b,c)
> > A
>  a   b   c
>  [1,] "1" "4" "a"
>  [2,] "1" "5" "b"
>  [3,] "1" "6" "c"
>  [4,] "1" "5" "d"
>  [5,] "1" "4" "a"
>  [6,] "2" "7" "b"
>  [7,] "2" "8" "b"
>  [8,] "2" "9" "a"
>  [9,] "2" "8" "d"
> [10,] "2" "7" "d"
> > test <- ftable(a,b,c)
> > test
>c a b c d
> a b
> 1 4   2 0 0 0
>  5   0 1 0 1
>  6   0 0 1 0
>  7   0 0 0 0
>  8   0 0 0 0
>  9   0 0 0 0
> 2 4   0 0 0 0
>  5   0 0 0 0
>  6   0 0 0 0
>  7   0 1 0 1
>  8   0 1 0 1
>  9   1 0 0 0
>
> my question: how can i "delete" the zero rows and preserve the structure
> and attributes of the original table?
>  simply doing something like:
> test2 <- test[test>0]
> obviously only returns the non-zero values, but not the nice structure and
> attributes of the original table.
>
> to do it by hand is not an option as the original table has like 2000 rows,
> 1500 of which are zero...
>
> thanks in advance
> marc
> __
> "Hostage" mit Bruce Willis kostenlos anschauen!
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

[[alternative HTML version deleted]]

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


Re: [R] How is the binding for a super assignment made visible?

2008-09-02 Thread Gabor Grothendieck
You could try using an assign statement instead:

assign("Globals", list(), .GlobalEnv)

and see if it complains about that or not.

On Mon, Sep 1, 2008 at 10:01 PM,  <[EMAIL PROTECTED]> wrote:
> The statement Globals <<- list() in the body of a function in a package was 
> intended to write an empty list to the R workspace to collect results during 
> the computations of the function.  A package name space has not been 
> specified.
>
> The package appears to function correctly, but
>
> during the R CMD check of the package while "checking R code for possible 
> problems ... NOTE",
>
> no visible binding for '<<-' assignment to 'Globals' is displayed.
>
> Can you tell in this case why the binding needs to be visible?  What 
> statement might do that?  A specific reference in the R manuals would be 
> appreciated.
>
> Bill Morphet
> ATK Space Systems, Launch Systems, Nozzle Structural Analysis
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] how to increase row number

2008-09-02 Thread ram basnet
Hi
 
I think this may be simple task but creating trouble to me. 
In my calculation, i have got large number of rows (5546) and column 182 but i 
am not able to see all the rows. May be somebody have idea how i can see all 
the rows ? I will be greatful if any body help me.
Thanks in advance
 
Sincerely,
Ram Kumar Basnet
Graduate student
Wageningen University
Netherlands
 


  
[[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] More help with stl?

2008-09-02 Thread rkevinburton
There was a typo. I wnated to form an array so it should be:

y <- numeric(365)

Now you should be able to reproduce it.

Kevin

 stephen sefick <[EMAIL PROTECTED]> wrote: 
> I can't reproduce this because the data has two points 0 and one at
> the ends of the data set, and I get an na.fail error.  There is no
> periodic part to this data- it doesn't seem because there are only two
> points.
> 
> stephen
> 
> On Tue, Sep 2, 2008 at 11:38 AM,  <[EMAIL PROTECTED]> wrote:
> > I don't understand the output of stl. As a simple example:
> >
> > y <- numeric(1:365)
> > y[250] = 1
> >
> > stl <- stl(ts(y, frequency=7), s.window="periodic")
> >
> > This returns without error but the results are puzzling to me. If you plot 
> > the results it is probably easiest to visualize what I mean.
> >
> > plot(stl)
> >
> > This shows the original data (a single spike at 250). A trend (which also 
> > shows a bump at 250). It is the rest that I have a question on. For the 
> > "seasonal" component it seems to show a sinusoid like wave with a period 
> > roughly a week (7 days) long all with the same amplitude. I can't see how a 
> > single spike can generate a "seasonal" component that is periodic for every 
> > period in the data. Finally the "remainder" portion of the data generated 
> > seems to show just what I want, a representation of the input. But if this 
> > is ruly the remainder (data - (trend + seasonal)) then shouldn't it have 
> > all entries close to zero? Please help me with my misunderstanding if you 
> > have any experience with stl.
> >
> > Finally it has been suggested that in order to find an overall formula to 
> > represent the data a model will need to be constructed. I unfortunately 
> > don't have any experience in developing a model. Any hints on where to 
> > start?
> >
> > Thank you.
> >
> > Kevin
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> 
> 
> 
> -- 
> Stephen Sefick
> Research Scientist
> Southeastern Natural Sciences Academy
> 
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods. We are mammals, and have not exhausted the
> annoying little problems of being mammals.
> 
>   -K. Mullis

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


Re: [R] subsetting the gene list

2008-09-02 Thread Adaikalavan Ramasamy
Have you tried reading some of the material from the BioConductor 
workshop http://bioconductor.org/workshops/ ?


Here is a simplistic way of proceeding:

 ## Calculate pvalues from t-test
 p <- apply( mat, function(x) t.test( x ~ cl )$p.value )

 ## Subset
 mat.sub <- mat[ p, ]

 ## Cluster
 heatmap(m)

Regards, Adai




Abhilash Venu wrote:

Hi all,

I am working on a single color expression data using limma. I would like to
perform a cluster analysis after selecting the differentially genes based on
the P value (say 0.001). As far as my knowledge is concerned I have to do
the sub setting of these selected genes on the normalized data (MA), to
retrieve the distribution across the samples.
But I am wondering whether I can perform using the R script?
I would appreciate any help.



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


Re: [R] How is the binding for a super assignment made visible?

2008-09-02 Thread Thomas Lumley

On Tue, 2 Sep 2008 [EMAIL PROTECTED] wrote:


The statement Globals <<- list() in the body of a function in a package was
intended to write an empty list to the R workspace to collect results during 
the computations of the function.  A package name space has not beenspecified.


The package appears to function correctly, but

during the R CMD check of the package while "checking R code for possible
problems ... NOTE",

no visible binding for '<<-' assignment to 'Globals' is displayed.

Can you tell in this case why the binding needs to be visible?  What
statement might do that?  A specific reference in the R manuals would be
appreciated.


The message really parses as "no binding (as far as we can see)" - it's a 
request for a binding, not a request for visibility. It's phrased that way because it is 
possible to have false positives in this code -- variables that really have been 
previously defined, but don't look as if they have.

You are using <<- on a variable that doesn't appear to have been previously 
defined in your code.  In fact, as you tell us, it wasn't previously defined, so the 
note is correct.  This isn't an error, but is often a bad idea: the assignment will 
overwrite any existing variable called Globals that happens to be in the workspace.

One alternative (since you don't have a namespace) would be to define a variable
  Globals <- list()
at the top level in your package.  The superassignments would then modify that 
variable.

   -thomas



Thomas Lumley   Assoc. Professor, Biostatistics
[EMAIL PROTECTED]   University of Washington, Seattle

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


[R] ftables package, zero rows

2008-09-02 Thread Marc Flockerzi
dear all,

i'm just about to do some straightforward contingency tables using ftables (and 
ctab() for percents).

the problem:
factor "a" are regions, factor "b" are subregions. 
every region "a" consists of some subregions "b", but obviously not every 
subregion "b" is part of every region "a".
if i use the ftable() function, the table contains a lot of zero rows which i 
don't want in my output.

minimal example:

a <- c(1,1,1,1,1,2,2,2,2,2)
> b <- c(4,5,6,5,4,7,8,9,8,7)
> c <- c("a","b","c","d","a","b","b","a","d","d")
> A <- cbind(a,b,c)
> A
  a   b   c  
 [1,] "1" "4" "a"
 [2,] "1" "5" "b"
 [3,] "1" "6" "c"
 [4,] "1" "5" "d"
 [5,] "1" "4" "a"
 [6,] "2" "7" "b"
 [7,] "2" "8" "b"
 [8,] "2" "9" "a"
 [9,] "2" "8" "d"
[10,] "2" "7" "d"
> test <- ftable(a,b,c)
> test
c a b c d
a b  
1 4   2 0 0 0
  5   0 1 0 1
  6   0 0 1 0
  7   0 0 0 0
  8   0 0 0 0
  9   0 0 0 0
2 4   0 0 0 0
  5   0 0 0 0
  6   0 0 0 0
  7   0 1 0 1
  8   0 1 0 1
  9   1 0 0 0

my question: how can i "delete" the zero rows and preserve the structure and 
attributes of the original table?
 simply doing something like:
test2 <- test[test>0] 
obviously only returns the non-zero values, but not the nice structure and 
attributes of the original table.

to do it by hand is not an option as the original table has like 2000 rows, 
1500 of which are zero...

thanks in advance
marc
__
"Hostage" mit Bruce Willis kostenlos anschauen!

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

2008-09-02 Thread stephen sefick
I can't reproduce this because the data has two points 0 and one at
the ends of the data set, and I get an na.fail error.  There is no
periodic part to this data- it doesn't seem because there are only two
points.

stephen

On Tue, Sep 2, 2008 at 11:38 AM,  <[EMAIL PROTECTED]> wrote:
> I don't understand the output of stl. As a simple example:
>
> y <- numeric(1:365)
> y[250] = 1
>
> stl <- stl(ts(y, frequency=7), s.window="periodic")
>
> This returns without error but the results are puzzling to me. If you plot 
> the results it is probably easiest to visualize what I mean.
>
> plot(stl)
>
> This shows the original data (a single spike at 250). A trend (which also 
> shows a bump at 250). It is the rest that I have a question on. For the 
> "seasonal" component it seems to show a sinusoid like wave with a period 
> roughly a week (7 days) long all with the same amplitude. I can't see how a 
> single spike can generate a "seasonal" component that is periodic for every 
> period in the data. Finally the "remainder" portion of the data generated 
> seems to show just what I want, a representation of the input. But if this is 
> ruly the remainder (data - (trend + seasonal)) then shouldn't it have all 
> entries close to zero? Please help me with my misunderstanding if you have 
> any experience with stl.
>
> Finally it has been suggested that in order to find an overall formula to 
> represent the data a model will need to be constructed. I unfortunately don't 
> have any experience in developing a model. Any hints on where to start?
>
> Thank you.
>
> Kevin
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Stephen Sefick
Research Scientist
Southeastern Natural Sciences Academy

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

-K. Mullis

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


Re: [R] Interpolation Problems

2008-09-02 Thread Pedro.Rodriguez
Hi Steve,

It could be the case that you are trying to find values that are not in
the range of values you are providing.

For example,
 
x <- c(1,2,3,4,5)
y <- c(10,11,12,13,14)
xout <- c(0.01,0.02)
approx(x,y,xout,method="linear")

R's output:
$x
[1] 0.01 0.02

$y
[1] NA NA

If you want to see the value of 10 when you Xs are below 1 and 14 when
the Xs are above 5, then code below may help.

Regards,

Pedro


interpolation_test  <- function(data,cum_prob,xout)
{
y   <- vector(length=length(xout))
for(i in 1:length(xout))
{
ValueToCheck <- xout[i]
j   <-1
while(cum_prob[j] < ValueToCheck && j < length(cum_prob) -2)
{
j <- j + 1
}

y0  <- data[j]
x0  <- cum_prob[j]

y1  <- data[j+1]
x1  <- cum_prob[j+1]

if(x0==ValueToCheck)
{
y[i]<- y0 
} else {

y[i]<- y0 +  (ValueToCheck-x0)*(y1-y0)/(x1-x0)
}
}
return(y)
}




-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On Behalf Of Steve Murray
Sent: Monday, September 01, 2008 6:17 PM
To: r-help@r-project.org
Subject: [R] Interpolation Problems


Dear all,

I'm trying to interpolate a dataset to give it twice as many values (I'm
giving the dataset a finer resolution by interpolating from 1 degree to
0.5 degrees) to match that of a corresponding dataset.

I have the data in both a data frame format (longitude column header
values along the top with latitude row header values down the side) or
column format (in the format latitude, longitude, value).

I have used Google to determine 'approxfun' the most appropriate command
to use for this purpose - I may well be wrong here though! Nevertheless,
I've tried using it with the default arguments for the data frame (i.e.
interp <- approxfun(dataset) ) but encounter the following errors:

> interp <- approxfun(JanAv)
Error in approxfun(JanAv) : 
  need at least two non-NA values to interpolate
In addition: Warning message:
In approxfun(JanAv) : collapsing to unique 'x' values


However, there are no NA values! And to double-check this, I did the
following:

> JanAv[is.na(JanAv)] <- 0

...to ensure that there really are no NAs, but receive the same error
message each time.

With regard to the latter 'collapsing to unique 'x' values', I'm not
sure what this means exactly, or how to deal with it.


Any words of wisdom on how I should go about this, or whether I should
use an alternative command (I want to perform a simple (e.g. linear)
interpolation), would be much appreciated.


Many thanks for any advice offered,

Steve

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] graphics: per mil plus delta notation

2008-09-02 Thread Nora Hanson
Hello All,

I am looking for some help with was is probably very simple:

I want the y-axis title '$B&D(B13C$B"s(B VPDB' and similarly for nitrogen 
on a
secondary axis. I am able to get the delta notation and per mil notation
separately but one messes up the other with the inclusion of quotations:

ylab=expression("{delta}^13*C \211 VPDB")
gives the per mil just fine but not the delta
ylab=expression({delta}^13*C \211 VPDB)
gives an error message for the inclusion of the backslash, and
ylab=expression({delta}^13*C "\211 VPDB")
gives an error message for the inclustion of the quotation marks.

Apologies for the simple question, any help is very much appreciated!

Thank you in advance,

Nora



-- 
Nora Hanson
Gatty Marine Institute
Sea Mammal Research Unit
University of St. Andrews
St. Andrews
Fife KY16 9AL
Scotland
Mobile: 07846140350
[EMAIL PROTECTED]

[[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] aov or lme effect size calculation

2008-09-02 Thread Doran, Harold
Greg

You haven't really explained what your problem is. If it is conceptual
(i.e., how do I do it) this is not really the right place for in-depth
statistical advice, but it is often given. OTOH, if your problem is
computational, please explain what that is? For example, maybe you know
how to compute eta-squared, but you want to extract the variance
component and you can't figure that out.

Without more info, it is hard to help. Now, with that said, with lme (or
mixed models) you have multiple variance components, so how would you go
about computing eta-squared anyhow?

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Greg Trafton
> Sent: Tuesday, September 02, 2008 10:25 AM
> To: r-help@r-project.org
> Subject: [R] aov or lme effect size calculation
> 
> (A repost of this request with a bit more detail)
> 
> Hi, All.  I'd like to calculate effect sizes for aov or lme 
> and seem to have a bit of a problem.  partial-eta squared 
> would be my first choice, but I'm open to suggestions.
> 
> I have a completely within design with 2 conditions 
> (condition and palette).
> 
> Here is the aov version:
> 
>  > fit.aov <- (aov(correct ~ cond * palette + Error(subject),
> data=data))
>  > summary(fit.aov)
> 
> Error: subject
>   Df  Sum Sq Mean Sq F value Pr(>F) Residuals 15 
> 0.17326 0.01155
> 
> Error: Within
>  Df  Sum Sq Mean Sq F valuePr(>F)
> cond  1 0.32890 0.32890  52.047 4.906e-09 ***
> palette   1 0.21971 0.21971  34.768 4.447e-07 ***
> cond:palette  1 0.50387 0.50387  79.735 1.594e-11 ***
> Residuals45 0.28437 0.00632
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> and here is the lme version:
> 
>  > fm1 <- lme(correct ~ cond * palette, random=~1 | subject, 
> data=data)  > anova(fm1)
>  numDF denDF  F-value p-value
> (Intercept)  145 4031.042  <.0001
> cond 145   52.047  <.0001
> palette  145   34.768  <.0001
> cond:palette 145   79.735  <.0001
> 
> Thanks so much!
> Greg
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] More help with stl?

2008-09-02 Thread rkevinburton
I don't understand the output of stl. As a simple example:

y <- numeric(1:365)
y[250] = 1

stl <- stl(ts(y, frequency=7), s.window="periodic")

This returns without error but the results are puzzling to me. If you plot the 
results it is probably easiest to visualize what I mean.

plot(stl)

This shows the original data (a single spike at 250). A trend (which also shows 
a bump at 250). It is the rest that I have a question on. For the "seasonal" 
component it seems to show a sinusoid like wave with a period roughly a week (7 
days) long all with the same amplitude. I can't see how a single spike can 
generate a "seasonal" component that is periodic for every period in the data. 
Finally the "remainder" portion of the data generated seems to show just what I 
want, a representation of the input. But if this is ruly the remainder (data - 
(trend + seasonal)) then shouldn't it have all entries close to zero? Please 
help me with my misunderstanding if you have any experience with stl. 

Finally it has been suggested that in order to find an overall formula to 
represent the data a model will need to be constructed. I unfortunately don't 
have any experience in developing a model. Any hints on where to start?

Thank you.

Kevin

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

2008-09-02 Thread LuriFax

Dear All,

I am doing experiments in live plant tissue using a laser confocal
microscope. The method is called "fluorescence recovery after
photo-bleaching"  (FRAP) and here follows a short summary:

1. Record/ measure fluorescence intensity in a defined, round region of
interest (ROI, in this case a small spot) to determine the initial intensity
value before the bleaching. This pre-bleach value is also used for
normalising the curve (pre-bleach is then set to 1).

2. Bleach this ROI (with high laser intensity).

3. Record/ measure the recovery of fluorescence over time in the ROI until
it reaches a steady state (a plateau).
.
n. Fit the measured intensity for each time point and mesure the half time
(the timepoint which the curve has reached half the plateau), and more...

The recovery of fluorescence in the ROI is used as a measurement of protein
diffusion in the time range of the experiment. A steep curve means that the
molecules has diffused rapidly into the observed ROI and vice versa.



When I do a regressional curve fit without any constraints I get a huge
deviation from the measured value and the fitted curve at the first data
point in the curve (se the bottom picture).

My question is simply: can I constrain the fitting so that the first point
used in fitting is equal to the measured first point? Also, is this method
of fitting statistically justified / a correct way of doing it when it comes
to statistical error? 

Since the first point in the curve is critical for the calculation of the
halftime I get a substantial deviation when I compare the halftime from a
"automatically" fitted curve (let software decide) and a fitting with a
constrained first-point (y0).

I assume that all measured values have the same amount of noise and
therefore it seems strange that the first residual deviates that strongly
(the curve fit is even not in the range of the standard deviation of the
first point). 


I will greatly appreciate some feedback. Thank you.

---
http://www.nabble.com/file/p19268931/CurveFit_SigmaPlot.png 
-- 
View this message in context: 
http://www.nabble.com/Help-with-nonlinear-regressional-tp19268931p19268931.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] aov or lme effect size calculation

2008-09-02 Thread Greg Trafton

(A repost of this request with a bit more detail)

Hi, All.  I'd like to calculate effect sizes for aov or lme and seem  
to have a bit of a problem.  partial-eta squared would be my first  
choice, but I'm open to suggestions.


I have a completely within design with 2 conditions (condition and  
palette).


Here is the aov version:

> fit.aov <- (aov(correct ~ cond * palette + Error(subject),  
data=data))

> summary(fit.aov)

Error: subject
 Df  Sum Sq Mean Sq F value Pr(>F)
Residuals 15 0.17326 0.01155

Error: Within
Df  Sum Sq Mean Sq F valuePr(>F)
cond  1 0.32890 0.32890  52.047 4.906e-09 ***
palette   1 0.21971 0.21971  34.768 4.447e-07 ***
cond:palette  1 0.50387 0.50387  79.735 1.594e-11 ***
Residuals45 0.28437 0.00632
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

and here is the lme version:

> fm1 <- lme(correct ~ cond * palette, random=~1 | subject, data=data)
> anova(fm1)
numDF denDF  F-value p-value
(Intercept)  145 4031.042  <.0001
cond 145   52.047  <.0001
palette  145   34.768  <.0001
cond:palette 145   79.735  <.0001

Thanks so much!
Greg

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] two lattice graphs in one object

2008-09-02 Thread Andreas Krause

When I create a lattice/Trellis type graph, I typically write a function that 
returns the graph, as in
do.graph <- function(x, y, ...)
  {
require(lattice)
return(xyplot(y~x, ...))
  }

My question today is this: 
If I want two graphs on one page, one way of achieving it is to print the 
objects into defined areas, as in

gr1 <- xyplot(rnorm(111) ~ runif(111))
gr2 <- xyplot(runif(111) ~ runif(111))
print(gr1, pos=c(0, 0, 1, 0.5), more=T)
print(gr2, pos=c(0, 0.5, 1, 1), more=F)

Instead of using the print method, can I create a single trellis object that 
contains those two "sub-graphs"?
I do not think so, given what I know about the design of these objects. 
I am hoping for a pleasant surprise though.

   Andreas

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] receiving "Error: cannot allocate vector of size 1.5 Gb"

2008-09-02 Thread Hayes, Daniel
Dear all,

 

In my attempt to run the below modelling command in R 2.7.0 under windows XP 
(4GB RAM with /3GB switch set) I receive the following error:

Error: cannot allocate vector of size 1.5 Gb

 

I have searched a bit and have tried adding: --max-mem-size=3071M  to the 
command line (when set to 3G I get the error that 3072M is too much)

I also run:

> memory.size()

[1] 11.26125

> memory.size(max=T)

[1] 13.4375 

 

Modelling script:

model.females <- quote(gamlss(WAZ11~cs(sqrtage,df=12)+country, 
sigma.formula=~cs(sqrtage,df=3)+country, 

  nu.formula=~cs(sqrtage,df=1), tau.formula=~cs(sqrtage,df=1), 

data=females, family=BCPE, control=con))

fit.females <- eval(model.females)

 

the females  (1,654KB) that is being modelled by the GAMLSS package contains 
158,533 observations

I have further installed various memory optimization programs under XP but to 
no avail.

I believe that I perhaps need to set the Vcells and Ncells but am not sure 
which nor to what limits.

Any other help in maximizing my RAM usage in R would be great

 

I am quite a novice so please excuse any obvious mistakes or omissions.

Thank you in advance for your help

 

Dr. Daniel Hayes


[[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] give all combinations

2008-09-02 Thread Martin Maechler
> "CW" == Carl Witthoft <[EMAIL PROTECTED]>
> on Mon, 01 Sep 2008 12:19:07 -0400 writes:

CW> I seem to be missing something here:
CW> given a set X:{a,b,c,whatever...}
CW> the mathematical definition of 'permutation' is the set of all possible 
CW> sequences of the elements of X.
CW> The definition of 'combination' is all elements of 'permutation' which 
CW> cannot be re-ordered to become a different element.

CW> example:  X:{a,b,c}

CW> perm(X) = ab, ab, bc, ba, ca, cb
CW> comb(X) = ab, ac, bc


CW> So maybe a better question for this mailing list is:  Are there 
CW> functions available in any R package which produce perm() and comb() 
CW> (perhaps as well as other standard combinatoric functions) ?

combn()  has been "standard R"  for quite a while now.

As others have mentioned in this thread (different branch of
same thread-tree), for permutations, one can typically easily
get what one wants  via  expand.grid() or  outer()  [or
mapply() .. ].
For that reason, those knowing R well haven't seen a big need
for a permutation() function.

Martin Maechler, ETH Zurich and R Core Team

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 in .local(object, ...) : test vector does not match model !

2008-09-02 Thread Nair, Murlidharan T
I am getting a really strange error when I am using predict on an ksvm model.  
The error is "Error in .local(object, ...) : test vector does not match model 
!". I do understand that this happens when the test vectors do not match the
Model.  But in this case it is not so. I am attaching a portion of both the 
test data used for prediction and the data used to build the model.  I could 
find any difference in the data structure. May be I am missing something can 
someone point that.  Sorry that the data is slightly large even though there 
are only ten points.  Is there a reason why my model will work on 
sequence.data.train and give the above mentioned error on sequence.data.test.  
I looked at the attributes and type and everything seems the same.  Any 
comments are greatly appreciated.
Thanks ../Murli


#===
sequence.data.train<- structure(list(Class = structure(c(2L, 2L, 2L, 2L, 2L, 
2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("-",
"+"), class = "factor"), V1 = structure(c(2L, 1L, 1L, 1L, 1L,
1L, 1L, 4L, 2L, 4L, 2L, 1L, 2L, 4L, 3L, 4L, 4L, 4L, 3L, 3L), .Label = c("a",
"c", "g", "t"), class = "factor"), V2 = structure(c(4L, 3L, 1L,
2L, 4L, 1L, 4L, 3L, 1L, 3L, 4L, 4L, 2L, 4L, 2L, 2L, 4L, 3L, 3L,
1L), .Label = c("a", "c", "g", "t"), class = "factor"), V3 = structure(c(4L,
4L, 3L, 4L, 3L, 4L, 2L, 2L, 3L, 1L, 4L, 1L, 4L, 4L, 4L, 4L, 1L,
4L, 1L, 4L), .Label = c("a", "c", "g", "t"), class = "factor"),
V4 = structure(c(4L, 2L, 1L, 4L, 2L, 4L, 4L, 4L, 1L, 1L,
4L, 1L, 1L, 3L, 3L, 4L, 2L, 4L, 2L, 3L), .Label = c("a",
"c", "g", "t"), class = "factor"), V5 = structure(c(2L, 4L,
4L, 2L, 4L, 1L, 3L, 2L, 4L, 1L, 1L, 1L, 2L, 2L, 1L, 4L, 4L,
2L, 1L, 2L), .Label = c("a", "c", "g", "t"), class = "factor"),
V6 = structure(c(4L, 4L, 3L, 1L, 3L, 3L, 3L, 2L, 2L, 3L,
3L, 2L, 2L, 4L, 1L, 3L, 2L, 4L, 4L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V7 = structure(c(3L, 1L,
1L, 2L, 4L, 1L, 2L, 4L, 1L, 2L, 2L, 1L, 3L, 2L, 3L, 2L, 2L,
2L, 2L, 4L), .Label = c("a", "c", "g", "t"), class = "factor"),
V8 = structure(c(2L, 1L, 1L, 1L, 2L, 1L, 4L, 1L, 3L, 1L,
1L, 1L, 1L, 4L, 2L, 2L, 4L, 1L, 4L, 2L), .Label = c("a",
"c", "g", "t"), class = "factor"), V9 = structure(c(2L, 2L,
4L, 3L, 1L, 4L, 4L, 1L, 2L, 1L, 1L, 4L, 4L, 2L, 4L, 3L, 4L,
3L, 3L, 1L), .Label = c("a", "c", "g", "t"), class = "factor"),
V10 = structure(c(2L, 4L, 2L, 2L, 1L, 4L, 2L, 4L, 2L, 1L,
3L, 4L, 1L, 1L, 3L, 3L, 4L, 3L, 1L, 2L), .Label = c("a",
"c", "g", "t"), class = "factor"), V11 = structure(c(4L,
4L, 4L, 4L, 2L, 3L, 1L, 3L, 4L, 4L, 4L, 3L, 4L, 4L, 1L, 2L,
4L, 1L, 1L, 2L), .Label = c("a", "c", "g", "t"), class = "factor"),
V12 = structure(c(1L, 3L, 4L, 3L, 3L, 2L, 3L, 4L, 3L, 3L,
3L, 1L, 2L, 1L, 3L, 2L, 4L, 1L, 3L, 2L), .Label = c("a",
"c", "g", "t"), class = "factor"), V13 = structure(c(2L,
4L, 2L, 2L, 3L, 4L, 1L, 2L, 1L, 2L, 1L, 4L, 4L, 2L, 4L, 4L,
4L, 4L, 1L, 4L), .Label = c("a", "c", "g", "t"), class = "factor"),
V14 = structure(c(1L, 3L, 3L, 2L, 3L, 2L, 3L, 1L, 3L, 4L,
1L, 1L, 4L, 1L, 4L, 3L, 1L, 3L, 3L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V15 = structure(c(3L,
3L, 3L, 3L, 3L, 2L, 1L, 1L, 1L, 4L, 2L, 4L, 1L, 4L, 1L, 3L,
1L, 1L, 1L, 4L), .Label = c("a", "c", "g", "t"), class = "factor"),
V16 = structure(c(1L, 1L, 1L, 4L, 4L, 4L, 1L, 2L, 1L, 4L,
3L, 4L, 4L, 2L, 4L, 4L, 4L, 1L, 1L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V17 = structure(c(4L,
1L, 3L, 3L, 3L, 1L, 3L, 2L, 3L, 3L, 4L, 4L, 4L, 1L, 3L, 3L,
3L, 1L, 3L, 2L), .Label = c("a", "c", "g", "t"), class = "factor"),
V18 = structure(c(2L, 3L, 4L, 2L, 1L, 1L, 3L, 3L, 2L, 4L,
3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V19 = structure(c(2L,
2L, 4L, 1L, 1L, 4L, 4L, 1L, 2L, 3L, 1L, 2L, 4L, 1L, 1L, 1L,
2L, 4L, 1L, 3L), .Label = c("a", "c", "g", "t"), class = "factor"),
V20 = structure(c(2L, 4L, 4L, 3L, 1L, 3L, 1L, 3L, 1L, 1L,
4L, 3L, 4L, 1L, 3L, 1L, 3L, 3L, 1L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V21 = structure(c(4L,
3L, 4L, 3L, 2L, 4L, 4L, 1L, 2L, 2L, 2L, 1L, 3L, 3L, 3L, 4L,
2L, 1L, 3L, 4L), .Label = c("a", "c", "g", "t"), class = "factor"),
V22 = structure(c(3L, 4L, 2L, 2L, 1L, 2L, 1L, 1L, 4L, 4L,
3L, 2L, 2L, 3L, 2L, 1L, 4L, 1L, 2L, 3L), .Label = c("a",
"c", "g", "t"), class = "factor"), V23 = structure(c(2L,
1L, 1L, 2L, 3L, 1L, 4L, 4L, 2L, 1L, 4L, 4L, 4L, 2L, 3L, 2L,
3L, 1L, 2L, 1L), .Label = c("a", "c", "g", "t"), class = "factor"),
V24 = structure(c(2L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 4L,
4L, 3L, 2L, 1L, 1L, 1L, 1L, 2L, 4L, 1L), .Label = c("a",
"c", "g", "t"), class = "factor"), V25 = structure(c(1L,
1L, 3L, 4L, 4L, 2L, 2L, 1L, 1L, 4L, 4L, 1L, 4L, 1L, 1L, 4L,
1L, 1L, 3L, 1L), .Label = c("a

Re: [R] Help with nonlinear regressional

2008-09-02 Thread Dieter Menne


LuriFax wrote:
> 
> 
> When I do a regressional curve fit without any constraints I get a huge
> deviation from the measured value and the fitted curve at the first data
> point in the curve (se the bottom picture).
> 
Note that this is a text-only list; most people cannot see your figure, I
read it on Nabble where it is possible to view the data


LuriFax wrote:
> 
> My question is simply: can I constrain the fitting so that the first point
> used in fitting is equal to the measured first point? 
> 
Yes, you can. Normalize all points by dividing through the first data point,
and fit the model with a fixed initial value set to 1.


LuriFax wrote:
> 
> Also, is this method of fitting statistically justified / a correct way of
> doing it when it comes to statistical error? 
> 
No, no, no, don't do it. We have similar curves from gastric emptying
(http://www.menne-biomed.de/gastempt/index.html), and 20 years ago some
authority recommended that method of clamping to the first data point. With
the effect that American medical journals, who usually do not have
statistical referees (British do), simply refuse to publish anything that
does not follow that advice. If you look carefully, it is not the first
point that is wrong, but there is to much tension in the fit so that the
first part as a whole is off.

So you have two choices: Either accept the slight deviation of the fit; or
add another parameter. If you have many sets of data, and can use nlme to
give "borrowing strength", the latter approach could work. If you have only
one curve, be careful when adding another parameter. With nls, it is not
that dangerous because this brutal function simply refuses to converge when
there is too much correlation between coefficients. With SigmaPlot, you can
end up with seeminglingly good fits; only when you look at the coefficient
StdDevs, you may note that these are 3.0 plus/minus 4000 or so!

Dieter











-- 
View this message in context: 
http://www.nabble.com/Help-with-nonlinear-regressional-tp19268931p19270899.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] give all combinations

2008-09-02 Thread Adaikalavan Ramasamy
Yuan Jian, sending 9 emails within the span of few seconds all with 
similar text is very confusing to say the least!


Carl, look up combinations() and permutations() in the gtools package.


For two case scenario, you can use combinations()

   v <- c("a","b","c")

   library(gtools)
   tmp <- combinations(3, 2, v,repeats=TRUE)
   apply( tmp, 1, paste, collapse="" )
   [1] "aa" "ab" "ac" "bb" "bc" "cc"


For more than two cases, I don't know of an elegant way except to 
generate all possible permutations and then eliminate those with the 
same ingredients. This function will be slow for large numbers!



   multiple.combinations <- function( vec, times ){

  input <- vector( mode="list", times )
  for(i in 1:times) input[[i]] <- vec

  out <- expand.grid( input )
  out <- apply( out, 1, function(x) paste( sort(x), collapse="" ) )
  unique.out <- unique(out)
  return(unique.out)
   }

  multiple.combinations( v, 3 )
  [1] "aaa" "aab" "aac" "abb" "abc" "acc" "bbb" "bbc" "bcc" "ccc"

  multiple.combinations( v, 6 )
  "aa" "ab" "ac" "bb" "bc" "cc" "aaabbb"
  "aaabbc" "aaabcc" "aaaccc" "aa" "aabbbc" "aabbcc" "aabccc"
  "aa" "ab" "ac" "abbbcc" "abbccc" "ab" "ac"
  "bb" "bc" "cc" "bbbccc" "bb" "bc" "cc"

Regards, Adai


Carl Witthoft wrote:

I seem to be missing something here:

given a set X:{a,b,c,whatever...}
the mathematical definition of 'permutation' is the set of all possible 
sequences of the elements of X.
The definition of 'combination' is all elements of 'permutation' which 
cannot be re-ordered to become a different element.


example:  X:{a,b,c}

perm(X) = ab, ab, bc, ba, ca, cb
comb(X) = ab, ac, bc


So maybe a better question for this mailing list is:  Are there 
functions available in any R package which produce perm() and comb() 
(perhaps as well as other standard combinatoric functions) ?


Carl

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

2008-09-02 Thread Chibisi Chima-Okereke
Hi all,

could someone tell me which is the best user-friendly free/open source sql
database system to use with R, particularly as a back-end for a data-hungry
web page?

Thanks in advance.

Kind Regards

Chibisi

[[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] Non-constant variance and non-Gaussian errors with gnls

2008-09-02 Thread Paul Suckling
I have been using the nls function to fit some simple non-linear
regression models for properties of graphite bricks to historical
datasets. I have then been using these fits to obtain mean predictions
for the properties of the bricks a short time into the future. I have
also been calculating approximate prediction intervals.

The information I have suggests that the assumption of a normal
distribution with constant variance is not necessarily the most
appropriate. I would like to see if I can obtain improved fits and
hence more accurate predictions and prediction intervals by
experimenting with a) a non-constant (time dependent) variance and b)
a non-normal
error distribution.

It looks to me like the gnls function from the nlme R package is
probably the appropriate one to use for both these situations.
However, I have looked at the gnls help files/documentation and am
still left unsure as to how to specify the arguments of the gnls
function in order to achieve what I want. In particular, I am unsure
how to use the params argument.

Is anyone here able to help me out or point me to some documentation
that is likely to help me achieve this?

Thank you.

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

2008-09-02 Thread Abhilash Venu
Hi all,

I am working on a single color expression data using limma. I would like to
perform a cluster analysis after selecting the differentially genes based on
the P value (say 0.001). As far as my knowledge is concerned I have to do
the sub setting of these selected genes on the normalized data (MA), to
retrieve the distribution across the samples.
But I am wondering whether I can perform using the R script?
I would appreciate any help.

-- 

Regards,
Abhilash

[[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] programming

2008-09-02 Thread stephen sefick
smultron is my favorite editor on mac.  It is more or less a text
editor with some text highlighting.  I can see the use of Emacs, but
my thinking isn't quite there yet.
for what its worth

Stephen

On Tue, Sep 2, 2008 at 8:41 AM, Michael Lawrence <[EMAIL PROTECTED]> wrote:
> Am I missing something or does that list not include Emacs/ESS? It's also
> missing TextMate (for the Mac people). There's probably a bunch more stuff
> for Eclipse than it mentions.
>
> Michael
>
> On Mon, Sep 1, 2008 at 6:17 PM, Gabor Grothendieck
> <[EMAIL PROTECTED]>wrote:
>
>> Check out:
>> http://www.sciviews.org/_rgui/projects/Editors.html
>>
>> On Mon, Sep 1, 2008 at 8:52 PM, Yuan Jian <[EMAIL PROTECTED]> wrote:
>> > Hi,
>> >
>> > I am looking for R editor. does anyone know good editor? which tells you
>> > syntax error and it has function to beautify format (insert TAB etc.).
>> >
>> > Yu
>> >
>> >
>> >
>> >
>> >[[alternative HTML version deleted]]
>> >
>> >
>> > __
>> > R-help@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>> >
>> >
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>[[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.
>



-- 
Stephen Sefick
Research Scientist
Southeastern Natural Sciences Academy

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

-K. Mullis

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


Re: [R] programming

2008-09-02 Thread Michael Lawrence
Am I missing something or does that list not include Emacs/ESS? It's also
missing TextMate (for the Mac people). There's probably a bunch more stuff
for Eclipse than it mentions.

Michael

On Mon, Sep 1, 2008 at 6:17 PM, Gabor Grothendieck
<[EMAIL PROTECTED]>wrote:

> Check out:
> http://www.sciviews.org/_rgui/projects/Editors.html
>
> On Mon, Sep 1, 2008 at 8:52 PM, Yuan Jian <[EMAIL PROTECTED]> wrote:
> > Hi,
> >
> > I am looking for R editor. does anyone know good editor? which tells you
> > syntax error and it has function to beautify format (insert TAB etc.).
> >
> > Yu
> >
> >
> >
> >
> >[[alternative HTML version deleted]]
> >
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[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] Interpolation Problems

2008-09-02 Thread Roger Bivand
Steve Murray  hotmail.com> writes:

> 
> 
> Thanks Duncan - a couple of extra points... I should have perhaps pointed 
> out that the data are on a *regular*
> 'box' grid (with each value currently spaced at 1 degree intervals). Also,  
> I'm looking for something
> fairly simple, like a bilinear interpolation (where each new point is 
> created based on the values of the
> four points surrounding it).
> 
> In answer to your question, JanAv is simply the data frame of values. 
> And yes, you're right, I think I'll need
> a 2D interpolation as it's a grid with latitude and longitude values 
> (which as an aside, I guess these need
> to be interpolated differently? In a 1D format??). I think you're also 
> right in that the 'akima' package
> isn't suitable for this job, as it's designed for irregular grids.

Yes, interp() is not intended to points on grid lines. My suggestion would 
be to review the Spatial task view on CRAN, and possibly to post to the 
R-sig-geo list. Note that the metric may need to be spherical if you are
close to the Poles. If you convert JanAv into a SpatialPointsDataFrame, with
an appropriate CRS("+proj=longlat") coordinate reference system, and
define a SpatialPixels object for the new grid, you may be able to use
idw() in gstat with a very limited search radius to do a simple interpolation
on the sphere.

Roger
> 
> Do you, or does anyone, have any suggestions as to what my best option 
> should be?
> 
> Thanks again,
> 
> Steve
>

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


Re: [R] How remove the small outline around R plots in LaTeX?

2008-09-02 Thread Prof Brian Ripley
I think this is already fixed: please try the current version of R (see 
the posting guide).  It is specific to the R.app GUI, so R-sig-mac would 
be the appropriate list.  (I think this is somewhere in that list's 
archives.)


An alternative is to use R's pdf() device, including using dev.copy2pdf.

On Mon, 1 Sep 2008, Kenneth Benoit wrote:


Dear R List,

I noticed something odd when I started saving R graphs in pdf format and 
including them in LaTeX documents using (e.g.) 
\includegraphics{mygraphic.pdf}.  A very thin dotted border is included in 
the LaTeX document around the graphic, and I cannot seem to make it go away. 
I searched Google as well as the R list archives but could not find anything 
on this.  Note that this only seemed to start when I upgraded R to 2.6 or 2.7 
-- I did not have the poblem before this.


I am using R 2.7.0 on Mac OS/X 10.5.4.

Any help greatly appreciated!

Ken

Kenneth Benoit
Professor of Quantitative Social Sciences
Head, Department of Political Science
Trinity College
Dublin 2, Ireland
http://kenbenoit.net
Tel: 353-1-896-2491

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


--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

2008-09-02 Thread Prof Brian Ripley

On Tue, 2 Sep 2008, Keith Jewell wrote:



"Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message
news:[EMAIL PROTECTED]

On Tue, 2 Sep 2008, Keith Jewell wrote:


As possible help to others, and also as a request for comments on how I
might do things better, I describe how I've recently altered my system to
handle this.

I'm on a Windows Server 2003 network and the R installation is accessible
to
many others. Everyone has "read" access to all installation files, but
only
I have write access. I do _not_ have Administrator privileges on the
server,
so I cannot make/change registry entries.

R versions are in an R folder tree, which also holds the version
independent
library folder.

//Server02/stats/R
//Server02/stats/R/R-2.7.0
//Server02/stats/R/R-2.7.1
//Server02/stats/R/R-2.7.1pat
//Server02/stats/R/R-2.7.2
  :
//Server02/stats/R/library

In each version I have edited /etc/Rprofile.site to include the line
.libPaths("//Server02/stats/R/library")

The "default" libraries (base, boot, class...) are installed into the
relevant version specific .../R-n.n.n/library/ folder by the windows
installer program (e.g. R-2.7.2-win32.exe).
Occasionaly, and after installing a new R-version,  I update all the
downloaded libraries in the version independent
//Server02/stats/R/library/
folder with a simple update.packages().

HTH. Comments welcome.


Please read the rw-FAQ for more details and insight, especially in using
Renviron.site (here you want to set R_LIBS_SITE: see ?libPaths) and using
update.packages(checkBuilt=TRUE, ask=FALSE)

For example, my sysadmins have Windows R installed locally (via the MSI
installer), and then have on each machine in etc/Renviron.site have

R_LIBS_SITE=N:/R/library/2.7
R_LIBS_USER=P:/R/win-library/2.7

This allows both common packages and user-specific packages to be stored
on SMB-mounted drives.  (N: is common to all machines, P: is 'personal' --
we use mapped drives rather than shares to allow quick changes of server)

There are two reasons we install locally.  The first is performance: if 20
machines in a lab start R up simultaneously the network load is high.  The
second is that our security settings (and I believe the defaults these
days) disallow use of CHM help on network drives -- we remove the chtml
directories from packages installed on R_LIBS_SITE, so R defaults to text
help for those packages.



Thanks Professor Ripley. I'd looked at the R Windows FAQ but  not carefully
enough.Using Renviron.site does seem much more elegant than Rprofile.site.
I've now included a line in Renviron.site
R_LIBS_SITE=//Server02/stats/R/library

I also note your use of site and personal libraries specific to "first
decimal"
versions of R. Would you recommend others to follow this practice?


Yes: that is the level of compatibility.  (It is done automatically for 
personal libraries created by R: we set R_LIBS_USER because HOME is 
set incorrectly on some of our systems.)


--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

2008-09-02 Thread Gabor Grothendieck
One possibility is to only keep one version when the third
digit in the version number changes and only create new fresh versions
when the first or second digit in the version number changes. Thus when
going from 2.7.1 to 2.7.2 you could just overwrite the 2.7.1 installation with
the new 2.7.2 installation but when going from 2.7.2 to 2.8.0 you would
create a new parallel installation of R.  AFAIK that has always worked in
the past and is what I normally do.

To overwrite your 2.7.1 with 2.7.2, when installing 2.7.2 just enter
the old folder
path in the "Select Destination Location" screen of the installer.
All your old
settings and libraries will be preserved.  (You will no longer have a
2.7.1 version.)

If you do want to create a new parallel installation when 2.8.0 comes out, say,
there are links in
http://batchfiles.googlecode.com
discussing different methods.  One of these
is provided by movedir.bat or copydir.bat in the batchfiles collection which
can be used to move or copy libraries safely as they will never overwrite.

On Mon, Sep 1, 2008 at 10:04 PM,  <[EMAIL PROTECTED]> wrote:
> More and more I am getting warnings from packages that I install that the 
> package was built with 2.7.2 (I am running 2.7.1). I would like to upgrade 
> but don't want to loose all of the packages that I have installed and the 
> settings. Is there a way to just "upgrade" without uninstalling and 
> reinstalling 'R'?
>
> Thank you.
>
> Kevin
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Upgrade 'R'

2008-09-02 Thread Keith Jewell

"Prof Brian Ripley" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> On Tue, 2 Sep 2008, Keith Jewell wrote:
>
>> As possible help to others, and also as a request for comments on how I
>> might do things better, I describe how I've recently altered my system to
>> handle this.
>>
>> I'm on a Windows Server 2003 network and the R installation is accessible 
>> to
>> many others. Everyone has "read" access to all installation files, but 
>> only
>> I have write access. I do _not_ have Administrator privileges on the 
>> server,
>> so I cannot make/change registry entries.
>>
>> R versions are in an R folder tree, which also holds the version 
>> independent
>> library folder.
>>
>> //Server02/stats/R
>> //Server02/stats/R/R-2.7.0
>> //Server02/stats/R/R-2.7.1
>> //Server02/stats/R/R-2.7.1pat
>> //Server02/stats/R/R-2.7.2
>>   :
>> //Server02/stats/R/library
>>
>> In each version I have edited /etc/Rprofile.site to include the line
>> .libPaths("//Server02/stats/R/library")
>>
>> The "default" libraries (base, boot, class...) are installed into the
>> relevant version specific .../R-n.n.n/library/ folder by the windows
>> installer program (e.g. R-2.7.2-win32.exe).
>> Occasionaly, and after installing a new R-version,  I update all the
>> downloaded libraries in the version independent 
>> //Server02/stats/R/library/
>> folder with a simple update.packages().
>>
>> HTH. Comments welcome.
>
> Please read the rw-FAQ for more details and insight, especially in using 
> Renviron.site (here you want to set R_LIBS_SITE: see ?libPaths) and using 
> update.packages(checkBuilt=TRUE, ask=FALSE)
>
> For example, my sysadmins have Windows R installed locally (via the MSI 
> installer), and then have on each machine in etc/Renviron.site have
>
> R_LIBS_SITE=N:/R/library/2.7
> R_LIBS_USER=P:/R/win-library/2.7
>
> This allows both common packages and user-specific packages to be stored 
> on SMB-mounted drives.  (N: is common to all machines, P: is 'personal' --  
> we use mapped drives rather than shares to allow quick changes of server)
>
> There are two reasons we install locally.  The first is performance: if 20 
> machines in a lab start R up simultaneously the network load is high.  The 
> second is that our security settings (and I believe the defaults these 
> days) disallow use of CHM help on network drives -- we remove the chtml 
> directories from packages installed on R_LIBS_SITE, so R defaults to text 
> help for those packages.
>

Thanks Professor Ripley. I'd looked at the R Windows FAQ but  not carefully
enough.Using Renviron.site does seem much more elegant than Rprofile.site.
I've now included a line in Renviron.site
R_LIBS_SITE=//Server02/stats/R/library

I also note your use of site and personal libraries specific to "first 
decimal"
versions of R. Would you recommend others to follow this practice?

Thanks again,

Keith Jewell.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] intercept of 3D line? (Orthogonal regression)

2008-09-02 Thread William Simpson
OK thanks Moshe I will think about your answer.

Cheers
Bill

On Tue, Sep 2, 2008 at 5:55 AM, Moshe Olshansky <[EMAIL PROTECTED]> wrote:
> I do not see why you can not use regression even in this case.
>
> To make things more simple suppose that the exact model is:
>
> y = a + b*x, i.e.
> y1 = a + b*x1
> ...
> yn = a + b*xn
>
> But you can not observe y and x. Instead you observe
> ui = xi + ei (i=1,...,n) and
> vi = yi + di (i=1,...,n)
>
> Now you have
>
> vi = yi + di = a + b*xi + di = a + b*(ui - ei) + di
> = a + b*ui + (di - b*ei)
>
> and under regular assumptions about ei's end di's we get a standard 
> regression problem (note that b is unknown to you but is constant).
>
>
> --- On Tue, 2/9/08, William Simpson <[EMAIL PROTECTED]> wrote:
>
>> From: William Simpson <[EMAIL PROTECTED]>
>> Subject: [R] intercept of 3D line? (Orthogonal regression)
>> To: r-help@r-project.org
>> Received: Tuesday, 2 September, 2008, 4:53 AM
>> I posted before recently about fitting 3D data x, y, z where
>> all have
>> error attached.
>> I want to predict z from x and y; something like
>> z = b0 + b1*x + b2*y
>> But multiple regression is not suitable because all of x,
>> y, and z have errors.
>>
>> I have plotted a 3D scatterplot of some data using rgl. I
>> see that the
>> data form a cigar-shaped cloud. I think multiple regression
>> is only
>> suitable when the points fall on a plane (forgetting about
>> the error
>> in x and y).
>>
>> I now know the right way how to find the best fitting plane
>> to x,y,z
>> data using princomp.
>> But a new problem is how to get the best fitting *line*. I
>> actually
>> know how to do that too using princomp. But there is a
>> mathematical
>> problem: there's no way to specify a line in 3D space
>> in the form
>> z=f(x,y) or in other words with an intercept and slopes.
>> Instead, one way to deal with the problem is to use a
>> parametric
>> version of the line: you use an arbitrary starting point
>> x0, y0, z0
>> and the direction vector of your line (I know how to get
>> the direction
>> vector).
>>
>> BUT how do I get the intercept??? At this point my lines
>> just go
>> through the origin.
>> Do I just use $center from the princomp output modified in
>> some way?
>>
>> Thanks for any help!
>>
>> Cheers
>> Bill
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained,
>> reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Newbie question - how to do it 'right'

2008-09-02 Thread JohnDoe365

I read the documentation carefully, yet I miss the experience to do it right

Suppose I have a data.frame DS with this data:

Homepage?   Willhaveone?When?
HP  WHP WHEN
levels='Yes', 'No'  levels='Yes', 'No', 'NV'levels='Now', 'Half 
Year',
'Year', 'Later'
=
yes NV  NV
no  no  NV
no  yes half year
.
.
.

1300 rows.

This is data I obtained from an online questionnaire.

The participant was  asked whether he/she has a homepage (varaible HP). In
that case zhe question about "Will you have one?" was not applicable nor the
question about the time when he/she will have one.

If the answer was No, the participant could either say 'yes' to the question
if he/she will ever have a homepage (variable WHP). If the answer was yes,
the third question was applicable (variable WHEN) with the values 'Now' ,
'Half Year' 'Year' or 'Later'

So all the three coloumns are of type factor. However I wonder what would be
the best way to analyse and plot that data.

The value 'NV' means not valid, as the question as whole is not useful to
answer when the participant will not have a homepage, he need not be asked
about when he will have one.

What I did so far were things like table(DS(factor(DS['HP'])['yes'])) to get
the amount of answers for for: 'yes I will  have a homepage' and so on for
those who will never have a homepage or if they have not yet a homepage,
when 'Now', 'Half year' , ... they will have one. The 'NV' value ist not of
my interesst and as such left out.

As a result I pasted the sub-results together c(X, Y, Z) and made a barplot.

I found this tedious and cumbersome so I wonder what would be the right
approach to analyse such data.

Regards,

Johann
-- 
View this message in context: 
http://www.nabble.com/Newbie-question---how-to-do-it-%27right%27-tp19266789p19266789.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] plotting glmmPQL function

2008-09-02 Thread rhelpless

hello all,

i'm an R newbie struggling a bit with the glmmPQL function (from the nlme
pack).  i think i've managed to run the model successfully, but can't seem
to plot the resulting function. plot(glmmPQL(...)) plots the residuals
rather than the model...  i know this should be basic, but i can't seem to
figure it out.  many thanks in advance.

j
-- 
View this message in context: 
http://www.nabble.com/plotting-glmmPQL-function-tp19266698p19266698.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] How to create functions with new "defaults"?

2008-09-02 Thread Prof Brian Ripley

On Tue, 2 Sep 2008, David Hajage wrote:


I can't test it right now, but is this ok ?

t.test.new <- function(x, y = NULL,
  alternative = c("two.sided", "less", "greater"),
  mu = 0, paired = TRUE, var.equal = FALSE,
  conf.level = 0.95, ...) {
 t.test(x, y,
  alternative,
  mu, paired, var.equal,
  conf.level, ...)
}


You should name the arguments in the t.test call for safety.  But just

t.test.new <- function(..., paired = TRUE) t.test(..., paired=paired)

is all you need, and for example allows formulae to be used.



2008/9/2 Prabhanjan Tattar <[EMAIL PROTECTED]>


Hi R!
I attempt to clear my question through this example. We know that the
default functions for "t.test" is
 t.test(x, y = NULL,
   alternative = c("two.sided", "less", "greater"),
   mu = 0, paired = FALSE, var.equal = FALSE,
   conf.level = 0.95, ...)
Suppose, I want to create a new function where "paired=TRUE". What is the
simpler way of creating such functions? I did try, but without much
success.
Thanks in advance
--
Prabhanjan N. Tattar
Lead Statistician
CustomerXPS
Bangalore

   [[alternative HTML version deleted]]

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



[[alternative HTML version deleted]]

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



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

2008-09-02 Thread Prof Brian Ripley

On Tue, 2 Sep 2008, Keith Jewell wrote:


As possible help to others, and also as a request for comments on how I
might do things better, I describe how I've recently altered my system to
handle this.

I'm on a Windows Server 2003 network and the R installation is accessible to
many others. Everyone has "read" access to all installation files, but only
I have write access. I do _not_ have Administrator privileges on the server,
so I cannot make/change registry entries.

R versions are in an R folder tree, which also holds the version independent
library folder.

//Server02/stats/R
//Server02/stats/R/R-2.7.0
//Server02/stats/R/R-2.7.1
//Server02/stats/R/R-2.7.1pat
//Server02/stats/R/R-2.7.2
  :
//Server02/stats/R/library

In each version I have edited /etc/Rprofile.site to include the line
.libPaths("//Server02/stats/R/library")

The "default" libraries (base, boot, class...) are installed into the
relevant version specific .../R-n.n.n/library/ folder by the windows
installer program (e.g. R-2.7.2-win32.exe).
Occasionaly, and after installing a new R-version,  I update all the
downloaded libraries in the version independent //Server02/stats/R/library/
folder with a simple update.packages().

HTH. Comments welcome.


Please read the rw-FAQ for more details and insight, especially in using 
Renviron.site (here you want to set R_LIBS_SITE: see ?libPaths) and using 
update.packages(checkBuilt=TRUE, ask=FALSE)


For example, my sysadmins have Windows R installed locally (via the MSI 
installer), and then have on each machine in etc/Renviron.site have


R_LIBS_SITE=N:/R/library/2.7
R_LIBS_USER=P:/R/win-library/2.7

This allows both common packages and user-specific packages to be stored 
on SMB-mounted drives.  (N: is common to all machines, P: is 'personal' 
-- we use mapped drives rather than shares to allow quick changes of 
server)


There are two reasons we install locally.  The first is performance: if 20 
machines in a lab start R up simultaneously the network load is high.  The 
second is that our security settings (and I believe the defaults these 
days) disallow use of CHM help on network drives -- we remove the chtml 
directories from packages installed on R_LIBS_SITE, so R defaults to text 
help for those packages.




Keith Jewell
---
"Leon Yee" <[EMAIL PROTECTED]> wrote in message
news:[EMAIL PROTECTED]

Hello, Kevin

   You can get some hints by browsing in this mailist with the subject of
" Upgrading R means I lose my packages", which were posted several days
ago.

HTH

Leon


[EMAIL PROTECTED] wrote:

More and more I am getting warnings from packages that I install that the
package was built with 2.7.2 (I am running 2.7.1). I would like to
upgrade but don't want to loose all of the packages that I have installed
and the settings. Is there a way to just "upgrade" without uninstalling
and reinstalling 'R'?

Thank you.

Kevin


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



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

2008-09-02 Thread Dieter Menne


rhelpless wrote:
> 
> 'm an R newbie struggling a bit with the glmmPQL function (from the nlme
> pack).  i think i've managed to run the model successfully, but can't seem
> to plot the resulting function. plot(glmmPQL(...)) plots the residuals
> rather than the model
> 

Use predict on the fit; you should probably convert the plus-minus infinity
range by using logistic transformation for plotting.

glmm = glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
family = binomial, data = bacteria)
bacteria$pred = predict(glmm)


Creating a special rhelpless mail account is not considered good etiquette
here, so better use your real name if you want a reply.


Dieter



-- 
View this message in context: 
http://www.nabble.com/plotting-glmmPQL-function-tp19266698p19267495.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] Exporting frequency distributions

2008-09-02 Thread Jim Lemon

Luz Milena Zea Fernandez wrote:

Dear support, we are trying to export a frequency distributions obtanined by 
stament freq(x,variable.labels=NULL,display.na=TRUE,levels=NULL)
of prettyR.
How can I do it?
  

Hi Luz,
If you want to export (save to a file) the output of freq:

sink("myfile.out")
freq(x,variable.labels=NULL,display.na=TRUE,levels=NULL)
sink()

If you want to export the object that freq returns (a list with as many 
components as there are frequency tables), you could try some variant of 
lapply, although I haven't figured out how to get it to print the 
components without also printing the object passed.


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] cluster a distance(analogue)-object using agnes(cluster)

2008-09-02 Thread Birgitle

I try to perform a clustering using an existing dissimilarity matrix that I
calculated using distance (analogue)
I tried two different things. One of them worked and one not and I don`t
understand why.
Here the code:

not working example

library(cluster)
library(analogue)

iris2<-as.data.frame(iris)
str(iris2)
'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 ...

Test.Gower <- distance(iris2, method ="mixed")
Test.Gower.agnes<-agnes(Test.Gower, diss=T)
Fehler in agnes(Test.Gower, diss = T) : 
  (list) Objekt kann nicht nach 'logical' umgewandelt werden
Error in agnes(Test.Gower, diss=T).
(list) object can`t be transformed to "logical"

working example only numerics used:

library(cluster)
library(analogue)

irisPart<-subset(iris, select= Sepal.Length:Petal.Width)
Dist.Gower <- distance(irisPart, method ="mixed")
AgnesA <- agnes(Dist.Gower, method="average", diss=TRUE) 

Would be great if somebody could help me.
The dataset that I would like to use for the clustering also contains
factors.
and gives me the same Error message as in the not working example.

Thanks in advance

B.






-
The art of living is more like wrestling than dancing.
(Marcus Aurelius)
-- 
View this message in context: 
http://www.nabble.com/cluster-a-distance%28analogue%29-object-using-agnes%28cluster%29-tp19267349p19267349.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] How to create functions with new "defaults"?

2008-09-02 Thread David Hajage
I can't test it right now, but is this ok ?

t.test.new <- function(x, y = NULL,
   alternative = c("two.sided", "less", "greater"),
   mu = 0, paired = TRUE, var.equal = FALSE,
   conf.level = 0.95, ...) {
  t.test(x, y,
   alternative,
   mu, paired, var.equal,
   conf.level, ...)
}

2008/9/2 Prabhanjan Tattar <[EMAIL PROTECTED]>

> Hi R!
> I attempt to clear my question through this example. We know that the
> default functions for "t.test" is
>  t.test(x, y = NULL,
>alternative = c("two.sided", "less", "greater"),
>mu = 0, paired = FALSE, var.equal = FALSE,
>conf.level = 0.95, ...)
> Suppose, I want to create a new function where "paired=TRUE". What is the
> simpler way of creating such functions? I did try, but without much
> success.
> Thanks in advance
> --
> Prabhanjan N. Tattar
> Lead Statistician
> CustomerXPS
> Bangalore
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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


[R] How to create functions with new "defaults"?

2008-09-02 Thread Prabhanjan Tattar
Hi R!
I attempt to clear my question through this example. We know that the
default functions for "t.test" is
 t.test(x, y = NULL,
alternative = c("two.sided", "less", "greater"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95, ...)
Suppose, I want to create a new function where "paired=TRUE". What is the
simpler way of creating such functions? I did try, but without much success.
Thanks in advance
-- 
Prabhanjan N. Tattar
Lead Statistician
CustomerXPS
Bangalore

[[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] Interpolation Problems

2008-09-02 Thread Steve Murray

Thanks Duncan - a couple of extra points... I should have perhaps pointed out 
that the data are on a *regular* 'box' grid (with each value currently spaced 
at 1 degree intervals). Also, I'm looking for something fairly simple, like a 
bilinear interpolation (where each new point is created based on the values of 
the four points surrounding it).

In answer to your question, JanAv is simply the data frame of values. And yes, 
you're right, I think I'll need a 2D interpolation as it's a grid with latitude 
and longitude values (which as an aside, I guess these need to be interpolated 
differently? In a 1D format??). I think you're also right in that the 'akima' 
package isn't suitable for this job, as it's designed for irregular grids.

Do you, or does anyone, have any suggestions as to what my best option should 
be?

Thanks again,

Steve



> Date: Mon, 1 Sep 2008 18:45:35 -0400
> From: [EMAIL PROTECTED]
> To: [EMAIL PROTECTED]
> CC: r-help@r-project.org
> Subject: Re: [R] Interpolation Problems
>
> On 01/09/2008 6:17 PM, Steve Murray wrote:
>> Dear all,
>>
>> I'm trying to interpolate a dataset to give it twice as many values (I'm 
>> giving the dataset a finer resolution by interpolating from 1 degree to 0.5 
>> degrees) to match that of a corresponding dataset.
>>
>> I have the data in both a data frame format (longitude column header values 
>> along the top with latitude row header values down the side) or column 
>> format (in the format latitude, longitude, value).
>>
>> I have used Google to determine 'approxfun' the most appropriate command to 
>> use for this purpose - I may well be wrong here though! Nevertheless, I've 
>> tried using it with the default arguments for the data frame (i.e. interp <- 
>> approxfun(dataset) ) but encounter the following errors:
>>
>>> interp <- approxfun(JanAv)
>> Error in approxfun(JanAv) :
>> need at least two non-NA values to interpolate
>> In addition: Warning message:
>> In approxfun(JanAv) : collapsing to unique 'x' values
>>
>>
>> However, there are no NA values! And to double-check this, I did the 
>> following:
>>
>>> JanAv[is.na(JanAv)] <- 0
>>
>> ...to ensure that there really are no NAs, but receive the same error 
>> message each time.
>>
>> With regard to the latter 'collapsing to unique 'x' values', I'm not sure 
>> what this means exactly, or how to deal with it.
>>
>>
>> Any words of wisdom on how I should go about this, or whether I should use 
>> an alternative command (I want to perform a simple (e.g. linear) 
>> interpolation), would be much appreciated.
>
> What is JanAv? approxfun needs to be able to construct x and y values
> to interpolate; it may be that your JanAv object doesn't allow it to do
> that. (The general idea is that it will consider y to be a function of
> x, and will construct a function that takes arbitrary x values and
> returns y values matching those in the dataset, with some sort of
> interpolation between values.)
>
> If you really have longitude and latitude on some sort of grid, you
> probably want a two-dimensional interpolation, not a 1-d interpolation
> as done by approxfun. The interp() function in the akima() package does
> this, but maybe not in the format you need.
>
> Duncan Murdoch

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


Re: [R] convenient way to calculate specificity, sensitivity and accuracy from raw data

2008-09-02 Thread drflxms
Hello Dimitris, Hello Gabor,

absolutely incredible! I can't tell you how happy I am about your code
which worked out of the box and saved me from days of boring and stupid
Excel-handwork. Thank you a thousand times!

Just for other newbies, that might be faced with a similar problem, I'd
like to make a few closing remarks to the way I calculate now:

The read.table command is not necessary in my case, because there are
already ready-to-use data.frames I created with the "reshape" package.
So I started with the line:

pairs<-data.frame(pred=factor(unlist(input.frame[2:21])),ref=factor(input.frame[,22]))
# explanation for other newbies: creates a data.frame named pairs, with
two columns. In the column pred(iction) you have the values from the
columns 2-21 of the original input-data.frame "input.frame" which
corresponds to all the observations the medical doctors made in my
specific case. In the column ref(erence) you have the observations from
gold-standard, which are assumed to be the truth.

pred<-pairs$pred
#saves column "pred" of "pairs" data.frame as vector named "pred"

lab <- pairs$ref
#saves column "ref" of "pairs" data.frame as vector named "lab"

library(caret)
#loads library "caret"

confusionMatrix(pred, ref, positive=1) 
#creates a confusion matrix with sensitivity, specificity, accuracy,
kappa and much more; please see documentation (?confusionMatrix) for
details.

Example output for the data.frame I sent with my original question:

Confusion Matrix and Statistics

  Reference
Prediction   0   1
 0 656 122
 1  24  38
  
Accuracy : 0.8262 
  95% CI : (0.7988, 0.8512)
 No Information Rate : 0.8095 
 P-Value [Acc > NIR] : 0.117  
  
   Kappa : 0.264  
  
 Sensitivity : 0.2375 
 Specificity : 0.9647 
  Pos Pred Value : 0.6129 
  Neg Pred Value : 0.8432

This works not only for input data that consists of 2 result-classes
like true or false, but for data with multiple categories/result-classes
as well! See example output:

Confusion Matrix and Statistics

  Reference
Prediction   0   1  10  11 100 101 110
   0   349  31  60  40  66   1  15
   125  80   1  22   3  17   3
   100   1  24   8   3   1  10
   111   6   5   3   1   0   2
   100   3   1   6   7  24   0   5
   101   0   0   0   0   0   1   0
   110   2   1   4   0   3   0   5

Overall Statistics
  
Accuracy : 0.5786 
  95% CI : (0.5444, 0.6122)
 No Information Rate : 0.4524 
 P-Value [Acc > NIR] : 1.506e-13  
  
   Kappa : 0.3571 

Statistics by Class:

   Sensitivity Specificity Pos Pred Value Neg Pred Value
Class: 00.9184  0.5370 0.6210 0.8885
Class: 10.6667  0.9014 0.5298 0.9419
Class: 10   0.2400  0.9689 0.5106 0.9042
Class: 11   0.0375  0.9803 0.1667 0.9063
Class: 100  0.2400  0.9703 0.5217 0.9043
Class: 101  0.0500  1. 1. 0.9774
Class: 110  0.1250  0.9875 0. 0.9576

This is much more than I ever had expected! (Thank you to Max Kuhn, the
creator of "caret"-package!)

The code from Dimitris (see below) perfectly re-samples the way  I did
the calculation in Excel by hand. Wow! This is very instructive for me.
I never had thought about real programming, cause I always believed this
is much too high for me. But as I now try to understand the code, that
solves "my problem", I'll re-think this. It is still "magic" to me, but
magic, one can learn ;-). I definitely like to become a  so(u)rcerer's
apprentice :-).

So again thank you for your quick and efficient help! Great software,
great community. I am really happy, that I decided "against all odds"
and advice from colleagues not to use SPSS or SAS, but to learn R. I
never had thought, that I might succeed in evaluating the results of our
small study in just a few weeks by my own using R.

Cheers,
Felix.

Dimitris Rizopoulos wrote:
> try something like this:
>
>
> dat <- read.table(textConnection("video 1 2 3 4 5 6 7 8 9 10 11 12 13
> 14 15 16 17 18 19 20 21
> 1  1 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0
> 2  2 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  1
> 3  3 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0
> 4  4 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0
> 5  5 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  1  0
> 6  6 0 0 0 0 0 0 0 0 0  0  0  0  0  1  0  0  0  0  0  0  0
> 7  7 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0
> 8  8 0 0 0 0 0 0 0 0 0  0  0  0  0  0  1  0  0  0  0  0  0
> 9 

Re: [R] Upgrade 'R'

2008-09-02 Thread Keith Jewell
As possible help to others, and also as a request for comments on how I 
might do things better, I describe how I've recently altered my system to 
handle this.

I'm on a Windows Server 2003 network and the R installation is accessible to 
many others. Everyone has "read" access to all installation files, but only 
I have write access. I do _not_ have Administrator privileges on the server, 
so I cannot make/change registry entries.

R versions are in an R folder tree, which also holds the version independent 
library folder.

//Server02/stats/R
//Server02/stats/R/R-2.7.0
//Server02/stats/R/R-2.7.1
//Server02/stats/R/R-2.7.1pat
//Server02/stats/R/R-2.7.2
   :
//Server02/stats/R/library

In each version I have edited /etc/Rprofile.site to include the line
 .libPaths("//Server02/stats/R/library")

The "default" libraries (base, boot, class...) are installed into the 
relevant version specific .../R-n.n.n/library/ folder by the windows 
installer program (e.g. R-2.7.2-win32.exe).
Occasionaly, and after installing a new R-version,  I update all the 
downloaded libraries in the version independent //Server02/stats/R/library/ 
folder with a simple update.packages().

HTH. Comments welcome.

Keith Jewell
---
"Leon Yee" <[EMAIL PROTECTED]> wrote in message 
news:[EMAIL PROTECTED]
> Hello, Kevin
>
>You can get some hints by browsing in this mailist with the subject of 
> " Upgrading R means I lose my packages", which were posted several days 
> ago.
>
> HTH
>
> Leon
>
>
> [EMAIL PROTECTED] wrote:
>> More and more I am getting warnings from packages that I install that the 
>> package was built with 2.7.2 (I am running 2.7.1). I would like to 
>> upgrade but don't want to loose all of the packages that I have installed 
>> and the settings. Is there a way to just "upgrade" without uninstalling 
>> and reinstalling 'R'?
>>
>> Thank you.
>>
>> Kevin

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

2008-09-02 Thread Prof Brian Ripley

On Tue, 2 Sep 2008, [EMAIL PROTECTED] wrote:


Hi Professor Ripley,
(Thank you very much)^infinite for your response. I still have few questions
1. R shows that I have two column in nlaw data.

ncol(nlaw)

[1] 2

corr(nlaw)

[1] 0.7763745


Since we don't have law, we don't know what nlaw is: you recomputed it 
after listing it.


2. Why can't I use boot::corr? I check corr in boot package and it seems to 
be used for computing correlation coefficient,right?


You can, but you need to use it correctly.  I gave you one suggestion that 
most people find simpler to understand.


3. The last question, cor is for computed the correlation variance and 
covariance matrix. Would you give a explanation for  cor(data[ind,])[1,2]

I do not get it.


Try the pieces for yourself.  What does cor(nlaw) give? cor(nlaw)[1,2]?



Appreciate,
Chunhao


Quoting Prof Brian Ripley <[EMAIL PROTECTED]>:


Please look at the help for boot, in particular what it says the form
of 'statistic' should be.  Possibly what you want is

corr.t <- function(data, ind) cor(data[ind,])[1,2]

nlaw appears to be a single-column matrix: that will not work either.

On Mon, 1 Sep 2008, [EMAIL PROTECTED] wrote:


Hi R users,
I have one simple question but I really don't know why I can't get it 
work.

The data was adopted from Efron's An introduction to the bootstrap book.

nlaw

  LSAT  GPA
[1,]  576 3.39
[2,]  635 3.30
[3,]  558 2.81
[4,]  578 3.03
[5,]  666 3.44
[6,]  580 3.07
[7,]  555 3.00
[8,]  661 3.43
[9,]  651 3.36
[10,]  605 3.13
[11,]  653 3.12
[12,]  575 2.74
[13,]  545 2.76
[14,]  572 2.88
[15,]  594 2.96

nlaw<-as.matrix(law[,-1])
corr.t<-function(data){

+return(boot::corr(data))
+}

law.boot<-boot(data=nlaw,statistic=corr.t,R=49)

Error in statistic(data, original, ...) : unused argument(s) (1:15)

many thanks for the help
Chunhao

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

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


--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595






--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] installation problem: package 'mgcv' could not be loaded

2008-09-02 Thread Simon Zhu
Hello all, i'm a newbie of R trying to make some statistical work in R 
environment. Now i have to laptops, one is Thinkpad X40 with Debian 
Lenny and the other is Thinkpad T43 with Ubuntu 8.10. Recently i met 
such problem and am wondering if anybody can do some help.


After upgrading my /etc/apt/sources.list , i install R by apt-get 
install command. It works fine in both laptops. Then i tried to install 
packages such as JGR and TSA. Problem occured when i installed TSA in my 
Ubuntu T43.


> install.packages('TSA')

and i get the following error message:

   Loading required package: leaps
   Loading required package: locfit
   Loading required package: akima
   Loading required package: lattice
   locfit 1.5-4  2007-11-27
   Loading required package: mgcv
   Error in dyn.load(file, DLLpath = DLLpath, ...) :
 unable to load shared library '/usr/lib/R/library/mgcv/libs/mgcv.so':
  libRlapack.so: cannot open shared object file: No such file or 
directory

Error: package 'mgcv' could not be loaded
Execution halted
ERROR: lazy loading failed for package 'TSA'
** Removing '/usr/local/lib/R/site-library/TSA'

The downloaded packages are in
   /tmp/RtmpClJMOJ/downloaded_packages
   Warning message:
   In install.packages("TSA") :
   installation of package 'TSA' had non-zero exit status

The problems is, i used exact the same commands and configuration of R 
in both my Debian-X40 and Ubuntu-T43. It works in Debian but not in T43. 
Anybody who can help me about how to fix it would be a lot appreciated. 
thanks !


--
---
Yanyuan Zhu
Doctoral Candidate
School of Economics & Management, Tongji University
MSN messenger: [EMAIL PROTECTED]
E-mail: [EMAIL PROTECTED]
BLOG http://yyzhu.net

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