Re: [R] Stratifying level-1 variance with lmer()

2008-02-06 Thread Dieter Menne
David Afshartous  med.miami.edu> writes:

> I've fit some models via lme() and now I'm trying to fit similar models with
> lmer() for some simulations I'm running.
> 
> The model below (fm1) has an intercept variance that depends on treatment
> group.  How would one accomplish a similar stratification for the level-1
> variance, i.e., the within-group or patient variance?
> 
> With lme() I was able to do this via:
> update(fm1,
> weights = varIdent(form = ~ 1 | treatment.ind) )
> 
> However, this does not work for lmer().  Is there any way around this?

This is not yet implemented. See Douglas Bates'

http://article.gmane.org/gmane.comp.lang.r.lme4.devel/518

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] How to generate table output of a series of lm's

2008-02-06 Thread Dieter Menne
Martin Elff  sowi.uni-mannheim.de> writes:

> 
> modco <- list(
>   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz),
>   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz, 
> subset=sector!="X"),
>   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz, 
> subset=sector!="A"),
>   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz, 
> subset=sector=="A"),
>   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz, 
> subset=sector=="M")
>   )
> 
> sapply(modco,function(x) coef(summary(x)))

This works, but consider if 

lm(normskvop ~ I(nts^0.5)-1+sector, data = colo,weights=wtz),

is a better approach anyway.

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] ARCH LM test for univariant time series

2008-02-06 Thread Pfaff, Bernhard Dr.
Hello Spencer,

splendid. Please go ahead. I am wondering if one should return the lm-object 
too and not only the htest-object. The benefit would be, that 
summary(lm-object) would return the mentioned F-test in the R-Help thread too, 
or one can return just the F-test result as a separate list element. If so, a 
more appropriate function name would be archtest().

What do you think?

Best,
Bernhard

Dr. Bernhard Pfaff
International Structured Products Group
Director

Invesco Asset Management Deutschland GmbH
Bleichstrasse 60-62
D-60313 Frankfurt am Main

Tel: +49(0)69 29807 230
Fax: +49(0)69 29807 178
Email: [EMAIL PROTECTED]

Geschäftsführer: Karl Georg Bayer, Bernhard Langer, Dr. Jens Langewand, 
Alexander Lehmann, Christian Puschmann
Handelsregister: Frankfurt am Main, HRB 28469
Sitz der Gesellschaft: Frankfurt am Main
  

>-Ursprüngliche Nachricht-
>Von: Spencer Graves [mailto:[EMAIL PROTECTED] 
>Gesendet: Mittwoch, 6. Februar 2008 05:02
>An: Pfaff, Bernhard Dr.
>Cc: tom soyer; r-help@r-project.org
>Betreff: Re: AW: [R] ARCH LM test for univariant time series
>
>Dear Bernhard: 
>
>  Thanks very much.  Unless you object, I shall add it to the 
>'FinTS' library as "ArchTest" (comparable to the S-PLUS Finmetrics 
>'archTest' function) -- with a worked example in '\scripts\ch03.R'. 
>
>  Best Wishes,
>  Spencer
>
>Pfaff, Bernhard Dr. wrote:
>> Dear All,
>>
>>
>> one can visually inspect ARCH-effects by plotting acf/pacf of the
>> squared residuals from an OLS-estimation. This can be as simple as a
>> demeaned series. Further one can run an auxiliary regression by
>> regressing q lagged squared values and a constant on the 
>squared series
>> itself. This test statistic (N-q)*R^2 is distributed as chisq with q
>> degrees of freedom.  
>>
>> Something along the lines:
>>
>> archlmtest <- function (x, lags, demean = FALSE) 
>> {
>>   x <- as.vector(x)
>>   if(demean) x <- scale(x, center = TRUE, scale = FALSE)
>> lags <- lags + 1
>> mat <- embed(x^2, lags)
>> arch.lm <- summary(lm(mat[, 1] ~ mat[, -1]))
>> STATISTIC <- arch.lm$r.squared * length(resid(arch.lm))
>> names(STATISTIC) <- "Chi-squared"
>> PARAMETER <- lags - 1
>> names(PARAMETER) <- "df"
>> PVAL <- 1 - pchisq(STATISTIC, df = PARAMETER)
>> METHOD <- "ARCH LM-test"
>> result <- list(statistic = STATISTIC, parameter = PARAMETER, 
>> p.value = PVAL, method = METHOD, data.name =
>> deparse(substitute(x)))
>> class(result) <- "htest"
>> return(result)
>> }
>>
>> should work and yield equal results as mentioned earlier in 
>this thread.
>>
>> Best,
>> Bernhard
>>
>>
>>   
>>> Spencer,
>>>
>>> The warning message is sent from VAR, it basically lets you 
>>> know that the
>>> data it used had no column names and it had to supply them 
>>> using y1, y2, y3,
>>> etc. It can be suppressed by including options(warn=-1) in the 
>>> function.
>>>
>>> Anyway, it seems that the p value from my function does not match
>>> FinMetrics'. I guess the function doesn't work... hmm...
>>>
>>>
>>> On 2/2/08, Spencer Graves <[EMAIL PROTECTED]> wrote:
>>> 
 Dear Tom:

  Your revised function eliminates the discrepancy in the 
   
>>> degrees of
>>> 
 freedom but is still very different from the numbers reports 
   
>>> on Tsay, p.
>>> 
 102:

 archTest(log(1+as.numeric(m.intc7303)), lag=12)

ARCH test (univariate)

 data:  Residual of y1 equation
 Chi-squared = 13.1483, df = 12, p-value = 0.3584

 Warning message:
 In VAR(s, p = 1, type = "const") :
 No column names supplied in y, using: y1, y2, y3, y4, y5, 
>y6, y7, y8,
 y9, y10, y11, y12 , instead.


  TOM:  What can you tell me about the warning message?

  Thanks for your help with this.
  Spencer Graves

 tom soyer wrote:
   
> Spencer,
>
> Sorry, I forgot that the default lag in arch is 16. Here 
> 
>>> is the fix. Can
>>> 
 you
   
> try it again and see if it gives the correct (or at least similar
> 
 compared
   
> to a true LM test) result?
>
> archTest=function(x, lags=12){
>  #x is a vector
>  require(vars)
>  s=embed(x,lags)
>  y=VAR(s,p=1,type="const")
>  result=arch(y,lags.single=lags,multi=F)$arch.uni[[1]]
>  return(result)
> }
>
> Thanks and sorry about the bug.
>
>
> On 2/2/08, Spencer Graves <[EMAIL PROTECTED]> wrote:
>
> 
>> Dear Tom, Bernhard, Ruey:
>>
>>  I can't get that to match Tsay's example, but I have other
>> questions about that.
>>
>>  1.  I got the following using Tom's 'archTest' 
>>   
>>> function (below):
>>> 
>>   
>>> archTest(log(1+as.numeric(m.intc7303)), lags=12)
>>>
>>> 
>>ARCH test (univariate)
>>
>> data:  Residu

Re: [R] maps and lattice

2008-02-06 Thread Paul Hiemstra
Hi,

I use the spplot function from the sp package, it uses lattice to make 
spatial plots. The sp-package provides spatial classes for R. See 
http://r-spatial.sourceforge.net/ for more information on the use of 
spatial data in R, http://r-spatial.sourceforge.net/gallery/ gives a 
gallery of spatial plots with the R code it was created by. For 
questions on spatial data in R the r-sig-geo mailing list is also a good 
option (https://stat.ethz.ch/mailman/listinfo/r-sig-geo).

hope this helps!

Paul

Jon Loehrke wrote:
> Is it possible to place maps onto lattice plots?
>
> With basic plotting you can add a map to a plot
>
> library(lattice)
> long<-c(-69.2, -69.5, -70.1, -70.3)
> lat<-c(41, 41.5, 43.2, 42.8)
> plot(long, lat)
> map('state', c("massachusetts"),add=TRUE)
>
> but is it possible with lattice?
>
>
> library(lattice)
> factor<-c(1,1,2,2)
> xyplot(lat~long|fact)
> ...now what?
>
> I have looked at panel and found through google an extravagant shape  
> file export/import to R.
> Is there a simpler fix?
>
> Thank you very much for your help.
>
> Jon
>
> School for Marine Science and Technology
> UMASS-Dartmouth
>
>
>   [[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.
>   


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +31302535773
Fax:+31302531145
http://intamap.geo.uu.nl/~paul

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 generate table output of a series of lm's

2008-02-06 Thread Martin Elff
On Wednesday 06 February 2008 (09:34:53), Dieter Menne wrote:
> Martin Elff  sowi.uni-mannheim.de> writes:
> > modco <- list(
> > lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz),
> > lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz,
> > subset=sector!="X"), lm(normskvop ~ I(nts^0.5)-1, data =
> > colo,weights=wtz, subset=sector!="A"), lm(normskvop ~ I(nts^0.5)-1, data
> > = colo,weights=wtz, subset=sector=="A"), lm(normskvop ~ I(nts^0.5)-1,
> > data = colo,weights=wtz, subset=sector=="M") )
> >
> > sapply(modco,function(x) coef(summary(x)))
>
> This works, but consider if
>
> lm(normskvop ~ I(nts^0.5)-1+sector, data = colo,weights=wtz),
>
> is a better approach anyway.

Hmm. Maybe if Daniel intended to do:

# Note now it is sector=="A" instead of sector!="A"

 modco <- list( 
   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz,
   subset=sector=="X"),
   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz,
   subset=sector=="A"),
   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz,
   subset=sector=="A"),
   lm(normskvop ~ I(nts^0.5)-1, data = colo,weights=wtz,
   subset=sector=="M")
   )

But in that case 
lm(normskvop ~ sector/I(nts^0.5), data = colo,weights=wtz)

or

lm(normskvop ~ sector/I(nts^0.5) - sector, data = colo,weights=wtz)

would be the "better" approach.

Martin

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

2008-02-06 Thread Martin Elff
On Wednesday 06 February 2008 (01:35:15), [EMAIL PROTECTED] wrote:
> how to generate correlated data which is correlated in three variables??

# Your correlation matrix
S <- rbind(
  c(1,.3,.3),
  c(.3,1,.3),
  c(.3,.3,1)
  )

# Three independent normal variates
x1 <- rnorm(1000)
x2 <- rnorm(1000)
x3 <- rnorm(1000)

# Using the cholesky decomposition of S
Y <- cbind(x1,x2,x3)%*%chol(S)

# Three correlated normal variates
y1 <- Y[,1]
y2 <- Y[,2]
y3 <- Y[,3]

# Check
cor(Y)

There may be more elegant and general solutions and
you can also use package "mvtnorm" to get correlated
normal (or t) variates. But the principle comes down
to this.

Hope that helps,

Martin

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

2008-02-06 Thread Senthil Kumar M
Hi,

I have the following data:

> Myvalues
 Gene  ES   MEF Embryo  ESHyp
1   GeneA   -0.38509507  0.00 1.6250  1.7039921
2   GeneB0.06262914  0.00 1.6250 -0.272033
and so on...

I want to plot the expression values of GeneA and GeneB in the
different cell/embryo/conditions (columns 2:5 above). Now, if I do:

>library(ggplot2)
> qplot(x=Gene, Embryo, geom = "bar")

I get a plot of GeneA, B...so on only for the Embryo (ie column 4).

How do I get to plot multiple instances of Y for the same value of X ?

I have tried to find this out myself, but it is a bit confusing and
so, I am writing to the list as a last resort.

Any help or guidance will be much appreciated.

TIA,

Senthil

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


Re: [R] inserting text lines in a dat frame

2008-02-06 Thread Richard . Cotton
>  I am trying to prepare a bed file to load as accustom track on the 
> UCSC genome browser.
> I have a data frame that looks like the one below.
> > x
>  V1V2 V3
> 1 chr1 11255 55
> 2 chr1 11320 29
> 3 chr1 11400 45
> 4 chr2 21680 35
> 5 chr2 21750 84
> 6 chr2 21820 29
> 7 chr2 31890 46
> 8 chr3 32100 29
> 9 chr3 52380 29
> 10 chr3 66450 46
> I would like to insert the following 4 lines at the beginning:
> browser position chr1:1-1
> browser hide all
> track type=wiggle_0 name=sample description=chr1_sample visibility=full
> variableStep chrom=chr1 span=1
> and then insert 2 lines before each chromosome:
> track type=wiggle_0 name=sample description=chr2_sample visibility=full
> vriableStep chrom=chr2 span=1
> The final result should be tab delimited file that looks like this:
> browser position chr1:1-1
> browser hide all
> track type=wiggle_0 name=sample description=chr1_sample visibility=full
> variableStep chrom=chr1 span=1
> chr1 11255 55
> chr1 11320 29
> chr1 11400 45
> track type=wiggle_0 name=sample description=chr2_sample visibility=full
> variableStep chrom=chr2 span=1
> chr2 21680 35
> chr2 21750 84
> chr2 21820 29
> track type=wiggle_0 name=sample description=chr3_sample visibility=full
> variableStep chrom=chr3 span=1
> chr3 32100 29
> chr3 32170 45
> chr3 32240 45

#First write your preamble text to a file
preamble <- c("browser position chr1:1-1",
  "browser hide all",
  "track type=wiggle_0 name=sample description=chr1_sample 
visibility=full",
  "variableStep chrom=chr1 span=1")
write(preamble, "myfile.txt")

#Now split your data frame up by the values in V1
x <- data.frame(V1=rep(c("chr1", "chr2", "chr3"),times=c(3,4,3)), 
V2=c(11255,11320,11400,21680,21750,21820,21890,32100,52380,66450),V3=c(55,29,45,35,84,29,46,29,29,46))
spx <- split(x, x$V1)

#Create lines of text to go before each piece of data frame
lV1 <- levels(x$V1)
maintext <- paste("track type=wiggle_0 name=sample description=", lV1,
   "_sample visibility=full\nvariableStep chrom=", lV1,
   " span=1", sep="")

#Use a loop to write the pieces to the file
for(i in 1: nlevels(x$V1))
  {
write(maintext[i], "myfile.txt", append=TRUE)
write.table(spx[[i]], "myfile.txt", append=TRUE, sep="\t", 
quote=FALSE,
col.names=FALSE, row.names=FALSE)
  }

Regards,
Richie.

Mathematical Sciences Unit
HSL



ATTENTION:

This message contains privileged and confidential inform...{{dropped:20}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Incomplete ouput with sink and split=TRUE

2008-02-06 Thread Alex Brown
you could use the unix function 'script' before invoking the R  
interpreter.

example session:

$ script
Script started, file is typescript
[x86_64|[EMAIL PROTECTED]:~]
$ R --quiet --vanilla
 > 1:10
[1]  1  2  3  4  5  6  7  8  9 10
 > q()
[x86_64|[EMAIL PROTECTED]:~]
$ exit
exit
Script done, file is typescript


contents of file typescript:

===
Script started on Tue Feb  5 19:01:32 2008
[x86_64|[EMAIL PROTECTED]:~]
$ R --quiet --vanilla
 > 1:10
[1]  1  2  3  4  5  6  7  8  9 10
 > q()
[x86_64|[EMAIL PROTECTED]:~]
$ exit
exit

Script done on Tue Feb  5 19:01:45 2008
[x86_64|[EMAIL PROTECTED]:~]


-Alex

On 5 Feb 2008, at 16:12, jiho wrote:

> Dear List,
>
> I am trying to get R's terminal output to a file and to the terminal
> at the same time, so that I can walk through some tests and keep a log
> concurrently. The function 'sink' with the option split=TRUE seems to
> do just that. It works fine for most output but for objects of class
> htest, the terminal output is incomplete (the lines are there but
> empty). Here is an example session which shows the problem:
>
>> sink("textout.txt", type="output", split=T)
>> b=bartlett.test(runif(10),c(1,1,1,1,2,2,2,2,2,2))
>> class(b)
> [1] "htest"
>> b
>
>
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
>
>> t=t.test(runif(10),c(1,1,1,1,2,2,2,2,2,2))
>> t
>
>
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
> -1.5807338 -0.7316803
> sample estimates:
> mean of x mean of y
> 0.4437929 1.600
>
>> sink()   # output in the file is complete
>> b
>
>   Bartlett test of homogeneity of variances
>
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> Bartlett's K-squared = 0.9959, df = 1, p-value = 0.3183
>
>> t
>
>   Welch Two Sample t-test
>
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> t = -5.7659, df = 16.267, p-value = 2.712e-05
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
> -1.5807338 -0.7316803
> sample estimates:
> mean of x mean of y
> 0.4437929 1.600
>
>>
>
> Is this a known bug (I'm using R 2.6.1 on OS X and Linux - FC8)? Is
> there an inherent reason why some portions of this output are not
> redirected?
>
> Thank you in advance for your help.
>
> JiHO
> ---
> http://jo.irisson.free.fr/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Histogram/Bar plot graph

2008-02-06 Thread ONKELINX, Thierry
You'll need to transform your dataset in a long format first.

library(ggplot2)
n <- 5
MyValues <- data.frame(Gene = factor(LETTERS[seq_len(n)]), ES =
rnorm(n), MEF = rnorm(n), Embrio = rnorm(n), EShyp = rnorm(n))
MyValuesMelt <- melt(MyValues, id.var = "Gene")
ggplot(MyValuesMelt, aes(x = Gene, y = value, fill = variable)) +
geom_bar(position = "dodge") 
ggplot(MyValuesMelt, aes(x = Gene, y = value)) + geom_bar(position =
"dodge") + facet_grid(. ~ variable)

HTH,

Thierry


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

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney

-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens Senthil Kumar M
Verzonden: woensdag 6 februari 2008 11:19
Aan: r-help@r-project.org
Onderwerp: [R] Histogram/Bar plot graph

Hi,

I have the following data:

> Myvalues
 Gene  ES   MEF Embryo  ESHyp
1   GeneA   -0.38509507  0.00 1.6250  1.7039921
2   GeneB0.06262914  0.00 1.6250 -0.272033
and so on...

I want to plot the expression values of GeneA and GeneB in the different
cell/embryo/conditions (columns 2:5 above). Now, if I do:

>library(ggplot2)
> qplot(x=Gene, Embryo, geom = "bar")

I get a plot of GeneA, B...so on only for the Embryo (ie column 4).

How do I get to plot multiple instances of Y for the same value of X ?

I have tried to find this out myself, but it is a bit confusing and so,
I am writing to the list as a last resort.

Any help or guidance will be much appreciated.

TIA,

Senthil

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

2008-02-06 Thread Senthil Kumar M
On Feb 6, 2008 2:35 AM, ONKELINX, Thierry <[EMAIL PROTECTED]> wrote:
> You'll need to transform your dataset in a long format first.
>
> library(ggplot2)
> n <- 5
> MyValues <- data.frame(Gene = factor(LETTERS[seq_len(n)]), ES =
> rnorm(n), MEF = rnorm(n), Embrio = rnorm(n), EShyp = rnorm(n))
> MyValuesMelt <- melt(MyValues, id.var = "Gene")
> ggplot(MyValuesMelt, aes(x = Gene, y = value, fill = variable)) +
> geom_bar(position = "dodge")
> ggplot(MyValuesMelt, aes(x = Gene, y = value)) + geom_bar(position =
> "dodge") + facet_grid(. ~ variable)
>

Hi Thierry,
Splendid! It is exactly what I wanted.
Now, I am actually studying your reply to understand what those
commands *actually* do.

Thanks a lot,
Sincerely Yours,
Senthil

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


Re: [R] inserting text lines in a dat frame

2008-02-06 Thread jim holtman
Try this and see if it is what you want:

x <- read.table(textConnection(" V1V2 V3
1 chr1 11255 55
2 chr1 11320 29
3 chr1 11400 45
4 chr2 21680 35
5 chr2 21750 84
6 chr2 21820 29
7 chr2 31890 46
8 chr3 32100 29
9 chr3 52380 29
10 chr3 66450 46" ), header=TRUE)
cat("browser position chr1:1-1\nrowser hide all\n", file='tempxx.txt')
result <- lapply(split(x, x$V1), function(.chro){
cat(sprintf("track type=wiggle_0 name=sample description=%s_sample
visibility=full\nvariableStep chrom=%s span=1\n",
as.character(.chro$V1[1]), as.character(.chro$V1[1])),
file="tempxx.txt", append=TRUE)
write.table(.chro, sep="\t", file="tempxx.txt", append=TRUE,
col.names=FALSE, row.names=FALSE)
})



On Feb 5, 2008 11:22 PM, joseph <[EMAIL PROTECTED]> wrote:
>
>
>
>
>
> Hi Jim
>  I am trying to prepare a bed file to load as accustom track on the UCSC
> genome browser.
> I have a data frame that looks like the one below.
> > x
>  V1V2 V3
> 1 chr1 11255 55
> 2 chr1 11320 29
> 3 chr1 11400 45
> 4 chr2 21680 35
> 5 chr2 21750 84
> 6 chr2 21820 29
> 7 chr2 31890 46
> 8 chr3 32100 29
> 9 chr3 52380
>  29
> 10 chr3 66450 46
> I would like to insert the following 4 lines at the beginning:
> browser position chr1:1-1
> browser hide all
> track type=wiggle_0 name=sample description=chr1_sample visibility=full
> variableStep chrom=chr1 span=1
> and then insert 2 lines before each chromosome:
> track type=wiggle_0 name=sample description=chr2_sample visibility=full
> vriableStep chrom=chr2 span=1
> The final result should be tab delimited file that looks like this:
> browser position chr1:1-1
> browser hide all
> track type=wiggle_0 name=sample description=chr1_sample visibility=full
> variableStep chrom=chr1 span=1
> chr1 11255 55
> chr1 11320 29
> chr1 11400 45
> track type=wiggle_0 name=sample description=chr2_sample visibility=full
> variableStep chrom=chr2 span=1
> chr2 21680 35
> chr2 21750 84
> chr2 21820 29
> track type=wiggle_0 name=sample description=chr3_sample visibility=full
> variableStep chrom=chr3
>  span=1
> chr3 32100 29
> chr3 32170 45
> chr3 32240 45
> Any kind of help or guidance will be much appreciated.
> Joseph
>
> 
> Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it
> now.



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

What is the problem 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] error message from apply()

2008-02-06 Thread jim holtman
Is 'thr' supposed to be the mean and sd of all the values in data2_1?
If so, then

thr <- mean(data2_1, na.rm=TRUE) + sd(data2_1,na.rm=TRUE)

I am not exactly sure of "what is the problem that you are trying to
solve".  You just have to make sure that the object you are creating
by precomputing has the right structure to do what you want.

On Feb 6, 2008 12:56 AM, Stanley Ng <[EMAIL PROTECTED]> wrote:
> Now I understand why 3 by 3 data2_1 works and not the 3x10 data2_1.
>
> How can I precompute thr and pass it safely to function(x) for the column
> operation ?
>
>
> -Original Message-
> From: jim holtman [mailto:[EMAIL PROTECTED]
> Sent: Wednesday, February 06, 2008 11:33
> To: Ng Stanley
> Cc: r-help
> Subject: Re: [R] error message from apply()
>
> You matrix only has 3 rows, so when you do 'apply(data2_1,2,...)' you are
> extracting columns which only have a length of 3 while thr has a length of
> 10
>
> > str(data2_1)
>  num [1:3, 1:10]  0.958  0.271 -0.950 -0.130 -0.754 ...
> > str(thr)
>  num [1:10]  1.060  0.528  0.104  0.925 -0.256 ...
> >
>  That is why you get the error message of a size mismatch.
>
> On Feb 5, 2008 10:21 PM, Ng Stanley <[EMAIL PROTECTED]> wrote:
> > Replacing colMeans by mean removed the warning messages. Thanks
> >
> > However, when I precompute thr, and pass it to function(x), the error
> > returns. Using the shorter data2_1, doesn't give any warnings. What is
> > happening ?
> >
> > data2_1 <- matrix(c(0.9584190, 0.2710325, -0.9495618, -0.1301772,
> > -0.7539687, 0.5344464, -0.8205933, 0.1581723, -0.5351588, 0.04448065,
> > 0.9936430, 0.2278786, -0.8160700, -0.3314779, -0.4047975, 0.1168152,
> > -0.7458182, - 0.2231588, -0.5051651, -0.74871174, 0.9450363,
> > 0.4797723, -0.9033313, - 0.5825065, 0.8523742, 0.7402795, -0.7134312,
> > -0.8162558, 0.6345438, - 0.05704138), 3,10) # data2_1 <-
> > matrix(c(0.9584190, 0.2710325, -0.9495618, -0.1301772, - 0.7539687,
> > 0.5344464, -0.8205933, 0.1581723, -0.5351588), 3,3)
> >
> > thr <- colMeans(data2_1, na.rm = TRUE) + sd(data2_1, na.rm = TRUE)
> >
> > num <- apply(data2_1, 2, function(x) {
> >sum(x > (thr), na.rm = TRUE)
> > })
> >
> >
> >
> > On 2/6/08, jim holtman <[EMAIL PROTECTED]> wrote:
> > >
> > > The error message was coming from the call to colMeans where 'x' was
> > > not a matrix; it was a vector that resulted from the 'apply' call.
> > > Did you intend to use 'mean' instead like this example:
> > >
> > > > data2_1 <- matrix(c(0.9584190, 0.2710325, -0.9495618, -0.1301772,
> > > > -
> > > 0.7539687,
> > > + 0.5344464, -0.8205933, 0.1581723, -0.5351588, 0.04448065,
> > > + 0.9936430, 0.2278786, -0.8160700, -0.3314779, -0.4047975,
> > > + 0.1168152, -0.7458182, - 0.2231588, -0.5051651, -0.74871174,
> > > + 0.9450363, 0.4797723, -0.9033313, - 0.5825065, 0.8523742,
> > > + 0.7402795, -0.7134312, -0.8162558, 0.6345438, - 0.05704138), 3,10)
> > > >
> > > > num <- apply(data2_1, 2, function(x) {sum(x > (mean(x, na.rm =
> > > > TRUE) +
> > > + 1*sd(x, na.rm = TRUE)), na.rm = TRUE)})
> > > > num
> > > [1] 0 1 1 1 0 0 1 1 0 0
> > > >
> > >
> > >
> > > On Feb 5, 2008 8:43 PM, Ng Stanley <[EMAIL PROTECTED]> wrote:
> > > > Hi,
> > > >
> > > > I keep getting the error message. Please help.
> > > >
> > > > Error in colMeans(x, na.rm = TRUE) :   'x' must be an array of at
> least
> > > two
> > > > dimensions
> > > >
> > > > The codes are:
> > > >
> > > > data2_1 <- matrix(c(0.9584190, 0.2710325, -0.9495618, -0.1301772,
> > > > -
> > > 0.7539687,
> > > > 0.5344464, -0.8205933, 0.1581723, -0.5351588, 0.04448065,
> > > > 0.9936430, 0.2278786, -0.8160700, -0.3314779, -0.4047975,
> > > > 0.1168152, -0.7458182, - 0.2231588, -0.5051651, -0.74871174,
> > > > 0.9450363, 0.4797723, -0.9033313, - 0.5825065, 0.8523742,
> > > > 0.7402795, -0.7134312, -0.8162558, 0.6345438, - 0.05704138), 3,10)
> > > >
> > > > num <- apply(data2_1, 2, function(x) {sum(x > (colMeans(x, na.rm =
> > > > TRUE)
> > > +
> > > > 1*sd(x, na.rm = TRUE)), na.rm = TRUE)})
> > > >
> > > >[[alternative HTML version deleted]]
> > > >
> > > > __
> > > > R-help@r-project.org mailing list
> > > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > > PLEASE do read the posting guide
> > > http://www.R-project.org/posting-guide.html
> > > > and provide commented, minimal, self-contained, reproducible code.
> > > >
> > >
> > >
> > >
> > > --
> > > Jim Holtman
> > > Cincinnati, OH
> > > +1 513 646 9390
> > >
> > > What is the problem you are trying to solve?
> > >
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
>
> --
> Jim Holtman
> Cincinnati, OH
> +1 513 646 9390
>
> What is the problem you are trying to solve?
>
> _

[R] Multivariate Maximum Likelihood Estimation

2008-02-06 Thread Konrad BLOCHER
Hi,

I am trying to perform Maximum Likelihood estimation of a Multivariate
model (2 independent variables + intercept) with autocorrelated errors of
1st order (ar(1)).

Does R have a function for that? I could only find an univariate option
(ar.mle function) and when writing my own I find that it is pretty
memory-consuming (and sometimes wrong) so there must be a better way.

Thanks,

KB

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 search for packages - wrap up!

2008-02-06 Thread Neil Shephard



Charilaos Skiadas-3 wrote:
> 
> On Feb 5, 2008, at 10:37 AM, Monica Pisica wrote:
> 
> But perhaps I am missing something very obvious?
> 
> 

I thought the task views were located where they are (linked from the page
that lists packages) as they summarise the available packages for the given
topic.

Neil
-- 
View this message in context: 
http://www.nabble.com/Re%3A-How-to-search-for-packages---wrap-up%21-tp15297545p15306481.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] Incomplete ouput with sink and split=TRUE

2008-02-06 Thread jiho

On 2008-February-06  , at 11:25 , Alex Brown wrote:
> you could use the unix function 'script' before invoking the R  
> interpreter.


Thanks for the suggestion. It would work for some cases and I did not  
know about this utility. But most of the time I just put a list of  
commands in a .R file and source it to execute it. In this case the  
transcript would only contain "source("whatever.R", print.eval=T)"  
which won't be very informative :/. And in the case of sourcing, the  
problem with disappearing text stays the same.

Thanks for you answer anyway.

JiHO
---
http://jo.irisson.free.fr/

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

2008-02-06 Thread Konrad BLOCHER
I get this message:

Error: could not find function "gls" (and also)
Error: could not find function "lm.gls"

Which package is that in?

Thanks,

KB

> Have you tried gls()?
>
>
> 
> 
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
> tel. + 32 54/436 185
> [EMAIL PROTECTED]
> www.inbo.be
>
> Do not put your faith in what statistics say until you have carefully
> considered what they do not say.  ~William W. Watt
> A statistical analysis, properly conducted, is a delicate dissection of
> uncertainties, a surgery of suppositions. ~M.J.Moroney
>
> -Oorspronkelijk bericht-
> Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> Namens Konrad BLOCHER
> Verzonden: woensdag 6 februari 2008 12:12
> Aan: r-help@r-project.org
> Onderwerp: [R] Multivariate Maximum Likelihood Estimation
>
> Hi,
>
> I am trying to perform Maximum Likelihood estimation of a Multivariate
> model (2 independent variables + intercept) with autocorrelated errors
> of
> 1st order (ar(1)).
>
> Does R have a function for that? I could only find an univariate option
> (ar.mle function) and when writing my own I find that it is pretty
> memory-consuming (and sometimes wrong) so there must be a better way.
>
> Thanks,
>
> KB
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 search for packages - wrap up!

2008-02-06 Thread Gavin Simpson
hits=-2.6 tests=BAYES_00
X-USF-Spam-Flag: NO

On Wed, 2008-02-06 at 03:23 -0800, Neil Shephard wrote:
> 
> 
> Charilaos Skiadas-3 wrote:
> > 
> > On Feb 5, 2008, at 10:37 AM, Monica Pisica wrote:
> > 
> > But perhaps I am missing something very obvious?
> > 
> > 
> 
> I thought the task views were located where they are (linked from the page
> that lists packages) as they summarise the available packages for the given
> topic.
> 
> Neil

Originally, that might have been the thinking behind locating the link
on CRAN and not the R homepage. Another reason might have been to keep
the number of menu options on the R homepage to a manageable number.

A new user will come to the R homepage, go to CRAN via the link under
download and from there see Packages and then be swamped by the huge
number available. Having Task Views as a link on the R homepage would
make these more visible.

However, getting this working so that users don't swamp the CRAN master
will be interesting. The Task Views themselves could be on r-project.org
but looking at the packages linked from the Task Views should be done on
one of the CRAN mirrors and how that is arranged could be more hassle
than it is worth, just to give Task Views higher visibility on R
homepage.

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

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


Re: [R] How to search for packages - wrap up!

2008-02-06 Thread Charilaos Skiadas

On Feb 6, 2008, at 6:23 AM, Neil Shephard wrote:

> Charilaos Skiadas-3 wrote:
>>
>> On Feb 5, 2008, at 10:37 AM, Monica Pisica wrote:
>>
>> But perhaps I am missing something very obvious?
>
> I thought the task views were located where they are (linked from  
> the page
> that lists packages) as they summarise the available packages for  
> the given
> topic.

That would make sense, and brings up the other point, that the  
package directory is nontrivial to find as well. I'm thinking here of  
someone who has not used R yet, but is considering it. One approach  
those people would take, it seems to me, is to go to the main page  
and look around at the links to see if there is any information that  
would help them decide if R is right for what they want it. They have  
not heard of packages, and have no idea that a lot of the  
functionality is in the packages.

So they would look on the list on the left looking for something that  
"clicks", and none of the items there would. If they eventually  
decide to click on the CRAN link,  they are now faced with a long  
list of mirrors, and no further explanation of what is behind that.  
If they have seen mirror systems elsewhere, they would immediately  
come to the conclusion that this page will simply lead to downloading  
the software, and dismiss it as not useful.


Just reading Gavin's reply, and he makes a good point about the  
difficulties related to the CRAN master. But couldn't we simply have  
the links to packages from the task views send someone through the  
mirror list? I would imagine something like that should be doable.

Or we could simply have the "Task Views" link on the main page send  
you to a mirrors list, with perhaps a one-line explanation on the top  
of why this part is necessary. But at least at that point the user  
has selected "Task Views" and is more certain that they are on the  
right track. Even better, this step could be even one step deeper,  
after the user has selected which task view they want to see.

> Neil

Haris Skiadas
Department of Mathematics and Computer Science
Hanover College

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Effect size of comparison of two levels of a factor in multiple linear regression

2008-02-06 Thread Christoph Mathys
On Feb 4, 08:49 AM, Chuck Cleland wrote:
> On 2/3/2008 10:09 AM, Christoph Mathys wrote:
>> Dear R users,
>> I have a linear model of the kind
>> outcome ~ treatment + covariate
>> where 'treatment' is a factor with three levels ("0", "1", and "2"),
>> and the covariate is continuous. Treatments "1" and "2" both have
>> regression coefficients significantly different from 0 when using
>> treatment contrasts with treatment "0" as the baseline. I would now like
>> to determine effect sizes (akin to Cohen's d in a two-sample comparison)
>> for the comparison to baseline of treatments "1" and "2". I have
>> illustrated a way to do this in the reproducible example below and am
>> grateful for any comments on the soundness of what I'm doing. I have not
>> yet found a way to determine confidence intervals for the effect sizes
>> derived below and would appreciate tips on that.
>
>   How about scaling the outcome by the residual standard error from the 
> unstandardized model?  For example:

Yes, that also works. It's actually what my simulation does implicitly
and what I should have done in the first place. It has the added benefit
that a confidence interval can be given. Thanks.

>
> library(MASS)
>
> cmat <- diag(3)
> diag(cmat) <- c(25,25,25)
>
> df <- data.frame(Y = c(mvrnorm(100, mu=c(10,30,40), Sigma=cmat, 
> empirical=TRUE)), TX = factor(rep(c(0,1,2), each=100)))
>
> fm1 <- lm(Y ~ TX, data = df)
> fm2 <- lm(scale(Y, scale=summary(fm1)$sigma) ~ TX, data = df)
>
> summary(fm1)
>
> Call:
> lm(formula = Y ~ TX, data = df)
>
> Residuals:
>  Min   1Q   Median   3Q  Max
> -12.7260  -3.5215  -0.1982   3.4517  12.1690
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept)  10. 0.5000   20.00   <2e-16 ***
> TX1  20. 0.7071   28.28   <2e-16 ***
> TX2  30. 0.7071   42.43   <2e-16 ***
> ---
> Signif. codes:  0 ???***??? 0.001 ???**??? 0.01 ???*??? 0.05 ???.??? 0.1 
> ??? ??? 1
>
> Residual standard error: 5 on 297 degrees of freedom
> Multiple R-squared: 0.8627, Adjusted R-squared: 0.8618
> F-statistic: 933.3 on 2 and 297 DF,  p-value: < 2.2e-16
>
> summary(fm2)
>
> Call:
> lm(formula = scale(Y, scale = summary(fm1)$sigma) ~ TX, data = df)
>
> Residuals:
>  Min   1Q   Median   3Q  Max
> -2.54521 -0.70431 -0.03965  0.69034  2.43380
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept)  -3. 0.1000  -33.33   <2e-16 ***
> TX1   4. 0.1414   28.28   <2e-16 ***
> TX2   6. 0.1414   42.43   <2e-16 ***
> ---
> Signif. codes:  0 ???***??? 0.001 ???**??? 0.01 ???*??? 0.05 ???.??? 0.1 
> ??? ??? 1
>
> Residual standard error: 1 on 297 degrees of freedom
> Multiple R-squared: 0.8627, Adjusted R-squared: 0.8618
> F-statistic: 933.3 on 2 and 297 DF,  p-value: < 2.2e-16
>
> confint(fm2)
> 2.5 %97.5 %
> (Intercept) -3.530132 -3.136535
> TX1  3.721685  4.278315
> TX2  5.721685  6.278315
>
>   I've never seen this approach before, and I'm not how it would work when 
> the variances within groups are heterogeneous or one or more covariates are 
> added to the model.

Strictly speaking, the approach breaks down when we don't have
homogeneity of variances across groups. But luckily, I have it in the
real data I'm looking at.

When covariates are added, one has to decide case by case what to do
about them, namely which model's sigma to standardize against.
Standardizing against the sigma of a model that includes a given
covariate implies looking at effect sizes in a population where that
covariate has a fixed value. On the other hand, standardizing against the
sigma of a model that does not include a given covariate implies looking
at effect sizes in a population where that covariate varies.

> hope this helps,

Sure did. Thanks again,

Christoph

--
Christoph Mathys, M.S.
Music and Neuroimaging Laboratory
Beth Israel Deaconess Medical Center
and Harvard Medical School

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

2008-02-06 Thread Gavin Simpson
hits=-2.6 tests=BAYES_00
X-USF-Spam-Flag: NO

On Wed, 2008-02-06 at 12:45 +0100, Konrad BLOCHER wrote:
> I get this message:
> 
> Error: could not find function "gls" (and also)
> Error: could not find function "lm.gls"
> 
> Which package is that in?

RSiteSearch("gls", restrict = "functions")

Tells you the answer.

Now might be a good time to review R's search/help tools - details of
which are in a section in the Posting Guide. You'll save yourself a lot
of time in the long run if you do.

HTH

G

> 
> Thanks,
> 
> KB
> 
> > Have you tried gls()?
> >
> >
> > 
> > 
> > ir. Thierry Onkelinx
> > Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> > and Forest
> > Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> > methodology and quality assurance
> > Gaverstraat 4
> > 9500 Geraardsbergen
> > Belgium
> > tel. + 32 54/436 185
> > [EMAIL PROTECTED]
> > www.inbo.be
> >
> > Do not put your faith in what statistics say until you have carefully
> > considered what they do not say.  ~William W. Watt
> > A statistical analysis, properly conducted, is a delicate dissection of
> > uncertainties, a surgery of suppositions. ~M.J.Moroney
> >
> > -Oorspronkelijk bericht-
> > Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> > Namens Konrad BLOCHER
> > Verzonden: woensdag 6 februari 2008 12:12
> > Aan: r-help@r-project.org
> > Onderwerp: [R] Multivariate Maximum Likelihood Estimation
> >
> > Hi,
> >
> > I am trying to perform Maximum Likelihood estimation of a Multivariate
> > model (2 independent variables + intercept) with autocorrelated errors
> > of
> > 1st order (ar(1)).
> >
> > Does R have a function for that? I could only find an univariate option
> > (ar.mle function) and when writing my own I find that it is pretty
> > memory-consuming (and sometimes wrong) so there must be a better way.
> >
> > Thanks,
> >
> > KB
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/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.
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

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


Re: [R] advice requested re: building "good" system (R, SQL db) for handling large datasets

2008-02-06 Thread Richard Pearson
Hi Thomas

I'm certainly no expert but thought I'd reply as I'm likely to be in a 
similar position soon.

With regards versions of R I think you should always have the latest 
release version. This will mean upgrading at least every 6 months, but 
this shouldn't be too much of a problem.

With OSs, you need to be aware that there is an upper limit to the 
amount of RAM than be handled (2-4 GB) with many. I think if you plan to 
use more than 4GB RAM, you should definitely consider 64-bit linux. I 
have no information or opinions as to which flavour of linux.

With databases, one issue that might be relevant is whether you want to 
store data in tables (e.g. one table to store one data.frame) that can 
subsequently be manipulated in the DB, or to store R objects as R 
objects (e.g. as BLOBs). My situation is likely to be the later case, 
and one of my concerns is that many DBs have an upper limit of 2GB on 
BLOBs, and I might potentially have objects that are larger than this.

Finally, you might get more response on database issues from R-sig-db 
than R-help.

Best wishes

Richard.


Thomas Pujol wrote:
> R-community,
> Sometime during the next 12-months, I plan on configuring a new computer 
> system on which I will primarily run "R" and a SQL database (Microsoft SQL 
> Server, MySQL, Oracle, etc).  My primary goal is to "optimize" the system for 
> R, and for passing data to and from R and the database.
>
> I work with large datasets, and therefore I "think" one of my most important 
> goals should be to maximize the amount of RAM that R can utilize effectively.
>
> I am seeking advice concerning the version of R, OS,  processor, 
> hard-drive/storage configuration, database, etc. that I should consider. (I 
> am guessing that I should build a system with lots of RAM, and a Linux OS, 
> but am seeking advice from the R community.) If I choose Linux, does it 
> matter which version I use? Any opinion regarding  implementing a 
> commercially supported version from a vendor such as Red Hat, Sun, etc? Is 
> any database particularly better at "exchanging" data with R?
>
> While cost is of course a consideration, it is probably a secondary 
> consideration to overall performance, reliability, and ease of ongoing 
> maintenance/support.
>
> Thanks!
>
>
> -
>
>   [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 search for packages - wrap up!

2008-02-06 Thread Neil Shephard
I'm not sure that this would make any difference to someone considering using R.

Would they know what CRAN stands for? Probably not unless they've used
CPAN or equivalent in the past.

Would they know what a 'Task View' is? Again probably not as its not
patently obvious what it is, it doesn't "click" (at least not to my
mind if I take a step back from already knowing what a Task View is).

I think the current set up with mirrors is fine, it only takes a few
seconds to click on a mirror and find that it takes you to a page that
isn't simply for downloading software.

Discussions of how the R web-site can be improved/altered seem to crop
up periodically on the list.  At the end of the day someone (who?) has
to make the changes to the site.  R developers already sacrifice their
own time to improving the software.  If you want to see changes, then
volunteer to (I presume) the R-core team and enter into dialogue with
them.

Neil

On Feb 6, 2008 12:02 PM, Charilaos Skiadas <[EMAIL PROTECTED]> wrote:
>
> On Feb 6, 2008, at 6:23 AM, Neil Shephard wrote:
>
> > Charilaos Skiadas-3 wrote:
> >>
> >> On Feb 5, 2008, at 10:37 AM, Monica Pisica wrote:
> >>
> >> But perhaps I am missing something very obvious?
> >
> > I thought the task views were located where they are (linked from
> > the page
> > that lists packages) as they summarise the available packages for
> > the given
> > topic.
>
> That would make sense, and brings up the other point, that the
> package directory is nontrivial to find as well. I'm thinking here of
> someone who has not used R yet, but is considering it. One approach
> those people would take, it seems to me, is to go to the main page
> and look around at the links to see if there is any information that
> would help them decide if R is right for what they want it. They have
> not heard of packages, and have no idea that a lot of the
> functionality is in the packages.
>
> So they would look on the list on the left looking for something that
> "clicks", and none of the items there would. If they eventually
> decide to click on the CRAN link,  they are now faced with a long
> list of mirrors, and no further explanation of what is behind that.
> If they have seen mirror systems elsewhere, they would immediately
> come to the conclusion that this page will simply lead to downloading
> the software, and dismiss it as not useful.
>
>
> Just reading Gavin's reply, and he makes a good point about the
> difficulties related to the CRAN master. But couldn't we simply have
> the links to packages from the task views send someone through the
> mirror list? I would imagine something like that should be doable.
>
> Or we could simply have the "Task Views" link on the main page send
> you to a mirrors list, with perhaps a one-line explanation on the top
> of why this part is necessary. But at least at that point the user
> has selected "Task Views" and is more certain that they are on the
> right track. Even better, this step could be even one step deeper,
> after the user has selected which task view they want to see.
>
> > Neil
>
> Haris Skiadas
> Department of Mathematics and Computer Science
> Hanover College
>
>
>
>
>



-- 
"Do you really need to send me the email I just sent to you?"  - Me

Email - [EMAIL PROTECTED] / [EMAIL PROTECTED]
Website - http://slack.ser.man.ac.uk/
Photos - http://www.flickr.com/photos/slackline/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] advice requested re: building "good" system (R, SQL db) for handling large datasets

2008-02-06 Thread Rainer M Krug
With databases, one issue that might be relevant is whether you want to
> store data in tables (e.g. one table to store one data.frame) that can
> subsequently be manipulated in the DB, or to store R objects as R
> objects (e.g. as BLOBs). My situation is likely to be the later case,
> and one of my concerns is that many DBs have an upper limit of 2GB on
> BLOBs, and I might potentially have objects that are larger than this.


R objects in blobs - I never thought about that. Could you elaborate on how
to do something like that (I am using RMySQL)?

Thanks

Rainer


Finally, you might get more response on database issues from R-sig-db
> than R-help.
>
> Best wishes
>
> Richard.
>
>
> Thomas Pujol wrote:
> > R-community,
> > Sometime during the next 12-months, I plan on configuring a new computer
> system on which I will primarily run "R" and a SQL database (Microsoft SQL
> Server, MySQL, Oracle, etc).  My primary goal is to "optimize" the system
> for R, and for passing data to and from R and the database.
> >
> > I work with large datasets, and therefore I "think" one of my most
> important goals should be to maximize the amount of RAM that R can utilize
> effectively.
> >
> > I am seeking advice concerning the version of R, OS,  processor,
> hard-drive/storage configuration, database, etc. that I should consider. (I
> am guessing that I should build a system with lots of RAM, and a Linux OS,
> but am seeking advice from the R community.) If I choose Linux, does it
> matter which version I use? Any opinion regarding  implementing a
> commercially supported version from a vendor such as Red Hat, Sun, etc? Is
> any database particularly better at "exchanging" data with R?
> >
> > While cost is of course a consideration, it is probably a secondary
> consideration to overall performance, reliability, and ease of ongoing
> maintenance/support.
> >
> > Thanks!
> >
> >
> > -
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 

-- 
Rainer M. Krug, Dipl. Phys. (Germany), MSc Conservation Biology (UCT)

Plant Conservation Unit Department of Botany
University of Cape Town
Rondebosch 7701
South Africa

[[alternative HTML version deleted]]

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


Re: [R] How to form a simple R package

2008-02-06 Thread Doran, Harold
Well, not exactly. package.skeleton() is very useful as a first step,
but it does *not* create a package entirely. It turns out that in
Windows, creating a package is very simple once you have downloaded all
programs needed (e.g., perl) and you have your path configured exaclty,
(exactly, exactly) as described in the document found at the link below.

After running package.skeleton, you must still run Rcmd build, Rcmd
check, and Rcmd install.

http://cran.r-project.org/doc/manuals/R-admin.html#The-Windows-toolset

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Benilton Carvalho
> Sent: Wednesday, February 06, 2008 2:33 AM
> To: [EMAIL PROTECTED]
> Cc: r-help@r-project.org
> Subject: Re: [R] How to form a simple R package
> 
> ?package.skeleton
> 
> b
> 
> On Feb 5, 2008, at 1:48 PM, [EMAIL PROTECTED] wrote:
> 
> > Is there a function to form in one step (configure files 
> and install) 
> > a simple R package of consisting of one script file of R functions?
> >
> > For example in Windows:
> >
> > form.package(name="mypkg", rcodefile ="c:\misc\mypkg.r" )
> >
> > Thank you for any comments.
> > Giles
> 
> 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 search for packages - wrap up!

2008-02-06 Thread hadley wickham
> A new user will come to the R homepage, go to CRAN via the link under
> download and from there see Packages and then be swamped by the huge
> number available. Having Task Views as a link on the R homepage would
> make these more visible.

I would think that a new user would see the download heading, then
think - I don't want to download a cran (whatever that is) I want to
downloaded R, and then continue to be confused for another 10 minutes
until they ask their colleague down the hall for the sequence of 5
(cran -> mirror -> windows -> base -> R-2.6.1-win32.exe) clicks that
they need to find the installer!

It would be nice to automatically provide a link to the current
version of R for the platform that the user is browsing from (see e.g.
getfirefox.com).  Automatic selection of a mirror would be just as
desirable, although I still maintain it would be better to dump the
mirror system entirely and move to a content delivery network (CDN).

Hadley

-- 
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] Mixed models quantile regression

2008-02-06 Thread roger koenker
There is no generally agreed upon notion of random effects for quantile
regression applications.   Insofar as one is willing to accept the  
idea that
random effects are just "shrunken fixed effects" one can consider  
similar
schemes in the QR context; one such is described in

“Quantile Regression for Longitudinal Data,” J. of Multivariate  
Analysis,
(2004),  91, 74-89.

But there are many open questions,  so this is not for the faint hearted
and there is certainly no "official" code.


url:www.econ.uiuc.edu/~rogerRoger Koenker
email[EMAIL PROTECTED]Department of Economics
vox: 217-333-4558University of Illinois
fax:   217-244-6678Champaign, IL 61820


On Feb 5, 2008, at 9:40 PM, Beth Strain wrote:

>
> Dear R,
> I have a question concerning quantile regression models.
>
> I am using the quantile regression model to test the relationship  
> between
> abalone and the percentage cover of algae etc at different sites and  
> depths.
>
> An example is
> fit<-rq(abalone~%coversessileinvertebrates+factor(Depth) 
> +factor(Site),tau=0.7)
>
> In this model depth is fixed and site is random. My question is, is it
> possible specify the fixed and random effects in this model. If so  
> could
> someone please give me an example of how to write the code.
>
> I can;t seem to find any information in the R help.
> Thanks in advance
> Beth Strain
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Incomplete ouput with sink and split=TRUE

2008-02-06 Thread Duncan Murdoch
On 2/5/2008 11:12 AM, jiho wrote:
> Dear List,
> 
> I am trying to get R's terminal output to a file and to the terminal  
> at the same time, so that I can walk through some tests and keep a log  
> concurrently. The function 'sink' with the option split=TRUE seems to  
> do just that. It works fine for most output but for objects of class  
> htest, the terminal output is incomplete (the lines are there but  
> empty). Here is an example session which shows the problem:

stats:::print.htest() uses writeLines to write some of its output to 
stdout(), and it looks as though sink(split=T) misses those bits.

I'll change print.htest to use cat(), but it is probably a sign of a 
bigger problem in sink(), and it's too late in the schedule to touch 
that for 2.6.2.

Duncan Murdoch

> 
>  > sink("textout.txt", type="output", split=T)
>  > b=bartlett.test(runif(10),c(1,1,1,1,2,2,2,2,2,2))
>  > class(b)
> [1] "htest"
>  > b
> 
> 
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> 
>  > t=t.test(runif(10),c(1,1,1,1,2,2,2,2,2,2))
>  > t
> 
> 
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
>   -1.5807338 -0.7316803
> sample estimates:
> mean of x mean of y
> 0.4437929 1.600
> 
>  > sink() # output in the file is complete
>  > b
> 
>   Bartlett test of homogeneity of variances
> 
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> Bartlett's K-squared = 0.9959, df = 1, p-value = 0.3183
> 
>  > t
> 
>   Welch Two Sample t-test
> 
> data:  runif(10) and c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2)
> t = -5.7659, df = 16.267, p-value = 2.712e-05
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
>   -1.5807338 -0.7316803
> sample estimates:
> mean of x mean of y
> 0.4437929 1.600
> 
>  >
> 
> Is this a known bug (I'm using R 2.6.1 on OS X and Linux - FC8)? Is  
> there an inherent reason why some portions of this output are not  
> redirected?
> 
> Thank you in advance for your help.
> 
> JiHO
> ---
> http://jo.irisson.free.fr/
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] counting row repetitions without loop

2008-02-06 Thread Waterman, DG (David)
Hi,
 
I have a data frame consisting of coordinates on a 10*10 grid, i.e.
 
> example
x  y
1   4  5
2   6  7
3   6  6
4   7  5
5   5  7
6   6  7
7   4  5
8   6  7
9   7  6
10  5  6

What I would like to do is return an 10*10 matrix consisting of counts
at each position, so in the above example I would have a matrix where,
for example, cell [4,5] contains 2 and [6,7] contains 3. At the moment I
have implemented this using a for loop over the rows of the data frame,
however the data frames I want to process are very long so the loop
takes many minutes to complete. Can I do this in a more efficient way?
 
Cheers,
David
This e-mail and any attachments may contain 
confidential, copyright and or privileged material, and are for the use of the 
intended addressee only. If you are not the intended addressee or an authorised 
recipient of the addressee please notify us of receipt by returning the e-mail 
and do not use, copy, retain, distribute or disclose the information in or 
attached to the e-mail.
Any opinions expressed within this e-mail are those of the individual and not 
necessarily of Diamond Light Source Ltd. 
Diamond Light Source Ltd. cannot guarantee that this e-mail or any attachments 
are free from viruses and we cannot accept liability for any damage which you 
may sustain as a result of software viruses which may be transmitted in or with 
the message.
Diamond Light Source Limited (company no. 4375679). Registered in England and 
Wales with its registered office at Diamond House, Harwell Science and 
Innovation Campus, Didcot, Oxfordshire, OX11 0DE, United Kingdom
 

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

2008-02-06 Thread [EMAIL PROTECTED]


Hi, i'm having trouble with my x and y axis. The commands i'm using are
   below. The problem is that the y axis starts at coordinate 0,1 and the x
   axis starts at coordinate 0,0. As far as I know the y axis can't  start at 0
   (because it's log scaled) ,so I would like to position the x axis at 0,1 but
   don't know how to do this. Could anyone advise?

   Thanks

   Rich

   df3<-data.frame(x=c(10,11,115,12,13,14,16,17,18,21,22,23,24,26,27,28,29,3,30
   ,32,33,34,35,4,41,45,5,50,52,56,58,6,67,6738,68,7,8,9,),fq=c(8,11,1,2,4,4,2,
   2,6,3,4,2,2,1,1,1,4,51,3,1,1,1,1,35,1,1,19,2,1,1,1,14,1,1,1,10,13,5,),fqcvd=
   c(5,8,1,1,3,3,2,2,5,3,4,2,2,0,1,1,3,13,2,1,1,1,1,17,1,0,11,2,1,1,1,7,1,1,1,7
   ,7,1,),fqcan=c(1,1,0,2,1,1,1,0,3,0,2,0,1,0,1,0,1,4,2,1,1,0,0,4,1,1,2,2,0,1,0
   ,2,0,1,1,2,3,1,),fqnondis=c(8,11,1,2,4,4,2,2,6,3,4,2,2,1,1,1,4,50,3,1,1,1,1,
   34,1,1,19,2,1,1,1,14,1,1,1,10,12,5,))
   k3<-with(df3,rep(x,times=fq))
   kcvd3<-with(df3,rep(x,times=fqcvd))
   kcvd3<-c(kcvd3,rep(NA,times=length(k3)-length(kcvd3)))
   kcan3<-with(df3,rep(x,times=fqcan))
   kcan3<-c(kcan3,rep(NA,times=length(k3)-length(kcan3)))
   knondis3<-with(df3,rep(x,times=fqnondis))
   knondis3<-c(knondis3,rep(NA,times=length(k3)-length(knondis3)))

   boxplot(kcvd3,kcan3,knondis3,log='y',
   ylim=c(1,4000),at=c(0.5,0.6,0.7),boxwex=   0.05,axes=FALSE,col   =
   c("yellow","orange","red"))
   axis(1,at=c(0,0.6,1.1,1.6,2.1,2.6,3.1,3.6),labels=c(0,3,4,5,6,7,8,NA))
   axis(2)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 generate table output of t-test

2008-02-06 Thread John Kane
There may be an easier way but you can extract the
desired values from the list values in t
str(t) to see the elements in t.

test <- matrix(c(1, 1,2,2), 2,2)
tt <- apply(test, 1, t.test)

ttable <- function(tlist) {
  tframe <- data.frame(NULL)
  for(i in 1:length(tlist)){
  t.value  <- tlist[[i]][1]
  dfs <- tlist[[i]][2]
  conf.int1 <- tlist[[i]]$conf.int[1]
  conf.int2 <- tlist[[i]]$conf.int[2]
  p.value <- tlist[[i]][3]
  t.vect <- cbind(t.value, dfs, conf.int1, conf.int2,
p.value)
  tframe <- rbind(tframe, t.vect)
  } 
}
  
myts <- ttable(tt) ; myts
  

--- Ng Stanley <[EMAIL PROTECTED]> wrote:

> Hi,
> 
> Given
> 
> test <- matrix(c(1, 1,2,2), 2,2)
> t <- apply(test, 1, t.test)
> 
> How can I obtain a table of p-values, confidence
> interval etc, instead of
> 
> 
> [[1]]
> 
> One Sample t-test
> 
> data:  newX[, i]
> t = 3, df = 1, p-value = 0.2048
> alternative hypothesis: true mean is not equal to 0
> 95 percent confidence interval:
>  -4.853102  7.853102
> sample estimates:
> mean of x
>   1.5
> 
> 
> [[2]]
> 
> One Sample t-test
> 
> data:  newX[, i]
> t = 3, df = 1, p-value = 0.2048
> alternative hypothesis: true mean is not equal to 0
> 95 percent confidence interval:
>  -4.853102  7.853102
> sample estimates:
> mean of x
>   1.5
> 
>   [[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] wilderSum

2008-02-06 Thread Shubha Vishwanath Karanth
Hi,

 

Can somebody tell me the formula for "?wilderSum" in TTR package? I mean how 
are these calculated?

 

BR, Shubha

Shubha Karanth | Amba Research

Ph +91 80 3980 8031 | Mob +91 94 4886 4510 

Bangalore * Colombo * London * New York * San José * Singapore * 
www.ambaresearch.com

 

This e-mail may contain confidential and/or privileged i...{{dropped:13}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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-Commander - pie charts menu blinded out

2008-02-06 Thread Peter Dalgaard
John Fox wrote:
> Dear Peter and Iksmax,
>
> To elaborate slightly, the Rcmdr tries to figure out which menu items are
> appropriate in a given context, and as Peter says, requires that you have at
> least one factor in the active dataset before activating the pie chart menu
> item; only factors will be included in the pie-chart dialog variable list. I
> suppose that this approach sometimes makes it more difficult than necessary
> to do some things that people legitimately want to do, but the idea is to
> protect students in introductory statistics courses -- the target audience
> for the Rcmdr -- from doing foolish things. (If I were really "holier than
> the prophet" I would have omitted pie charts entirely!)
>
>   
Perhaps stating the obvious: That remark was intended as tongue-in-cheek

Apart from protecting people from doing foolish things (there's a famous
quote about that, though), there's also the aspect that an interface
gets simpler if you reduce the number of choices. So of course it is a
valid design decision.

-p

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

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


Re: [R] help with oop in R - class structure and syntex

2008-02-06 Thread tom soyer
Thanks Gabor for illustrating the basics oop in R using S3. Maybe I didn't
have the right documents, but you example taught me more about oop in R than
everything else I read combined! Thanks for the tip on R.oo, I plan to check
it out later.

I have a few followup questions...

(1) how do I encapsulate the generics? i.e., if a class has 100 methods,
then does it mean 100 generics would be dumped in the global environment?
Or, is it possible to define a local environment and restrict the generics
from one class to a particular environment? But then how do you call the
generics from a special environment? Also, is it possible to inherit classes
across different environment?

(2) it seems that an R specific IDE would improve productivity dramatically
(maybe even necessary?) with respect to oop. Is there such an IDE and does
it work for oop? I recall a group (in Germany?) was working on it but I
can't remember where I read about it.

Thanks!

On 2/5/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
>
> This illustration uses S3.  Note that functions do not modify their
> arguments
> so to modify an object we have to pass it to the method and then pass the
> object back.  There is also another system called S4 which involves typing
> of arguments and there are packages proto and R.oo which provide different
> OO models.
>
>
> # define constructors for bicycle and mountainBike classes
> Bicycle <- function(cadence, gear, speed) structure(list(cadence =
> cadence, gear = gear, speed = speed), class = "bicycle")
> MountainBike <- function(cadence, gear, speed, seatHeight) {
>   x <- Bicycle(cadence, gear, speed)
>   x$seatHeight <- seatHeight
>   class(x) <- c(class(x), "mountainBike")
>   x
> }
>
> # define generic setCadence and then methods for each class
>
> setCadence <- function(x, cadence) UseMethod("setCadence")
> setCadence.bicycle <- function(x, cadence) { x$cadence <- cadence; x }
>
> # mountBike's setCadence method overrides bicycle's setCadence method
> setCadence.mountBike <- function(x, cadence) { x$cadence <- cadence + 1; x
> }
>
> # list the setCadence methods avialable
> methods(setCadence)
>
> # define a generic applyBrake and a "bicycle" method
> # mountainBike will inherit the bicycle method by default
> applyBrake <- function(x, decrement) UseMethod("applyBrake")
> applyBrake.bicycle <- function(x, decrement) { x$speed <- x$speed -
> decrement }
>
> # list the applyBrake methods available
> methods(applyBrake)
>
> b <- Bicycle(1, 2, 3)
> m <- MountainBike(4, 5, 6, 7)
> m <- applyBrake(m, 1)
>
>
> On Feb 5, 2008 8:21 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > Hi,
> >
> > I read section 5, oop, of the R lang doc, and I am still not sure I
> > understand how to build a class in R for oop. I thought that since I
> > understand the oop syntex of Java and VB, I am wondering if the R
> programmig
> > experts could help me out by comparing and contrasting the oop syntex in
> R
> > with that of Java. For example, the basic class structure in Java is
> like
> > this:
> >
> > public class Bicycle {
> >
> >// *the Bicycle class has three fields*
> >public int cadence;
> >public int gear;
> >public int speed;
> >
> >// *the Bicycle class has one constructor*
> >public Bicycle(int startCadence, int startSpeed, int startGear) {
> >gear = startGear;
> >cadence = startCadence;
> >speed = startSpeed;
> >}
> >
> >// *the Bicycle class has four methods*
> >public void setCadence(int newValue) {
> >cadence = newValue;
> >}
> >
> >public void setGear(int newValue) {
> >gear = newValue;
> >}
> >
> >public void applyBrake(int decrement) {
> >speed -= decrement;
> >}
> >
> >public void speedUp(int increment) {
> >speed += increment;
> >}
> >
> > }
> >
> > Could one of the R experts please illustrate the R class syntex for
> writing
> > the R equivalent of the Java Bicycle class above?
> >
> > Also, in Java, inheritance is done like this:
> >
> > public class MountainBike extends Bicycle {
> >
> >// *the MountainBike subclass has one field*
> >public int seatHeight;
> >
> >// *the MountainBike subclass has one constructor*
> >public MountainBike(int startHeight, int startCadence, int
> startSpeed,
> > int startGear) {
> >super(startCadence, startSpeed, startGear);
> >seatHeight = startHeight;
> >}
> >
> >// *the MountainBike subclass has one method*
> >public void setHeight(int newValue) {
> >seatHeight = newValue;
> >}
> >
> > }
> >
> > What would be the R oop syntex for inheritance in the case of the
> > MontainBike class?
> >
> > Thanks!
> >
> > --
> > Tom
> >
> >[[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 pro

Re: [R] Incomplete ouput with sink and split=TRUE

2008-02-06 Thread jiho
On 2008-February-06  , at 14:45 , Duncan Murdoch wrote:
> On 2/5/2008 11:12 AM, jiho wrote:
>> Dear List,
>> I am trying to get R's terminal output to a file and to the  
>> terminal  at the same time, so that I can walk through some tests  
>> and keep a log  concurrently. The function 'sink' with the option  
>> split=TRUE seems to  do just that. It works fine for most output  
>> but for objects of class  htest, the terminal output is incomplete  
>> (the lines are there but  empty). Here is an example session which  
>> shows the problem:
>
> stats:::print.htest() uses writeLines to write some of its output to  
> stdout(), and it looks as though sink(split=T) misses those bits.
>
> I'll change print.htest to use cat(), but it is probably a sign of a  
> bigger problem in sink(), and it's too late in the schedule to touch  
> that for 2.6.2.
>
> Duncan Murdoch


Thank you for your attention on this. Given your comment I filed on a  
bug report on this (same subject that this email, no bug ID yet)

JiHO
---
http://jo.irisson.free.fr/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] counting row repetitions without loop

2008-02-06 Thread James Foadi
On Wednesday 06 February 2008 14:08, Waterman, DG (David) wrote:
> Hi,
>
> I have a data frame consisting of coordinates on a 10*10 grid, i.e.
>
> > example
>
> x  y
> 1   4  5
> 2   6  7
> 3   6  6
> 4   7  5
> 5   5  7
> 6   6  7
> 7   4  5
> 8   6  7
> 9   7  6
> 10  5  6
>
> What I would like to do is return an 10*10 matrix consisting of counts
> at each position, so in the above example I would have a matrix where,
> for example, cell [4,5] contains 2 and [6,7] contains 3. At the moment I
> have implemented this using a for loop over the rows of the data frame,
> however the data frames I want to process are very long so the loop
> takes many minutes to complete. Can I do this in a more efficient way?
>
> Cheers,


David,
have a look at "mapply" (?mapply). This does what you need very quickly.

J

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] counting row repetitions without loop

2008-02-06 Thread Doran, Harold
I think this does what you want, but there may be a more efficient way

x  y
4  5
6  7
6  6
7  5
5  7
6  7
4  5
6  7
7  6
5  6
dat <- read.table('clipboard', header=TRUE) # copy sample data above
dat$patt <- paste(dat$x,dat$y, sep='')
mm <- as.data.frame(with(dat, table(patt)))
dat <- merge(dat, mm, by='patt')
mat <- matrix(0, ncol=10, nrow=10)
gg <- matrix(c(dat$x, dat$y), ncol=2)
mat[gg] <- dat$Freq 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Waterman, 
> DG (David)
> Sent: Wednesday, February 06, 2008 9:08 AM
> To: r-help@r-project.org
> Subject: [R] counting row repetitions without loop
> 
> Hi,
>  
> I have a data frame consisting of coordinates on a 10*10 grid, i.e.
>  
> > example
> x  y
> 1   4  5
> 2   6  7
> 3   6  6
> 4   7  5
> 5   5  7
> 6   6  7
> 7   4  5
> 8   6  7
> 9   7  6
> 10  5  6
> 
> What I would like to do is return an 10*10 matrix consisting 
> of counts at each position, so in the above example I would 
> have a matrix where, for example, cell [4,5] contains 2 and 
> [6,7] contains 3. At the moment I have implemented this using 
> a for loop over the rows of the data frame, however the data 
> frames I want to process are very long so the loop takes many 
> minutes to complete. Can I do this in a more efficient way?
>  
> Cheers,
> David
> This e-mail and any 
> attachments may contain confidential, copyright and or 
> privileged material, and are for the use of the intended 
> addressee only. If you are not the intended addressee or an 
> authorised recipient of the addressee please notify us of 
> receipt by returning the e-mail and do not use, copy, retain, 
> distribute or disclose the information in or attached to the e-mail.
> Any opinions expressed within this e-mail are those of the 
> individual and not necessarily of Diamond Light Source Ltd. 
> Diamond Light Source Ltd. cannot guarantee that this e-mail 
> or any attachments are free from viruses and we cannot accept 
> liability for any damage which you may sustain as a result of 
> software viruses which may be transmitted in or with the message.
> Diamond Light Source Limited (company no. 4375679). 
> Registered in England and Wales with its registered office at 
> Diamond House, Harwell Science and Innovation Campus, Didcot, 
> Oxfordshire, OX11 0DE, United Kingdom  
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] counting row repetitions without loop

2008-02-06 Thread Doran, Harold
Sorry, word wrap made that incomprehensible, I think

x  y
4  5
6  7
6  6
7  5
5  7
6  7
4  5
6  7
7  6
5  6
dat <- read.table('clipboard', header=TRUE)

dat$patt <- paste(dat$x,dat$y, sep='')

mm <- as.data.frame(with(dat, table(patt)))

dat <- merge(dat, mm, by='patt')

mat <- matrix(0, ncol=10, nrow=10)

gg <- matrix(c(dat$x, dat$y), ncol=2)

mat[gg] <- dat$Freq 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of Doran, Harold
> Sent: Wednesday, February 06, 2008 9:56 AM
> To: Waterman, DG (David); r-help@r-project.org
> Subject: Re: [R] counting row repetitions without loop
> 
> I think this does what you want, but there may be a more efficient way
> 
> x  y
> 4  5
> 6  7
> 6  6
> 7  5
> 5  7
> 6  7
> 4  5
> 6  7
> 7  6
> 5  6
> dat <- read.table('clipboard', header=TRUE) # copy sample 
> data above dat$patt <- paste(dat$x,dat$y, sep='') mm <- 
> as.data.frame(with(dat, table(patt))) dat <- merge(dat, mm, 
> by='patt') mat <- matrix(0, ncol=10, nrow=10) gg <- 
> matrix(c(dat$x, dat$y), ncol=2) mat[gg] <- dat$Freq 
> 
> > -Original Message-
> > From: [EMAIL PROTECTED]
> > [mailto:[EMAIL PROTECTED] On Behalf Of Waterman, DG 
> > (David)
> > Sent: Wednesday, February 06, 2008 9:08 AM
> > To: r-help@r-project.org
> > Subject: [R] counting row repetitions without loop
> > 
> > Hi,
> >  
> > I have a data frame consisting of coordinates on a 10*10 grid, i.e.
> >  
> > > example
> > x  y
> > 1   4  5
> > 2   6  7
> > 3   6  6
> > 4   7  5
> > 5   5  7
> > 6   6  7
> > 7   4  5
> > 8   6  7
> > 9   7  6
> > 10  5  6
> > 
> > What I would like to do is return an 10*10 matrix 
> consisting of counts 
> > at each position, so in the above example I would have a 
> matrix where, 
> > for example, cell [4,5] contains 2 and [6,7] contains 3. At 
> the moment 
> > I have implemented this using a for loop over the rows of the data 
> > frame, however the data frames I want to process are very 
> long so the 
> > loop takes many minutes to complete. Can I do this in a 
> more efficient 
> > way?
> >  
> > Cheers,
> > David
> > This e-mail and any 
> attachments may 
> > contain confidential, copyright and or privileged material, and are 
> > for the use of the intended addressee only. If you are not the 
> > intended addressee or an authorised recipient of the 
> addressee please 
> > notify us of receipt by returning the e-mail and do not use, copy, 
> > retain, distribute or disclose the information in or 
> attached to the 
> > e-mail.
> > Any opinions expressed within this e-mail are those of the 
> individual 
> > and not necessarily of Diamond Light Source Ltd.
> > Diamond Light Source Ltd. cannot guarantee that this e-mail or any 
> > attachments are free from viruses and we cannot accept 
> liability for 
> > any damage which you may sustain as a result of software 
> viruses which 
> > may be transmitted in or with the message.
> > Diamond Light Source Limited (company no. 4375679). 
> > Registered in England and Wales with its registered office 
> at Diamond 
> > House, Harwell Science and Innovation Campus, Didcot, Oxfordshire, 
> > OX11 0DE, United Kingdom 
> > 
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/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.


[R] GLM coefficients

2008-02-06 Thread Munyandorero, Joseph
Dear all,

After running a glm, I use the summary ( ) function to extract its
coefficients and related statistics for further use. Unfortunately, the
screen only displays a small (last) part of the results. I tried to
overcome the problem by creating/saving an object "coef" for
coefficients of the model and export/save it e.g. as a cvs document.
While I succed with this operatiion, I do not manage to pick the
standard deviations and other statistics associated with individual
coefficients. What may I do to get the whole set of coefficients and
associated statistics for my GLM? Thanks for your help.

J. Munyandorero

[[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] scatterplot3d + box lines

2008-02-06 Thread Pedro Mardones
Dear all;
I've been trying to change the type of line used to draw the box
around the 3d scatterplot (package scatterplot3d) from lty=1 to lty=2
without sucess. I would appreciate suggestions of how to do it.
Thanks
PM

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


Re: [R] scatterplot3d + box lines

2008-02-06 Thread Duncan Murdoch
On 2/6/2008 10:07 AM, Pedro Mardones wrote:
> Dear all;
> I've been trying to change the type of line used to draw the box
> around the 3d scatterplot (package scatterplot3d) from lty=1 to lty=2
> without sucess. I would appreciate suggestions of how to do it.

Use lty.axis=2.  (This is mentioned on the help page, but there are a 
lot of parameters there, and I can see how you might have missed it.)

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] GLM coefficients

2008-02-06 Thread Ioannis Kosmidis
If I understand correctly.. using the example in ?glm

 > utils::data(anorexia, package="MASS")
 > anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
 family = gaussian, data = anorexia)
 > summ <- summary(anorex.1)

Then check

 > names(summ)

for the available objects.

For coefficients standard errors and t-tests

 > summ$coefficients

Again do this for any object other than "coefficients" you would want  
from names(summ).

Best regards,

Ioannis Kosmidis


On 6 Feb 2008, at 15:06, Munyandorero, Joseph wrote:

> Dear all,
>
> After running a glm, I use the summary ( ) function to extract its
> coefficients and related statistics for further use. Unfortunately,  
> the
> screen only displays a small (last) part of the results. I tried to
> overcome the problem by creating/saving an object "coef" for
> coefficients of the model and export/save it e.g. as a cvs document.
> While I succed with this operatiion, I do not manage to pick the
> standard deviations and other statistics associated with individual
> coefficients. What may I do to get the whole set of coefficients and
> associated statistics for my GLM? Thanks for your help.
>
> J. Munyandorero
>
>   [[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.

===
Dr Ioannis Kosmidis
Department of Statistics ,
University of Warwick,
Coventry, CV4 7AL,
United Kingdom

Email   :   [EMAIL PROTECTED]
Home:   http://go.warwick.ac.uk/kosmidis
Voice   :   +44(0)2476 575754
Fax :   +44(0)2476 524532

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

2008-02-06 Thread Henrique Dallazuanna
See this example:

 coef(summary(glm(Postwt ~ Prewt + Treat +
offset(Prewt),family=gaussian, data=anorexia)))

Is this you want?

On 06/02/2008, Munyandorero, Joseph <[EMAIL PROTECTED]> wrote:
> Dear all,
>
> After running a glm, I use the summary ( ) function to extract its
> coefficients and related statistics for further use. Unfortunately, the
> screen only displays a small (last) part of the results. I tried to
> overcome the problem by creating/saving an object "coef" for
> coefficients of the model and export/save it e.g. as a cvs document.
> While I succed with this operatiion, I do not manage to pick the
> standard deviations and other statistics associated with individual
> coefficients. What may I do to get the whole set of coefficients and
> associated statistics for my GLM? Thanks for your help.
>
> J. Munyandorero
>
> [[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

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

2008-02-06 Thread Christopher Chizinski
Is there anyway to produce a 3D scatterplot with error bars in the x,y,z 
directions?  I have searched around and know of scatterplot3d but did 
not see any way to put error bars on the points.  Any ideas?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 oop in R - class structure and syntex

2008-02-06 Thread Gabor Grothendieck
On Feb 6, 2008 9:45 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> Thanks Gabor for illustrating the basics oop in R using S3. Maybe I didn't
> have the right documents, but you example taught me more about oop in R than
> everything else I read combined! Thanks for the tip on R.oo, I plan to check
> it out later.
>
> I have a few followup questions...
>
> (1) how do I encapsulate the generics? i.e., if a class has 100 methods,
> then does it mean 100 generics would be dumped in the global environment?
> Or, is it possible to define a local environment and restrict the generics
> from one class to a particular environment? But then how do you call the
> generics from a special environment? Also, is it possible to inherit classes
> across different environment?

Both ofthe following return 3.  If that is not what you are asking
then I suggest you experiment with a few small examples of
this nature:

f <- function() {
F <- function(x) UseMethod("F")
F.default <- function(x) x
F(3)
}
F.default <- function(x) NULL
f()


F <- function(x) UseMethod("F")
f <- function() {
F.default <- function(x) x
F(3)
}
F.default <- function(x) NULL
f()

>
> (2) it seems that an R specific IDE would improve productivity dramatically
> (maybe even necessary?) with respect to oop. Is there such an IDE and does
> it work for oop? I recall a group (in Germany?) was working on it but I
> can't remember where I read about it.
>
> Thanks!
>
>
> On 2/5/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> > This illustration uses S3.  Note that functions do not modify their
> arguments
> > so to modify an object we have to pass it to the method and then pass the
> > object back.  There is also another system called S4 which involves typing
> > of arguments and there are packages proto and R.oo which provide different
> > OO models.
> >
> >
> > # define constructors for bicycle and mountainBike classes
> > Bicycle <- function(cadence, gear, speed) structure(list(cadence =
> > cadence, gear = gear, speed = speed), class = "bicycle")
> > MountainBike <- function(cadence, gear, speed, seatHeight) {
> >   x <- Bicycle(cadence, gear, speed)
> >   x$seatHeight <- seatHeight
> >   class(x) <- c(class(x), "mountainBike")
> >   x
> > }
> >
> > # define generic setCadence and then methods for each class
> >
> > setCadence <- function(x, cadence) UseMethod("setCadence")
> > setCadence.bicycle <- function(x, cadence) { x$cadence <- cadence; x }
> >
> > # mountBike's setCadence method overrides bicycle's setCadence method
> > setCadence.mountBike <- function(x, cadence) { x$cadence <- cadence + 1; x
> }
> >
> > # list the setCadence methods avialable
> > methods(setCadence)
> >
> > # define a generic applyBrake and a "bicycle" method
> > # mountainBike will inherit the bicycle method by default
> > applyBrake <- function(x, decrement) UseMethod("applyBrake")
> > applyBrake.bicycle <- function(x, decrement) { x$speed <- x$speed -
> decrement }
> >
> > # list the applyBrake methods available
> > methods(applyBrake)
> >
> > b <- Bicycle(1, 2, 3)
> > m <- MountainBike(4, 5, 6, 7)
> > m <- applyBrake(m, 1)
> >
> >
> > On Feb 5, 2008 8:21 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > Hi,
> > >
> > > I read section 5, oop, of the R lang doc, and I am still not sure I
> > > understand how to build a class in R for oop. I thought that since I
> > > understand the oop syntex of Java and VB, I am wondering if the R
> programmig
> > > experts could help me out by comparing and contrasting the oop syntex in
> R
> > > with that of Java. For example, the basic class structure in Java is
> like
> > > this:
> > >
> > > public class Bicycle {
> > >
> > >// *the Bicycle class has three fields*
> > >public int cadence;
> > >public int gear;
> > >public int speed;
> > >
> > >// *the Bicycle class has one constructor*
> > >public Bicycle(int startCadence, int startSpeed, int startGear) {
> > >gear = startGear;
> > >cadence = startCadence;
> > >speed = startSpeed;
> > >}
> > >
> > >// *the Bicycle class has four methods*
> > >public void setCadence(int newValue) {
> > >cadence = newValue;
> > >}
> > >
> > >public void setGear(int newValue) {
> > >gear = newValue;
> > >}
> > >
> > >public void applyBrake(int decrement) {
> > >speed -= decrement;
> > >}
> > >
> > >public void speedUp(int increment) {
> > >speed += increment;
> > >}
> > >
> > > }
> > >
> > > Could one of the R experts please illustrate the R class syntex for
> writing
> > > the R equivalent of the Java Bicycle class above?
> > >
> > > Also, in Java, inheritance is done like this:
> > >
> > > public class MountainBike extends Bicycle {
> > >
> > >// *the MountainBike subclass has one field*
> > >public int seatHeight;
> > >
> > >// *the MountainBike subclass has one constructor*
> > >public MountainBike(int startHeight, int startCadence, int
> startSpeed,
> > >

Re: [R] Inconsistent lattice scales$x$at, label behaviour for POSIXct

2008-02-06 Thread Alex Brown
Here's a related problem that works for numeric but not for POSIXct

I am seeing it where a panel has no at labels, but others do.

This simple example only has one panel with no at labels.

 > baseval = 0;
 > xyplot(1:10 ~ (baseval + c(1:10)), scales=list(x=list
+ (at=list(c()), labels=list(c()), relation="free")), type="l")

(works)

baseval = as.POSIXct("2007-01-01");
 > xyplot(1:10 ~ (baseval + c(1:10)), scales=list(x=list
+ (at=list(c()), labels=list(c()), relation="free")), type="l")
Error in as.POSIXct.default(at) :
  do not know how to convert 'at' to class "POSIXlt"

(fails)

-Alex


On 5 Feb 2008, at 20:59, Deepayan Sarkar wrote:

> On 2/5/08, Alex Brown <[EMAIL PROTECTED]> wrote:
>> I have encountered the following behaviour in lattice in 2.6.1 (and
>> 2.4.0) which differs depending upon the type you use.  I believe the
>> numeric behaviour to be correct, and the POSIXct behaviour to be in
>> error.
>>
>> When the x data and x axis in a lattice graph are POSIXct, and when
>> using scales$x$at and scales$x$labels to add custom labels:
>>
>> If the first visible at value is not the first specified at value  
>> (due
>> to x limit settings), the first visible label nonetheless receives  
>> the
>> first specified label, instead of the label corresponding to the  
>> first
>> visible at.
>
> Yes, the relevant "POSIXct" method was subsetting 'at' to lie within
> the specified limits but not 'labels'. Should be fixed in the next
> update. Thanks,
>
> -Deepayan
>
>
>> In the following example, at = baseval + 1:4, label = letters[1:4]
>>
>> I have set it up so that baseval + 1:2 are not visible in the graph
>>
>> for numeric and Date types, the visible labels are letters[3:4] - ie
>> "c" "d"
>>
>> for POSIXct, the visible labels are letters[1:3] - "a" "b".
>>
>>
>> Simplest to show this by example:
>>
>> # numeric
>>
>> baseval=0;
>> xyplot(1:10 ~ (baseval + 1:10) , scales=list(x=list
>> (at=baseval+1:4, labels=letters[1:4])),xlim=baseval+c(3,10))
>>
>> # Date
>>
>> baseval = as.Date("2007-01-01");
>> xyplot(1:10 ~ (baseval + 1:10) , scales=list(x=list
>> (at=baseval+1:4, labels=letters[1:4])),xlim=baseval+c(3,10))
>>
>> # POSIXct
>>
>> baseval = as.POSIXct("2007-01-01");
>> xyplot(1:10 ~ (baseval + 1:10) , scales=list(x=list
>> (at=baseval+1:4, labels=letters[1:4])),xlim=baseval+c(3,10))
>>
>> in particular, compare the Date and POSIXct versions
>>
>> -Alex Brown
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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.


[R] kinship package: drawing pedigree error

2008-02-06 Thread Iris Kolder
Hi

Im using the kinship package to draw a pedigree. On my data set this works fine 
but when i add indivudals to the pedigree i keep getting an error i hope 
someone can help me!

This is the code im using:

Data<-read.table("Tree.txt", header=T, sep=",")
attach(Data)
ped<-pedigree(id, dadid, momid, sex, aff)
par(xpd=T)
plot.pedigree(ped)

This is my data looks like without the last 3 individuals it works perfect,when 
i add them i get the following error:
Error in if (min(pos2) < 0) { : missing value where TRUE/FALSE needed

famid,id,dadid,momid,sex,aff
1,8860,9972,8856,2,0
1,8855,9972,8856,2,0
1,8859,9976,8860,2,0
1,8854,9971,8855,2,0
1,8863,9975,8859,2,0
1,8858,9975,8859,2,0
1,8865,9975,8859,2,0
1,9969,9970,8854,1,0
1,8864,9980,8863,2,0
1,8862,9980,8863,2,0
1,23834,9981,8865,2,0
1,9968,9969,8853,1,0
1,21141,9974,8858,1,0
1,21142,9974,8858,2,0
1,23172,265,21142,2,0
1,25409,265,21142,1,0
1,23171,265,21142,2,0
1,265,NA,NA,1,0
1,9981,NA,NA,1,0
1,8853,NA,NA,2,0
1,9974,NA,NA,1,0
1,9980,NA,NA,1,0
1,9972,NA,NA,1,0
1,8856,NA,NA,2,0
1,9976,NA,NA,1,0
1,9971,NA,NA,1,0
1,9975,NA,NA,1,0
1,9970,NA,NA,1,0

1,9979,NA,NA,1,0
1,20644,9979,8862,2,0
1,20670,9979,8862,1,0

I tried putting nuclear families closer together in de pedigree but this 
doesn't seem to help. All individuals have both parents in the pedigree or have 
2 missing parents.


I hope someone can help me with this!

Iris


  

Never miss a thing.  Make Yahoo your home page. 

[[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] Multivariate Maximum Likelihood Estimation

2008-02-06 Thread Konrad BLOCHER
OK I got gls running but this is the error message that I got:

> gls_results <- gls(LnWRPK ~ LnWGDP+LnWYIELD,correlation = corAR1)

Error in switch(mode(x), "NULL" = structure(NULL, class = "formula"),  :
  invalid formula

Google didn't give any results unfortunately.

Help! :)

Thanks,

KB

> Have you tried gls()?
>
>
> 

> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
> tel. + 32 54/436 185
> [EMAIL PROTECTED]
> www.inbo.be
>
> Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
> A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney
>
> -Oorspronkelijk bericht-
> Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens Konrad BLOCHER
> Verzonden: woensdag 6 februari 2008 12:12
> Aan: r-help@r-project.org
> Onderwerp: [R] Multivariate Maximum Likelihood Estimation
>
> Hi,
>
> I am trying to perform Maximum Likelihood estimation of a Multivariate
model (2 independent variables + intercept) with autocorrelated errors
of
> 1st order (ar(1)).
>
> Does R have a function for that? I could only find an univariate option
(ar.mle function) and when writing my own I find that it is pretty
memory-consuming (and sometimes wrong) so there must be a better way.
>
> Thanks,
>
> KB
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Regression with time-dependent coefficients

2008-02-06 Thread achronal
Hi,

I was wondering if someone might be willing to indulge a question about
R and the estimation of a linear regression with time-varying coefficients.

The model I am trying to estimate is of the form:

y(t) = beta(t) * x(t) + v(t)
beta(t) = gamma * beta(t-1) + w(t)

where gamma is a constant, v(t) and w(t) are Gaussian innovations,
and where y(t) and x(t) are univariate time series that are known.
I would like to estimate gamma and the variance of v(t) and w(t).

I have tried using the following packages with no success:  sspir,
dse (specifically dse1) and dlm.  In each case, any attempt to
specify the above model meets with some sort of error.  Hence, I
have not been able to estimate the model.  I am using "manufactured"
data, i.e., I assume a value for gamma, beta(0), and values
for the variances of v(t) and w(t), so I know what result I
should attain.

I would really appreciate any help on how to appropriately specify
the above model in any of the packages I mentioned.

Best regards,
-Paul

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

2008-02-06 Thread Gang Chen
Normally I can run an R script in batch mode with a command like this

R CMD BATCH MyScript.R MyOutput &

However I prefer to write another script containing something like

R --no-restore --save --no-readline < $1 >$2

so that I could run the original script simply on the prompt as

MyScript.R MyOutput &

The second method is almost the same as the R CMD BATCH method except  
MyOutput does not contain the output from the OS calls such as system 
(). So what option(s) should add in the second approach to make it  
identical to R CMD BATCH?

Thanks,
Gang

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 oop in R - class structure and syntex

2008-02-06 Thread Jeffrey J. Hallman
"tom soyer" <[EMAIL PROTECTED]> writes:
> (1) how do I encapsulate the generics? i.e., if a class has 100 methods,
> then does it mean 100 generics would be dumped in the global environment?
> Or, is it possible to define a local environment and restrict the generics
> from one class to a particular environment? But then how do you call the
> generics from a special environment? Also, is it possible to inherit classes
> across different environment?

I think the answer is: you can't.  There is no sensible way to assign
multi-methods to a particular class.  So you have to have some sort of
database of methods and signatures that gets consulted by the dispatcher.  How 
you do a sensible IDE with this setup is beyond me.  This is my biggest
gripe with S4 methods: there doesn't seem to be any good way to browse the
environment and discover whats already in there.  S3 methods are much easier
to follow. You can do multiple dispatch with S3 methods, but hardly anyone
ever does.  My view is that multiple dispatch is not worth the complexity it
adds to a language. 

> (2) it seems that an R specific IDE would improve productivity dramatically
> (maybe even necessary?) with respect to oop. Is there such an IDE and does
> it work for oop? I recall a group (in Germany?) was working on it but I
> can't remember where I read about it.

-- 
Jeff

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 oop in R - class structure and syntex

2008-02-06 Thread tom soyer
Thanks Gabor. I guess true oo encapsulation is not possible in R.

It seems that there is an IDE for S+  in Eclipse...


On 2/6/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
>
> On Feb 6, 2008 9:45 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > Thanks Gabor for illustrating the basics oop in R using S3. Maybe I
> didn't
> > have the right documents, but you example taught me more about oop in R
> than
> > everything else I read combined! Thanks for the tip on R.oo, I plan to
> check
> > it out later.
> >
> > I have a few followup questions...
> >
> > (1) how do I encapsulate the generics? i.e., if a class has 100 methods,
> > then does it mean 100 generics would be dumped in the global
> environment?
> > Or, is it possible to define a local environment and restrict the
> generics
> > from one class to a particular environment? But then how do you call the
> > generics from a special environment? Also, is it possible to inherit
> classes
> > across different environment?
>
> Both ofthe following return 3.  If that is not what you are asking
> then I suggest you experiment with a few small examples of
> this nature:
>
> f <- function() {
>F <- function(x) UseMethod("F")
>F.default <- function(x) x
>F(3)
> }
> F.default <- function(x) NULL
> f()
>
>
> F <- function(x) UseMethod("F")
> f <- function() {
>F.default <- function(x) x
>F(3)
> }
> F.default <- function(x) NULL
> f()
>
> >
> > (2) it seems that an R specific IDE would improve productivity
> dramatically
> > (maybe even necessary?) with respect to oop. Is there such an IDE and
> does
> > it work for oop? I recall a group (in Germany?) was working on it but I
> > can't remember where I read about it.
> >
> > Thanks!
> >
> >
> > On 2/5/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> > > This illustration uses S3.  Note that functions do not modify their
> > arguments
> > > so to modify an object we have to pass it to the method and then pass
> the
> > > object back.  There is also another system called S4 which involves
> typing
> > > of arguments and there are packages proto and R.oo which provide
> different
> > > OO models.
> > >
> > >
> > > # define constructors for bicycle and mountainBike classes
> > > Bicycle <- function(cadence, gear, speed) structure(list(cadence =
> > > cadence, gear = gear, speed = speed), class = "bicycle")
> > > MountainBike <- function(cadence, gear, speed, seatHeight) {
> > >   x <- Bicycle(cadence, gear, speed)
> > >   x$seatHeight <- seatHeight
> > >   class(x) <- c(class(x), "mountainBike")
> > >   x
> > > }
> > >
> > > # define generic setCadence and then methods for each class
> > >
> > > setCadence <- function(x, cadence) UseMethod("setCadence")
> > > setCadence.bicycle <- function(x, cadence) { x$cadence <- cadence; x }
> > >
> > > # mountBike's setCadence method overrides bicycle's setCadence method
> > > setCadence.mountBike <- function(x, cadence) { x$cadence <- cadence +
> 1; x
> > }
> > >
> > > # list the setCadence methods avialable
> > > methods(setCadence)
> > >
> > > # define a generic applyBrake and a "bicycle" method
> > > # mountainBike will inherit the bicycle method by default
> > > applyBrake <- function(x, decrement) UseMethod("applyBrake")
> > > applyBrake.bicycle <- function(x, decrement) { x$speed <- x$speed -
> > decrement }
> > >
> > > # list the applyBrake methods available
> > > methods(applyBrake)
> > >
> > > b <- Bicycle(1, 2, 3)
> > > m <- MountainBike(4, 5, 6, 7)
> > > m <- applyBrake(m, 1)
> > >
> > >
> > > On Feb 5, 2008 8:21 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > > Hi,
> > > >
> > > > I read section 5, oop, of the R lang doc, and I am still not sure I
> > > > understand how to build a class in R for oop. I thought that since I
> > > > understand the oop syntex of Java and VB, I am wondering if the R
> > programmig
> > > > experts could help me out by comparing and contrasting the oop
> syntex in
> > R
> > > > with that of Java. For example, the basic class structure in Java is
> > like
> > > > this:
> > > >
> > > > public class Bicycle {
> > > >
> > > >// *the Bicycle class has three fields*
> > > >public int cadence;
> > > >public int gear;
> > > >public int speed;
> > > >
> > > >// *the Bicycle class has one constructor*
> > > >public Bicycle(int startCadence, int startSpeed, int startGear) {
> > > >gear = startGear;
> > > >cadence = startCadence;
> > > >speed = startSpeed;
> > > >}
> > > >
> > > >// *the Bicycle class has four methods*
> > > >public void setCadence(int newValue) {
> > > >cadence = newValue;
> > > >}
> > > >
> > > >public void setGear(int newValue) {
> > > >gear = newValue;
> > > >}
> > > >
> > > >public void applyBrake(int decrement) {
> > > >speed -= decrement;
> > > >}
> > > >
> > > >public void speedUp(int increment) {
> > > >speed += increment;
> > > >}
> > > >
> > > > }
> > > >
> > > > Could one of the R expe

[R] Discriminant function analysis

2008-02-06 Thread Birgit Lemcke
Hello R-Cracks,

I am using R 2.6.1 on a PowerBook G4.
I would like to perform a discriminant function analysis. I found lda  
in MASS but as far as I understood, is it only working with  
explanatory variables of the class factor. My dataset contains  
variables of the classes factor and numeric.
Is there another function that is able to handle this?

Thank you all in advance.

Greetings

Birgit

Birgit Lemcke
Institut für Systematische Botanik
Zollikerstrasse 107
CH-8008 Zürich
Switzerland
Ph: +41 (0)44 634 8351
[EMAIL PROTECTED]

175 Jahre UZH
«staunen.erleben.begreifen. Naturwissenschaft zum Anfassen.»
MNF-Jubiläumsevent für gross und klein.
19. April 2008, 10.00 Uhr bis 02.00 Uhr
Campus Irchel, Winterthurerstrasse 190, 8057 Zürich
Weitere Informationen www.175jahre.uzh.ch

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

2008-02-06 Thread ONKELINX, Thierry
It seems like you didn't look at the examples in the helpfiles. See ?gls
and ?corAR1.




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

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney

-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens Konrad BLOCHER
Verzonden: woensdag 6 februari 2008 16:42
Aan: r-help@r-project.org
Onderwerp: Re: [R] Multivariate Maximum Likelihood Estimation

OK I got gls running but this is the error message that I got:

> gls_results <- gls(LnWRPK ~ LnWGDP+LnWYIELD,correlation = corAR1)

Error in switch(mode(x), "NULL" = structure(NULL, class = "formula"),  :
  invalid formula

Google didn't give any results unfortunately.

Help! :)

Thanks,

KB

> Have you tried gls()?
>
>
>


> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
> tel. + 32 54/436 185
> [EMAIL PROTECTED]
> www.inbo.be
>
> Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
> A statistical analysis, properly conducted, is a delicate dissection
of
uncertainties, a surgery of suppositions. ~M.J.Moroney
>
> -Oorspronkelijk bericht-
> Van: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
Namens Konrad BLOCHER
> Verzonden: woensdag 6 februari 2008 12:12
> Aan: r-help@r-project.org
> Onderwerp: [R] Multivariate Maximum Likelihood Estimation
>
> Hi,
>
> I am trying to perform Maximum Likelihood estimation of a Multivariate
model (2 independent variables + intercept) with autocorrelated errors
of
> 1st order (ar(1)).
>
> Does R have a function for that? I could only find an univariate
option
(ar.mle function) and when writing my own I find that it is pretty
memory-consuming (and sometimes wrong) so there must be a better way.
>
> Thanks,
>
> KB
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] unload & reload a (new version of a) package

2008-02-06 Thread Harte, Thomas P
i returned to this problem and found that the solution was a disarmingly

straight-forward oversight on my part: 'detach' requires the version
number of the library to be specified! 

thus, if you install packages '--with-package-versions' switched on,
and if you load the version of choice in an R session with:

> library(foo, version="1.0")

then you must 'detach' it in the same R session using

> detach(foo, version="1.0")

in the *same* R session you can load another version:

> library(foo, version="1.2")

and happily detach it with

> detach(foo, version="1.2")

and then go back to version 1.0 again (all in the same R session):

> library(foo, version="1.0")

thus obviating the behavior that i had previously reported.



-Original Message-
From: Harte, Thomas P 
Sent: Tuesday, January 15, 2008 10:47 AM
To: r-help@r-project.org
Subject: RE: [R] unload & reload a (new version of a) package


update: this phenomenon disappears if i do a straight install *without*
package versioning, i.e. 

C:\packages>c:\R\R-2.6.1\bin\R CMD BUILD foo_1.0
[creates foo_1.0.tar.gz]
C:\packages>c:\R\R-2.6.1\bin\R CMD INSTALL -l c:/library foo_1.0.tar.gz

i can update the package, re-build and re-install and load it into the
existing R session with a call to 'library(foo, lib.loc="c:/library")'

when i posted, i had been doing things *with versioning* switched on:

C:\packages>c:\R\R-2.6.1\bin\R CMD BUILD foo_1.0
[creates foo_1.0.tar.gz]
C:\packages>c:\R\R-2.6.1\bin\R CMD INSTALL --with-package-versions -l
c:/library foo_1.0.tar.gz

and loading it into the existing R session with 'library(foo,
lib.loc="c:/library", version="1.0")'.

so the problem seemes to exist if i install a newer version of the
package this way, thereby clobbering the old version foo_1.0 in the
process. however, if I copy the source tree for foo_1.0 and make a
*completely new version* foo_1.1, say, and install
--with-package-versions:

C:\packages>c:\R\R-2.6.1\bin\R CMD BUILD foo_1.1 
C:\packages>c:\R\R-2.6.1\bin\R CMD INSTALL --with-package-versions -l
c:/library foo_1.1.zip 
library(foo, lib.loc="c:/library", version="1.1")

then i can successfully load version 1.1 

something is possibly wrong somewhere with versioning, or with my
package. either way, i can live without versioning switched on and just
clobber versions of the library foo with the updated package foo.


-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On Behalf Of Harte, Thomas P
Sent: Monday, January 14, 2008 6:49 PM
To: r-help@r-project.org
Subject: [R] unload & reload a (new version of a) package

i'm putting the final touches on a package that i'm developing and i
noticed that if i detach the package, and then re-build & re-install it
(using R
CMD)
then I can't get the newer version of the package to load in the
existing R session (i have to close it out and start a new session, then
the newer 
version of the package is loaded). 

looking through the source of 'detach'  i see :

  .Call("R_lazyLoadDBflush", paste(libpath, "/R/", pkgname, 
".rdb", sep = ""), PACKAGE = "base")

is there some absolute way similar to the above to flush the package db 
and ensure that a newer version of the package can be loaded into the
existing R session? detach calls .Last.lib and seems to go through the
motions of purging the loaded package; why, then, is the package still
lurking around in the existing R session?

it's not a big deal; it's only a minor pain having to re-start an R
session. i'm more interesting in why this is happening.

cheers,

thomas.   






This message, including any attachments, contains\ confi...{{dropped:26}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Levelplot of percentages always using 0 to 100 in the color scheme

2008-02-06 Thread Deepayan Sarkar
On 2/5/08, Kelvin <[EMAIL PROTECTED]> wrote:
> I am trying to create levelplot's of cpu usage for systems.
> print(levelplot(util.mean ~ x.hour * x.day, colorkey=T, cut=20,
> scales=list(x=list(at=seq(0,96,length=25),
> labels=ifelse(seq(0,24) %% 4 == 0, seq(0,24), ''))),   # add
> tick marks at the hour
> main="CPU Utilization During November on idnprod",
>
> col.regions=colorRampPalette(c('white','green','yellow','red'))(101), #
> 101 colors
> xlab="Hours (15 minute intervals)", ylab="Day of the Month")
> )
> util.mean has the utilizations by timestamps
> If the values in util.mean are only 10-50, the plot codes 10 as white
> and low and 50 as red and high.
> What I would desire is the 0 is always the low and white and 100 always
> the high and red.

Specify a suitable 'at' argument; e.g.,

at = seq(0, 100, length = 51)

-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] advice requested re: building "good" system (R, SQL db) for handling large datasets

2008-02-06 Thread Jeffrey Horner
Richard Pearson wrote on 02/06/2008 06:25 AM:
> Hi Thomas
[...]
> With databases, one issue that might be relevant is whether you want to 
> store data in tables (e.g. one table to store one data.frame) that can 
> subsequently be manipulated in the DB, or to store R objects as R 
> objects (e.g. as BLOBs). My situation is likely to be the later case, 
> and one of my concerns is that many DBs have an upper limit of 2GB on 
> BLOBs, and I might potentially have objects that are larger than this.
[...]

I'd be curious as to why you'd want to store and retrieve R objects from 
a BLOB column in a table. I've often thought about this, but 
unfortunately neither the DBI package nor the RODBC package support this.

Jeff
-- 
http://biostat.mc.vanderbilt.edu/JeffreyHorner

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

2008-02-06 Thread Josh Ulrich
Shubha,

The calculation is:

wilderSum[1] <- x[1]
for(i in 2:NROW(x)) {
  wilderSum[i] <- x[i] + wilderSum[i-1] * (n-1)/n
}

Where 'x' is the price series and 'n' is the number of periods.

You can see the calculations for all functions in TTR's source code,
which is on CRAN and r-forge
(http://r-forge.r-project.org/projects/ttr/).

Please feel free to contact me directly regarding questions about TTR.

Best,
Josh

On Feb 6, 2008 8:22 AM, Shubha Vishwanath Karanth
<[EMAIL PROTECTED]> wrote:
>
> Hi,
>
> Can somebody tell me the formula for "?wilderSum" in TTR package? I mean how
> are these calculated?
>
> BR, Shubha
>
> Shubha Karanth | Amba Research
>
> Ph +91 80 3980 8031 | Mob +91 94 4886 4510
>
> Bangalore • Colombo • London • New York • San José • Singapore •
> www.ambaresearch.com
>
>   This e-mail may contain confidential and/or privileged information. If you
> are not the intended recipient (or have received this
> e-mail in error) please notify the sender immediately and destroy this
> e-mail. Any unauthorized copying, disclosure or distribution of
> the material in this e-mail is strictly forbidden. Any views or opinions
> presented are solely those of the author and do not
> necessarily represent those of Amba Holdings Inc., and/or its affiliates.
> Important additional terms relating to this email can be obtained
> at http://www.ambaresearch.com/disclaimer
>
>

-- 
http://quantemplation.blogspot.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] 3d scatterplot with error bars

2008-02-06 Thread Uwe Ligges


Christopher Chizinski wrote:
> Is there anyway to produce a 3D scatterplot with error bars in the x,y,z 
> directions?  I have searched around and know of scatterplot3d but did 
> not see any way to put error bars on the points.  Any ideas?


Yes, type

vignette("s3d", package="scatterplot3d")

and see page 22.

Uwe Ligges




> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] advice requested re: building "good" system (R, SQL db) for handling large datasets

2008-02-06 Thread Dirk Eddelbuettel
On Wed, Feb 06, 2008 at 02:34:52PM +0200, Rainer M Krug wrote:
> R objects in blobs - I never thought about that. Could you elaborate on how
> to do something like that (I am using RMySQL)?

Look at help(serialize) -- any R object can be turned into a suitable
representation, either binary (more efficient) or ascii (possibly
'safer'). 

Store that, retrieve it later, reconstruct the object and be merry :) 

> tmpDf <- data.frame(a=1:10, b=LETTERS[1:10], c=rnorm(10))
> serDf <- serialize(tmpDf, NULL, ascii=TRUE)
> rm(tmpDf)
> head(unserialize(serDf))
  a b  c
  1 1 A -0.6945820
  2 2 B -0.2960084
  3 3 C -0.2514302
  4 4 D -0.7318635
  5 5 E -0.1698489
  6 6 F  0.4331521
  > 

Dirk

-- 
Three out of two people have difficulties with fractions.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Inconsistent lattice scales$x$at, label behaviour for POSIXct

2008-02-06 Thread Deepayan Sarkar
On 2/6/08, Alex Brown <[EMAIL PROTECTED]> wrote:
> Here's a related problem that works for numeric but not for POSIXct
>
> I am seeing it where a panel has no at labels, but others do.
>
> This simple example only has one panel with no at labels.
>
>  > baseval = 0;
>  > xyplot(1:10 ~ (baseval + c(1:10)), scales=list(x=list
> + (at=list(c()), labels=list(c()), relation="free")), type="l")
>
> (works)
>
> baseval = as.POSIXct("2007-01-01");
>  > xyplot(1:10 ~ (baseval + c(1:10)), scales=list(x=list
> + (at=list(c()), labels=list(c()), relation="free")), type="l")
> Error in as.POSIXct.default(at) :
>   do not know how to convert 'at' to class "POSIXlt"

There may be a better long-term solution, but does this work for now?:

> xyplot(1:10 ~ (baseval + c(1:10)), scales=list(x=list
+ (at=list(baseval + c()), labels=list(c()), relation="free")), type="l")

-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] 3d scatterplot with error bars

2008-02-06 Thread hadley wickham
> Is there anyway to produce a 3D scatterplot with error bars in the x,y,z
> directions?  I have searched around and know of scatterplot3d but did
> not see any way to put error bars on the points.  Any ideas?

It may be possible, but do you think that's the best solution for your
problem?  Given that any 3d plot must be projected into two
dimensions, how on earth will people be able to compare the relative
magnitude in each direction or between points?

One could imagine drawing error ellispoids (since you might as well
include covariance information as well) but I would imagine that the
tight connection between the view of the data and the two-D projection
of the ellipsoid (representing a particular linear combination of the
x, y and z errors) would make it very difficult to see what's going
on.  Not to mention the difficulties with near points obscuring far
points, etc etc.

Hadley

-- 
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] 3d scatterplot with error bars

2008-02-06 Thread hadley wickham
> Yes, type
>
> vignette("s3d", package="scatterplot3d")
>
> and see page 22.

which of course won't work unless you have scatterplot3d installed,
but I wonder if vignette() could be modified to try to find the
vignette on the web if the package isn't installed.  It would make it
a little easy to read about a package before installing it.

Hadley

-- 
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] Multivariate Maximum Likelihood Estimation

2008-02-06 Thread Konrad BLOCHER
Thanks,

I managed to run gls, but my problem isn't solved :)

1. I do not know the autocorrelation between errors and its estimates
change for different methods. Therefore I wasn't including it in the
corAR1 function.

2. GLS yields counteraintuitive results with one of the variables (its
sign) and I need to make sure that it is either only the case for Least
Squares, or for the data (disregarding the influence of different
estimation techniques) - thus a need for multivariate MLE.

Cheers,

KB


> It seems like you didn't look at the examples in the helpfiles. See ?gls
> and ?corAR1.
>
>
> 
> 
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
> tel. + 32 54/436 185
> [EMAIL PROTECTED]
> www.inbo.be
>
> Do not put your faith in what statistics say until you have carefully
> considered what they do not say.  ~William W. Watt
> A statistical analysis, properly conducted, is a delicate dissection of
> uncertainties, a surgery of suppositions. ~M.J.Moroney
>
> -Oorspronkelijk bericht-
> Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> Namens Konrad BLOCHER
> Verzonden: woensdag 6 februari 2008 16:42
> Aan: r-help@r-project.org
> Onderwerp: Re: [R] Multivariate Maximum Likelihood Estimation
>
> OK I got gls running but this is the error message that I got:
>
>> gls_results <- gls(LnWRPK ~ LnWGDP+LnWYIELD,correlation = corAR1)
>
> Error in switch(mode(x), "NULL" = structure(NULL, class = "formula"),  :
>   invalid formula
>
> Google didn't give any results unfortunately.
>
> Help! :)
>
> Thanks,
>
> KB
>
>> Have you tried gls()?
>>
>>
>>
> 
> 
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
>> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
>> Gaverstraat 4
>> 9500 Geraardsbergen
>> Belgium
>> tel. + 32 54/436 185
>> [EMAIL PROTECTED]
>> www.inbo.be
>>
>> Do not put your faith in what statistics say until you have carefully
> considered what they do not say.  ~William W. Watt
>> A statistical analysis, properly conducted, is a delicate dissection
> of
> uncertainties, a surgery of suppositions. ~M.J.Moroney
>>
>> -Oorspronkelijk bericht-
>> Van: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED]
> Namens Konrad BLOCHER
>> Verzonden: woensdag 6 februari 2008 12:12
>> Aan: r-help@r-project.org
>> Onderwerp: [R] Multivariate Maximum Likelihood Estimation
>>
>> Hi,
>>
>> I am trying to perform Maximum Likelihood estimation of a Multivariate
> model (2 independent variables + intercept) with autocorrelated errors
> of
>> 1st order (ar(1)).
>>
>> Does R have a function for that? I could only find an univariate
> option
> (ar.mle function) and when writing my own I find that it is pretty
> memory-consuming (and sometimes wrong) so there must be a better way.
>>
>> Thanks,
>>
>> KB
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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] help with oop in R - class structure and syntex

2008-02-06 Thread hadley wickham
On Feb 6, 2008 10:13 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> Thanks Gabor. I guess true oo encapsulation is not possible in R.

Before making such a claim, I would encourage you to actually learn
what oo means.  A couple of good readings are:

Structure and Interpretation of Computer Programs.
http://mitpress.mit.edu/sicp/

Concepts, Techniques, and Models of Computer Programming.
http://www.info.ucl.ac.be/~pvr/book.html

Java style (message passing) oo is not the entirety of oo-based programming!

Hadley

-- 
http://had.co.nz/

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


[R] ci.pd() (Epi) and Newcombe method

2008-02-06 Thread Ted Harding
Greetings!

I suspect that there is an error in the code for the
function ci.pd() in the Epi package.

This function is for computing confidence intervals
for a difference of proportions between two independent
groups of 0/1 responses, and implements the Newcombe
("Nc") method and the Agrasti-Caffo "AC" method.

I think there is an error in the computation for the
Newcombe method.

Specifically, where the code says

# Put the data in the right form
  x1 <- aa; n1 <- aa+cc; p1 <- x1/n1
  x2 <- bb; n2 <- bb+dd; p2 <- x2/n2
  pd <- x1/n1 - x2/n2
  z <- qnorm(1-alpha/2)
  zz <- qchisq(1-alpha/2,1)

I believe zz (which is used in the Newcombe method,
while z is used for Agresti-Caffo) should be

  zz <- qchisq(1-alpha,1)

(though I think they could just as well have written
xx <- z^2), since z^2 seems to be what they want for the
Newcombe method, and qnorm(1-alpha/2)^2 = qchisq(1-alpha,1)

This was noticed as a result of comparing results from
ci.pd() with results from the RDCI module in STATA (run
by someone else, since I don't have access to STATA).
With the original code in ci.pd(), the Newcombe results
were different between ci.pd() and STATA. After making
the above change, they agreed.

I would be most grateful if anyone who has also been
down this path could confirm whether I am correct.
Since I don't have access to Newcombe's original
paper at the moment either, I am not able to check
the code against his own description of the method!

With thanks,
Ted.


E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 06-Feb-08   Time: 16:58:19
-- XFMail --

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] adding original axis for x and y to a plot

2008-02-06 Thread Lassana TOURE
Hi.

I have a AFTD plot. I want to add a original axis for x and y.

I did "help(plot)". But I did not see an option for that.

I need help.


[[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] ci.pd() (Epi) and Newcombe method

2008-02-06 Thread Peter Dalgaard
(Ted Harding) wrote:
> Greetings!
>
> I suspect that there is an error in the code for the
> function ci.pd() in the Epi package.
>   
Any particular reason not to include the maintainer among the 
recipients?? Not sure whether Bendix is subscribed.

-p

> This function is for computing confidence intervals
> for a difference of proportions between two independent
> groups of 0/1 responses, and implements the Newcombe
> ("Nc") method and the Agrasti-Caffo "AC" method.
>
> I think there is an error in the computation for the
> Newcombe method.
>
> Specifically, where the code says
>
> # Put the data in the right form
>   x1 <- aa; n1 <- aa+cc; p1 <- x1/n1
>   x2 <- bb; n2 <- bb+dd; p2 <- x2/n2
>   pd <- x1/n1 - x2/n2
>   z <- qnorm(1-alpha/2)
>   zz <- qchisq(1-alpha/2,1)
>
> I believe zz (which is used in the Newcombe method,
> while z is used for Agresti-Caffo) should be
>
>   zz <- qchisq(1-alpha,1)
>
> (though I think they could just as well have written
> xx <- z^2), since z^2 seems to be what they want for the
> Newcombe method, and qnorm(1-alpha/2)^2 = qchisq(1-alpha,1)
>
> This was noticed as a result of comparing results from
> ci.pd() with results from the RDCI module in STATA (run
> by someone else, since I don't have access to STATA).
> With the original code in ci.pd(), the Newcombe results
> were different between ci.pd() and STATA. After making
> the above change, they agreed.
>
> I would be most grateful if anyone who has also been
> down this path could confirm whether I am correct.
> Since I don't have access to Newcombe's original
> paper at the moment either, I am not able to check
> the code against his own description of the method!
>
> With thanks,
> Ted.
>
> 
> E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
> Fax-to-email: +44 (0)870 094 0861
> Date: 06-Feb-08   Time: 16:58:19
> -- XFMail --
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>   


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

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


Re: [R] help with oop in R - class structure and syntex

2008-02-06 Thread tom soyer
Thanks Hardley. I see what you mean. You are right, I am not an expert in
oop AND I don't really know how R oo works, so certainly I shouldn't be
making any sweeping statement. I was just thinking about the issue of local
vs. global, i.e. changes intended for the local environment shouldn't affect
the global environment. That was all I meant.


On 2/6/08, hadley wickham <[EMAIL PROTECTED]> wrote:
>
> On Feb 6, 2008 10:13 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > Thanks Gabor. I guess true oo encapsulation is not possible in R.
>
> Before making such a claim, I would encourage you to actually learn
> what oo means.  A couple of good readings are:
>
> Structure and Interpretation of Computer Programs.
> http://mitpress.mit.edu/sicp/
>
> Concepts, Techniques, and Models of Computer Programming.
> http://www.info.ucl.ac.be/~pvr/book.html
>
> Java style (message passing) oo is not the entirety of oo-based
> programming!
>
> Hadley
>
> --
> http://had.co.nz/
>



-- 
Tom

[[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] Maximum number of variables allowed in a multiple linearregression model

2008-02-06 Thread Tony Plate
Bert Gunter wrote:
> I strongly suggest you collaborate with a local statistician. I can think of
> no circumstance where multiple regression on "hundreds of thousands of
> variables" is anything more than a fancy random number generator. 

That sounds like a challenge!  What is the largest regression problem (in 
terms of numbers of variables) that people have encountered where it made 
sense to do some sort of linear regression (and gave useful results)? 
(Including multilevel and Bayesian techniques.)

However, the original poster did say "hundreds to thousands", which is 
smaller than "hundreds of thousands".  When I try a regression problem with 
3,000 coefficients in R running under Windows XP 64 bit with 8Gb of memory 
on the machine and the /3Gb option active (i.e., R can get up to 3Gb), R 
2.6.1 runs out of memory (apparently trying to duplicate the model matrix):

R version 2.6.1 (2007-11-26)
Copyright (C) 2007 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

 > m <- 3000
 > n <- m * 10
 > x <- matrix(rnorm(n*m), ncol=m, nrow=n, 
dimnames=list(paste("C",1:n,sep=""), paste("X",1:m,sep="")))
 > dim(x)
[1] 3  3000
 > k <- sample(m, 10)
 > y <- rowSums(x[,k]) + 10 * rnorm(n)
 > fit <- lm.fit(y=y, x=x)
Error: cannot allocate vector of size 686.6 Mb
 > object.size(x)/2^20
[1] 687.7787
 > memory.size()
[1] -2022.552
 >
and the Windows process monitor shows the peak memory usage for Rgui.exe at 
2,137,923K.  But in a 64 bit version of R, I would be surprised if it was 
not possible to run this (given sufficient memory).

However, R easily handles a slightly smaller problem:
 > m <- 1000 # of variables
 > n <- m * 10 # of rows
 > k <- sample(m, 10)
 > x <- matrix(rnorm(n*m), ncol=m, nrow=n, 
dimnames=list(paste("C",1:n,sep=""), paste("X",1:m,sep="")))
 > y <- rowSums(x[,k]) + 10 * rnorm(n)
 > fit <- lm.fit(y=y, x=x)
 > # distribution of coefs that should be one vs zero
 > round(rbind(one=quantile(fit$coefficients[k]), 
zero=quantile(fit$coefficients[-k])), digits=2)
 0%   25%   50%  75% 100%
one   0.94  0.98  1.04 1.10 1.18
zero -0.30 -0.08 -0.01 0.06 0.29
 >

To echo Bert Gunter's cautions, one must be careful doing ordinary linear 
regression with large numbers of coefficients.  It does seem a little 
unlikely that there is sufficient data to get useful estimates of three 
thousand coefficients using linear regression in data managed in Excel 
(though I guess it could be possible using Excel 12.0, which can handle up 
to 1 million rows - recent versions prior to 2008 could handle on 64K rows 
- see http://en.wikipedia.org/wiki/Microsoft_Excel#Versions ).  So, the 
suggestion to consult a local statistician is good advice - there may be 
other more suitable approaches, and if some form of linear regression is an 
appropriate approach, there are things to do to gain confidence that the 
results of the linear regression convey useful information.

-- Tony Plate

> 
> -- Bert Gunter
> Genentech Nonclinical Statistics
> 
> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
> Behalf Of Michelle Chu
> Sent: Tuesday, February 05, 2008 9:00 AM
> To: R-help@r-project.org
> Subject: [R] Maximum number of variables allowed in a multiple
> linearregression model
> 
> Hi,
> 
> I appreciate it if someone can confirm the maximum number of variables
> allowed in a multiple linear regression model.  Currently, I am looking for
> a software with the capacity of handling approximately 3,000 variables.  I
> am using Excel to process the results.  Any information for processing a
> matrix from Excel with hundreds to thousands of variables will helpful.
> 
> Best Regards,
> Michelle
> 
>   [[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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ci.pd() (Epi) and Newcombe method

2008-02-06 Thread Ted Harding
On 06-Feb-08 17:09:59, Peter Dalgaard wrote:
> (Ted Harding) wrote:
>> Greetings!
>>
>> I suspect that there is an error in the code for the
>> function ci.pd() in the Epi package.
>>   
> Any particular reason not to include the maintainer among the 
> recipients??

Oversight, and haste! Thank you for copying your reply
(with my original message) to him. And apologies to Bendix.

Ted.

> Not sure whether Bendix is subscribed.
> 
> -p
> 
>> This function is for computing confidence intervals
>> for a difference of proportions between two independent
>> groups of 0/1 responses, and implements the Newcombe
>> ("Nc") method and the Agrasti-Caffo "AC" method.
>>
>> I think there is an error in the computation for the
>> Newcombe method.
>>
>> Specifically, where the code says
>>
>> # Put the data in the right form
>>   x1 <- aa; n1 <- aa+cc; p1 <- x1/n1
>>   x2 <- bb; n2 <- bb+dd; p2 <- x2/n2
>>   pd <- x1/n1 - x2/n2
>>   z <- qnorm(1-alpha/2)
>>   zz <- qchisq(1-alpha/2,1)
>>
>> I believe zz (which is used in the Newcombe method,
>> while z is used for Agresti-Caffo) should be
>>
>>   zz <- qchisq(1-alpha,1)
>>
>> (though I think they could just as well have written
>> xx <- z^2), since z^2 seems to be what they want for the
>> Newcombe method, and qnorm(1-alpha/2)^2 = qchisq(1-alpha,1)
>>
>> This was noticed as a result of comparing results from
>> ci.pd() with results from the RDCI module in STATA (run
>> by someone else, since I don't have access to STATA).
>> With the original code in ci.pd(), the Newcombe results
>> were different between ci.pd() and STATA. After making
>> the above change, they agreed.
>>
>> I would be most grateful if anyone who has also been
>> down this path could confirm whether I am correct.
>> Since I don't have access to Newcombe's original
>> paper at the moment either, I am not able to check
>> the code against his own description of the method!
>>
>> With thanks,
>> Ted.
>>
>> 
>> E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
>> Fax-to-email: +44 (0)870 094 0861
>> Date: 06-Feb-08   Time: 16:58:19
>> -- XFMail --
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>   
> 
> 
> -- 
>O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
>   c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
>  (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45)
> 35327918
> ~~ - ([EMAIL PROTECTED])  FAX: (+45)
> 35327907
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 06-Feb-08   Time: 17:29:34
-- XFMail --

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


Re: [R] inserting text lines in a dat frame

2008-02-06 Thread joseph
Hi Jim

yes, this exactly want I want. However,  I need one more step to get
rid of the first column so that the final file looks like this:
browser 
position 
chr1:1-1
browser 
hide 
all
track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full
variableStep 
chrom=chr1 
span=1
11255  55
11320  
29
 11400  
45
track 
type=wiggle_0 
name=sample 
description=chr2_sample 
visibility=full
variableStep 
chrom=chr2 
span=1
21680 
35
21750 
84
21820 
29
track 
type=wiggle_0 
name=sample 
description=chr3_sample 
visibility=full
variableStep 
chrom=chr3 
span=1
32100 
29
32170 
45
32240 
45
 Thank you very much
- Original Message 
From: jim holtman <[EMAIL PROTECTED]>
To: joseph <[EMAIL PROTECTED]>
Cc: r-help@r-project.org
Sent: Wednesday, February 6, 2008 2:37:38 AM
Subject: Re: inserting text lines in a dat frame


Try 
this 
and 
see 
if 
it 
is 
what 
you 
want:

x 
<- 
read.table(textConnection("  
  
 
V1  
  
V2 
V3
1 
chr1 
11255 
55
2 
chr1 
11320 
29
3 
chr1 
11400 
45
4 
chr2 
21680 
35
5 
chr2 
21750 
84
6 
chr2 
21820 
29
7 
chr2 
31890 
46
8 
chr3 
32100 
29
9 
chr3 
52380 
29
10 
chr3 
66450 
46" 
), 
header=TRUE)
cat("browser 
position 
chr1:1-1\nrowser 
hide 
all\n", 
file='tempxx.txt')
result 
<- 
lapply(split(x, 
x$V1), 
function(.chro){
  
  
cat(sprintf("track 
type=wiggle_0 
name=sample 
description=%s_sample
visibility=full\nvariableStep 
chrom=%s 
span=1\n",
  
  
  
  
as.character(.chro$V1[1]), 
as.character(.chro$V1[1])),
file="tempxx.txt", 
append=TRUE)
  
  
write.table(.chro, 
sep="\t", 
file="tempxx.txt", 
append=TRUE,
col.names=FALSE, 
row.names=FALSE)
})



On 
Feb 
5, 
2008 
11:22 
PM, 
joseph 
<[EMAIL PROTECTED]> 
wrote:
>
>
>
>
>
> 
Hi 
Jim
>  
I 
am 
trying 
to 
prepare 
a 
bed 
file 
to 
load 
as 
accustom 
track 
on 
the 
UCSC
> 
genome 
browser.
> 
I 
have 
a 
data 
frame 
that 
looks 
like 
the 
one 
below.
> 
> 
x
>  
  
  
V1  
  
V2 
V3
> 
1 
chr1 
11255 
55
> 
2 
chr1 
11320 
29
> 
3 
chr1 
11400 
45
> 
4 
chr2 
21680 
35
> 
5 
chr2 
21750 
84
> 
6 
chr2 
21820 
29
> 
7 
chr2 
31890 
46
> 
8 
chr3 
32100 
29
> 
9 
chr3 
52380
>  
29
> 
10 
chr3 
66450 
46
> 
I 
would 
like 
to 
insert 
the 
following 
4 
lines 
at 
the 
beginning:
> 
browser 
position 
chr1:1-1
> 
browser 
hide 
all
> 
track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full
> 
variableStep 
chrom=chr1 
span=1
> 
and 
then 
insert 
2 
lines 
before 
each 
chromosome:
> 
track 
type=wiggle_0 
name=sample 
description=chr2_sample 
visibility=full
> 
vriableStep 
chrom=chr2 
span=1
> 
The 
final 
result 
should 
be 
tab 
delimited 
file 
that 
looks 
like 
this:
> 
browser 
position 
chr1:1-1
> 
browser 
hide 
all
> 
track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full
> 
variableStep 
chrom=chr1 
span=1
> 
chr1 
11255 
55
> 
chr1 
11320 
29
> 
chr1 
11400 
45
> 
track 
type=wiggle_0 
name=sample 
description=chr2_sample 
visibility=full
> 
variableStep 
chrom=chr2 
span=1
> 
chr2 
21680 
35
> 
chr2 
21750 
84
> 
chr2 
21820 
29
> 
track 
type=wiggle_0 
name=sample 
description=chr3_sample 
visibility=full
> 
variableStep 
chrom=chr3
>  
span=1
> 
chr3 
32100 
29
> 
chr3 
32170 
45
> 
chr3 
32240 
45
> 
Any 
kind 
of 
help 
or 
guidance 
will 
be 
much 
appreciated.
> 
Joseph
>
> 

> 
Be 
a 
better 
friend, 
newshound, 
and 
know-it-all 
with 

Mobile. 
Try 
it
> 
now.



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

What 
is 
the 
problem 
you 
are 
trying 
to 
solve?






  

Be a better friend, newshound, and 


[[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] kruskal's MONANOVA algorithm

2008-02-06 Thread Jaime Brugueras
I am trying to obtain the Kruskal (1964) secondary least-squares monotonic
transformation of a rank variable given 4 categorical variables in order to
obtain optimal transformation for regression. The academic problem assigned
is to compare R, SPSS (Conjoint Analysis), and SAS' proc transreg in speed
and accuracy. Currently, SAS and SPSS are giving similar results, but R's
are quite different. There is something I am misunderstanding about
acepackage and/or isoMDS.

The data looks like this:

Brand  Price   Life Hazard Rank
1  Goodstone $69.99 60,000Yes3
2  Goodstone $69.99 70,000 No2
...
7 Pirogi $69.99 50,000 No7
8 Pirogi $69.99 70,000 No1
9 Pirogi $74.99 50,000Yes8

The ace and avas functions transform the y values into very small values of
rank, like this:
$ty
 [1] -1.3552125 -1.6732919  0.8859707

and hence the estimates are quite different.

The R-squared is .93 while SAS and SPSS give .99.

The isoMDS from MASS package gives weird results when i choose k=4.

Here is my acepackage code and isoMDS function:

X <- cbind(Brand, Price, Life, Hazard)  # independent variables
Y <-  Rank # response variable
cate <- as.vector(c(1,2,3,4))   # categorical variables(columns) in
X

mycon <- avas(x=X, y=Y, cat=cate)

mymatrix <- as.matrix((X)
row.names(mymatrix) <- Rank

Any help is well appreciated.

Thanks.



myc <- isoMDS(dist(mymatrix), k=?)

[[alternative HTML version deleted]]

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


Re: [R] help with oop in R - class structure and syntex

2008-02-06 Thread Gabor Grothendieck
If you believe that certain changes intended only to affect the local
environment nevertheless affect the global environment please give a
code example.

On Feb 6, 2008 12:11 PM, tom soyer <[EMAIL PROTECTED]> wrote:
> Thanks Hardley. I see what you mean. You are right, I am not an expert in
> oop AND I don't really know how R oo works, so certainly I shouldn't be
> making any sweeping statement. I was just thinking about the issue of local
> vs. global, i.e. changes intended for the local environment shouldn't affect
> the global environment. That was all I meant.
>
>
> On 2/6/08, hadley wickham <[EMAIL PROTECTED]> wrote:
> >
> > On Feb 6, 2008 10:13 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > Thanks Gabor. I guess true oo encapsulation is not possible in R.
> >
> > Before making such a claim, I would encourage you to actually learn
> > what oo means.  A couple of good readings are:
> >
> > Structure and Interpretation of Computer Programs.
> > http://mitpress.mit.edu/sicp/
> >
> > Concepts, Techniques, and Models of Computer Programming.
> > http://www.info.ucl.ac.be/~pvr/book.html
> >
> > Java style (message passing) oo is not the entirety of oo-based
> > programming!
> >
> > Hadley
> >
> > --
> > http://had.co.nz/
> >
>
>
>
> --
> Tom
>
>[[alternative HTML version deleted]]
>
> __
>
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] inserting text lines in a dat frame

2008-02-06 Thread joseph
thanks Richie. It woks perfectly.

- Original Message 
From: "[EMAIL PROTECTED]" <[EMAIL PROTECTED]>
To: joseph <[EMAIL PROTECTED]>
Cc: r-help@r-project.org; [EMAIL PROTECTED]
Sent: Wednesday, February 6, 2008 2:22:16 AM
Subject: Re: [R] inserting text lines in a dat  frame


>  
I 
am 
trying 
to 
prepare 
a 
bed 
file 
to 
load 
as 
accustom 
track 
on 
the 
> 
UCSC 
genome 
browser.
> 
I 
have 
a 
data 
frame 
that 
looks 
like 
the 
one 
below.
> 
> 
x
>  
  
  
V1  
  
V2 
V3
> 
1 
chr1 
11255 
55
> 
2 
chr1 
11320 
29
> 
3 
chr1 
11400 
45
> 
4 
chr2 
21680 
35
> 
5 
chr2 
21750 
84
> 
6 
chr2 
21820 
29
> 
7 
chr2 
31890 
46
> 
8 
chr3 
32100 
29
> 
9 
chr3 
52380 
29
> 
10 
chr3 
66450 
46
> 
I 
would 
like 
to 
insert 
the 
following 
4 
lines 
at 
the 
beginning:
> 
browser 
position 
chr1:1-1
> 
browser 
hide 
all
> 
track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full
> 
variableStep 
chrom=chr1 
span=1
> 
and 
then 
insert 
2 
lines 
before 
each 
chromosome:
> 
track 
type=wiggle_0 
name=sample 
description=chr2_sample 
visibility=full
> 
vriableStep 
chrom=chr2 
span=1
> 
The 
final 
result 
should 
be 
tab 
delimited 
file 
that 
looks 
like 
this:
> 
browser 
position 
chr1:1-1
> 
browser 
hide 
all
> 
track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full
> 
variableStep 
chrom=chr1 
span=1
> 
chr1 
11255 
55
> 
chr1 
11320 
29
> 
chr1 
11400 
45
> 
track 
type=wiggle_0 
name=sample 
description=chr2_sample 
visibility=full
> 
variableStep 
chrom=chr2 
span=1
> 
chr2 
21680 
35
> 
chr2 
21750 
84
> 
chr2 
21820 
29
> 
track 
type=wiggle_0 
name=sample 
description=chr3_sample 
visibility=full
> 
variableStep 
chrom=chr3 
span=1
> 
chr3 
32100 
29
> 
chr3 
32170 
45
> 
chr3 
32240 
45

#First 
write 
your 
preamble 
text 
to 
a 
file
preamble 
<- 
c("browser 
position 
chr1:1-1",
  
  
  
  
  
  
  
"browser 
hide 
all",
  
  
  
  
  
  
  
"track 
type=wiggle_0 
name=sample 
description=chr1_sample 
visibility=full",
  
  
  
  
  
  
  
"variableStep 
chrom=chr1 
span=1")
write(preamble, 
"myfile.txt")

#Now 
split 
your 
data 
frame 
up 
by 
the 
values 
in 
V1
x 
<- 
data.frame(V1=rep(c("chr1", 
"chr2", 
"chr3"),times=c(3,4,3)), 
V2=c(11255,11320,11400,21680,21750,21820,21890,32100,52380,66450),V3=c(55,29,45,35,84,29,46,29,29,46))
spx 
<- 
split(x, 
x$V1)

#Create 
lines 
of 
text 
to 
go 
before 
each 
piece 
of 
data 
frame
lV1 
<- 
levels(x$V1)
maintext 
<- 
paste("track 
type=wiggle_0 
name=sample 
description=", 
lV1,
  
  
  
  
  
  
  
 
"_sample 
visibility=full\nvariableStep 
chrom=", 
lV1,
  
  
  
  
  
  
  
 
" 
span=1", 
sep="")

#Use 
a 
loop 
to 
write 
the 
pieces 
to 
the 
file
for(i 
in 
1: 
nlevels(x$V1))
  
{
  
  
write(maintext[i], 
"myfile.txt", 
append=TRUE)
  
  
write.table(spx[[i]], 
"myfile.txt", 
append=TRUE, 
sep="\t", 
quote=FALSE,
  
  
  
  
  
  
  
  
col.names=FALSE, 
row.names=FALSE)
  
}

Regards,
Richie.

Mathematical 
Sciences 
Unit
HSL



ATTENTION:

This 
message 
contains 
privileged 
and 
confidential 
information 
intended
for 
the 
addressee(s) 
only. 
If 
this 
message 
was 
sent 
to 
you 
in 
error,
you 
must 
not 
disseminate, 
copy 
or 
take 
any 
action 
in 
reliance 
on 
it 
and
we 
request 
that 
you 
notify 
the 
sender 
immediately 
by 
return 
email.

Opinions 
expressed 
in 
this 
message 
and 
any 
attachments 
are 
not
necessarily 
those 
held 
by 
the 
Health 
and 
Safety 
Laboratory 
or 
any 
person
connected 
with 
the 
organisation, 
save 
those 
by 
whom 
the 
opinions 
were
expressed.

Please 
note 
that 
any 
messages 
sent 
or 
received 
by 
the 
Health 
and 
Safety
Laboratory 
email 
system 
may 
be 
monitored 
and 
stored 
in 
an 
information
retrieval 
system.



This 
e-mail 
message 
has 
been 
scanned 
for 
Viruses 
and 
Content 
and 
cleared 
by 
NetIQ 
MailMarshal







  

Be a better friend, newshound, and 


[[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] box.Cox.powers() warning

2008-02-06 Thread cmhcordei
Dear Rlist,

Using an example in box.cox.powers() help, I have the following warning message.

example:
library(car)
>attach(Prestige)

> box.cox.powers(income)
Box-Cox Transformation to Normality 

 Est.Power Std.Err. Wald(Power=0) Wald(Power=1)
0.1793   0.11081.6179   -7.4062

L.R. test, power = 0:  2.7103   df = 1   p = 0.0997
L.R. test, power = 1:  47.261   df = 1   p = 0 
Warning message:
In optimize(f = function(lambda) univ.neg.kernel.logL(x = X[, j],  :
  NA/Inf replaced by maximum positive value

I do not know the reason for such warning message. 
Can someone help me on this matter?

Best regards, 
Clara Cordeiro

> sessionInfo()
R version 2.6.1 (2007-11-26) 
i386-pc-mingw32 

locale:
LC_COLLATE=Portuguese_Portugal.1252;LC_CTYPE=Portuguese_Portugal.1252;LC_MONETARY=Portuguese_Portugal.1252;LC_NUMERIC=C;LC_TIME=Portuguese_Portugal.1252

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

other attached packages:
[1] nnet_7.2-38car_1.2-7  MASS_7.2-38FitAR_1.0  leaps_2.7 
[6] lattice_0.17-2

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 oop in R - class structure and syntex

2008-02-06 Thread tom soyer
Gabor, maybe I am not understanding it right. I was thinking for example
something like how a call to plot actually calls plot.zoo when zoo is
loaded. And a specific call to plot.zoo would not work, etc. One would think
that plot and plot.zoo are separate and that one could call the parent as
well as the offsprings independently. But if this kind of setup is by
design, than it's fine. Does that make sense?


On 2/6/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
>
> If you believe that certain changes intended only to affect the local
> environment nevertheless affect the global environment please give a
> code example.
>
> On Feb 6, 2008 12:11 PM, tom soyer <[EMAIL PROTECTED]> wrote:
> > Thanks Hardley. I see what you mean. You are right, I am not an expert
> in
> > oop AND I don't really know how R oo works, so certainly I shouldn't be
> > making any sweeping statement. I was just thinking about the issue of
> local
> > vs. global, i.e. changes intended for the local environment shouldn't
> affect
> > the global environment. That was all I meant.
> >
> >
> > On 2/6/08, hadley wickham <[EMAIL PROTECTED]> wrote:
> > >
> > > On Feb 6, 2008 10:13 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > > Thanks Gabor. I guess true oo encapsulation is not possible in R.
> > >
> > > Before making such a claim, I would encourage you to actually learn
> > > what oo means.  A couple of good readings are:
> > >
> > > Structure and Interpretation of Computer Programs.
> > > http://mitpress.mit.edu/sicp/
> > >
> > > Concepts, Techniques, and Models of Computer Programming.
> > > http://www.info.ucl.ac.be/~pvr/book.html
> > >
> > > Java style (message passing) oo is not the entirety of oo-based
> > > programming!
> > >
> > > Hadley
> > >
> > > --
> > > http://had.co.nz/
> > >
> >
> >
> >
> > --
> > Tom
> >
> >[[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.
> >
>



-- 
Tom

[[alternative HTML version deleted]]

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


Re: [R] help with oop in R - class structure and syntex

2008-02-06 Thread Gabor Grothendieck
zoo has a NAMESPACE which restricts access to the
internals of the package except for explicitly exported items.
zoo exports plot.zoo as an S3 method so that the plot generic
can find it.  It is also possible to simply export a function
not specifically as an S3 method in which case it will be
visible more generally.  For example see seq.dates in
chron.  seq.dates is a seq method but can be called directly
too since it was exported directly.

# plot.zoo exported as S3 method
library(zoo)
z <- zoo(1:10)
plot(z) # plot.zoo has been exported as S3 method
zoo:::plot.zoo(z) # same output - not recommended
# plot.zoo(z) # won't work

# seq.dates exported
library(chron)
x <- chron(1:10)
seq(x[1], x[10])
seq.dates(x[1], x[10]) # direct access - same output

If an object is not exported at all then it cannot be accessed
(other than via ::: notation).

The other type of package is one with no NAMESPACE.
In that case everything in the package is visible.

On Feb 6, 2008 1:16 PM, tom soyer <[EMAIL PROTECTED]> wrote:
> Gabor, maybe I am not understanding it right. I was thinking for example
> something like how a call to plot actually calls plot.zoo when zoo is
> loaded. And a specific call to plot.zoo would not work, etc. One would think
> that plot and plot.zoo are separate and that one could call the parent as
> well as the offsprings independently. But if this kind of setup is by
> design, than it's fine. Does that make sense?
>
>
> On 2/6/08, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
>
> > If you believe that certain changes intended only to affect the local
> > environment nevertheless affect the global environment please give a
> > code example.
> >
> > On Feb 6, 2008 12:11 PM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > Thanks Hardley. I see what you mean. You are right, I am not an expert
> in
> > > oop AND I don't really know how R oo works, so certainly I shouldn't be
> > > making any sweeping statement. I was just thinking about the issue of
> local
> > > vs. global, i.e. changes intended for the local environment shouldn't
> affect
> > > the global environment. That was all I meant.
> > >
> > >
> > > On 2/6/08, hadley wickham <[EMAIL PROTECTED]> wrote:
> > > >
> > > > On Feb 6, 2008 10:13 AM, tom soyer <[EMAIL PROTECTED]> wrote:
> > > > > Thanks Gabor. I guess true oo encapsulation is not possible in R.
> > > >
> > > > Before making such a claim, I would encourage you to actually learn
> > > > what oo means.  A couple of good readings are:
> > > >
> > > > Structure and Interpretation of Computer Programs.
> > > > http://mitpress.mit.edu/sicp/
> > > >
> > > > Concepts, Techniques, and Models of Computer Programming.
> > > > http://www.info.ucl.ac.be/~pvr/book.html
> > > >
> > > > Java style (message passing) oo is not the entirety of oo-based
> > > > programming!
> > > >
> > > > Hadley
> > > >
> > > > --
> > > > http://had.co.nz/
> > > >
> > >
> > >
> > >
> > > --
> > > Tom
> > >
> > >[[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.
> > >
> >
>
>
>
> --
> Tom

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

2008-02-06 Thread Tim Hesterberg
>   I want to generate different samples using the
>followindg code:
>
>g<-sample(LETTERS[1:2], 24, replace=T)
>
>   How can I specify that I need 12 "A"s and 12 "B"s?

I introduced the concept of "sampling with minimal replacement" into the
S-PLUS version of sample to handle things like this:
sample(LETTERS[1:2], 24, minimal = T)

This is very useful in variance reduction applications, to approximately
stratify but with introducing bias.  I'd like to see this in R.


I'll raise a related issue - sampling with unequal probabilities,
without replacement.  R does the wrong thing, in my opinion:

> values <- sapply(1:1000, function(i) sample(1:3, size=2, prob = c(.5, .25, 
> .25)))
> table(values)
values
  1   2   3 
834 574 592 

The selection probabilities are not proportional to the specified
probabilities.  

In contrast, in S-PLUS:
> values <- sapply(1:1000, function(i) sample(1:3, size=2, prob = c(.5, .25, 
> .25)))
> table(values)
1   2   3 
 1000 501 499

You can specify minimal = FALSE to get the same behavior as R:
> values <- sapply(1:1000, function(i) sample(1:3, size=2, prob = c(.5, .25, 
> .25), minimal = F))
> table(values)
   1   2   3 
 844 592 564

There is a reason this is associated with the concept of sampling with
minimal replacement.  Consider for example:
sample(1:4, size = 3, prob = 1:4/10)
The expected frequencies of (1,2,3,4) should be proportional
to size*prob = c(.3,.6,.9,1.2).  That isn't possible when sampling
without replacement.  Sampling with minimal replacement allows this;
observation 4 is included in every sample, and is included twice in
20% of the samples.

Tim Hesterberg

Disclaimer - these are my opinions, not those of my employer.


| Tim Hesterberg   Senior Research Scientist   |
| [EMAIL PROTECTED]  Insightful Corp.|
| (206)802-23191700 Westlake Ave. N, Suite 500 |
| (206)283-8691 (fax)  Seattle, WA 98109-3044, U.S.A.  |
|  www.insightful.com/Hesterberg   |

I'll teach short courses:
Advanced Programming in S-PLUS: San Antonio TX, March 26-27, 2008.
Bootstrap Methods and Permutation Tests: San Antonio, March 28, 2008.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] counting row repetitions without loop

2008-02-06 Thread Douglas Bates
On Feb 6, 2008 8:08 AM, Waterman, DG (David)
<[EMAIL PROTECTED]> wrote:
> Hi,

> I have a data frame consisting of coordinates on a 10*10 grid, i.e.

> > example
> x  y
> 1   4  5
> 2   6  7
> 3   6  6
> 4   7  5
> 5   5  7
> 6   6  7
> 7   4  5
> 8   6  7
> 9   7  6
> 10  5  6

> What I would like to do is return an 10*10 matrix consisting of counts
> at each position, so in the above example I would have a matrix where,
> for example, cell [4,5] contains 2 and [6,7] contains 3. At the moment I
> have implemented this using a for loop over the rows of the data frame,
> however the data frames I want to process are very long so the loop
> takes many minutes to complete. Can I do this in a more efficient way?

What you are describing is essentially a cross-tabulation so you could use

> examp
   x y
1  4 5
2  6 7
3  6 6
4  7 5
5  5 7
6  6 7
7  4 5
8  6 7
9  7 6
10 5 6
> xtabs(~ x + y, examp)
   y
x   5 6 7
  4 2 0 0
  5 0 1 1
  6 0 1 3
  7 1 1 0

This omits the rows and columns which are completely empty but you can
work around that.

If you have a very large collection of such pairs to summarize you
could consider the version of xtabs in the Matrix package that allows
for the argument sparse = TRUE.  That uses conversion of the "triplet"
form of a sparse matrix to the compressed column for to do the
counting.

If you want to do this without converting the integers in 'x' and 'y'
to factors you can use a distinctly unobvious function like

library(Matrix)
sparsetab <- function(x, y)
{
x <- as.integer(x)
y <- as.integer(y)
stopifnot(length(x) == length(y))
lx <- length(x)
mx <- max(x)
my <- max(y)
as(new("dgTMatrix", i = x - 1L, j = y - 1L,
   x = rep(1, length(x)), Dim = c(mx, my),
   Dimnames = list(1:mx,1:my)), "dgCMatrix")
}

which produces

> with(examp, sparsetab(x, y))
7 x 7 sparse Matrix of class "dgCMatrix"
  1 2 3 4 5 6 7
1 . . . . . . .
2 . . . . . . .
3 . . . . . . .
4 . . . . 2 . .
5 . . . . . 1 1
6 . . . . . 1 3
7 . . . . 1 1 .

One reason to use such a function instead of xtabs is because xtabs
will convert 'x' and 'y' to factors and the default ordering of the
levels is lexicographic so '11' occurs before '2'.  Again, you can get
around that but the function shown above is more direct and should be
fast enough for most any application.

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


Re: [R] inserting text lines in a dat frame

2008-02-06 Thread jim holtman
This should do it for you:

x <- read.table(textConnection("V1V2 V3
1 chr1 11255 55
2 chr1 11320 29
3 chr1 11400 45
4 chr2 21680 35
5 chr2 21750 84
6 chr2 21820 29
7 chr2 31890 46
8 chr3 32100 29
9 chr3 52380 29
10 chr3 66450 46" ), header=TRUE)
outfile <- '/tempxx.txt'
cat("browser position chr1:1-1\nrowser hide all\n", file=outfile)
result <- lapply(split(x, x$V1), function(.chro){
cat(sprintf("track type=wiggle_0 name=sample description=%s_sample
visibility=full\nvariableStep chrom=%s span=1\n",
as.character(.chro$V1[1]), as.character(.chro$V1[1])),
file=outfile, append=TRUE)
write.table(.chro[, -1], sep="\t", file=outfile, append=TRUE,
quote=FALSE, col.names=FALSE, row.names=FALSE)
})


On 2/6/08, joseph <[EMAIL PROTECTED]> wrote:
>
> Hi Jim
> yes, this exactly want I want. However,  I need one more step to get rid of
> the first column so that the final file looks like this:
> browser position chr1:1-1
> browser hide all
> track type=wiggle_0 name=sample description=chr1_sample visibility=full
> variableStep chrom=chr1 span=1
> 11255  55
> 11320  29
>  11400  45
> track type=wiggle_0 name=sample description=chr2_sample visibility=full
> variableStep chrom=chr2 span=1
> 21680 35
> 21750 84
> 21820 29
> track type=wiggle_0 name=sample description=chr3_sample visibility=full
> variableStep chrom=chr3 span=1
> 32100 29
> 32170 45
> 32240 45
>  Thank you very much
> - Original Message 
> From: jim holtman <[EMAIL PROTECTED]>
> To: joseph <[EMAIL PROTECTED]>
> Cc: r-help@r-project.org
> Sent: Wednesday, February 6, 2008 2:37:38 AM
> Subject: Re: inserting text lines in a dat frame
>
> Try this and see if it is what you want:
>
> x <- read.table(textConnection("V1V2 V3
> 1 chr1 11255 55
> 2 chr1 11320 29
> 3 chr1 11400 45
> 4 chr2 21680 35
> 5 chr2 21750 84
> 6 chr2 21820 29
> 7 chr2 31890 46
> 8 chr3 32100 29
> 9 chr3 52380 29
> 10 chr3 66450 46" ), header=TRUE)
> cat("browser position chr1:1-1\nrowser hide all\n", file='tempxx.txt')
> result <- lapply(split(x, x$V1), function(.chro){
> cat(sprintf("track type=wiggle_0 name=sample description=%s_sample
> visibility=full\nvariableStep chrom=%s span=1\n",
> as.character(.chro$V1[1]), as.character(.chro$V1[1])),
> file="tempxx.txt", append=TRUE)
> write.table(.chro, sep="\t", file="tempxx.txt", append=TRUE,
> col.names=FALSE, row.names=FALSE)
> })
>
>
>
> On Feb 5, 2008 11:22 PM, joseph <[EMAIL PROTECTED]> wrote:
> >
> >
> >
> >
> >
> > Hi Jim
> >  I am trying to prepare a bed file to load as accustom track on the UCSC
> > genome browser.
> > I have a data frame that looks like the one below.
> > > x
> >  V1V2 V3
> > 1 chr1 11255 55
> > 2 chr1 11320 29
> > 3 chr1 11400 45
> > 4 chr2 21680 35
> > 5 chr2 21750 84
> > 6 chr2 21820 29
> > 7 chr2 31890 46
> > 8 chr3 32100 29
> > 9 chr3 52380
> >  29
> > 10 chr3 66450 46
> > I would like to insert the following 4 lines at the beginning:
> > browser position chr1:1-1
> > browser hide all
> > track type=wiggle_0 name=sample description=chr1_sample visibility=full
> > variableStep chrom=chr1 span=1
> > and then insert 2 lines before each chromosome:
> > track type=wiggle_0 name=sample description=chr2_sample visibility=full
> > vriableStep chrom=chr2 span=1
> > The final result should be tab delimited file that looks like this:
> > browser position chr1:1-1
> > browser hide all
> > track type=wiggle_0 name=sample description=chr1_sample visibility=full
> > variableStep chrom=chr1 span=1
> > chr1 11255 55
> > chr1 11320 29
> > chr1 11400 45
> > track type=wiggle_0 name=sample description=chr2_sample visibility=full
> > variableStep chrom=chr2 span=1
> > chr2 21680 35
> > chr2 21750 84
> > chr2 21820 29
> > track type=wiggle_0 name=sample description=chr3_sample visibility=full
> > variableStep chrom=chr3
> >  span=1
> > chr3 32100 29
> > chr3 32170 45
> > chr3 32240 45
> > Any kind of help or guidance will be much appreciated.
> > Joseph
> >
> > 
> > Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it
> > now.
>
>
>
> --
> Jim Holtman
> Cincinnati, OH
> +1 513 646 9390
>
> What is the problem you are trying to solve?
>
>
> 
> Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it
> now.


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

What is the problem you are trying to solve?

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


[R] RSF

2008-02-06 Thread Luca Brugnaro











how can i compare a cox model and a rsf?
someone can tell me how to do that in R?
Thanks
[[alternative HTML version deleted]]

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


Re: [R] axis help

2008-02-06 Thread Martin Elff
On Wednesday 06 February 2008 (15:14:28), [EMAIL PROTECTED] wrote:
> Hi, i'm having trouble with my x and y axis. The commands i'm using are
>below. The problem is that the y axis starts at coordinate 0,1 and the x
>axis starts at coordinate 0,0. As far as I know the y axis can't  start
> at 0 (because it's log scaled) ,so I would like to position the x axis at
> 0,1 but don't know how to do this. Could anyone advise?
>axis(1,at=c(0,0.6,1.1,1.6,2.1,2.6,3.1,3.6),labels=c(0,3,4,5,6,7,8,NA))
>axis(2)

Try:

axis(1,at=c(0,0.6,1.1,1.6,2.1,2.6,3.1,3.6),labels=c(0,3,4,5,6,7,8,NA),pos=1)

The 'pos=' argument does the trick.

Best,
Martin

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

2008-02-06 Thread Peter Dalgaard
Tim Hesterberg wrote:
>> values <- sapply(1:1000, function(i) sample(1:3, size=2, prob = c(.5, .25, 
>> .25)))
>> table(values)
>> 
> values
>   1   2   3 
> 834 574 592 
>
> The selection probabilities are not proportional to the specified
> probabilities.  
>
> In contrast, in S-PLUS:
>   
>> values <- sapply(1:1000, function(i) sample(1:3, size=2, prob = c(.5, .25, 
>> .25)))
>> table(values)
>> 
> 1   2   3 
>  1000 501 499
>
>   
But is that the right thing? If you can use prob=c(.6, .2, .2) and get 
1200 - 400 - 400, then I'm not going to play poker with you

The interpretation in R is that you are dealing with "fat cards", i.e. 
card 1 is twice as thick as the other two, so there is 50% chance of 
getting the _first_ card as a 1 and additionally, (.25+.25)*2/3 to get 
the 2nd card as a 1 for a total of .8333. And since the two cards are 
different the expected number of occurrences of card 1 in 1000 samples 
is 833.  What is the interpretation in S-PLUS?

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

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


[R] filling data into objects

2008-02-06 Thread Nair, Murlidharan T
I am trying to generate artificial data for feature selection.  Basically 
trying to generate a total of 1000 features with 100 that are informative and 
rest are uninformative.
Informative.data.class1<-rnorm(100,0.25,1)
Uninformative.data.class1<-rnorm(900,0,1)
Informative.data.class2<-rnorm(100, -0.25,1)
Uninformative.data.class2<-rnorm(900,0,1)

The above will give me one set of data for the two classes. I am interested in 
generating about a 100 set for each class.  What is a neat way to write it in R?
Thanks ../Murli


[[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] Discriminant function analysis

2008-02-06 Thread Tyler Smith
On 2008-02-06, Birgit Lemcke <[EMAIL PROTECTED]> wrote:
>
> I am using R 2.6.1 on a PowerBook G4.
> I would like to perform a discriminant function analysis. I found lda  
> in MASS but as far as I understood, is it only working with  
> explanatory variables of the class factor. 

I think you are mistaken. If you reread the help for lda() you'll see
that the grouping has to be a factor, but the explanatory variables
should be numeric.

> My dataset contains variables of the classes factor and numeric. Is
> there another function that is able to handle this?

The numeric variables are fine. The factor variables may have to be
recoded into dummy binary variables, I'm not sure if lda() will deal
with them properly otherwise.

HTH,

Tyler

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

2008-02-06 Thread John Fox
Dear Clara,

The warning is nothing to worry about. If start values aren't specified,
box.cox.powers() uses optimize() to find start values for each
transformation parameter (in this case, there's only one) prior to using
optim() to maximize the (possibly multivariate) likelihood. In the course of
finding start values, optimize() apparently tried an unreasonable value. In
this case, box.cox.powers(income, start=1) gives you the same result minus
the warning.

Regards,
 John


John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox


> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> project.org] On Behalf Of [EMAIL PROTECTED]
> Sent: February-06-08 1:05 PM
> To: R-help@r-project.org
> Subject: [R] box.Cox.powers() warning
> 
> Dear Rlist,
> 
> Using an example in box.cox.powers() help, I have the following warning
> message.
> 
> example:
> library(car)
> >attach(Prestige)
> 
> > box.cox.powers(income)
> Box-Cox Transformation to Normality
> 
>  Est.Power Std.Err. Wald(Power=0) Wald(Power=1)
> 0.1793   0.11081.6179   -7.4062
> 
> L.R. test, power = 0:  2.7103   df = 1   p = 0.0997
> L.R. test, power = 1:  47.261   df = 1   p = 0
> Warning message:
> In optimize(f = function(lambda) univ.neg.kernel.logL(x = X[, j],  :
>   NA/Inf replaced by maximum positive value
> 
> I do not know the reason for such warning message.
> Can someone help me on this matter?
> 
> Best regards,
> Clara Cordeiro
> 
> > sessionInfo()
> R version 2.6.1 (2007-11-26)
> i386-pc-mingw32
> 
> locale:
> LC_COLLATE=Portuguese_Portugal.1252;LC_CTYPE=Portuguese_Portugal.1252;L
> C_MONETARY=Portuguese_Portugal.1252;LC_NUMERIC=C;LC_TIME=Portuguese_Por
> tugal.1252
> 
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
> 
> other attached packages:
> [1] nnet_7.2-38car_1.2-7  MASS_7.2-38FitAR_1.0
> leaps_2.7
> [6] lattice_0.17-2
> 
> loaded via a namespace (and not attached):
> [1] grid_2.6.1  tools_2.6.1
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] filling data into objects

2008-02-06 Thread Henrique Dallazuanna
Perhaps:

for(i in 1:100)assign(sprintf("Informative.data.class%s", i),
rnorm(100, 0.25,1))
for(i in 1:100)assign(sprintf("Uninformative.data.class%s", i), rnorm(900))

Or working with a list:

Informative <- replicate(100, rnorm(100, 0.25,1))
Uninformative <- replicate(100, rnorm(900))

On 06/02/2008, Nair, Murlidharan T <[EMAIL PROTECTED]> wrote:
> I am trying to generate artificial data for feature selection.  Basically 
> trying to generate a total of 1000 features with 100 that are informative and 
> rest are uninformative.
> Informative.data.class1<-rnorm(100,0.25,1)
> Uninformative.data.class1<-rnorm(900,0,1)
> Informative.data.class2<-rnorm(100, -0.25,1)
> Uninformative.data.class2<-rnorm(900,0,1)
>
> The above will give me one set of data for the two classes. I am interested 
> in generating about a 100 set for each class.  What is a neat way to write it 
> in R?
> Thanks ../Murli
>
>
> [[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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lme predicted value confidence intervals

2008-02-06 Thread Colin
Dear R users,

Does anyone know of a way to obtain approximate 95% confidence intervals 
for predicted values for factor levels of fixed effects from lme?  Our 
goal is to use these intervals to interpret patterns across our 
predicted values for certain factor levels.

Our mixed model has the following form with 7 levels of mtDNA, 2 levels 
of autosome, 2 levels of brood and 2 levels of block,

 > lme(fitness ~ mtDNA*autosome + brood, random = ~1 | block)

We have used the predict.lme function to obtain predicted values, but 
are unsure how to obtain appropriate standard errors on these predicted 
values.

Using predict.lme to predict "fitness" across a subset of our factor 
levels (2 mtDNA, 2 autosome) generates the following output,

autosome   mtDNA brood block predict.fixed predict.block
1   ore ore A A 0.4977047 0.5016255
2   ore simw501 A A 0.4278287 0.4317495
3   ore ore B A 0.5042857 0.5082065
4   ore simw501 B A 0.4344098 0.4383306
5   ore ore A B 0.4977047 0.4937839
6   ore simw501 A B 0.4278287 0.4239079
7   ore ore B B 0.5042857 0.5003649
8   ore simw501 B B 0.4344098 0.4304890
9   aut ore A A 0.5321071 0.5360279
10  aut simw501 A A 0.4866497 0.4905705
11  aut ore B A 0.5386882 0.5426090
12  aut simw501 B A 0.4932308 0.4971516
13  aut ore A B 0.5321071 0.5281863
14  aut simw501 A B 0.4866497 0.4827289
15  aut ore B B 0.5386882 0.5347674
16  aut simw501 B B 0.4932308 0.4893099

We would like to calculate, for example, the appropriate 95% confidence 
intervals for the predicted values of autosome=ore + mtDNA=ore, 
autosome=ore + mtDNA=simw501, etc.

Sincerely,
Kristi Montooth and Colin Meiklejohn
Ecology and Evolutionary Biology
Brown University

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


Re: [R] filling data into objects

2008-02-06 Thread Nair, Murlidharan T
Oh yes, that's what I was looking for. I can then do a
Class1<-rbind(Informative.data, Uninformative.data)

Thanks ../Murli


-Original Message-
From: Henrique Dallazuanna [mailto:[EMAIL PROTECTED]
Sent: Wednesday, February 06, 2008 3:06 PM
To: Nair, Murlidharan T
Cc: [EMAIL PROTECTED]
Subject: Re: [R] filling data into objects

Perhaps:

for(i in 1:100)assign(sprintf("Informative.data.class%s", i),
rnorm(100, 0.25,1))
for(i in 1:100)assign(sprintf("Uninformative.data.class%s", i), rnorm(900))

Or working with a list:

Informative <- replicate(100, rnorm(100, 0.25,1))
Uninformative <- replicate(100, rnorm(900))

On 06/02/2008, Nair, Murlidharan T <[EMAIL PROTECTED]> wrote:
> I am trying to generate artificial data for feature selection.  Basically 
> trying to generate a total of 1000 features with 100 that are informative and 
> rest are uninformative.
> Informative.data.class1<-rnorm(100,0.25,1)
> Uninformative.data.class1<-rnorm(900,0,1)
> Informative.data.class2<-rnorm(100, -0.25,1)
> Uninformative.data.class2<-rnorm(900,0,1)
>
> The above will give me one set of data for the two classes. I am interested 
> in generating about a 100 set for each class.  What is a neat way to write it 
> in R?
> Thanks ../Murli
>
>
> [[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

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

2008-02-06 Thread Tim Hesterberg
>Tim Hesterberg wrote:
>>I'll raise a related issue - sampling with unequal probabilities,
>>without replacement.  R does the wrong thing, in my opinion:
>>...
>Peter Dalgaard wrote:
>But is that the right thing? ...
(See bottom for more of the previous messages.)


First, consider the common case, where size * max(prob) < 1 --
sampling with unequal probabilities without replacement.

Why do people do sampling with unequal probabilities, without
replacement?  A typical application would be sampling with probability
proportional to size, or more generally where the desire is that
selection probabilities match some criterion.

The default S-PLUS algorithm does that.  The selection probabilities
at each of step 1, 2, ..., size are all equal to prob, and the overall
probabilities of selection are size*prob.

The R algorithm (and old S+ algorithm, prior to S+6.2) did not;
at step 1 the selection probabilities are correct, but not in
subsequent steps.  With size = 2 and prob = c(.5, .25, .25),
for example, the selection probabilities at step 2 are (.33, .33, .33).

The R/old S+ algorithm is the result of programming "sampling with
probability proportional to prob without replacement" literally,
sampling greedily at each step, ignoring dependence, without thinking
about that fact that that doesn't give the desired unconditional
probabilities after step 1.  In practice one often needs to know the
actual selection probabilities, and they are difficult to compute
except in toy problems with small size.


Now turn to the less common case -- possible overselection -- sampling
with unequal probabilities with minimal replacement.  One example of
this would be in multistage sampling.  At stage one you select cities
with probability proportional to prob.  If size*max(prob) > 1, then a
city may be selected more than once.  If selected more than once, the
second stage selects that many individuals from the city.


Another example of both cases is an improvement on SIR
(sampling-importance-resampling).  In Bayesian simulation, if
importance sampling is used, each observation ends up with a value and
a weight.  Often one would prefer the output to be a sample without
weights.  In SIR this is done by sampling from the observations with
probability proportional to the weight, <>.  As a
result some observations may be sampled multiple times.

It would be better to do this with minimal replacement.  If
size*max(prob) <= 1 then no observation is sampled multiple times.
If size*max(prob) > 1 then some observations may be sampled
multiple times, but there will tend to be less of this (and more
unique values sampled) than if sampling with replacement.


Here's a final example of the less common case.  This one is more
colorful, but a bit long; at the end I'll relate this to another
improvement on SIR.

"Splitting" and "Russian Roulette" are complementary techniques that
date back to the early days of Monte Carlo simulation, developed for
estimating the probability of neutrons or other particles penetrating
a shield.

The simplest procedure is to generate a bunch of random neutrons, and
follow each on its path through the shield - it flies a random
distance until it hits an atom, there is a certain probability of
absorption, otherwise it is reflected in a random direction, etc.  In
the end you count the proportion of neutrons that penetrate.  The
problem is that the penetration probability is so small that a
you have to start simulating with an enormous number of neutrons, to
accurately estimate the small probability.

The next improvement is to bias the sampling, and assign each
neutron a weight.  The weight W starts at 1.  When the neutron hits
an atom, instead of being absorbed with probability p, multiply W by p.
You can also bias the direction, in favor of of the direction of
the shield boundary, and counteract that sampling bias by multiplying
W by the relative likelihood (likelihood of observed direction under
random sampling / likelihood under biased sampling).  As a neutron
proceeds its W keeps getting up-weighted and down-weighted.  The final
estimate of penetration is based on the average W.  The problem with
this is that the variance of W across neutrons can be huge; most
W's are very small, and a few are much larger than the others.

That is where Russian roulette and splitting come in.  In the middle
of a simulated path, if a neutron has a very small W, it undergoes
Russian roulette - there is a probability P that it is killed; if not
killed its weight is divided by (1-P).  Conversely, if the W is
relatively large, the particle is split into K simulated particles,
each with initial weight W/K, which are subsequently simulated
independently.  The two techniques together reduce the variance of W
across neutrons.

Sampling with unequal probabilities and minimal replacement can be
used to similar effect.  Sample from the current population of
particles with probabilities proportional to 'prob', then divide each
particle's W 

[R] "Inverting" a list

2008-02-06 Thread hadley wickham
Is there a built in function to invert a list?  i.e. to go from

list(
  a = c("1", "2", "3"),
  b = c("1"),
  d = c("2", "4")
)

to

list(
 "1" = c("a", "b"),
 "2" = c("a", "d"),
 "3" = "a",
 "4" = "2"
)

Hadley

-- 
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] "Inverting" a list

2008-02-06 Thread Gabor Grothendieck
This isn't a single command but its pretty short:

unstack(stack(L)[2:1])

On Feb 6, 2008 5:51 PM, hadley wickham <[EMAIL PROTECTED]> wrote:
> Is there a built in function to invert a list?  i.e. to go from
>
> list(
>  a = c("1", "2", "3"),
>  b = c("1"),
>  d = c("2", "4")
> )
>
> to
>
> list(
>  "1" = c("a", "b"),
>  "2" = c("a", "d"),
>  "3" = "a",
>  "4" = "2"
> )
>
> Hadley
>
> --
> http://had.co.nz/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

2008-02-06 Thread Richard M. Heiberger
A more primitive method is about 5 times faster than Gabor's.

L <- list(
a = c("1", "2", "3"),
b = c("1"),
d = c("2", "4")
)

system.time(
for (i in 1:100)
{t1 <- unlist(L)
names(t1) <- rep(names(L), lapply(L, length))
tapply(names(t1), t1, c)
}
)

system.time(
for (i in 1:100)
unstack(stack(L)[2:1])
)

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

2008-02-06 Thread John Smith
Thank all of you for your helps. They are very helpful.

But I have a further question. Suppose I have  the following mixed effect
model

thetaMixed <- function(tau, i)
  {
w <- 1 / (s^2 + tau^2)
mu <- sum(theta * w) / sum(w)
b <- s[i]^2 / (s[i]^2 + tau^2)
theta[i]*(1-b) + mu*b)
  }

I want draw all the mixed effects in a single figure using
for (i in 1:10)
{
   plot(Vectorize(thetaMixed), 0, 2, i=i)
}
and hope plot will recognize that i is the argument for function thetaMixed,
and 0, 2 are the range for tau. But it fails.

Could anyone kindly help me on this issue?

Thanks




On Feb 5, 2008 10:45 PM, Duncan Murdoch <[EMAIL PROTECTED]> wrote:

> jim holtman wrote:
> > Your function 'll' only returns a single value when passed a vector:
> >
> >
> >> x <- seq(0,2,.1)
> >> ll(x)
> >>
> > [1] -7.571559
> >
> >
> > 'plot' expects to pass a vector to the function and have it return a
> > vector of the same length; e.g.,
> >
> >
> >> sin(x)
> >>
> >  [1] 0. 0.09983342 0.19866933 0.29552021 0.38941834 0.47942554
> > 0.56464247 0.64421769 0.71735609
> > [10] 0.78332691 0.84147098 0.89120736 0.93203909 0.96355819 0.98544973
> > 0.99749499 0.99957360 0.99166481
> > [19] 0.97384763 0.94630009 0.90929743
> >
> >
> > So you either have to rewrite your function, or have a loop that will
> > evaluate the function at each individual point and then plot it.
> >
> Or use Vectorize, e.g.
>
> plot(Vectorize(ll), 0, 2)
>
> Duncan Murdoch
> > On Feb 5, 2008 7:06 PM, John Smith <[EMAIL PROTECTED]> wrote:
> >
> >> Dear R-users,
> >>
> >> Suppose I have defined a likelihood function as ll(tau), how can I plot
> this
> >> likelihood function by calling it by plot?
> >>
> >> I want to do it like this:
> >>
> >> ll <- function(tau)
> >>  {
> >>w <- 1 / (s^2 + tau^2)
> >>mu <- sum(theta * w) / sum(w)
> >>-1/2*sum((theta-mu)^2*w -log(w))
> >>  }
> >> plot(ll, 0, 2)
> >>
> >>
> >>
> >> But have the following error:
> >> Error in xy.coords(x, y, xlabel, ylabel, log) :
> >>  'x' and 'y' lengths differ
> >> In addition: Warning messages:
> >> 1: In s^2 + tau^2 :
> >>  longer object length is not a multiple of shorter object length
> >> 2: In theta * w :
> >>  longer object length is not a multiple of shorter object length
> >> 3: In (theta - mu)^2 * w :
> >>  longer object length is not a multiple of shorter object length
> >>
> >>
> >> Thanks
> >>
> >>[[alternative HTML version deleted]]
> >>
> >> __
> >> R-help@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >>
> >>
> >
> >
> >
> >
>
>

[[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] time series plot

2008-02-06 Thread Ana Quit�rio
Dear all.

I want to add a vertical line in my time series series plot.
If i had an anual data I do:
plot(data,v="2008"), but if I want to add a vertical line in fisrt quarter
of 2008, how can i do?

Thanks a lot

Ana
--

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

2008-02-06 Thread mohamed nur anisah
Dear list,
   
  I'm trying to make a loop of a (5x10) matrix and below are my codes. Could 
anybody help me figure out why my loop is not working. Thanks in advance!!
   
   
  m<-1:5
n<-1:10
for(i in 1:length(m))
{ for(j in 1:length(n))
 { 
  y[i,j]=sum(i,j)
  y<-as.matrix(y[i,j]) 
 }
  }
  cheers,
  Anisah

   
-

[[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] time series plot

2008-02-06 Thread jim holtman
 PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html and provide commented,
minimal, self-contained, reproducible code.
>

It all depends on how you constructed the x-axis.  You can try abline(v=2008.25)

On Feb 6, 2008 6:41 PM, Ana Quitério <[EMAIL PROTECTED]> wrote:
> Dear all.
>
> I want to add a vertical line in my time series series plot.
> If i had an anual data I do:
> plot(data,v="2008"), but if I want to add a vertical line in fisrt quarter
> of 2008, how can i do?
>
> Thanks a lot
>
> Ana
> --
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 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] matrix loop

2008-02-06 Thread jim holtman
What exactly are you intending the loop to do?  Why do you have the
'as.matrix' in the middle of the loop?  Where was 'y' defined?   Does
this do what you want?

> outer(1:5, 1:10, "+")
 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]23456789   1011
[2,]3456789   10   1112
[3,]456789   10   11   1213
[4,]56789   10   11   12   1314
[5,]6789   10   11   12   13   1415



On Feb 6, 2008 7:52 PM, mohamed nur anisah <[EMAIL PROTECTED]> wrote:
> Dear list,
>
>  I'm trying to make a loop of a (5x10) matrix and below are my codes. Could 
> anybody help me figure out why my loop is not working. Thanks in advance!!
>
>
>  m<-1:5
> n<-1:10
> for(i in 1:length(m))
> { for(j in 1:length(n))
>  {
>  y[i,j]=sum(i,j)
>  y<-as.matrix(y[i,j])
>  }
>  }
>  cheers,
>  Anisah
>
>
> -
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

What is the problem you are trying to solve?

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


[R] Res: gc() and memory efficiency

2008-02-06 Thread Milton Cezar Ribeiro
Dear Harold,

I had the same problem some times ago. I noticed that after I run a set 
commands (cleaning all non-usefull variables) for 5 times, the system 
broken-down. I solved it building several scritpsNN.R and call them in a .BAT 
DOS file. It worked so fine, almost in my case, and the computer runned for 
several days without stop :-)

If you need more info on this solutions, fill free to write me again.

By other side, if you find a better solutions (I also unfortunatelly run 
windows), share with us.

Kind regards

Miltinho
Brazil




- Mensagem original 
De: Prof Brian Ripley <[EMAIL PROTECTED]>
Para: "Doran, Harold" <[EMAIL PROTECTED]>
Cc: r-help@r-project.org
Enviadas: Terça-feira, 5 de Fevereiro de 2008 3:06:51
Assunto: Re: [R] gc() and memory efficiency

1) See ?"Memory-limits": it is almost certainly memory fragmentation.
You don't need to give the memory back to the OS (and few OSes actually do 
so).

2) I've never seen this running a 64-bit version of R.

3) You can easily write a script to do this.  Indeed, you could write an R 
script to run multiple R scripts in separate processes in turn (via
system("Rscript fileN.R") ).  For example. Uwe Ligges uses R to script 
building and testing of packages on Windows.

On Mon, 4 Feb 2008, Doran, Harold wrote:

> I have a program which reads in a very large data set, performs some 
> analyses, and then repeats this process with another data set. As soon 
> as the first set of analyses are complete, I remove the very large 
> object and clean up to try and make memory available in order to run the 
> second set of analyses. The process looks something like this:
>
> 1) read in data set 1 and perform analyses
> rm(list=ls())
> gc()
> 2) read in data set 2 and perform analyses
> rm(list=ls())
> gc()
> ...
>
> But, it appears that I am not making the memory that was consumed in 
> step 1 available back to the OS as R complains that it cannot allocate a 
> vector of size X as the process tries to repeat in step 2.
>
> So, I close and reopen R and then drop in the code to run the second 
> analysis. When this is done, I close and reopen R and run the third 
> analysis.
>
> This is terribly inefficient. Instead I would rather just source in the 
> R code and let the analyses run over night.
>
> Is there a way that I can use gc() or some other function more 
> efficiently rather than having to close and reopen R at each iteration?
>
> I'm using Windows XP and r 2.6.1
>
> Harold
>
> [[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.



 para armazenamento!

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


  1   2   >