[R] Welch Anova ?

[Hardi]

Hi,

I need to do factor analysis with non-constant variance. Is there a package 
that contains Welch ANOVA ?

Thanks,

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

[R] multiple error bars when NA present

[Junqian Gordon Xu]

package: psych
function: error.bars()
R version: 2.9.1
OS: both linux and windows

a<-c(1,2,3,4,5)
b<-c(1,2,3,4,4)
c<-c(1,2,3,NA,5)
data<-data.frame(a,b,c)

error.bars(data[,1:2],ylim=c(0,6))   looks fine

however

error.bars(data,ylim=c(0.6))

shows multiple error bars for each vector and the number of extra error 
bars equals to the number of vectors with NA value present.


Is this a bug?

Thanks
Gordon

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

Re: [R] How to find the min and max of two variables in a data frame

[Dimitris Rizopoulos]

try this:

set.seed(123)
dat <- data.frame(x = round(rnorm(10)), y = round(rnorm(10)))

dat$Min <- pmin(dat$x, dat$y)
dat$Max <- pmax(dat$x, dat$y)
dat

ind <- dat$Min != dat$Max
dat[ind, ]


I hope it helps.

Best,
Dimitris


kxk wrote:

I have two variables in a data frame, I want to generate two additional
variables. For every observations (i.e. every row), I want the first new
variable 'min' to carry the minimum of the two existing variables, and I
want the second new variable 'max' to carry the maximum of the two existing
variables.  I then want to sort the data frame by min and max, and delete
duplicated rows if both rows has the same min and same max. 


I am new to R so not sure how to fix my code.  Here is my attempt and it is
not working.  Thanks!

edge_dir$min=edge_dir[if (edge_dir$cid_dir1 < edeg_dir$cid_dir2)
edge_dir$cid_dir1 else edge_dir$cid_dir2]
edge_dir$max=edge_dir[if (edge_dir$cid_dir1 > edeg_dir$cid_dir2)
edge_dir$cid_dir1 else edge_dir$cid_dir2]



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

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

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

Re: [R] ggplot2 : commands on one line vs two lines.

[Kingsford Jones]

Hi John,

Thank you for the easily cut-and-pastable example.

You are providing syntactically complete code on the first line so the
second line is not seen as a continuation.  Try moving the '+' on the
second line to the first.  E.g,

pb1 <- p + geom_point(aes(Area, Pcode)) + geom_segment() +
 labs(x="Area", y="Povs")
pb1


hth,

Kingsford Jones

On Thu, Jul 23, 2009 at 5:13 PM, John Kane wrote:
>
> I have just started using ggplot2 and I seem to be doing something stupid
> in writing ggplot2 commands on more than one line.
>
> In the example below the commands on one line are working fine, but as
> soon as I put them on two lines I get an error.  Can any one point out
> what I am doing wrong? It must be something blindingly simple.
>
> Thanks
>
> Example
> #==
> # sample data file
> provs <- structure(list(Name = structure(c(8L, 11L, 6L, 2L, 9L, 1L, 12L,
> 3L, 13L, 5L, 4L, 7L, 10L), .Label = c("Alberta", "British Columbia",
> "Manitoba", "New Brunswick", "Newfoundland", "Northwest Territories",
> "Nova Scotia", "Nunavut", "Ontario", "Prince Edward Island",
> "Québec", "Saskatchewan", "Yukon"), class = "factor"), Pcode = structure(c(8L,
> 11L, 7L, 2L, 9L, 1L, 12L, 3L, 13L, 5L, 4L, 6L, 10L), .Label = c("AB",
> "BC", "MB", "NB", "NL", "NS", "NT", "NU", "ON", "PE", "QC", "SK",
> "YT"), class = "factor"), Area = c(1936113L, 1365128L, 1183085L,
> 925186L, 917741L, 642317L, 591670L, 553556L, 474391L, 373872L,
> 71450L, 53338L, 5660L)), .Names = c("Name", "Pcode", "Area"),
> class = "data.frame", row.names = c(NA,
> -13L))
>
> library(ggplot2)
>
> p <- ggplot(provs, aes(x = 0, xend = Area,  y = Pcode, yend = Pcode))
>
> # commmand below is on one line
> pa <- p + geom_point(aes(Area, Pcode)) + geom_segment() + xlab("Area")
>
> pa
>
> pa1 <- p + geom_point(aes(Area, Pcode)) + geom_segment()
>           + xlab("Area")
> # Error in +xlab("Area") : invalid argument to unary operator
>
>
> # commmand below is on one line
> pb <- p + geom_point(aes(Area, Pcode)) + geom_segment() + labs(x="Area", 
> y="Povs")
> pb
>
>
> pb1 <- p + geom_point(aes(Area, Pcode)) + geom_segment()
>        + labs(x="Area", y="Povs")
>
> # Error in +labs(x = "Area", y = "Povs") :
> #  invalid argument to unary operator
>
> #==
>
> R version 2.9.1 (2009-06-26)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_Canada.1252;LC_CTYPE=English_Canada.1252;
> LC_MONETARY=English_Canada.1252;LC_NUMERIC=C;LC_TIME=English_Canada.1252
>
> attached base packages:
> [1] grid      stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] ggplot2_0.8.3 reshape_0.8.3 plyr_0.1.9    proto_0.3-8
>
>
>
>
>      __
> Looking for the perfect gift? Give the gift of Flickr!
>
> http://www.flickr.com/gift/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] Simply matrix slice question

[Dimitris Rizopoulos]

yes, you need to use the 'drop' argument; have a look at the help file 
?"[" and then try this:


z <- matrix(rnorm(20), ncol = 5)
z[1, , drop = FALSE]
z[, 1, drop = FALSE]


I hope it helps.

Best,
Dimitris


Jim Nemesh wrote:
Is there any way to force a slice of a matrix to stay a matrix? R tends 
to convert a single row of a matrix into a vector.


Example:

z<-matrix (rnorm(20), ncol=5)
zz<-z[1,]
is.matrix(zz) #FALSE

I usually resort to:

zz<-matrix(z[1,], ncol=dim(z)[2], dimnames=list(rownames(z)[1], 
colnames(z)))


But that seems horribly kludgy.

Thanks!

-Jim

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

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



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

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

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

[R] How to find the min and max of two variables in a data frame

[kxk]

I have two variables in a data frame, I want to generate two additional
variables. For every observations (i.e. every row), I want the first new
variable 'min' to carry the minimum of the two existing variables, and I
want the second new variable 'max' to carry the maximum of the two existing
variables.  I then want to sort the data frame by min and max, and delete
duplicated rows if both rows has the same min and same max. 

I am new to R so not sure how to fix my code.  Here is my attempt and it is
not working.  Thanks!

edge_dir$min=edge_dir[if (edge_dir$cid_dir1 < edeg_dir$cid_dir2)
edge_dir$cid_dir1 else edge_dir$cid_dir2]
edge_dir$max=edge_dir[if (edge_dir$cid_dir1 > edeg_dir$cid_dir2)
edge_dir$cid_dir1 else edge_dir$cid_dir2]

-- 
View this message in context: 
http://www.nabble.com/How-to-find-the-min-and-max-of-two-variables-in-a-data-frame-tp24638856p24638856.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Simply matrix slice question

[Jim Nemesh]

Is there any way to force a slice of a matrix to stay a matrix? R  
tends to convert a single row of a matrix into a vector.


Example:

z<-matrix (rnorm(20), ncol=5)
zz<-z[1,]
is.matrix(zz) #FALSE

I usually resort to:

zz<-matrix(z[1,], ncol=dim(z)[2], dimnames=list(rownames(z)[1],  
colnames(z)))


But that seems horribly kludgy.

Thanks!

-Jim

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

Re: [R] Question about the lars package

[lulu9797]

Thanks a lot for the reply. 

1. How I compare the lars v.s. lm

Your understand is correct. 

2. Answer to my question. 

I guess in each step of lars, the coeff are not the same as the lm function
returns using the same "model", because the coeff of lars are the
"accumulated steps" towards the selected direction. 

It's actually my first time using this mailing list. Sorry that I wasn't
aware it's inappropriate to push for answers. But thanks a lot for the
recommendation of glmnet. I'll definitely look into it. 


Steve Lianoglou-6 wrote:
> 
> Hi,
> 
> On Jul 22, 2009, at 6:18 PM, lulu9797 wrote:
> 
>> The returned values of lars function include R squares along the  
>> variable
>> selection path.
> 
> Correct.
> 
>> However, such values are always slightly different from the
>> R squares returned by the regression function lm using the same  
>> models.
>> Anyone know the reasons?
> 
> How are you comparing the models from lars vs. lm?
> 
> Are you just using the non-zero coefs you get from lars in setting up  
> your formula for lm (or something)?
> 
> How different is "slightly different"?
> 
> No idea/speculation: if you're comparing the "dropped coefficients"  
> models to each other, perhaps the coefs on your predictors are  
> slightly different in lars vs. lm since you've got the l1-penalty on  
> them?
> 
> Perhaps a guru will respond with better insight.
> 
>> Very important, and needs quick answers.
> 
> As a point of etiquette, I reckon most people don't really respond too  
> well to requests like these since posting to this mailing list is  
> essentially asking for free help. If it's so urgent/important to you,  
> you can get some professional-for-hire to drop whatever s/he is doing  
> at the moment and work it all out for you.
> 
> By the by, you might want to look at the glmnet package if you find  
> yourself using lars often. Setting the glmnet alpha parameter to 1  
> essentially gives you the lasso regularization path and I've found it  
> to be quite a bit faster than lars.
> 
> -steve
> 
> --
> Steve Lianoglou
> Graduate Student: Physiology, Biophysics and Systems Biology
> Weill Medical College of Cornell University
> 
> Contact Info: http://cbio.mskcc.org/~lianos/contact
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Question-about-the-lars-package-tp24615778p24635487.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Making rq and bootcov play nice

[John Gardner]

I have a quick question, and I apologize in advance if, in asking, I
expose my woeful ignorance of R and its packages. I am trying to use
the bootcov function to estimate the standard errors for some
regression quantiles using a cluster bootstrap. However, it seems that
bootcov passes arguments that rq.fit doesn't like, preventing the
command from executing. Here is an example:

 e<-bootcov(rq(y~x),clust,B=10,fitter=rq.fit)

(where clust is my clustering variable) results in

Error in rq.fit.br(x, y, tau = tau, ...) :
  unused argument(s) (maxit = 15, penalty.matrix = NULL)

In contrast, the lm.fit function seems to just ignore these arguments,
resulting in the following warning:

10: In fitter(X[obs, , drop = FALSE], Y[obs, , drop = FALSE], maxit = maxit,  :
  extra arguments maxitpenalty.matrix are just disregarded.

Is there a way that I can either (a) modify bootcov so that it doesn't
pass these arguments or (b) modify rq so that it ignores them?

Thanks in advance,
John Gardner

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

Re: [R] alternative to rbind within a loop

[Denis Chabot]

Hi Greg,

Thanks, very encouraging: with my example, this is 10x more efficient  
than my loop:

utilisateur système  écoulé
 13.819   5.510  20.204

utilisateur système  écoulé
156.206  44.859 202.150


In real life, I did some work on each file before doing rbind. I'll  
see if this work can be put in a custom-built function that would go  
into the lapply call you suggested.


Denis

Le 09-07-23 à 17:27, Greg Snow a écrit :


Try something like (untested):


mylist <- lapply(all.files, function(i) read.csv(i) )
mydf <- do.call('rbind', mylist)


If all the csv files are conformable that rbind works on them (if  
the loop method works then that should be the case) then this will  
read in each file, store the data frames as a list, then rbind them  
all together.


It seems that this should be faster than the loop, but testing will  
be needed to be sure.


Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
project.org] On Behalf Of Denis Chabot
Sent: Thursday, July 23, 2009 1:54 PM
To: list R
Subject: [R] alternative to rbind within a loop

Hi,

I often have to do this:

select a folder (directory) containing a few hundred data files in  
csv

format (up to 1000 files, in fact)

open each file, transform some character variables in date-tiime  
format


make into a dataframe (involves getting rid of a few variables I  
don't

need

concatenate to the master dataframe that will eventually contain the
data from all the files in the folder.

I use a loop going from 1 to the number of files. I have added a
command to print an incrementing number to the R console each time  
the

loop completes one iteration, to judge the speed of the process.

At the beginning, 3-4 files are processed each second. After a few
hundred iterations it slows down to about 1 file per second. Before I
reach the last file (898 in the case at hand), it has become much
slower, about 1 file every 2-3 seconds.

This progressive slowing down suggests the problem is linked to the
size of the growing "master" dataframe that rbind combines with each
new file.

In fact, the small script below confirms this as nothing at all
happens within the loop but rbind. You can cut the size of this
example not to waste to much of your time:


# create a dummy data.frame and copy it in a large number of csv  
files


test  <- file.path("test")

a <- 1:350
b <- rnorm(350,100,10)
c <- runif(350, 0, 100)
d <- month.name[runif(350,1,12)]

the.data <- data.frame(a,b,c,d)

for(i in 1:850){
write.csv(the.data, file=paste(test, "/file_", i, ".csv",
sep=""))
}

# now lets make a single dataframe from all these csv files

all.files <- list.files(path=test,full.names=T,pattern=".csv")

new.data <- NULL

system.time({
for(i in all.files){
in.data <- read.csv(i)
if (is.null(new.data)) {new.data = in.data} else {new.data =
rbind(new.data, in.data)}
cat(paste(i, ", ", sep=""))
} # end for
}) # end system.time

utilisateur système  écoulé
156.206  44.859 202.150
This is with

sessionInfo()
R version 2.9.1 Patched (2009-07-16 r48939)
x86_64-apple-darwin9.7.0

locale:
fr_CA.UTF-8/fr_CA.UTF-8/C/C/fr_CA.UTF-8/fr_CA.UTF-8

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

other attached packages:
[1] doBy_3.7chron_2.3-30timeDate_290.84

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  Hmisc_3.5-2 lattice_0.17-25
tools_2.9.1


Would it be better to somehow save all 850 files in one dataframe
each, and then rbind them all in a single operation?

Can I combine all my files without using a loop? I've never quite
mastered the "apply" family of functions but have not seen examples  
to

read files.

Thanks in advance,

Denis Chabot

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

[Bert Gunter]

 

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Gabor Grothendieck
Sent: Thursday, July 23, 2009 7:33 PM
To: Nair, Murlidharan T
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] computing the radius of an arc

See ?draw.arc in the plotrix package to draw an arc.

-- Yes, but if I understand the query, the hard part is determining what arc
to draw. This in turn requires understanding how the curve is defined-- is
it defined in closed form mathematically, or does Murli seek a fit that is
an arc of a circle to a sample of data? If the former, this is a problem in
calculus or perhaps differential geometry, not statistics, and that's where
you'd look for methodology. "Best approximates" needs to be defined, too
(what's the distance function?). So all in all, at least as I understand it,
the question seems incoherent, incapable of being answered.

Bert Gunter
Genentech Nonclinical Statistics



On Thu, Jul 23, 2009 at 10:13 PM, Nair, Murlidharan T wrote:
> Hi!!
>
> I am interesting in computing the radius of an arc that best approximates
a curve. Is there an R function that I can use to draw an arc?
> Nothing useful came up when I searched help.search.  Does anyone have any
suggestion to do this?
> Thanks ../Murli

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

[wapita wapita]

Hello,
I have been trying to plot correctly a graph for 2 month now, with no success.I 
want to put a grid, adjusted on the X axis tickers.
Here is the way I build my X-Axis and my grid:


grid(11, NULL, col="grey40")
axis(2)  

# build the tickers on the beginning of each monthticks.at <- 
seq(ISOdatetime(lastYear, 01, 01, hour=0, min=0, sec=0, tz="GMT"), 
ISOdatetime(lastYear, 12, 31, hour=0, min=0, sec=0, tz="GMT"), by = "months")
ticks.lab <- format(ticks.at, format = "%b")
Axis(ISOdatetime(lastYear, substr(index(xts1),6,7), substr(index(xts1),9,10), 
hour=0, min=0, sec=0, tz="GMT"), at = ticks.at, side = 1, labels = ticks.lab, 
las = 2)
box()


If anyone can help me, I would be very grateful.
Xavier
_

[[elided Hotmail spam]]

[[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] Computer Modern

[Paul Murrell]

Hi


Mr Derik wrote:

Thank you for your help.

I've read all the documentation I can find and I still can't get this to
work.

postscriptFonts()

in my console produces a list of fonts already mapped yes?

one of which is:

$ComputerModernItalic
$family
[1] "ComputerModernItalic"

$metrics
[1] "CM_regular_10.afm" "CM_boldx_10.afm"   "cmti10.afm"
"cmbxti10.afm"  "CM_symbol_10.afm"

$encoding
[1] "TeXtext.enc"

attr(,"class")
[1] "Type1Font"

So modifying the code of Paul's website:

pdf("destructiontest.pdf")
par(family="ComputerModernItalic")
demo(plotmath)
dev.off()

Should give me what I want. Instead I get a file which can't be opened and
25 warnings which all look like this:

25: In text.default(2 * c - 0.5, -r, title) :
  font family not found in PostScript font database



But you are creating a PDF file, not a PostScript file, and if you look 
at names(pdfFonts()) you'll see that there is no ComputerModernItalic there.




I've tried the "Test full Adobe Symbol font" code off Paul's website with
the fonts loaded into my working dir as described. I get 35 warnings of the
type:

1: In grid.Call.graphics("L_text", as.graphicsAnnot(x$label),  ... :
  font width unknown for character 0x7f

and a pdf with dots in place of every character.

I can open the example pdf on Paul's website and see the correct CM font so
my viewer must be working.

Any other ideas would be welcome.



Is it possible that you have set the fonts up using PostScriptFonts() 
when you should really be using pdfFonts() ?


Paul



Paul Murrell wrote:

Hi

Also see http://www.stat.auckland.ac.nz/~paul/R/CM/CMR.html

Paul


(Ted Harding) wrote:

On 02-Jul-09 09:06:44, Mr Derik wrote:

I am trying to use computer modern fonts in postscript files
for a latex document. Ultimately I want to automate this through
sweave. I've read the documentation ans have tried the following
code to use lattice to produce a graph using computer modern:

library(lattice)
library(grid)
testPlot=(
xyplot(seq(1:10) ~ seq(1:10),
main="one to ten",
xlab="the quick fox",
ylab="jumped over the lazy brown dog",
xlim=c(0,1),
ylim=c(0,1),
col="black",
type="l" ,
lwd=2
)
)
setwd("C:\\R_folder\\CMtests")
postscript("cm_test.eps", width = 4.0, height = 3.0,
   horizontal = FALSE, onefile = FALSE, paper = "special",
   family = "ComputerModern", encoding = "TeXtext.enc")
   print(testPlot)
dev.off()

This produces a plot with courier.

I am using R 2.9.0 on a windows XP machine. I did manage to produce
one plot with CM as the font so I know it's possible with my set up.
I can't get back to that. Please help me with the code.
Thank You

I think you may need to also use the "fonts" pAramater to postscript().
See in '?postscript':

  fonts: a character vector specifying additional R graphics font
family names for font families whose declarations will be
included in the PostScript file and are available for use
with the device. See 'Families' below.  Defaults to 'NULL'.

Since the Computer Modern family is most probably not built in
to your printer, the PostScript file will need to include font
definitions for these fonts. If I understand aright, this is what
would be achieved by appropriate use of the "fonts" parameter.

If the font definitions are not included, the calls for them will
not be recognised by the printer which may then substitute a default
(likely to be Courier).

See also the section "TeX fonts" in '?postscript'.

Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 02-Jul-09   Time: 11:59:29
-- 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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




--
View this message in context: 
http://www.nabble.com/Computer-Modern-tp24303553p24424418.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.


--
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag

Re: [R] Duplicated date values aren't duplicates

[jim holtman]

What you are reporting is that there are two duplicated dates in your
data.  'duplicated' returns a logical vector that is TRUE for the
second and subsequent duplicates.  Notice what is returned:

> x <- c(1,2,2,3,4,4,5,6,4,7,3)
> x[duplicated(x)]
[1] 2 4 4 3
>



On Thu, Jul 23, 2009 at 8:50 PM, Tim Clark wrote:
>
> Dear list,
>
> I just had a function (as.ltraj in Adehabitat) give me the following error:
>
> "Error in as.ltraj(xy, id, date = da) : non unique dates for a given burst"
>
> I checked my dates and got the following:
>
>>   dupes<-mydata$DateTime[duplicated(mydata$DateTime)]
>> dupes
> [1] (07/30/02 00:00:00) (08/06/03 17:45:00)
>
> Is there a reason different dates would come up as duplicate values?  I would 
> prefer not to have to delete them if I don't have to.  Any suggestions on how 
> to get R to realize they are different?
>
> Thanks,
>
> Tim
>
>
>
> Tim Clark
> Department of Zoology
> University of Hawaii
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

What is the problem that you are trying to solve?

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

Re: [R] how to predict dynamic model in R

[Gabor Grothendieck]

If by "problem" you mean the problem of determining how, in general,
you should proceed perhaps you need an introductory guide on R
and time series such as Cowpertwait's book.

On Thu, Jul 23, 2009 at 10:37 PM, Hongwei Dong wrote:
> Hi, Gabor, it seems ARIMA model does not have that problem. For example:
> set.seed(123)
> y<-ts(c(1:20))
> x = ts(rnorm(20))
> z = ts(rnorm(20))
> tt<-ts(cbind(x, lag(x,-1),lag(x,-2),z))
> fit <- arima(y[1:15],order=c(1,0,0),xreg=tt[(1:15),])
> fit
> pred <- predict(fit, n.ahead=5,tt[(16:20),])
> pred
> What do you make of this? Thanks.
> Harry
>
> On Thu, Jul 23, 2009 at 6:36 PM, Gabor Grothendieck
>  wrote:
>>
>> Try this:
>>
>> library(dyn)
>> set.seed(123)
>> tz <- zoo(cbind(Y = 0, x = rnorm(10), z = rnorm(10)))
>>
>> # simulate values
>> for(i in 2:10) {
>>  tz$Y[i] <- with(tz, 2*Y[i-1] + 3*z[i] +4* x[i] + 5*x[i-1] + rnorm(1))
>> }
>>
>> # keep copy of tz to compare later to simulated Y's
>> tz.orig <- tz
>>
>> # NA out Y's that are to be predicted
>> tz[7:10, "Y"] <- NA
>>
>> L <- function(x, k = 1) lag(x, -k)
>>
>> # predict 1 ahead each iteration
>> for(i in 7:10) {
>>    # fit based on first i-1 values
>>    fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tz, subset = seq_len(i-1))
>>    # get prediction for ith value
>>    tz[i, "Y"] <- tail(predict(fit, tz[1:i,]), 1)
>> }
>> cbind(pred = tz[7:10, "Y"], act = tz.orig[7:10, "Y"])
>>
>>
>> On Thu, Jul 23, 2009 at 9:02 PM, Hongwei Dong wrote:
>> > What I want R to do is to use the estimated Y at t-1 to be the lag(Y,-1)
>> > in
>> > the forecast equation for time t. Is there anyway I can realize this
>> > with R?
>> > For example, when the Y value for year 18 is forecast, the estimated Y
>> > for
>> > year 17 is used, not the actual Y for year 17 already in the data.
>> > Thanks for you patience. I appreciate it.
>> > Harry
>> >
>> >
>> >
>> > On Thu, Jul 23, 2009 at 5:44 PM, Gabor Grothendieck
>> >  wrote:
>> >>
>> >> You can't remove Y since its in the rhs of your model.
>> >>
>> >> On Thu, Jul 23, 2009 at 8:25 PM, Hongwei Dong wrote:
>> >> > Thanks, Gabor. Here are the problems I'm trying to solve.
>> >> > FIRST, I run this to simulate a 20 years time series process. The
>> >> > data
>> >> > from
>> >> > 1-15 years are used to estimate the model, and this model is used to
>> >> > predict
>> >> > the year from 16-20. The following script works.
>> >> > set.seed(123)
>> >> > tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
>> >> > L <- function(x, k = 1) lag(x, -k)
>> >> > tt.zoo <- as.zoo(tt)
>> >> > fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
>> >> > fit
>> >> > pred <- predict(fit, tt.zoo[(16:20),])
>> >> > pred
>> >> > SECONDLY, I use similar script, but pretend that we do not know the Y
>> >> > data
>> >> > from year 16-20. We know x and z for year 16-20, and use them predict
>> >> > Y
>> >> > based on the model estimated from 1-15 years. So, in the "newdata"
>> >> > part,
>> >> > I
>> >> > use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
>> >> > Unfortunately, it does not work in that way.
>> >> > pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
>> >> > pred1
>> >> > It will be greatly appreciated if you can give me some guide on this.
>> >> > Thanks.
>> >> > Harry
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >> > On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck
>> >> >  wrote:
>> >> >>
>> >> >> Use dyn.predict like this:
>> >> >>
>> >> >> > library(dyn)
>> >> >> > x <- y <- zoo(1:5)
>> >> >> > mod <- dyn$lm(y ~ lag(x, -1))
>> >> >> > predict(mod, list(x = zoo(6:10, 6:10)))
>> >> >>  7  8  9 10
>> >> >>  7  8  9 10
>> >> >>
>> >> >>
>> >> >> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
>> >> >> wrote:
>> >> >> > I have a dynamic time series model like this:
>> >> >> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
>> >> >> >
>> >> >> > I need to do an out of sample forecast with this model. Is there
>> >> >> > any
>> >> >> > way
>> >> >> > I
>> >> >> > can do this with R?
>> >> >> > It would be greatly appreciated if some one can give me an
>> >> >> > example.
>> >> >> > Thanks.
>> >> >> >
>> >> >> >
>> >> >> > Harry
>> >> >> >
>> >> >> >        [[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 predict dynamic model in R

[Hongwei Dong]

Hi, Gabor, it seems ARIMA model does not have that problem. For example:
set.seed(123)
y<-ts(c(1:20))
x = ts(rnorm(20))
z = ts(rnorm(20))
tt<-ts(cbind(x, lag(x,-1),lag(x,-2),z))

fit <- arima(y[1:15],order=c(1,0,0),xreg=tt[(1:15),])
fit
pred <- predict(fit, n.ahead=5,tt[(16:20),])
pred

What do you make of this? Thanks.

Harry


On Thu, Jul 23, 2009 at 6:36 PM, Gabor Grothendieck  wrote:

> Try this:
>
> library(dyn)
> set.seed(123)
> tz <- zoo(cbind(Y = 0, x = rnorm(10), z = rnorm(10)))
>
> # simulate values
> for(i in 2:10) {
>   tz$Y[i] <- with(tz, 2*Y[i-1] + 3*z[i] +4* x[i] + 5*x[i-1] + rnorm(1))
> }
>
> # keep copy of tz to compare later to simulated Y's
> tz.orig <- tz
>
> # NA out Y's that are to be predicted
> tz[7:10, "Y"] <- NA
>
> L <- function(x, k = 1) lag(x, -k)
>
> # predict 1 ahead each iteration
> for(i in 7:10) {
># fit based on first i-1 values
>fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tz, subset = seq_len(i-1))
># get prediction for ith value
>tz[i, "Y"] <- tail(predict(fit, tz[1:i,]), 1)
> }
> cbind(pred = tz[7:10, "Y"], act = tz.orig[7:10, "Y"])
>
>
> On Thu, Jul 23, 2009 at 9:02 PM, Hongwei Dong wrote:
> > What I want R to do is to use the estimated Y at t-1 to be the lag(Y,-1)
> in
> > the forecast equation for time t. Is there anyway I can realize this with
> R?
> > For example, when the Y value for year 18 is forecast, the estimated Y
> for
> > year 17 is used, not the actual Y for year 17 already in the data.
> > Thanks for you patience. I appreciate it.
> > Harry
> >
> >
> >
> > On Thu, Jul 23, 2009 at 5:44 PM, Gabor Grothendieck
> >  wrote:
> >>
> >> You can't remove Y since its in the rhs of your model.
> >>
> >> On Thu, Jul 23, 2009 at 8:25 PM, Hongwei Dong wrote:
> >> > Thanks, Gabor. Here are the problems I'm trying to solve.
> >> > FIRST, I run this to simulate a 20 years time series process. The data
> >> > from
> >> > 1-15 years are used to estimate the model, and this model is used to
> >> > predict
> >> > the year from 16-20. The following script works.
> >> > set.seed(123)
> >> > tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
> >> > L <- function(x, k = 1) lag(x, -k)
> >> > tt.zoo <- as.zoo(tt)
> >> > fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
> >> > fit
> >> > pred <- predict(fit, tt.zoo[(16:20),])
> >> > pred
> >> > SECONDLY, I use similar script, but pretend that we do not know the Y
> >> > data
> >> > from year 16-20. We know x and z for year 16-20, and use them predict
> Y
> >> > based on the model estimated from 1-15 years. So, in the "newdata"
> part,
> >> > I
> >> > use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
> >> > Unfortunately, it does not work in that way.
> >> > pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
> >> > pred1
> >> > It will be greatly appreciated if you can give me some guide on this.
> >> > Thanks.
> >> > Harry
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> > On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck
> >> >  wrote:
> >> >>
> >> >> Use dyn.predict like this:
> >> >>
> >> >> > library(dyn)
> >> >> > x <- y <- zoo(1:5)
> >> >> > mod <- dyn$lm(y ~ lag(x, -1))
> >> >> > predict(mod, list(x = zoo(6:10, 6:10)))
> >> >>  7  8  9 10
> >> >>  7  8  9 10
> >> >>
> >> >>
> >> >> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
> >> >> wrote:
> >> >> > I have a dynamic time series model like this:
> >> >> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
> >> >> >
> >> >> > I need to do an out of sample forecast with this model. Is there
> any
> >> >> > way
> >> >> > I
> >> >> > can do this with R?
> >> >> > It would be greatly appreciated if some one can give me an example.
> >> >> > Thanks.
> >> >> >
> >> >> >
> >> >> > Harry
> >> >> >
> >> >> >[[alternative HTML version deleted]]
> >> >> >
> >> >> > __
> >> >> > R-help@r-project.org mailing list
> >> >> > https://stat.ethz.ch/mailman/listinfo/r-help
> >> >> > PLEASE do read the posting guide
> >> >> > http://www.R-project.org/posting-guide.html
> >> >> > and provide commented, minimal, self-contained, reproducible code.
> >> >> >
> >> >
> >> >
> >
> >
>

[[alternative HTML version deleted]]

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

Re: [R] computing the radius of an arc

[Gabor Grothendieck]

See ?draw.arc in the plotrix package to draw an arc.

On Thu, Jul 23, 2009 at 10:13 PM, Nair, Murlidharan T wrote:
> Hi!!
>
> I am interesting in computing the radius of an arc that best approximates a 
> curve. Is there an R function that I can use to draw an arc?
> Nothing useful came up when I searched help.search.  Does anyone have any 
> suggestion to do this?
> Thanks ../Murli

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

[Nair, Murlidharan T]

Hi!!

I am interesting in computing the radius of an arc that best approximates a 
curve. Is there an R function that I can use to draw an arc? 
Nothing useful came up when I searched help.search.  Does anyone have any 
suggestion to do this?
Thanks ../Murli

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

Re: [R] Sum with Conditions

[jimdare]

Thanks very much :)


Daniel Malter wrote:
> 
> length(datasetname$Events[datasetname$Length<28 &
> datasetname$Method%in%c(1,2,3) ...additional conditions...  , ]
> returns
> the number of rows
> sum(datasetname$Events[datasetname$Length<28 &
> datasetname$Method%in%c(1,2,3) ...additional conditions...  , ]
> returns
> the sum of the Events variable given your conditions
> 
> does that work for you?
> Daniel
> 
> -
> cuncta stricte discussurus
> -
> 
> -Ursprüngliche Nachricht-
> Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Im
> Auftrag von jimdare
> Gesendet: Thursday, July 23, 2009 8:07 PM
> An: r-help@r-project.org
> Betreff: [R] Sum with Conditions
> 
> 
> Dear R-Users,
> 
> I have a huge data set (172,696 rows) that is set out in the following
> way:
> 
> Area | Length | Species | Method | Date | Events  
> 
> I want to sum the number of Events if:
> 
> Length < 28AND
> Species is NOT "X" ,"Y", or "Z"   AND
> Method is either "1" ,"2", or "3"   AND
> Date is between 31/12/08 AND 1/02/09 (i.e. JAN 09) 
> 
> Is there an easy way to do this?
> 
> Your help is much appreciated.
> 
> James  
> 
>   
> --
> View this message in context:
> http://www.nabble.com/Sum-with-Conditions-tp24636551p24636551.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Sum-with-Conditions-tp24636551p24637520.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] Help with Loop!

[Steve Lianoglou]

Hi,

On Jul 23, 2009, at 7:30 PM, Lars Bishop wrote:


Dear experts,

I'm new in R and trying to learn by writing a version of the  
Perceptron
Algorithm. How can I tell in the code below to stop the iteration  
when the
condition in the "for loop" is not satisfied for all training  
examples?


Thanks in advance for your help!


## Generate a linearly separable data set in R2

sample <- as.data.frame(cbind(runif(n=100),runif(n=100)))
sample$Y <- ifelse(sample$V1>sample$V2,1,-1)
## Plot data;
attach(sample)
plot(V1, V2, type="p", pch=Y,main="Sample Data")

##Perceptron algorithm

sample_m <- as.matrix(sample)
w <- c(0,0); b <- 0; k <- 0; nu <- 1e-3
R <- max(sqrt(V1^2+V2^2))

repeat {
for (i in 1:nrow(sample_m)){
   if  (sample_m[i,3]*(t(w)%*%sample_m[i,1:2] + b) <= 0)
 w <- w + nu*sample_m[i,3]*sample_m[i,1:2]
 b <- b +nu*sample_m[i,3]*R^2
 k <- k+1
 cat(k, w, b, "\n") }
}



One thing to realize is that by using a variable named "sample",  
you're trampling over the "sample" function ... just don't trip over  
that later in your R career :-)


I whipped up two version of your repeat block (I'm using while  
instead, though). One uses the inner for loop like you have, the  
second one is doing the same thing but w/ no for loops:


set.seed(123)
examples <- as.data.frame(cbind(runif(n=100),runif(n=100)))
examples$Y <- ifelse(examples$V1 > examples$V2, 1, -1)
examples <- as.matrix(examples)

# plot(examples[,1], examples[,2], type="p", pch=examples[,3],  
main="Sample Data")


b <- 0; k <- 0; nu <- 1e-3;
R <- max(sqrt(examples[,1]^2 + examples[,2]^2))
w <- matrix(0, nrow=2, ncol=1)

misclassified <- examples[,3] * (examples[,1:2] %*% w) <= 0
while (any(misclassified)) {
  for (i in which(misclassified)) {
  w <- w + nu * examples[i,3] * examples[i,1:2]
  b <- b + nu * examples[i,3]*R^2
  k <- k + 1
  cat(k, w, b, "\n")
  }
  misclassified <- examples[,3] * (examples[,1:2] %*% w) <= 0
}


## Here's the while loop w/o for loops
misclassified <- examples[,3] * (examples[,1:2] %*% w) <= 0
while (any(misclassified)) {
  w <- w + colSums(nu * examples[misclassified,3] *  
examples[misclassified,1:2])

  b <- b + sum(nu * examples[misclassified,3] * R^2)
  k <- k + 1
  cat(k, w, b, "\n")
  misclassified <- examples[,3] * (examples[,1:2] %*% w) <= 0
}

HTH,
-steve

--
Steve Lianoglou
Graduate Student: Physiology, Biophysics and Systems Biology
Weill Medical College of Cornell University

Contact Info: http://cbio.mskcc.org/~lianos

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

[Gabor Grothendieck]

Try this:

library(dyn)
set.seed(123)
tz <- zoo(cbind(Y = 0, x = rnorm(10), z = rnorm(10)))

# simulate values
for(i in 2:10) {
  tz$Y[i] <- with(tz, 2*Y[i-1] + 3*z[i] +4* x[i] + 5*x[i-1] + rnorm(1))
}

# keep copy of tz to compare later to simulated Y's
tz.orig <- tz

# NA out Y's that are to be predicted
tz[7:10, "Y"] <- NA

L <- function(x, k = 1) lag(x, -k)

# predict 1 ahead each iteration
for(i in 7:10) {
# fit based on first i-1 values
fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tz, subset = seq_len(i-1))
# get prediction for ith value
tz[i, "Y"] <- tail(predict(fit, tz[1:i,]), 1)
}
cbind(pred = tz[7:10, "Y"], act = tz.orig[7:10, "Y"])


On Thu, Jul 23, 2009 at 9:02 PM, Hongwei Dong wrote:
> What I want R to do is to use the estimated Y at t-1 to be the lag(Y,-1) in
> the forecast equation for time t. Is there anyway I can realize this with R?
> For example, when the Y value for year 18 is forecast, the estimated Y for
> year 17 is used, not the actual Y for year 17 already in the data.
> Thanks for you patience. I appreciate it.
> Harry
>
>
>
> On Thu, Jul 23, 2009 at 5:44 PM, Gabor Grothendieck
>  wrote:
>>
>> You can't remove Y since its in the rhs of your model.
>>
>> On Thu, Jul 23, 2009 at 8:25 PM, Hongwei Dong wrote:
>> > Thanks, Gabor. Here are the problems I'm trying to solve.
>> > FIRST, I run this to simulate a 20 years time series process. The data
>> > from
>> > 1-15 years are used to estimate the model, and this model is used to
>> > predict
>> > the year from 16-20. The following script works.
>> > set.seed(123)
>> > tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
>> > L <- function(x, k = 1) lag(x, -k)
>> > tt.zoo <- as.zoo(tt)
>> > fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
>> > fit
>> > pred <- predict(fit, tt.zoo[(16:20),])
>> > pred
>> > SECONDLY, I use similar script, but pretend that we do not know the Y
>> > data
>> > from year 16-20. We know x and z for year 16-20, and use them predict Y
>> > based on the model estimated from 1-15 years. So, in the "newdata" part,
>> > I
>> > use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
>> > Unfortunately, it does not work in that way.
>> > pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
>> > pred1
>> > It will be greatly appreciated if you can give me some guide on this.
>> > Thanks.
>> > Harry
>> >
>> >
>> >
>> >
>> >
>> >
>> > On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck
>> >  wrote:
>> >>
>> >> Use dyn.predict like this:
>> >>
>> >> > library(dyn)
>> >> > x <- y <- zoo(1:5)
>> >> > mod <- dyn$lm(y ~ lag(x, -1))
>> >> > predict(mod, list(x = zoo(6:10, 6:10)))
>> >>  7  8  9 10
>> >>  7  8  9 10
>> >>
>> >>
>> >> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
>> >> wrote:
>> >> > I have a dynamic time series model like this:
>> >> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
>> >> >
>> >> > I need to do an out of sample forecast with this model. Is there any
>> >> > way
>> >> > I
>> >> > can do this with R?
>> >> > It would be greatly appreciated if some one can give me an example.
>> >> > Thanks.
>> >> >
>> >> >
>> >> > Harry
>> >> >
>> >> >        [[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 predict dynamic model in R

[Hongwei Dong]

What I want R to do is to use the estimated Y at t-1 to be the lag(Y,-1) in
the forecast equation for time t. Is there anyway I can realize this with R?
For example, when the Y value for year 18 is forecast, the estimated Y for
year 17 is used, not the actual Y for year 17 already in the data.
Thanks for you patience. I appreciate it.
Harry



On Thu, Jul 23, 2009 at 5:44 PM, Gabor Grothendieck  wrote:

> You can't remove Y since its in the rhs of your model.
>
> On Thu, Jul 23, 2009 at 8:25 PM, Hongwei Dong wrote:
> > Thanks, Gabor. Here are the problems I'm trying to solve.
> > FIRST, I run this to simulate a 20 years time series process. The data
> from
> > 1-15 years are used to estimate the model, and this model is used to
> predict
> > the year from 16-20. The following script works.
> > set.seed(123)
> > tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
> > L <- function(x, k = 1) lag(x, -k)
> > tt.zoo <- as.zoo(tt)
> > fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
> > fit
> > pred <- predict(fit, tt.zoo[(16:20),])
> > pred
> > SECONDLY, I use similar script, but pretend that we do not know the Y
> data
> > from year 16-20. We know x and z for year 16-20, and use them predict Y
> > based on the model estimated from 1-15 years. So, in the "newdata" part,
> I
> > use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
> > Unfortunately, it does not work in that way.
> > pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
> > pred1
> > It will be greatly appreciated if you can give me some guide on this.
> > Thanks.
> > Harry
> >
> >
> >
> >
> >
> >
> > On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck
> >  wrote:
> >>
> >> Use dyn.predict like this:
> >>
> >> > library(dyn)
> >> > x <- y <- zoo(1:5)
> >> > mod <- dyn$lm(y ~ lag(x, -1))
> >> > predict(mod, list(x = zoo(6:10, 6:10)))
> >>  7  8  9 10
> >>  7  8  9 10
> >>
> >>
> >> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
> wrote:
> >> > I have a dynamic time series model like this:
> >> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
> >> >
> >> > I need to do an out of sample forecast with this model. Is there any
> way
> >> > I
> >> > can do this with R?
> >> > It would be greatly appreciated if some one can give me an example.
> >> > Thanks.
> >> >
> >> >
> >> > Harry
> >> >
> >> >[[alternative HTML version deleted]]
> >> >
> >> > __
> >> > R-help@r-project.org mailing list
> >> > https://stat.ethz.ch/mailman/listinfo/r-help
> >> > PLEASE do read the posting guide
> >> > http://www.R-project.org/posting-guide.html
> >> > and provide commented, minimal, self-contained, reproducible code.
> >> >
> >
> >
>

[[alternative HTML version deleted]]

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

Re: [R] Sum with Conditions

[Daniel Malter]

length(datasetname$Events[datasetname$Length<28 &
datasetname$Method%in%c(1,2,3) ...additional conditions...  , ] returns
the number of rows
sum(datasetname$Events[datasetname$Length<28 &
datasetname$Method%in%c(1,2,3) ...additional conditions...  , ] returns
the sum of the Events variable given your conditions

does that work for you?
Daniel

-
cuncta stricte discussurus
-

-Ursprüngliche Nachricht-
Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Im
Auftrag von jimdare
Gesendet: Thursday, July 23, 2009 8:07 PM
An: r-help@r-project.org
Betreff: [R] Sum with Conditions


Dear R-Users,

I have a huge data set (172,696 rows) that is set out in the following way:

Area | Length | Species | Method | Date | Events

I want to sum the number of Events if:

Length < 28AND
Species is NOT "X" ,"Y", or "Z"   AND
Method is either "1" ,"2", or "3"   AND
Date is between 31/12/08 AND 1/02/09 (i.e. JAN 09) 

Is there an easy way to do this?

Your help is much appreciated.

James  

  
--
View this message in context:
http://www.nabble.com/Sum-with-Conditions-tp24636551p24636551.html
Sent from the R help mailing list archive at Nabble.com.

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

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

Re: [R] Assigning rank based on total count

[ws]

> 
> Here's a way to get to your solution, but it's not very pretty:
> 
> testdfr <- data.frame(POB=c("Oregon","Oregon","Oregon","New
> York","California","California"))
> 
> nstates <- length(unique(testdfr$POB))
> testdfr$ POBR <- c(nstates:1)[table(testdfr$POB)][testdfr$POB]

Hmm  I will have to think on that for a while  An ugly
solution is still better than no solution, but that doesn't quite 
make sense to me yet.

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

[Tim Clark]

Dear list,

I just had a function (as.ltraj in Adehabitat) give me the following error:

"Error in as.ltraj(xy, id, date = da) : non unique dates for a given burst"

I checked my dates and got the following:

>   dupes<-mydata$DateTime[duplicated(mydata$DateTime)]
> dupes
[1] (07/30/02 00:00:00) (08/06/03 17:45:00)

Is there a reason different dates would come up as duplicate values?  I would 
prefer not to have to delete them if I don't have to.  Any suggestions on how 
to get R to realize they are different?

Thanks,

Tim



Tim Clark
Department of Zoology 
University of Hawaii

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

[Gabor Grothendieck]

You can't remove Y since its in the rhs of your model.

On Thu, Jul 23, 2009 at 8:25 PM, Hongwei Dong wrote:
> Thanks, Gabor. Here are the problems I'm trying to solve.
> FIRST, I run this to simulate a 20 years time series process. The data from
> 1-15 years are used to estimate the model, and this model is used to predict
> the year from 16-20. The following script works.
> set.seed(123)
> tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
> L <- function(x, k = 1) lag(x, -k)
> tt.zoo <- as.zoo(tt)
> fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
> fit
> pred <- predict(fit, tt.zoo[(16:20),])
> pred
> SECONDLY, I use similar script, but pretend that we do not know the Y data
> from year 16-20. We know x and z for year 16-20, and use them predict Y
> based on the model estimated from 1-15 years. So, in the "newdata" part, I
> use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
> Unfortunately, it does not work in that way.
> pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
> pred1
> It will be greatly appreciated if you can give me some guide on this.
> Thanks.
> Harry
>
>
>
>
>
>
> On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck
>  wrote:
>>
>> Use dyn.predict like this:
>>
>> > library(dyn)
>> > x <- y <- zoo(1:5)
>> > mod <- dyn$lm(y ~ lag(x, -1))
>> > predict(mod, list(x = zoo(6:10, 6:10)))
>>  7  8  9 10
>>  7  8  9 10
>>
>>
>> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong wrote:
>> > I have a dynamic time series model like this:
>> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
>> >
>> > I need to do an out of sample forecast with this model. Is there any way
>> > I
>> > can do this with R?
>> > It would be greatly appreciated if some one can give me an example.
>> > Thanks.
>> >
>> >
>> > Harry
>> >
>> >        [[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 predict dynamic model in R

[Hongwei Dong]

Thanks, Gabor. Here are the problems I'm trying to solve.
*FIRST*, I run this to simulate a 20 years time series process. The data
from 1-15 years are used to estimate the model, and this model is used to
predict the year from 16-20. The following script works.

set.seed(123)
tt <- ts(cbind(Y = 1:20, x = rnorm(20), z = rnorm(20)))
L <- function(x, k = 1) lag(x, -k)
tt.zoo <- as.zoo(tt)
fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), tt.zoo[(1:15), ])
fit
pred <- predict(fit, tt.zoo[(16:20),])
pred

*SECONDLY*, I use similar script, but pretend that we do not know the Y data
from year 16-20. We know x and z for year 16-20, and use them predict Y
based on the model estimated from 1-15 years. So, in the "newdata" part, I
use tt.zoo[(16:20), (2:3)] to remove the Y out. here is the script.
Unfortunately, it does not work in that way.

pred1 <- predict(fit, tt.zoo[(16:20),(2:3)])
pred1

It will be greatly appreciated if you can give me some guide on this.
Thanks.

Harry






On Wed, Jul 22, 2009 at 10:04 PM, Gabor Grothendieck <
ggrothendi...@gmail.com> wrote:

> Use dyn.predict like this:
>
> > library(dyn)
> > x <- y <- zoo(1:5)
> > mod <- dyn$lm(y ~ lag(x, -1))
> > predict(mod, list(x = zoo(6:10, 6:10)))
>  7  8  9 10
>  7  8  9 10
>
>
> On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong wrote:
> > I have a dynamic time series model like this:
> > dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
> >
> > I need to do an out of sample forecast with this model. Is there any way
> I
> > can do this with R?
> > It would be greatly appreciated if some one can give me an example.
> Thanks.
> >
> >
> > Harry
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>

[[alternative HTML version deleted]]

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

Re: [R] rows missing after dataset loaded to R

[Gabor Grothendieck]

Another situation would be if you have comment characters in strings that
are intended to be content.

On Thu, Jul 23, 2009 at 5:35 PM, Greg Snow wrote:
> Some programs quote everything to be "safe", others only quote when needed.  
> The only case that I know of that read.table and friends require quotes for 
> is when a separator is inside of a string, for example if you are using 
> spaces as the separator and have some names with spaces in them (e.g. "North 
> Dakota"), without the quotes that would be seen as 2 fields, with the quotes 
> it is a single field.  If the coma (,) is the separator and you have names 
> (e.g. "Snow, Greg") then you would need the quotes.  If you don't have any 
> cases of the separators other than where they are separating fields, then the 
> quotes are probably not needed.
>
> Hope this helps,
>
> --
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> greg.s...@imail.org
> 801.408.8111
>
>
>> -Original Message-
>> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
>> project.org] On Behalf Of Rnewbie
>> Sent: Thursday, July 23, 2009 11:02 AM
>> To: r-help@r-project.org
>> Subject: Re: [R] rows missing after dataset loaded to R
>>
>>
>> Thank you very much for the reply. I checked the rows and it was the
>> unbalanced " quote marks in some of the rows that caused the problem.
>> Once I
>> disabled quoting altogether, the problem is solved.
>>
>> I have one more basic question.  I disabled quoting when loading the
>> file to
>> R, and all the columns consisting of characters with or without
>> quotation
>> marks displayed normally. What does quoting actually do in R?
>>
>> --
>> View this message in context: http://www.nabble.com/rows-missing-after-
>> dataset-loaded-to-R-tp24625882p24630154.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-
>> guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

[jimdare]

Dear R-Users,

I have a huge data set (172,696 rows) that is set out in the following way:

Area | Length | Species | Method | Date | Events

I want to sum the number of Events if:

Length < 28AND
Species is NOT "X" ,"Y", or "Z"   AND
Method is either "1" ,"2", or "3"   AND
Date is between 31/12/08 AND 1/02/09 (i.e. JAN 09) 

Is there an easy way to do this?

Your help is much appreciated.

James  

  
-- 
View this message in context: 
http://www.nabble.com/Sum-with-Conditions-tp24636551p24636551.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] Navigate to Index page of a package from R command prompt

[Gabor Grothendieck]

Try

enter at the R console: help.start()
and then when the help comes up in your browser click on Packages
and then click on the package you want
and then click on the help file you want

On Thu, Jul 23, 2009 at 3:30 PM, Steven McKinney wrote:
>
> Hi all,
>
> Is there a way to navigate directly to the "Index" page of help
> for a package?
>
> Here's my connundrum:
>
> I download and install package "foo".
> I don't know what functions are in package "foo",
> so I can't invoke the help for package "foo" via
>> ?someFunction
>
> help(package = "foo")
> pops up some non-hyperlinked information page, not
> package foo's help Index page.
>
> If the package author kindly made a "foo" object or function
> and put that in package "foo", then
>> ?foo
> works and yields a help page for package "foo".
> Now at the bottom of the help page is a hyperlink "Index"
> and I can click that to navigate to the main help Index page
> (the page I really want to get to straight from the R
> command line).
>
> I see that the link to "Index" for package "foo" appears always to be
> (on my Mac)
> file:///Library/Frameworks/R.framework/Resources/library/foo/html/00Index.html
>
> e.g.
> file:///Library/Frameworks/R.framework/Resources/library/cmprsk/html/00Index.html
>
> file:///Library/Frameworks/R.framework/Resources/library/utils/html/00Index.html
>
> Is there a command from the R listener that can take me directly to
> this "00Index.html" page of help for package "foo"?
>
> something like
>> help("00Index", package = "utils")
>
> (but this does not work)?
>
> Any info appreciated
>
> Best
>
> Steven McKinney
>
> Statistician
> Molecular Oncology and Breast Cancer Program
> British Columbia Cancer Research Centre
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Random # generator accuracy

[Jim Bouldin]

Perfectly explained Ted.  One might, at first reflection, consider that
simply repeating the values 7 through 12 and  sampling (w/o replacement)
from among the 18 resulting values, would be similar to just doubling the
selection probabilities for 7 through 12 and then sampling. That would
clearly not be true though.
Jim
> 
> Whereas, if you replace x = c(1,2,3,4,5,6,7,8,9,110,11,12)
> with the "weighted" equivalent, doubling up 7-12 as in your
>   x2 = c(1,2,3,4,5,6,7,7,8,8,9,9,10,10,11,11,12,12),
> each of the 18 items now has the same weight as the others,
> and the "unweighted" sampling
>   mean(replicate(100,(sample(x2, 3
> now gives the mean of the 18 values (7.5); whereas -- as my
> calculation showed -- the effect of the "sequential" weighting is
> to bias the result slightly downwards (in your example; generally,
> in favour of the items with lower weights), since the way weighting
> works in sample() is not equivalent to replicating each item "weight"
> times.
> 
> The general problem, of sampling without replacement in such a way
> that for each item the probability that it is included in the sample
> is proportional to a pre-assigned weight ("sampling with probability
> proportional to size") is quite tricky and, for certain choices
> of weights, impossible. For a glimpse of what's inside the can of
> worms, have a look at the reference manual for the 'sampfling'
> package, in particular the function samprop():
> 
> http://www.stats.bris.ac.uk/R/web/packages/sampfling/sampfling.pdf
> 
> Ted.
> 
> On 23-Jul-09 20:56:43, Jim Bouldin wrote:
> > 
> > You are absolutely correct Ted.  When no weights are applied it doesn't
> > matter if you sample with or without replacement, because the
> > probability
> > of choosing any particular value is equally distributed among all such.
> > But when they're weighted unequally that's not the case.
> > 
> > It is also interesting to note that if the problem is set up slightly
> > differently, by say defining the vector x as:
> > x = c(1,2,3,4,5,6,7,7,8,8,9,9,10,10,11,11,12,12), effectively giving
> > the
> > same probability of selection for the 12 integers as before, the same
> > problem does not arise, or at least not as severely:
> > 
> >> x2
> >  [1]  1  2  3  4  5  6  7  8  9 10 11 12  7  8  9 10 11 12
> > 
> >> d = mean(replicate(100,(sample(x2, 3;d  # (1 million samples
> >> from
> > x2, of size 3; the mean should be 7.50)
> > [1] 7.499233
> > 
> >> e = mean(replicate(100,(sample(x2, 3, replace = TRUE;e  # (1
> > million samples from x2, of size 3; with replacement this time the mean
> > should still be 7.50)
> > [1] 7.502085
> > 
> >> d/e
> > [1] 0.9996198
> > 
> > Jim
> >> 
> >> To obtain the result you expected, you would need to explicitly
> >> specify "replace=TRUE", since the default for "replace" is FALSE.
> >> (though probably what you really intended was sampling without
> >> replacement).
> >> 
> >>  -- when
> >> replace=FALSE, the probability of inclusion of an element is not
> >> proportional to its weight in 'prob'.
> >> 
> >> The reason is that elements with higher weights are more likely
> >> to be chosen early on. This then knocks that higher weight out
> >> of the contest, making it more likely that elements with smaller
> >> weights will be sampled subsequently. Hence the mean of the sample
> >> will be biased slightly downwards, relative to the weighted mean
> >> of the values in x.
> >> 
> >>   table(replicate(100,(sample(x, 3
> >>   #  1  2  3  4  5  6
> >>   # 250235 250743 249603 250561 249828 249777
> >>   #  7  8  9 10 11 12
> >>   # 249780 250478 249591 249182 249625 250597
> >> 
> >> (so all nice equal frequencies)
> >> 
> >>   table(replicate(100,(sample(x, 3,prob=weights
> >>   #  1  2  3  4  5  6
> >>   # 174873 175398 174196 174445 173240 174110
> >>   #  7  8  9 10 11 12
> >>   # 325820 326140 325289 325098 325475 325916
> >> 
> >> Note that the frequencies of the values with weight=2 are a bit
> >> less than twice the frequencies of the values with weight=1:
> >> 
> >>   (325820+326140+325289+325098+325475+325916)/
> >> (174873+175398+174196+174445+173240+174110)
> >>   # [1] 
> >> 
> >> 
> >> In fact this is fairly easily caluclated. The possible combinations
> >> (in order of sampling) of the two weights, with their probabilities,
> >> are:
> >>  1s  2s
> >> ---
> >> 1 1 1   P =  6/18 *  5/17 *  4/163   0
> >> 1 1 2   P =  6/18 *  5/17 * 12/162   1
> >> 1 2 1   P =  6/18 * 12/17 *  5/152   1
> >> 1 2 2   P =  6/18 * 12/17 * 10/151   2
> >> 2 1 1   P = 12/18 *  6/16 *  5/152   1
> >> 2 1 2   P = 12/18 *  6/16 * 10/151   2
> >> 2 2 1   P = 12/18 * 10/16 *  6/141   2
> >> 2 2 2   P = 12/18 * 10/16 *  8/140   3
> >> 
> >> So the expected number of "weight=1" in the sample is
> >> 

[R] Help with Loop!

[Lars Bishop]

Dear experts,

I'm new in R and trying to learn by writing a version of the Perceptron
Algorithm. How can I tell in the code below to stop the iteration when the
condition in the "for loop" is not satisfied for all training examples?

Thanks in advance for your help!


## Generate a linearly separable data set in R2

sample <- as.data.frame(cbind(runif(n=100),runif(n=100)))
sample$Y <- ifelse(sample$V1>sample$V2,1,-1)
## Plot data;
attach(sample)
plot(V1, V2, type="p", pch=Y,main="Sample Data")

##Perceptron algorithm

sample_m <- as.matrix(sample)
w <- c(0,0); b <- 0; k <- 0; nu <- 1e-3
R <- max(sqrt(V1^2+V2^2))

repeat {
 for (i in 1:nrow(sample_m)){
if  (sample_m[i,3]*(t(w)%*%sample_m[i,1:2] + b) <= 0)
  w <- w + nu*sample_m[i,3]*sample_m[i,1:2]
  b <- b +nu*sample_m[i,3]*R^2
  k <- k+1
  cat(k, w, b, "\n") }
 }

[[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] ggplot2 : commands on one line vs two lines.

[John Kane]

I have just started using ggplot2 and I seem to be doing something stupid
in writing ggplot2 commands on more than one line.

In the example below the commands on one line are working fine, but as
soon as I put them on two lines I get an error.  Can any one point out
what I am doing wrong? It must be something blindingly simple.

Thanks

Example
#==
# sample data file
provs <- structure(list(Name = structure(c(8L, 11L, 6L, 2L, 9L, 1L, 12L,
3L, 13L, 5L, 4L, 7L, 10L), .Label = c("Alberta", "British Columbia",
"Manitoba", "New Brunswick", "Newfoundland", "Northwest Territories",
"Nova Scotia", "Nunavut", "Ontario", "Prince Edward Island",
"Québec", "Saskatchewan", "Yukon"), class = "factor"), Pcode = structure(c(8L,
11L, 7L, 2L, 9L, 1L, 12L, 3L, 13L, 5L, 4L, 6L, 10L), .Label = c("AB",
"BC", "MB", "NB", "NL", "NS", "NT", "NU", "ON", "PE", "QC", "SK",
"YT"), class = "factor"), Area = c(1936113L, 1365128L, 1183085L,
925186L, 917741L, 642317L, 591670L, 553556L, 474391L, 373872L,
71450L, 53338L, 5660L)), .Names = c("Name", "Pcode", "Area"),
class = "data.frame", row.names = c(NA,
-13L))

library(ggplot2)

p <- ggplot(provs, aes(x = 0, xend = Area,  y = Pcode, yend = Pcode))

# commmand below is on one line
pa <- p + geom_point(aes(Area, Pcode)) + geom_segment() + xlab("Area")

pa

pa1 <- p + geom_point(aes(Area, Pcode)) + geom_segment()
   + xlab("Area")
# Error in +xlab("Area") : invalid argument to unary operator


# commmand below is on one line
pb <- p + geom_point(aes(Area, Pcode)) + geom_segment() + labs(x="Area", 
y="Povs")
pb


pb1 <- p + geom_point(aes(Area, Pcode)) + geom_segment()
+ labs(x="Area", y="Povs")

# Error in +labs(x = "Area", y = "Povs") :
#  invalid argument to unary operator

#==

R version 2.9.1 (2009-06-26)
i386-pc-mingw32

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

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

other attached packages:
[1] ggplot2_0.8.3 reshape_0.8.3 plyr_0.1.9proto_0.3-8




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

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

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

Re: [R] Random # generator accuracy

[Ted Harding]

On 23-Jul-09 22:16:39, Thomas Lumley wrote:
> On Thu, 23 Jul 2009 ted.hard...@manchester.ac.uk wrote:
> 
>> The general problem, of sampling without replacement in such a way
>> that for each item the probability that it is included in the sample
>> is proportional to a pre-assigned weight ("sampling with probability
>> proportional to size") is quite tricky and, for certain choices
>> of weights, impossible. For a glimpse of what's inside the can of
>> worms, have a look at the reference manual for the 'sampfling'
>> package, in particular the function samprop():
>>
> 
> The 'sampling' package also provides a nice selection of ways to draw
> PPS samples without replacement.
> 
> -thomas

Thanks, Thomas. There is of course also the package 'pps'.
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 24-Jul-09   Time: 00:00:35
-- 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] alternative to rbind within a loop

[Don MacQueen]

Another approach that might be worth trying is to 
create an empty data frame with lots and lots of 
rows before looping, and then replace rather than 
append. Of course, this requires knowing at least 
approximately how many rows total you will have. 
This suggestion comes from the help page for 
read.table(), which says;


Using 'nrows', even as a mild over-estimate, will help memory usage.

You may be doing a lot of unnecessary processing 
if you are allowing your character variables to 
be automatically converted to factors. This would 
especially be the case if each data frame has new 
character values not in the previous ones, since 
more levels would be added to the factor 
variables each time a data frame is appended.


Another approach would be to concatenate the 
files outside of R (in unix, this would be the 
"cat" command) and then read the single large 
file into R. This can be controlled from within 
R, i.e., using the system() command. It can even 
be done without actually writing the extra file, 
with something like


  read.csv( pipe( 'cat *.csv') )

Despite those ideas, I like Greg Snow's approach; 
I'd try it before any of these.


Finally, if you really want to find out where the 
cpu time is being spent, look into profiling; see 
?Rprof.


-Don

At 3:53 PM -0400 7/23/09, Denis Chabot wrote:

Hi,

I often have to do this:

select a folder (directory) containing a few 
hundred data files in csv format (up to 1000 
files, in fact)


open each file, transform some character variables in date-tiime format

make into a dataframe (involves getting rid of a few variables I don't need

concatenate to the master dataframe that will 
eventually contain the data from all the files 
in the folder.


I use a loop going from 1 to the number of 
files. I have added a command to print an 
incrementing number to the R console each time 
the loop completes one iteration, to judge the 
speed of the process.


At the beginning, 3-4 files are processed each 
second. After a few hundred iterations it slows 
down to about 1 file per second. Before I reach 
the last file (898 in the case at hand), it has 
become much slower, about 1 file every 2-3 
seconds.


This progressive slowing down suggests the 
problem is linked to the size of the growing 
"master" dataframe that rbind combines with each 
new file.


In fact, the small script below confirms this as 
nothing at all happens within the loop but 
rbind. You can cut the size of this example not 
to waste to much of your time:



# create a dummy data.frame and copy it in a large number of csv files

test  <- file.path("test")

a <- 1:350
b <- rnorm(350,100,10)
c <- runif(350, 0, 100)
d <- month.name[runif(350,1,12)]

the.data <- data.frame(a,b,c,d)

for(i in 1:850){
write.csv(the.data, file=paste(test, "/file_", i, ".csv", sep=""))
}

# now lets make a single dataframe from all these csv files

all.files <- list.files(path=test,full.names=T,pattern=".csv")

new.data <- NULL

system.time({
for(i in all.files){
in.data <- read.csv(i)
	if (is.null(new.data)) {new.data = 
in.data} else {new.data = rbind(new.data, 
in.data)}

cat(paste(i, ", ", sep=""))
} # end for
}) # end system.time

utilisateur système  écoulé
156.206  44.859 202.150
This is with

sessionInfo()
R version 2.9.1 Patched (2009-07-16 r48939)
x86_64-apple-darwin9.7.0

locale:
fr_CA.UTF-8/fr_CA.UTF-8/C/C/fr_CA.UTF-8/fr_CA.UTF-8

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

other attached packages:
[1] doBy_3.7chron_2.3-30timeDate_290.84

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  Hmisc_3.5-2 
lattice_0.17-25 tools_2.9.1



Would it be better to somehow save all 850 files 
in one dataframe each, and then rbind them all 
in a single operation?


Can I combine all my files without using a loop? 
I've never quite mastered the "apply" family of 
functions but have not seen examples to read 
files.


Thanks in advance,

Denis Chabot

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



--
--
Don MacQueen
Environmental Protection Department
Lawrence Livermore National Laboratory
Livermore, CA, USA
925-423-1062

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

Re: [R] Assigning rank based on total count

[Remko Duursma]

Here's a way to get to your solution, but it's not very pretty:

testdfr <- data.frame(POB=c("Oregon","Oregon","Oregon","New
York","California","California"))

nstates <- length(unique(testdfr$POB))
testdfr$ POBR <- c(nstates:1)[table(testdfr$POB)][testdfr$POB]


greetings,
Remko



-
Remko Duursma
Post-Doctoral Fellow

Centre for Plants and the Environment
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753

Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia

Mobile: +61 (0)422 096908



On Fri, Jul 24, 2009 at 8:45 AM, ws wrote:
> Hi all,
>
> I am using ACS micro data (PUMS) with one of the columns as a
> factor for the place of birth (POBPF).  I would like to create
> a column (POBR) containing a rank
> corresponding to the place of the observation
> in the POBPF rankings.  For example,
> if a person is from Oregon, Oregon is
> the most popular Place of Birth, so each
> person who came from Oregon should have a 1
> in their POBR.  Each row from
> California would have POBR 2, etc.
>
> This way I can subset on the most populous sending
> places when doing analyses.
>
> For the life of me I can't figure this out.
> I think some combination of "sort()",
> "xtabs()" and "merge()" should work,
> but my brain seems to be fried.
>
> Tx
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Assigning rank based on total count

[ws]

Hi all,

I am using ACS micro data (PUMS) with one of the columns as a 
factor for the place of birth (POBPF).  I would like to create 
a column (POBR) containing a rank 
corresponding to the place of the observation 
in the POBPF rankings.  For example, 
if a person is from Oregon, Oregon is 
the most popular Place of Birth, so each 
person who came from Oregon should have a 1 
in their POBR.  Each row from 
California would have POBR 2, etc.  

This way I can subset on the most populous sending 
places when doing analyses.

For the life of me I can't figure this out.  
I think some combination of "sort()", 
"xtabs()" and "merge()" should work, 
but my brain seems to be fried.

Tx

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

Re: [R] Calculate weighted mean for each group

[David Freedman]

After you fix your data frame and if you don't using 2 packages, you might
try something like:

lib(plyr) #for 'by' processing
lib(Hmisc) # for its wtd.mean function
d=data.frame(x=c(15,12,3,10,10),g=c(1,1,2,2,3),w=c(2,1,5,2,5)) ; d
ddply(d,~g,function(df) wtd.mean(df$x,df$w))



milton ruser wrote:
> 
> try your first "reproducible line" first :-)
> 
> 
> 
> On Thu, Jul 23, 2009 at 5:18 PM, Alexis Maluendas
> wrote:
> 
>> Hi R experts,
>>
>> I need know how calculate a weighted mean by group in a data frame. I
>> have
>> tried with aggragate() function:
>>
>> data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5)) -> d
>> aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)
>>
>> Generating the following error:
>>
>> Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length
>>
>> Thanks in advance
>>
>>[[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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Calculate-weighted-mean-for-each-group-tp24634873p24635837.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] Random # generator accuracy

[Thomas Lumley]

On Thu, 23 Jul 2009 ted.hard...@manchester.ac.uk wrote:



The general problem, of sampling without replacement in such a way
that for each item the probability that it is included in the sample
is proportional to a pre-assigned weight ("sampling with probability
proportional to size") is quite tricky and, for certain choices
of weights, impossible. For a glimpse of what's inside the can of
worms, have a look at the reference manual for the 'sampfling'
package, in particular the function samprop():



The 'sampling' package also provides a nice selection of ways to draw PPS 
samples without replacement.

   -thomas

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

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

Re: [R] Random # generator accuracy

[Ted Harding]

Indeed, Jim! And that's why I said to read carefully what is said about
"prob" in '?sample':

  If 'replace' is false, these probabilities are applied
  sequentially, that is the probability of choosing the
  next item is proportional to the weights amongst the
  remaining items.

Whereas, if you replace x = c(1,2,3,4,5,6,7,8,9,110,11,12)
with the "weighted" equivalent, doubling up 7-12 as in your
  x2 = c(1,2,3,4,5,6,7,7,8,8,9,9,10,10,11,11,12,12),
each of the 18 items now has the same weight as the others,
and the "unweighted" sampling
  mean(replicate(100,(sample(x2, 3
now gives the mean of the 18 values (7.5); whereas -- as my
calculation showed -- the effect of the "sequential" weighting is
to bias the result slightly downwards (in your example; generally,
in favour of the items with lower weights), since the way weighting
works in sample() is not equivalent to replicating each item "weight"
times.

The general problem, of sampling without replacement in such a way
that for each item the probability that it is included in the sample
is proportional to a pre-assigned weight ("sampling with probability
proportional to size") is quite tricky and, for certain choices
of weights, impossible. For a glimpse of what's inside the can of
worms, have a look at the reference manual for the 'sampfling'
package, in particular the function samprop():

http://www.stats.bris.ac.uk/R/web/packages/sampfling/sampfling.pdf

Ted.

On 23-Jul-09 20:56:43, Jim Bouldin wrote:
> 
> You are absolutely correct Ted.  When no weights are applied it doesn't
> matter if you sample with or without replacement, because the
> probability
> of choosing any particular value is equally distributed among all such.
> But when they're weighted unequally that's not the case.
> 
> It is also interesting to note that if the problem is set up slightly
> differently, by say defining the vector x as:
> x = c(1,2,3,4,5,6,7,7,8,8,9,9,10,10,11,11,12,12), effectively giving
> the
> same probability of selection for the 12 integers as before, the same
> problem does not arise, or at least not as severely:
> 
>> x2
>  [1]  1  2  3  4  5  6  7  8  9 10 11 12  7  8  9 10 11 12
> 
>> d = mean(replicate(100,(sample(x2, 3;d  # (1 million samples
>> from
> x2, of size 3; the mean should be 7.50)
> [1] 7.499233
> 
>> e = mean(replicate(100,(sample(x2, 3, replace = TRUE;e  # (1
> million samples from x2, of size 3; with replacement this time the mean
> should still be 7.50)
> [1] 7.502085
> 
>> d/e
> [1] 0.9996198
> 
> Jim
>> 
>> To obtain the result you expected, you would need to explicitly
>> specify "replace=TRUE", since the default for "replace" is FALSE.
>> (though probably what you really intended was sampling without
>> replacement).
>> 
>>  -- when
>> replace=FALSE, the probability of inclusion of an element is not
>> proportional to its weight in 'prob'.
>> 
>> The reason is that elements with higher weights are more likely
>> to be chosen early on. This then knocks that higher weight out
>> of the contest, making it more likely that elements with smaller
>> weights will be sampled subsequently. Hence the mean of the sample
>> will be biased slightly downwards, relative to the weighted mean
>> of the values in x.
>> 
>>   table(replicate(100,(sample(x, 3
>>   #  1  2  3  4  5  6
>>   # 250235 250743 249603 250561 249828 249777
>>   #  7  8  9 10 11 12
>>   # 249780 250478 249591 249182 249625 250597
>> 
>> (so all nice equal frequencies)
>> 
>>   table(replicate(100,(sample(x, 3,prob=weights
>>   #  1  2  3  4  5  6
>>   # 174873 175398 174196 174445 173240 174110
>>   #  7  8  9 10 11 12
>>   # 325820 326140 325289 325098 325475 325916
>> 
>> Note that the frequencies of the values with weight=2 are a bit
>> less than twice the frequencies of the values with weight=1:
>> 
>>   (325820+326140+325289+325098+325475+325916)/
>> (174873+175398+174196+174445+173240+174110)
>>   # [1] 
>> 
>> 
>> In fact this is fairly easily caluclated. The possible combinations
>> (in order of sampling) of the two weights, with their probabilities,
>> are:
>>  1s  2s
>> ---
>> 1 1 1   P =  6/18 *  5/17 *  4/163   0
>> 1 1 2   P =  6/18 *  5/17 * 12/162   1
>> 1 2 1   P =  6/18 * 12/17 *  5/152   1
>> 1 2 2   P =  6/18 * 12/17 * 10/151   2
>> 2 1 1   P = 12/18 *  6/16 *  5/152   1
>> 2 1 2   P = 12/18 *  6/16 * 10/151   2
>> 2 2 1   P = 12/18 * 10/16 *  6/141   2
>> 2 2 2   P = 12/18 * 10/16 *  8/140   3
>> 
>> So the expected number of "weight=1" in the sample is
>> 
>>   3*(6/18 *  5/17 *  4/16)  + 2*(6/18 *  5/17 * 12/16) +
>>   2*(6/18 * 12/17 *  5/15)  + 1*(6/18 * 12/17 * 10/15) +
>>   2*(12/18 *  6/16 *  5/15) + 1*(12/18 *  6/16 * 10/15) +
>>   1*(12/18 * 10/16 *  6/14) + 0
>>   = 1.046218
>> 
>> Hence the expected number of

Re: [R] Creating a loop to read 3D dataset

[Greg Snow]

Something like this (untested) may work for you:

> for (i in 1:132) {
+  a <- clim[ , , i]
+  nm <- sprint('clim%03d.txt',i)
+  write.table(a,nm)
+ }

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Andrew Aldersley
> Sent: Thursday, July 23, 2009 2:38 PM
> To: r-help@r-project.org
> Subject: [R] Creating a loop to read 3D dataset
> 
> 
> Dear all,
> 
> I have in my possession a netcdf from which I want to extract some data
> files. I have used the "ncdf" package to read the netcdf file and used
> the "get.var.ncdf" function to identify the variable i wish to use. The
> data is in the form of a time-series of geographical data points that
> relate to climatology variables. As such I have a large data frame,
> "clim", of the following dimensions:
> 
> [720, 360, 132]
> 
> The 720 and 360 values relate to longitude and latitude points on the
> globe and there are 132 different time-series grids. What I wish to do
> is extract each of the individual time-series data grids and save them
> separately. So something along the lines of...
> 
> a <- clim[1:720, 1:360, 1]
> b <- clim[1:720, 1:360, 2]
> c <- clim[1:720, 1:360, 3]
> 
> ...and so on, all the way up to 132.
> 
> To do this manually for each different z-coordinate of "clim" and save
> them is possible but requires me to write a fairly large amount of code
> for what seems like a simple operation. I was wondering if anyone could
> give me some advice on setting up some sort of loop function to step
> through each z-coordinate and extract the appropriate 720x360 grid and
> then save this grid in the form of a text file. It is neccessary that
> each grid has a different name. I had a read around and think that a
> "for" loop would do the job, i.e. something in the form of...
> 
> for (i in 1:132)
> a <- clim[1:720, 1:360, i]
> 
> ...and then use the "write.table" function to save my grid as a text
> file. However, my problem is that I'm not sure how to structure my loop
> so that each 720x360 grid will be given a different name. I'm also
> unaware as to how to incorporate my different names into "write.table".
> Is it possible to set up a loop for this function and again step
> through each file and save them separately? Or should I use a different
> function to save my grids?
> 
> Thanks in advance for your help.
> 
> Andy
> 
> 
> 
> _
> 
> icons.
> 
>   [[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] Calculate weighted mean for each group

[Marc Schwartz]

On Jul 23, 2009, at 4:18 PM, Alexis Maluendas wrote:


Hi R experts,

I need know how calculate a weighted mean by group in a data frame.  
I have

tried with aggragate() function:

data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5))  
-> d

aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)

Generating the following error:

Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length

Thanks in advance



Did you not notice the error message when creating the data frame:

> d <-  
data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5))

Error in data.frame(x = c(15, 12, 3, 10, 10), g = c(1, 1, 1, 2, 2, 3,  :
  arguments imply differing number of rows: 5, 7

You have 5 elements in 'x' and 7 in each of 'g' and 'w'...

In addition, you are passing all 7 elements in d$w to each of the  
subsets created by d$g, hence you are getting the aggregate() error  
message.


This is one of those cases where you may be better served by using  
split() directly to break up the data frame into groups and then use  
sapply() over the subsets:


# I am adding data here to create the data frame
d <-  
data 
.frame(x=c(15,12,3,10,10,12,12),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5))


> d
   x g w
1 15 1 2
2 12 1 3
3  3 1 1
4 10 2 5
5 10 2 5
6 12 3 2
7 12 3 5

> split(d, d$g)
$`1`
   x g w
1 15 1 2
2 12 1 3
3  3 1 1

$`2`
   x g w
4 10 2 5
5 10 2 5

$`3`
   x g w
6 12 3 2
7 12 3 5



> sapply(split(d, d$g), function(x) weighted.mean(x$x, w = x$w))
   123
11.5 10.0 12.0


See ?split, which is used by tapply(), which in turn is used in  
aggregate().


HTH,

Marc Schwartz

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

Re: [R] Calculate weighted mean for each group

[Chuck Cleland]

On 7/23/2009 5:18 PM, Alexis Maluendas wrote:
> Hi R experts,
> 
> I need know how calculate a weighted mean by group in a data frame. I have
> tried with aggragate() function:
> 
> data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5)) -> d
> aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)
> 
> Generating the following error:
> 
> Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length

DF <- data.frame(x=c(15,12, 3,10,10,14,12),
 g=c( 1, 1, 1, 2, 2, 3, 3),
 w=c( 2, 3, 1, 5, 5, 2, 5))

sapply(split(DF, DF$g), function(x){weighted.mean(x$x, x$w)})
   123
11.5 10.0 12.57143

> Thanks in advance
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code. 

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

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

Re: [R] Calculate weighted mean for each group

[milton ruser]

try your first "reproducible line" first :-)



On Thu, Jul 23, 2009 at 5:18 PM, Alexis Maluendas wrote:

> Hi R experts,
>
> I need know how calculate a weighted mean by group in a data frame. I have
> tried with aggragate() function:
>
> data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5)) -> d
> aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)
>
> Generating the following error:
>
> Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length
>
> Thanks in advance
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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

Re: [R] Java to R interface

[peterburzhec]

I already posted this on another forum
I kept having the same problem. Turns out the issue was - I set my
pathwrong.
Make sure you have both, the path to R.dll and jri.dll in your path
CORRECTLY, and make sure that you don't put any spaces inbetween values in
the path.



madhura wrote:
> 
> The path to R/bin is in the Windows PATH variable.  Yet I get this
> error.
> 
> On Jun 6, 10:37 am, "Dumblauskas, Jerry"  suisse.com> wrote:
>> Try and make sure that R is in your windows Path variable
>>
>> I got your message when I first did this, but when I did the about it
>> then worked...
>>
>> ===­===
>> Please access the attached hyperlink for an important electronic
>> communications disclaimer:
>>
>> http://www.credit-suisse.com/legal/en/disclaimer_email_ib.html
>> ===­===
>>
>>         [[alternative HTML version deleted]]
>>
>> __
>> r-h...@r-project.org mailing
>> listhttps://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting
>> guidehttp://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Java-to-R-interface-tp17672475p24634848.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] rows missing after dataset loaded to R

[Greg Snow]

Some programs quote everything to be "safe", others only quote when needed.  
The only case that I know of that read.table and friends require quotes for is 
when a separator is inside of a string, for example if you are using spaces as 
the separator and have some names with spaces in them (e.g. "North Dakota"), 
without the quotes that would be seen as 2 fields, with the quotes it is a 
single field.  If the coma (,) is the separator and you have names (e.g. "Snow, 
Greg") then you would need the quotes.  If you don't have any cases of the 
separators other than where they are separating fields, then the quotes are 
probably not needed.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Rnewbie
> Sent: Thursday, July 23, 2009 11:02 AM
> To: r-help@r-project.org
> Subject: Re: [R] rows missing after dataset loaded to R
> 
> 
> Thank you very much for the reply. I checked the rows and it was the
> unbalanced " quote marks in some of the rows that caused the problem.
> Once I
> disabled quoting altogether, the problem is solved.
> 
> I have one more basic question.  I disabled quoting when loading the
> file to
> R, and all the columns consisting of characters with or without
> quotation
> marks displayed normally. What does quoting actually do in R?
> 
> --
> View this message in context: http://www.nabble.com/rows-missing-after-
> dataset-loaded-to-R-tp24625882p24630154.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Navigate to Index page of a package from R command prompt

[Steven McKinney]

> -Original Message-
> From: Marc Schwartz [mailto:marc_schwa...@me.com]
> Sent: Thursday, July 23, 2009 2:04 PM
> To: Steven McKinney
> Cc: R-help@r-project.org
> Subject: Re: [R] Navigate to Index page of a package from R command
> prompt
> 
> On Jul 23, 2009, at 3:30 PM, Steven McKinney wrote:
> 
> >
> > Hi all,
> >
> > Is there a way to navigate directly to the "Index" page of help
> > for a package?
> >
> > Here's my connundrum:
> >
> > I download and install package "foo".
> > I don't know what functions are in package "foo",
> > so I can't invoke the help for package "foo" via
> >> ?someFunction
> >
> > help(package = "foo")
> > pops up some non-hyperlinked information page, not
> > package foo's help Index page.
> >
> > If the package author kindly made a "foo" object or function
> > and put that in package "foo", then
> >> ?foo
> > works and yields a help page for package "foo".
> > Now at the bottom of the help page is a hyperlink "Index"
> > and I can click that to navigate to the main help Index page
> > (the page I really want to get to straight from the R
> > command line).
> >
> > I see that the link to "Index" for package "foo" appears always to be
> > (on my Mac)
> >
> file:///Library/Frameworks/R.framework/Resources/library/foo/html/00Ind
> ex.html
> >
> > e.g.
> >
> file:///Library/Frameworks/R.framework/Resources/library/cmprsk/html/00
> Index.html
> >
> >
> file:///Library/Frameworks/R.framework/Resources/library/utils/html/00I
> ndex.html
> >
> > Is there a command from the R listener that can take me directly to
> > this "00Index.html" page of help for package "foo"?
> >
> > something like
> >> help("00Index", package = "utils")
> >
> > (but this does not work)?
> >
> > Any info appreciated
> >
> > Best
> >
> > Steven McKinney
> 
> 
> Steven,
> 
> Once a package has been loaded with 'library(PackageName)', you can
> use .path.package("PackageName") to get the path to the main package
> installation folder and then append the remainder:
> 
> library(survival)
> 
>  > .path.package("survival")
> [1] "/Library/Frameworks/R.framework/Resources/library/survival"
> 
>  > file.path(.path.package("survival"), "html", "00Index.html")
> [1] "/Library/Frameworks/R.framework/Resources/library/survival/html/
> 00Index.html"
> 
> See ?.path.package and ?file.path for more information.
> 
> You can also use .find.package() for unloaded packages and that is
> described on the same help page as .path.package.
> 
> Finally, you can use browseURL() to bring the index file up in your
> browser:
> 
>browseURL(file.path(.path.package("survival"), "html",
> "00Index.html"))
> 
> HTH,
> 
> Marc Schwartz

Thanks Marc,

Certainly useful information.

I'm just baffled and confused as to why
> ?survival
or
> ?survival:Index
or some such incantation
doesn't just take me right to the index page 
for package "survival" in the default help browser 
for the R session.

Why does it take so much work to get to the 
index page?

> library("survival")
> browseURL(file.path(.path.package("survival"), "html", "00Index.html"))
does conceptually do what I'm looking for, but 
not in the default browser.

I'd like to file a feature request on R-devel if this
is not too naïve of a request.

?survival:Index 
or some such syntax should do what you illustrate
above, but in the default help mechanism.
I don't yet know enough about the ? and ?? machinery
to develop such a capability.  Any ideas?




Best

Steve McKinney

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

Re: [R] alternative to rbind within a loop

[Greg Snow]

Try something like (untested):

> mylist <- lapply(all.files, function(i) read.csv(i) )
> mydf <- do.call('rbind', mylist)

If all the csv files are conformable that rbind works on them (if the loop 
method works then that should be the case) then this will read in each file, 
store the data frames as a list, then rbind them all together.

It seems that this should be faster than the loop, but testing will be needed 
to be sure.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Denis Chabot
> Sent: Thursday, July 23, 2009 1:54 PM
> To: list R
> Subject: [R] alternative to rbind within a loop
> 
> Hi,
> 
> I often have to do this:
> 
> select a folder (directory) containing a few hundred data files in csv
> format (up to 1000 files, in fact)
> 
> open each file, transform some character variables in date-tiime format
> 
> make into a dataframe (involves getting rid of a few variables I don't
> need
> 
> concatenate to the master dataframe that will eventually contain the
> data from all the files in the folder.
> 
> I use a loop going from 1 to the number of files. I have added a
> command to print an incrementing number to the R console each time the
> loop completes one iteration, to judge the speed of the process.
> 
> At the beginning, 3-4 files are processed each second. After a few
> hundred iterations it slows down to about 1 file per second. Before I
> reach the last file (898 in the case at hand), it has become much
> slower, about 1 file every 2-3 seconds.
> 
> This progressive slowing down suggests the problem is linked to the
> size of the growing "master" dataframe that rbind combines with each
> new file.
> 
> In fact, the small script below confirms this as nothing at all
> happens within the loop but rbind. You can cut the size of this
> example not to waste to much of your time:
> 
> 
> # create a dummy data.frame and copy it in a large number of csv files
> 
> test  <- file.path("test")
> 
> a <- 1:350
> b <- rnorm(350,100,10)
> c <- runif(350, 0, 100)
> d <- month.name[runif(350,1,12)]
> 
> the.data <- data.frame(a,b,c,d)
> 
> for(i in 1:850){
>   write.csv(the.data, file=paste(test, "/file_", i, ".csv",
> sep=""))
> }
> 
> # now lets make a single dataframe from all these csv files
> 
> all.files <- list.files(path=test,full.names=T,pattern=".csv")
> 
> new.data <- NULL
> 
> system.time({
>   for(i in all.files){
>   in.data <- read.csv(i)
>   if (is.null(new.data)) {new.data = in.data} else {new.data =
> rbind(new.data, in.data)}
>   cat(paste(i, ", ", sep=""))
> } # end for
> }) # end system.time
> 
> utilisateur système  écoulé
>  156.206  44.859 202.150
> This is with
> 
> sessionInfo()
> R version 2.9.1 Patched (2009-07-16 r48939)
> x86_64-apple-darwin9.7.0
> 
> locale:
> fr_CA.UTF-8/fr_CA.UTF-8/C/C/fr_CA.UTF-8/fr_CA.UTF-8
> 
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
> 
> other attached packages:
> [1] doBy_3.7chron_2.3-30timeDate_290.84
> 
> loaded via a namespace (and not attached):
> [1] cluster_1.12.0  grid_2.9.1  Hmisc_3.5-2 lattice_0.17-25
> tools_2.9.1
> 
> 
> Would it be better to somehow save all 850 files in one dataframe
> each, and then rbind them all in a single operation?
> 
> Can I combine all my files without using a loop? I've never quite
> mastered the "apply" family of functions but have not seen examples to
> read files.
> 
> Thanks in advance,
> 
> Denis Chabot
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Navigate to Index page of a package from R command prompt

[Marc Schwartz]

On Jul 23, 2009, at 3:30 PM, Steven McKinney wrote:



Hi all,

Is there a way to navigate directly to the "Index" page of help
for a package?

Here's my connundrum:

I download and install package "foo".
I don't know what functions are in package "foo",
so I can't invoke the help for package "foo" via

?someFunction


help(package = "foo")
pops up some non-hyperlinked information page, not
package foo's help Index page.

If the package author kindly made a "foo" object or function
and put that in package "foo", then

?foo

works and yields a help page for package "foo".
Now at the bottom of the help page is a hyperlink "Index"
and I can click that to navigate to the main help Index page
(the page I really want to get to straight from the R
command line).

I see that the link to "Index" for package "foo" appears always to be
(on my Mac)
file:///Library/Frameworks/R.framework/Resources/library/foo/html/00Index.html

e.g.
file:///Library/Frameworks/R.framework/Resources/library/cmprsk/html/00Index.html

file:///Library/Frameworks/R.framework/Resources/library/utils/html/00Index.html

Is there a command from the R listener that can take me directly to
this "00Index.html" page of help for package "foo"?

something like

help("00Index", package = "utils")


(but this does not work)?

Any info appreciated

Best

Steven McKinney



Steven,

Once a package has been loaded with 'library(PackageName)', you can  
use .path.package("PackageName") to get the path to the main package  
installation folder and then append the remainder:


library(survival)

> .path.package("survival")
[1] "/Library/Frameworks/R.framework/Resources/library/survival"

> file.path(.path.package("survival"), "html", "00Index.html")
[1] "/Library/Frameworks/R.framework/Resources/library/survival/html/ 
00Index.html"


See ?.path.package and ?file.path for more information.

You can also use .find.package() for unloaded packages and that is  
described on the same help page as .path.package.


Finally, you can use browseURL() to bring the index file up in your  
browser:


  browseURL(file.path(.path.package("survival"), "html",  
"00Index.html"))


HTH,

Marc Schwartz

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Calculate weighted mean for each group

[Alexis Maluendas]

Hi R experts,

I need know how calculate a weighted mean by group in a data frame. I have
tried with aggragate() function:

data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5)) -> d
aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)

Generating the following error:

Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length

Thanks in advance

[[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] Automatic differentiation in R

[Hans W. Borchers]

Having given a lecture on "Numerical Derivatives" just a short time ago, I
would like to mention the following:

Many functions, especially in engineering, are not available as formulas
built simply from arithmetical operators and elementary functions.  They are
provided as intricate procedures, as results from simulation runs, or as
return values from optimization routines.

In all these cases, symbolic differentiation is not a possibility at all. 
Still, for further processing there will be a need to differentiate these
functions.

As I understand the literature, Automated Differentiations (AD) is mostly
meant for these kinds of applications(*).  Another push for AD came from
increased interest in the "complex-step derivative approximation" where
first an analytic continuation -- often nontrivial -- has to be established,
before the derivation can be computed.

It may be different in statistical tasks.  But a general automated
differentiation in R may have to take these kinds of applications into
account.

Regards
Hans Werner

(*) See the "Community Portal for AD"  and the definition
of automated differentiation there as an example.


nashjc wrote:
> 
> [...]
> 
> The clear issue in my mind is that users who need gradients/Jacobians 
> for R want to be able to send a function X to some process that will 
> return another function gradX or JacX that computes analytic 
> derivatives. This has to be "easy", which implies a very simple command 
> or GUI interface. I am pretty certain the users have almost no interest 
> in the mechanism, as long as it works. Currently, most use numerical 
> derivatives, not realizing the very large time penalty and quite large 
> loss in accuracy that can compromise some optimization and differential 
> equation codes. I'll try to prepare a few examples to illustrate this 
> and post them somewhere in the next few weeks. Time, as always, ...
> 
> However, the topic does appear to be on the table.
> 
> JN
> 

-- 
View this message in context: 
http://www.nabble.com/Automatic-differentiation-in-R-tp24602805p24634481.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] Random # generator accuracy

[Jim Bouldin]

You are absolutely correct Ted.  When no weights are applied it doesn't
matter if you sample with or without replacement, because the probability
of choosing any particular value is equally distributed among all such. 
But when they're weighted unequally that's not the case.

It is also interesting to note that if the problem is set up slightly
differently, by say defining the vector x as:
x = c(1,2,3,4,5,6,7,7,8,8,9,9,10,10,11,11,12,12), effectively giving the
same probability of selection for the 12 integers as before, the same
problem does not arise, or at least not as severely:

> x2
 [1]  1  2  3  4  5  6  7  8  9 10 11 12  7  8  9 10 11 12

> d = mean(replicate(100,(sample(x2, 3;d  # (1 million samples from
x2, of size 3; the mean should be 7.50)
[1] 7.499233

> e = mean(replicate(100,(sample(x2, 3, replace = TRUE;e  # (1
million samples from x2, of size 3; with replacement this time the mean
should still be 7.50)
[1] 7.502085

> d/e
[1] 0.9996198

Jim
> 
> To obtain the result you expected, you would need to explicitly
> specify "replace=TRUE", since the default for "replace" is FALSE.
> (though probably what you really intended was sampling without
> replacement).
> 
> Read carefully what is said about "prob" in '?sample' -- when
> replace=FALSE, the probability of inclusion of an element is not
> proportional to its weight in 'prob'.
> 
> The reason is that elements with higher weights are more likely
> to be chosen early on. This then knocks that higher weight out
> of the contest, making it more likely that elements with smaller
> weights will be sampled subsequently. Hence the mean of the sample
> will be biased slightly downwards, relative to the weighted mean
> of the values in x.
> 
>   table(replicate(100,(sample(x, 3
>   #  1  2  3  4  5  6
>   # 250235 250743 249603 250561 249828 249777
>   #  7  8  9 10 11 12
>   # 249780 250478 249591 249182 249625 250597
> 
> (so all nice equal frequencies)
> 
>   table(replicate(100,(sample(x, 3,prob=weights
>   #  1  2  3  4  5  6
>   # 174873 175398 174196 174445 173240 174110
>   #  7  8  9 10 11 12
>   # 325820 326140 325289 325098 325475 325916
> 
> Note that the frequencies of the values with weight=2 are a bit
> less than twice the frequencies of the values with weight=1:
> 
>   (325820+326140+325289+325098+325475+325916)/
> (174873+175398+174196+174445+173240+174110)
>   # [1] 
> 
> 
> In fact this is fairly easily caluclated. The possible combinations
> (in order of sampling) of the two weights, with their probabilities,
> are:
>  1s  2s
> ---
> 1 1 1   P =  6/18 *  5/17 *  4/163   0
> 1 1 2   P =  6/18 *  5/17 * 12/162   1
> 1 2 1   P =  6/18 * 12/17 *  5/152   1
> 1 2 2   P =  6/18 * 12/17 * 10/151   2
> 2 1 1   P = 12/18 *  6/16 *  5/152   1
> 2 1 2   P = 12/18 *  6/16 * 10/151   2
> 2 2 1   P = 12/18 * 10/16 *  6/141   2
> 2 2 2   P = 12/18 * 10/16 *  8/140   3
> 
> So the expected number of "weight=1" in the sample is
> 
>   3*(6/18 *  5/17 *  4/16)  + 2*(6/18 *  5/17 * 12/16) +
>   2*(6/18 * 12/17 *  5/15)  + 1*(6/18 * 12/17 * 10/15) +
>   2*(12/18 *  6/16 *  5/15) + 1*(12/18 *  6/16 * 10/15) +
>   1*(12/18 * 10/16 *  6/14) + 0
>   = 1.046218
> 
> Hence the expected number of "weight=2" in the sample is
> 
>   3 - 1.046218 = 1.953782
> 
> and their ratio 1.953782/1.046218 = 1.867471
> 
> Compare this with the value 1.867351 (above) obtained by simulation!
> 
> Ted.
> 
> 
> E-Mail: (Ted Harding) 
> Fax-to-email: +44 (0)870 094 0861
> Date: 23-Jul-09   Time: 21:05:07
> -- XFMail --
> 

Jim Bouldin, PhD
Research Ecologist
Department of Plant Sciences, UC Davis
Davis CA, 95616
530-554-1740

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

Re: [R] mathematical notation in R

[Steven McKinney]

Hi Mary,

Something such as

> plot(1,1)
> title(expression(paste("LBAuo = -", infinity)))

Is this what you're after?

Best

Steve McKinney


From: Mary A. Marion [mms...@comcast.net]
Sent: July 23, 2009 4:23 PM
To: Steven McKinney
Subject: Re: [R] mathematical notation in R

Hi,

LBAuo <- -Inf was helpful but what I really would like to do is use the "Greek" 
expression for
infinity that is the number 8 on its side.

syntax   meaning


infinityinfinity symbol
was found using ?plotmath but I don't know how to use this information.  Any 
ideas?

Thank you.

Sincerely,
Mary A. Marion


Steven McKinney wrote:

Does this approach what you're looking for?
See
?is.finite

for more info



LBAuo<- -
LBAuo


[1] -


if (LBAuo <= -) LBAuo <- -Inf
LBAuo


[1] -Inf



HTH
Steven McKinney



-Original Message-
From: r-help-boun...@r-project.org 
[mailto:r-help-boun...@r-
project.org] On Behalf Of Mary A. Marion
Sent: Wednesday, July 22, 2009 8:00 PM
To: r-help@r-project.org
Subject: [R] mathematical notation in R

Hello,

I have an R function which includes the statement LBAuo<- -.
Rather
than printing out - I want it to print out - infinity.
As of yet I have not been able to do that.  I am not in a graphics
window just outputting a set of variables.  Can you help?
Thanks.

Sincerely,
Mary A. Marion

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

[Andrew Aldersley]

Dear all,

I have in my possession a netcdf from which I want to extract some data files. 
I have used the "ncdf" package to read the netcdf file and used the 
"get.var.ncdf" function to identify the variable i wish to use. The data is in 
the form of a time-series of geographical data points that relate to 
climatology variables. As such I have a large data frame, "clim", of the 
following dimensions:

[720, 360, 132]

The 720 and 360 values relate to longitude and latitude points on the globe and 
there are 132 different time-series grids. What I wish to do is extract each of 
the individual time-series data grids and save them separately. So something 
along the lines of...

a <- clim[1:720, 1:360, 1]
b <- clim[1:720, 1:360, 2]
c <- clim[1:720, 1:360, 3]

...and so on, all the way up to 132.

To do this manually for each different z-coordinate of "clim" and save them is 
possible but requires me to write a fairly large amount of code for what seems 
like a simple operation. I was wondering if anyone could give me some advice on 
setting up some sort of loop function to step through each z-coordinate and 
extract the appropriate 720x360 grid and then save this grid in the form of a 
text file. It is neccessary that each grid has a different name. I had a read 
around and think that a "for" loop would do the job, i.e. something in the form 
of...

for (i in 1:132)
a <- clim[1:720, 1:360, i]

...and then use the "write.table" function to save my grid as a text file. 
However, my problem is that I'm not sure how to structure my loop so that each 
720x360 grid will be given a different name. I'm also unaware as to how to 
incorporate my different names into "write.table". Is it possible to set up a 
loop for this function and again step through each file and save them 
separately? Or should I use a different function to save my grids?

Thanks in advance for your help.

Andy



_

icons.

[[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] Navigate to Index page of a package from R command prompt

[Steven McKinney]

Hi all,

Is there a way to navigate directly to the "Index" page of help
for a package?

Here's my connundrum:

I download and install package "foo".
I don't know what functions are in package "foo",
so I can't invoke the help for package "foo" via
> ?someFunction

help(package = "foo")
pops up some non-hyperlinked information page, not
package foo's help Index page.

If the package author kindly made a "foo" object or function
and put that in package "foo", then 
> ?foo
works and yields a help page for package "foo".
Now at the bottom of the help page is a hyperlink "Index"
and I can click that to navigate to the main help Index page
(the page I really want to get to straight from the R
command line).

I see that the link to "Index" for package "foo" appears always to be
(on my Mac)
file:///Library/Frameworks/R.framework/Resources/library/foo/html/00Index.html

e.g.
file:///Library/Frameworks/R.framework/Resources/library/cmprsk/html/00Index.html

file:///Library/Frameworks/R.framework/Resources/library/utils/html/00Index.html

Is there a command from the R listener that can take me directly to
this "00Index.html" page of help for package "foo"?

something like
> help("00Index", package = "utils")

(but this does not work)?

Any info appreciated

Best

Steven McKinney

Statistician
Molecular Oncology and Breast Cancer Program
British Columbia Cancer Research Centre

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

[Jorge Ivan Velez]

Dear Alberto,
Try this:

Res <- apply( bm[,colnames(bm) %in% d], 1, function(x) sum( as.numeric(x) )
)
cbind(bm, Res)

HTH,

Jorge


On Thu, Jul 23, 2009 at 3:31 PM, Alberto Lora M wrote:

>  Dear R project group
>
> I have the following problem
>
> Let suppose the following data:
>
>
>
> b<-c("0","1","1","1","0","1","0","0","1","0","1","1","1","0","1","1","0","1","1","0")
> bm<-matrix(b,ncol=4)
> colnames(bm)<-c("F1", "F2", "F3", "F4")
> d<-c("F1","F4")
>
> For the matrix bm i need to create a fifth column. 1This column should
> contain the sums  of values where the columname is equal to the vector d.
> 2. I need to get the maximum and minimum value of this last column
> The matrix result should be like this
>
> *F1*  F2  F3  *F4 Res*
> [1,] *"0"* "1" "1" "*1" "1"
> *[2,] "1" "0" "1" "0"  *"1"*
> [3,] "1" "0" "1" "1"  *"2"*
> [4,] "1" "1" "0" "1"  *"2"*
> [5,] "0" "0" "1" "0"  *"0"*
>
> Thxs again for your help
>
>
> --
> Alberto Lora Michiels
> Rue du Progrès,  6B
> 7860 Lessines
> GSM 32(0)496659457
>
>[[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

[[alternative HTML version deleted]]

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

Re: [R] Random # generator accuracy

[Ted Harding]

OOPS! The result of a calculation below somehow got omitted!

  (325820+326140+325289+325098+325475+325916)/
(174873+175398+174196+174445+173240+174110)
  # [1] 1.867351

to be compared (as at the end) with the ratio 1.867471 of the
expected number of "weight=2" to expected number of "weight=1".
Sorry.
Ted.

On 23-Jul-09 20:05:09, Ted Harding wrote:
> On 23-Jul-09 17:59:56, Jim Bouldin wrote:
>> Dan Nordlund wrote:
>> "It would be necessary to see the code for your 'brief test'
>> before anyone could meaningfully comment on your results.
>> But your results for a single test could have been a valid
>> "random" result."
>> 
>> I've re-created what I did below.  The problem appears to be
>> with the weighting process: the unweighted sample came out much
>> closer to the actual than the weighted sample (>1% error) did. 
>> Comments?
>> Jim
>> 
>>> x
>>  [1]  1  2  3  4  5  6  7  8  9 10 11 12
>>> weights
>>  [1] 1 1 1 1 1 1 2 2 2 2 2 2
>> 
>>> a = mean(replicate(100,(sample(x, 3, prob = weights;a  # (1
>> million samples from x, of size 3, weighted by "weights"; the mean
>> should
>> be 7.50)
>> [1] 7.406977
>>> 7.406977/7.5
>> [1] 0.987597
>> 
>>> b = mean(replicate(100,(sample(x, 3;b  # (1 million samples
>>> from
>> x, of size 3, not weighted this time; the mean should be 6.50)
>> [1] 6.501477
>>> 6.501477/6.5
>> [1] 1.000227
>> 
>> Jim Bouldin, PhD
> 
> To obtain the result you expected, you would need to explicitly
> specify "replace=TRUE", since the default for "replace" is FALSE.
> (though probably what you really intended was sampling without
> replacement).
> 
> Read carefully what is said about "prob" in '?sample' -- when
> replace=FALSE, the probability of inclusion of an element is not
> proportional to its weight in 'prob'.
> 
> The reason is that elements with higher weights are more likely
> to be chosen early on. This then knocks that higher weight out
> of the contest, making it more likely that elements with smaller
> weights will be sampled subsequently. Hence the mean of the sample
> will be biased slightly downwards, relative to the weighted mean
> of the values in x.
> 
>   table(replicate(100,(sample(x, 3
>   #  1  2  3  4  5  6
>   # 250235 250743 249603 250561 249828 249777
>   #  7  8  9 10 11 12
>   # 249780 250478 249591 249182 249625 250597
> 
> (so all nice equal frequencies)
> 
>   table(replicate(100,(sample(x, 3,prob=weights
>   #  1  2  3  4  5  6
>   # 174873 175398 174196 174445 173240 174110
>   #  7  8  9 10 11 12
>   # 325820 326140 325289 325098 325475 325916
> 
> Note that the frequencies of the values with weight=2 are a bit
> less than twice the frequencies of the values with weight=1:
> 
>   (325820+326140+325289+325098+325475+325916)/
> (174873+175398+174196+174445+173240+174110)
>   # [1] 
> 
> 
> In fact this is fairly easily caluclated. The possible combinations
> (in order of sampling) of the two weights, with their probabilities,
> are:
>  1s  2s
> ---
> 1 1 1   P =  6/18 *  5/17 *  4/163   0
> 1 1 2   P =  6/18 *  5/17 * 12/162   1
> 1 2 1   P =  6/18 * 12/17 *  5/152   1
> 1 2 2   P =  6/18 * 12/17 * 10/151   2
> 2 1 1   P = 12/18 *  6/16 *  5/152   1
> 2 1 2   P = 12/18 *  6/16 * 10/151   2
> 2 2 1   P = 12/18 * 10/16 *  6/141   2
> 2 2 2   P = 12/18 * 10/16 *  8/140   3
> 
> So the expected number of "weight=1" in the sample is
> 
>   3*(6/18 *  5/17 *  4/16)  + 2*(6/18 *  5/17 * 12/16) +
>   2*(6/18 * 12/17 *  5/15)  + 1*(6/18 * 12/17 * 10/15) +
>   2*(12/18 *  6/16 *  5/15) + 1*(12/18 *  6/16 * 10/15) +
>   1*(12/18 * 10/16 *  6/14) + 0
>   = 1.046218
> 
> Hence the expected number of "weight=2" in the sample is
> 
>   3 - 1.046218 = 1.953782
> 
> and their ratio 1.953782/1.046218 = 1.867471
> 
> Compare this with the value 1.867351 (above) obtained by simulation!
> 
> Ted.
> 
> 
> E-Mail: (Ted Harding) 
> Fax-to-email: +44 (0)870 094 0861
> Date: 23-Jul-09   Time: 21:05:07
> -- XFMail --


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 23-Jul-09   Time: 21:15:52
-- 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] rows missing after dataset loaded to R

[Rnewbie]

Thank you very much for the reply. I checked the rows and it was the
unbalanced " quote marks in some of the rows that caused the problem. Once I
disabled quoting altogether, the problem is solved.

I have one more basic question.  I disabled quoting when loading the file to
R, and all the columns consisting of characters with or without quotation
marks displayed normally. What does quoting actually do in R?

-- 
View this message in context: 
http://www.nabble.com/rows-missing-after-dataset-loaded-to-R-tp24625882p24630154.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] Question about the lars package

[Steve Lianoglou]

Hi,

On Jul 22, 2009, at 6:18 PM, lulu9797 wrote:

The returned values of lars function include R squares along the  
variable

selection path.


Correct.


However, such values are always slightly different from the
R squares returned by the regression function lm using the same  
models.

Anyone know the reasons?


How are you comparing the models from lars vs. lm?

Are you just using the non-zero coefs you get from lars in setting up  
your formula for lm (or something)?


How different is "slightly different"?

No idea/speculation: if you're comparing the "dropped coefficients"  
models to each other, perhaps the coefs on your predictors are  
slightly different in lars vs. lm since you've got the l1-penalty on  
them?


Perhaps a guru will respond with better insight.


Very important, and needs quick answers.


As a point of etiquette, I reckon most people don't really respond too  
well to requests like these since posting to this mailing list is  
essentially asking for free help. If it's so urgent/important to you,  
you can get some professional-for-hire to drop whatever s/he is doing  
at the moment and work it all out for you.


By the by, you might want to look at the glmnet package if you find  
yourself using lars often. Setting the glmnet alpha parameter to 1  
essentially gives you the lasso regularization path and I've found it  
to be quite a bit faster than lars.


-steve

--
Steve Lianoglou
Graduate Student: Physiology, Biophysics and Systems Biology
Weill Medical College of Cornell University

Contact Info: http://cbio.mskcc.org/~lianos/contact

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

Re: [R] Help

[Greg Snow]

Try:

> cbind( bm, res=apply(bm[,d],1 , sum) )

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Alberto Lora M
> Sent: Thursday, July 23, 2009 1:32 PM
> To: r-h...@stat.math.ethz.ch; R mailing list
> Subject: [R] Help
> 
>  Dear R project group
> 
> I have the following problem
> 
> Let suppose the following data:
> 
> 
> b<-
> c("0","1","1","1","0","1","0","0","1","0","1","1","1","0","1","1","0","
> 1","1","0")
> bm<-matrix(b,ncol=4)
> colnames(bm)<-c("F1", "F2", "F3", "F4")
> d<-c("F1","F4")
> 
> For the matrix bm i need to create a fifth column. 1This column should
> contain the sums  of values where the columname is equal to the vector
> d.
> 2. I need to get the maximum and minimum value of this last column
> The matrix result should be like this
> 
>  *F1*  F2  F3  *F4 Res*
> [1,] *"0"* "1" "1" "*1" "1"
> *[2,] "1" "0" "1" "0"  *"1"*
> [3,] "1" "0" "1" "1"  *"2"*
> [4,] "1" "1" "0" "1"  *"2"*
> [5,] "0" "0" "1" "0"  *"0"*
> 
> Thxs again for your help
> 
> 
> --
> Alberto Lora Michiels
> Rue du Progrès,  6B
> 7860 Lessines
> GSM 32(0)496659457
> 
>   [[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] alternative to rbind within a loop

[Denis Chabot]

Hi,

I often have to do this:

select a folder (directory) containing a few hundred data files in csv  
format (up to 1000 files, in fact)


open each file, transform some character variables in date-tiime format

make into a dataframe (involves getting rid of a few variables I don't  
need


concatenate to the master dataframe that will eventually contain the  
data from all the files in the folder.


I use a loop going from 1 to the number of files. I have added a  
command to print an incrementing number to the R console each time the  
loop completes one iteration, to judge the speed of the process.


At the beginning, 3-4 files are processed each second. After a few  
hundred iterations it slows down to about 1 file per second. Before I  
reach the last file (898 in the case at hand), it has become much  
slower, about 1 file every 2-3 seconds.


This progressive slowing down suggests the problem is linked to the  
size of the growing "master" dataframe that rbind combines with each  
new file.


In fact, the small script below confirms this as nothing at all  
happens within the loop but rbind. You can cut the size of this  
example not to waste to much of your time:



# create a dummy data.frame and copy it in a large number of csv files

test  <- file.path("test")

a <- 1:350
b <- rnorm(350,100,10)
c <- runif(350, 0, 100)
d <- month.name[runif(350,1,12)]

the.data <- data.frame(a,b,c,d)

for(i in 1:850){
write.csv(the.data, file=paste(test, "/file_", i, ".csv", sep=""))
}

# now lets make a single dataframe from all these csv files

all.files <- list.files(path=test,full.names=T,pattern=".csv")

new.data <- NULL

system.time({
for(i in all.files){
in.data <- read.csv(i)
	if (is.null(new.data)) {new.data = in.data} else {new.data =  
rbind(new.data, in.data)}

cat(paste(i, ", ", sep=""))
} # end for
}) # end system.time

utilisateur système  écoulé
156.206  44.859 202.150
This is with

sessionInfo()
R version 2.9.1 Patched (2009-07-16 r48939)
x86_64-apple-darwin9.7.0

locale:
fr_CA.UTF-8/fr_CA.UTF-8/C/C/fr_CA.UTF-8/fr_CA.UTF-8

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

other attached packages:
[1] doBy_3.7chron_2.3-30timeDate_290.84

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1  Hmisc_3.5-2 lattice_0.17-25  
tools_2.9.1



Would it be better to somehow save all 850 files in one dataframe  
each, and then rbind them all in a single operation?


Can I combine all my files without using a loop? I've never quite  
mastered the "apply" family of functions but have not seen examples to  
read files.


Thanks in advance,

Denis Chabot

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

[Alberto Lora M]

 Dear R project group

I have the following problem

Let suppose the following data:


b<-c("0","1","1","1","0","1","0","0","1","0","1","1","1","0","1","1","0","1","1","0")
bm<-matrix(b,ncol=4)
colnames(bm)<-c("F1", "F2", "F3", "F4")
d<-c("F1","F4")

For the matrix bm i need to create a fifth column. 1This column should
contain the sums  of values where the columname is equal to the vector d.
2. I need to get the maximum and minimum value of this last column
The matrix result should be like this

 *F1*  F2  F3  *F4 Res*
[1,] *"0"* "1" "1" "*1" "1"
*[2,] "1" "0" "1" "0"  *"1"*
[3,] "1" "0" "1" "1"  *"2"*
[4,] "1" "1" "0" "1"  *"2"*
[5,] "0" "0" "1" "0"  *"0"*

Thxs again for your help


-- 
Alberto Lora Michiels
Rue du Progrès,  6B
7860 Lessines
GSM 32(0)496659457

[[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] Random # generator accuracy

[Ted Harding]

On 23-Jul-09 17:59:56, Jim Bouldin wrote:
> Dan Nordlund wrote:
> "It would be necessary to see the code for your 'brief test'
> before anyone could meaningfully comment on your results.
> But your results for a single test could have been a valid
> "random" result."
> 
> I've re-created what I did below.  The problem appears to be
> with the weighting process: the unweighted sample came out much
> closer to the actual than the weighted sample (>1% error) did. 
> Comments?
> Jim
> 
>> x
>  [1]  1  2  3  4  5  6  7  8  9 10 11 12
>> weights
>  [1] 1 1 1 1 1 1 2 2 2 2 2 2
> 
>> a = mean(replicate(100,(sample(x, 3, prob = weights;a  # (1
> million samples from x, of size 3, weighted by "weights"; the mean
> should
> be 7.50)
> [1] 7.406977
>> 7.406977/7.5
> [1] 0.987597
> 
>> b = mean(replicate(100,(sample(x, 3;b  # (1 million samples
>> from
> x, of size 3, not weighted this time; the mean should be 6.50)
> [1] 6.501477
>> 6.501477/6.5
> [1] 1.000227
> 
> Jim Bouldin, PhD

To obtain the result you expected, you would need to explicitly
specify "replace=TRUE", since the default for "replace" is FALSE.
(though probably what you really intended was sampling without
replacement).

Read carefully what is said about "prob" in '?sample' -- when
replace=FALSE, the probability of inclusion of an element is not
proportional to its weight in 'prob'.

The reason is that elements with higher weights are more likely
to be chosen early on. This then knocks that higher weight out
of the contest, making it more likely that elements with smaller
weights will be sampled subsequently. Hence the mean of the sample
will be biased slightly downwards, relative to the weighted mean
of the values in x.

  table(replicate(100,(sample(x, 3
  #  1  2  3  4  5  6
  # 250235 250743 249603 250561 249828 249777
  #  7  8  9 10 11 12
  # 249780 250478 249591 249182 249625 250597

(so all nice equal frequencies)

  table(replicate(100,(sample(x, 3,prob=weights
  #  1  2  3  4  5  6
  # 174873 175398 174196 174445 173240 174110
  #  7  8  9 10 11 12
  # 325820 326140 325289 325098 325475 325916

Note that the frequencies of the values with weight=2 are a bit
less than twice the frequencies of the values with weight=1:

  (325820+326140+325289+325098+325475+325916)/
(174873+175398+174196+174445+173240+174110)
  # [1] 


In fact this is fairly easily caluclated. The possible combinations
(in order of sampling) of the two weights, with their probabilities,
are:
 1s  2s
---
1 1 1   P =  6/18 *  5/17 *  4/163   0
1 1 2   P =  6/18 *  5/17 * 12/162   1
1 2 1   P =  6/18 * 12/17 *  5/152   1
1 2 2   P =  6/18 * 12/17 * 10/151   2
2 1 1   P = 12/18 *  6/16 *  5/152   1
2 1 2   P = 12/18 *  6/16 * 10/151   2
2 2 1   P = 12/18 * 10/16 *  6/141   2
2 2 2   P = 12/18 * 10/16 *  8/140   3

So the expected number of "weight=1" in the sample is

  3*(6/18 *  5/17 *  4/16)  + 2*(6/18 *  5/17 * 12/16) +
  2*(6/18 * 12/17 *  5/15)  + 1*(6/18 * 12/17 * 10/15) +
  2*(12/18 *  6/16 *  5/15) + 1*(12/18 *  6/16 * 10/15) +
  1*(12/18 * 10/16 *  6/14) + 0
  = 1.046218

Hence the expected number of "weight=2" in the sample is

  3 - 1.046218 = 1.953782

and their ratio 1.953782/1.046218 = 1.867471

Compare this with the value 1.867351 (above) obtained by simulation!

Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 23-Jul-09   Time: 21:05:07
-- 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] I might be dumb : a simple question about "foreach"

[Daya Atapattu]

This problem is due to a bug in the iterators package.
It has been fixed in version 1.0.2 of iterators, which is
on CRAN now.

- Daya

Daya Atapattu
Revolution Computing


Olivier ETERRADOSSI wrote:

Many thanks David for making the connection betweeen the two reports.
If I understand, other french users may have faced the same problem...
However I'll wait until it's solved to have a look at "foreach".
Regards. Olivier

David M Smith a écrit :
A user in Japan reported a similar problem on the Revolutions blog ( 
http://bit.ly/FKP2I ), and my best guess is that it's an 
(unintended!) effect of using locales. The developers in New Haven 
are looking at it, and I expect they'll be able to post an update to 
CRAN soon.


# David Smith

On Mon, Jul 20, 2009 at 5:48 AM, Olivier ETERRADOSSI 
mailto:olivier.eterrado...@ema.fr>> wrote:


Hi list,
My attention was drawn to the foreach package by recent posts...I
decided to have a look...
I'm using R.2.9.1 on Windows, I have downloaded the foreach
package today (v 1.2.1), together with iterators (v. 1.0.1) and
codetools (v.0.2-2).
Full of hope I try the most simple thing of all out of the package
vignette :

 x <- foreach(i = 1:3) %do% sqrt(i)

and get :

Erreur dans sqrt(i) : indice hors limites ( i.e. "error in
sqrt(i) : index out of bounds")

but when trying :

x<-foreach(i = 1:3) %do% print(sqrt(i))

I get :

[1] 1
[1] 1.414214
[1] 1.732051
Erreur dans print(sqrt(i)) : indice hors limites

Probably I didn't drink enough coffee this morning and I'm still
asleep : it is obvious that  I miss a point... but I am unable to
see which one.
Any help appreciated ! Many thanks, and very best regards Olivier

-- Olivier ETERRADOSSI
Maître-Assistant
CMGD / Equipe "Propriétés Psycho-Sensorielles des Matériaux"
Ecole des Mines d'Alès
Hélioparc, 2 av. P. Angot, F-64053 PAU CEDEX 9
tel std: +33 (0)5.59.30.54.25
tel direct: +33 (0)5.59.30.90.35 fax: +33 (0)5.59.30.63.68
http://www.ema.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.




--
David M Smith >
Director of Community, REvolution Computing 
www.revolution-computing.com 

Tel: +1 (206) 577-4778 x3203 (San Francisco, USA)

Check out our upcoming events schedule at 
www.revolution-computing.com/events 







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

Re: [R] Constructing lists (yet, again)

[roger koenker]

Marc,

Thanks,

sapply(ls(pat = "^name"),get)

was exactly what I was after.  The default behavior for vectors of
equal length is nice too, but I was most interested in the ragged
case, which produces a list.


url:www.econ.uiuc.edu/~rogerRoger Koenker
emailrkoen...@uiuc.eduDepartment of Economics
vox: 217-333-4558University of Illinois
fax:   217-244-6678Urbana, IL 61801



On Jul 23, 2009, at 2:41 PM, Marc Schwartz wrote:


On Jul 23, 2009, at 9:19 AM, roger koenker wrote:

This is an attempt to rescue an old R-help question that apparently  
received
no response from the oblivion of collective silence, and besides  
I'm also

curious about the answer


From: Griffith Feeney (gfee...@hawaii.edu)
Date: Fri 28 Jan 2000 - 07:48:45 EST   wrote (to R-help)
Constructing lists with

list(name1=name1, name2=name2, ...)

is tedious when there are many objects and names are long. Is  
there an R
function that takes a character vector of object names as an  
argument and

returns a list with each objected tagged by its name?


The idiom

lapply(ls(pat = "^name"), function(x) eval(as.name(x)))

makes the list, but (ironically)  doesn't assign the names to the  
components.


Roger,

How about something like this:

name1 <- 1:3
name2 <- 1:5
name3 <- 1:9
name4 <- 1:7


> ls(pat = "^name")
[1] "name1" "name2" "name3" "name4"


> sapply(ls(pat = "^name"), get, simplify = FALSE)
$name1
[1] 1 2 3

$name2
[1] 1 2 3 4 5

$name3
[1] 1 2 3 4 5 6 7 8 9

$name4
[1] 1 2 3 4 5 6 7


Is that what you are after? With sapply(), you can take advantage of  
the USE.NAMES argument, which defaults to TRUE and then set simplify  
to FALSE to force the result to be a list rather than a matrix. Of  
course, in the case I have above, when there are uneven length  
elements, the result will be a list anyway.


HTH,

Marc Schwartz


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

[Gabor Grothendieck]

Best thing is to test it out on simulated data for which you
already know the answer.

library(dyn)
set.seed(123)
DF <- data.frame(x = rnorm(10), z = rnorm(10))
DF$Y <- 0
for(i in 2:10) {
  DF$Y[i] <- with(DF, 2*Y[i-1] + 3*z[i] +4* x[i] + 5*x[i-1] + rnorm(1))
}
DF.zoo <- do.call(merge, lapply(DF, zoo))
L <- function(x, k = 1) lag(x, -k)
fit <- dyn$lm(Y ~ L(Y) + z + L(x, 0:1), DF.zoo[-10, ])
fit
pred <- predict(fit, DF.zoo)
cbind(pred, DF.zoo$Y)
plot(cbind(pred, act = DF.zoo$Y), type = c("l", "p"), screen = 1)




On Thu, Jul 23, 2009 at 2:48 PM, Hongwei Dong wrote:
> Hi, Gabor, got it. Thanks a lot.
> I have one more question about how the "predict" function works here,
> especially for the lag(Y,-1) part.
> In my model, I assume I know predictors x and z in the next two years, and
> use them to predict Y. For each forecast step at time t, the lag(Y,-1) in
> the model should be the estimated result from the last forecast step. For
> example, the prediction of the Y in 2007 should be based on the estimated Y
> for 2006, not the actual number already in my dataset. Before I report my
> results, I want to make sure the "predict" function works in this way.
> Thanks.
> Harry
>
>
> On Thu, Jul 23, 2009 at 10:28 AM, Gabor Grothendieck
>  wrote:
>>
>> Please provide your code in a reproducible form (as requested previously).
>>
>> Here we fit with first 9 points and then add a point for prediction.  (Of
>> course your model can only predict the current value of Y so you may
>> have to rethink your model even aside from the implementation if you
>> really want to predict future values of Y.)
>>
>> > library(dyn)
>> > set.seed(123)
>> > tt <- ts(cbind(Y = 1:10, x = rnorm(10), z = rnorm(10)))
>> > L <- function(x, k = 1) lag(x, -k)
>> > tt.zoo <- as.zoo(tt)
>> > fit <- dyn$lm(Y ~ L(Y) + L(x, 0:2) + z, tt.zoo[-10, ])
>> > predict(fit, tt.zoo)
>>  1  2  3  4  5  6  7  8  9 10 11 12
>> NA NA  3  4  5  6  7  8  9 10 NA NA
>>
>>
>>
>>
>> On Thu, Jul 23, 2009 at 1:01 PM, Hongwei Dong wrote:
>> > Hi, Gabor and Other R users,
>> >   I'm re-posting my script and the results I got.
>> >
>> >   here is the dynamic model I used to estimate in-sample model
>> > (1996-2006)
>> > and it works:
>> >
>> >   fit<-dyn$lm(Y~lag(Y,-1)+z+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4))
>> >   Then I used this model to do out sample forecast with the following
>> > scripts, which do not work:
>> >
>> >  z<-ts(Z[41:52],start=2006,frequency=4)
>> >  x<-ts(X[41:52],start=2006,frequency=4)
>> > newdata<-data.frame(cbind(z,x))
>> > newdata<-ts(newdata,start=2006,frequency=4)
>> > pred<-predict(fit,newdata)
>> > Here is the results I got from R:
>> >          Qtr1     Qtr2     Qtr3     Qtr4
>> > 2006       NA       NA       NA       NA
>> > 2007 3083.362       NA       NA       NA
>> > 2008       NA       NA       NA       NA
>> > 2009       NA       NA       NA       NA
>> >
>> > I got only one prediction for the first quarter in 2007. Intuitively,
>> > there
>> > might be two problems: the definition of the newdata and how to define Y
>> > in
>> > newdata. But I just can't figure this out. It will greatly appreciated
>> > if
>> > someone can give me some help. Thanks.
>> > Harry
>> >
>> >
>> > On Thu, Jul 23, 2009 at 5:15 AM, Gabor Grothendieck
>> >  wrote:
>> >>
>> >> You have to use consistent classes.  You can't start out using
>> >> ts class and then suddenly switch to zoo class as if you had been
>> >> using zoo all along.   Either use ts everywhere or zoo everywhere.
>> >>
>> >> Also in the future please post reproducible examples as requested
>> >> at the bottom of every message to r-help.  That means include
>> >> a minimal amount of data so we can get exactly what you
>> >> are getting.
>> >>
>> >> On Thu, Jul 23, 2009 at 4:48 AM, Hongwei Dong wrote:
>> >> > Thanks, Gabor. This is really helpful.
>> >> > When the regressive part, lag(Y,-1), is not included, my sytax works
>> >> > well.
>> >> > However, when I include the lag(Y) in the model, the prediction part
>> >> > does
>> >> > not work. Here is my sytax for in-sample estimation and it works
>> >> > well:
>> >> > fit<-dyn$lm(Y~lag(Y,-1)+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4)+z).
>> >> > Then I use this model to do out of sample prediction:
>> >> > x<-ts(X[41:52],start=2006,frequency=4)
>> >> > z<-ts(Z[41:52],start=2006,frequency=4)
>> >> > newdata<-data.frame(cbind(x,z))
>> >> > newdata<-zooreg(newdata)
>> >> > pred<-predict(fit,newdata)
>> >> > With these, I got weird result, a prediction for each year from year
>> >> > 1
>> >> > to
>> >> > the first quarter of year 2007 (all "NA"). What I expect is a
>> >> > prediction
>> >> > for
>> >> > the 8 quarters from 2007 to 2008. Intuitively, I know there must be
>> >> > something wrong with my newdata definition. But I can't figure it
>> >> > out.
>> >> > I'll
>> >> > appreciate it very much if you can give some suggestions to modify my
>> >> > syntax. Thanks.
>> >> >
>> >> > Harry
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >

Re: [R] Constructing lists (yet, again)

[Greg Snow]

There are a couple of options:

The help page for lapply also includes the help for sapply and sapply has a 
USE.NAMES argument that may do what you want (specify simplify=FALSE to force 
the same behavior as lapply).

You can post specify the names like:

> names(mylist) <- vector.of.names

Do either of those do what you want?

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of roger koenker
> Sent: Thursday, July 23, 2009 8:20 AM
> To: r-help@r-project.org
> Subject: [R] Constructing lists (yet, again)
> 
> This is an attempt to rescue an old R-help question that apparently
> received
> no response from the oblivion of collective silence, and besides I'm
> also
> curious about the answer
> 
> > From: Griffith Feeney (gfee...@hawaii.edu)
> > Date: Fri 28 Jan 2000 - 07:48:45 EST   wrote (to R-help)
> > Constructing lists with
> >
> > list(name1=name1, name2=name2, ...)
> >
> > is tedious when there are many objects and names are long. Is there
> > an R
> > function that takes a character vector of object names as an
> > argument and
> > returns a list with each objected tagged by its name?
> >
> The idiom
> 
>   lapply(ls(pat = "^name"), function(x) eval(as.name(x)))
> 
> makes the list, but (ironically)  doesn't assign the names to the
> components.
> 
> 
> 
> url:www.econ.uiuc.edu/~rogerRoger Koenker
> emailrkoen...@uiuc.eduDepartment of Economics
> vox: 217-333-4558University of Illinois
> fax:   217-244-6678Urbana, IL 61801
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Constructing lists (yet, again)

[Marc Schwartz]

On Jul 23, 2009, at 9:19 AM, roger koenker wrote:

This is an attempt to rescue an old R-help question that apparently  
received
no response from the oblivion of collective silence, and besides I'm  
also

curious about the answer


From: Griffith Feeney (gfee...@hawaii.edu)
Date: Fri 28 Jan 2000 - 07:48:45 EST   wrote (to R-help)
Constructing lists with

list(name1=name1, name2=name2, ...)

is tedious when there are many objects and names are long. Is there  
an R
function that takes a character vector of object names as an  
argument and

returns a list with each objected tagged by its name?


The idiom

lapply(ls(pat = "^name"), function(x) eval(as.name(x)))

makes the list, but (ironically)  doesn't assign the names to the  
components.


Roger,

How about something like this:

name1 <- 1:3
name2 <- 1:5
name3 <- 1:9
name4 <- 1:7


> ls(pat = "^name")
[1] "name1" "name2" "name3" "name4"


> sapply(ls(pat = "^name"), get, simplify = FALSE)
$name1
[1] 1 2 3

$name2
[1] 1 2 3 4 5

$name3
[1] 1 2 3 4 5 6 7 8 9

$name4
[1] 1 2 3 4 5 6 7


Is that what you are after? With sapply(), you can take advantage of  
the USE.NAMES argument, which defaults to TRUE and then set simplify  
to FALSE to force the result to be a list rather than a matrix. Of  
course, in the case I have above, when there are uneven length  
elements, the result will be a list anyway.


HTH,

Marc Schwartz

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

Re: [R] Random # generator accuracy

[Greg Snow]

Well one quick way (for non-generics) is the 'args' function:

> args(sample)
function (x, size, replace = FALSE, prob = NULL) 
NULL

A similar line appears near the top of the help page when you do '?sample'.  
The "replace = FALSE" in the line above means that false is the default (with 
the assumption that FALSE => NO Replacement).  The context of the first example 
and part of the details section supports that 'without replacement' is the 
default, but you are correct that it could be clearer (some help pages list the 
defaults for each of the arguments).  Some functions don't show the default, or 
use NULL as the default so that you need to read closer in the argument list or 
details about what will happen when that argument is left blank.  

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: Jim Bouldin [mailto:jrboul...@ucdavis.edu]
> Sent: Thursday, July 23, 2009 12:49 PM
> To: Greg Snow; r-help@r-project.org
> Subject: RE: [R] Random # generator accuracy
> 
> 
> Thanks Greg, that most definitely was it.  So apparently the default is
> sampling without replacement.  Fine, but this brings up a question I've
> had
> for a bit now, which is, how do you know what the default settings are
> for
> the arguments of any given function?  The HTML help files don't seem to
> indicate in many (most) cases.  Thanks.
> 
> > Try adding replace=TRUE to your call to sample, then you will get
> numbers
> > closer to what you are expecting.
> >
> > --
> > Gregory (Greg) L. Snow Ph.D.
> > Statistical Data Center
> > Intermountain Healthcare
> > greg.s...@imail.org
> > 801.408.8111
> >
> >
> > > -Original Message-
> > > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> > > project.org] On Behalf Of Jim Bouldin
> > > Sent: Thursday, July 23, 2009 12:00 PM
> > > To: r-help@r-project.org
> > > Subject: [R] Random # generator accuracy
> > >
> > >
> > > Dan Nordlund wrote:
> > >
> > > "It would be necessary to see the code for your 'brief test' before
> > > anyone
> > > could meaningfully comment on your results.  But your results for a
> > > single
> > > test could have been a valid "random" result."
> > >
> > > I've re-created what I did below.  The problem appears to be with
> the
> > > weighting process: the unweighted sample came out much closer to
> the
> > > actual
> > > than the weighted sample (>1% error) did.  Comments?
> > > Jim
> > >
> > > > x
> > >  [1]  1  2  3  4  5  6  7  8  9 10 11 12
> > > > weights
> > >  [1] 1 1 1 1 1 1 2 2 2 2 2 2
> > >
> > > > a = mean(replicate(100,(sample(x, 3, prob = weights;a  #
> (1
> > > million samples from x, of size 3, weighted by "weights"; the mean
> > > should
> > > be 7.50)
> > > [1] 7.406977
> > > > 7.406977/7.5
> > > [1] 0.987597
> > >
> > > > b = mean(replicate(100,(sample(x, 3;b  # (1 million
> samples
> > > from
> > > x, of size 3, not weighted this time; the mean should be 6.50)
> > > [1] 6.501477
> > > > 6.501477/6.5
> > > [1] 1.000227
> > >
> > >
> > > Jim Bouldin, PhD
> > > Research Ecologist
> > > Department of Plant Sciences, UC Davis
> > > Davis CA, 95616
> > > 530-554-1740
> > >
> > > __
> > > R-help@r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/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 Bouldin, PhD
> Research Ecologist
> Department of Plant Sciences, UC Davis
> Davis CA, 95616
> 530-554-1740

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

Re: [R] howto create a list row-by-row as input to function call?

[Mark Knecht]

On Thu, Jul 23, 2009 at 11:04 AM, Mark Knecht wrote:
> Hi,
>   I'm having trouble within my function CalcPos to get it to call
> CalcHorz with values from each row. I *think* it's calling CalcHorz
> with the final values of the inputs and not the values from each row.
> How can I do this properly in R?
>
>   The values aa,bb,cc,dd are inputs. CalcPos first calculates V1 and
> V2 vertically, and then I attempt to call CalcHorz to handle H1, H2 &
> H3 using a list L1 comprised of df[aa:V2]. I wanted that list to be
> the values on that row and disposed of each pass, but I think all
> calls are probably using the final values. Finally I use the returned
> values H1:H3 for more vertical calculations in V3.
>
>   How can I get list L1 to be the unique row values for each call to
> CalcHorz? L1 would just be tossed when the call to CalcPos is
> complete.
>
>   PLEASE - do not combine any of the inputs or return values. I need
> all of them aa through V3 and want to keep them separate. Thanks!
>
> Cheers,
> Mark
>

Doing the CalcHorz as a data.frame solved the problem for me.

Sorry for the noise.

- Mark

CalcHorz = function (L1) {
# Intended to calculate 3 values for each row
H1 <- with(L1, V1*cc)
H2 <- with(L1, V2*dd)
H3 <- with(L1, H1 - H2)
return(data.frame(H1, H2, H3))
}

CalcPos = function (df) {

#calculate vertical initially
df$V1 <- with(df, aa*bb)
df$V2 <- with(df, aa*cc-dd)

#At this point do horizontal calculations
#Goal is this is done for each row using
#the values aa,bb,cc,dd,V1,V2 for that row

L1 = df[1:dim(df)[2]]

R1 = CalcHorz(L1)
df$H1 <- R1[1]
df$H2 <- R1[2]
df$H3 <- R1[3]

#Switch back to vertical

df$V3 = with (df, cumsum(H3))

return(df)
}


DF <- data.frame(aa=1:12, bb=12:15, cc=1:3, dd=1:6)
DF

DF <- CalcPos(DF)
DF

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

[roger koenker]

This is an attempt to rescue an old R-help question that apparently  
received
no response from the oblivion of collective silence, and besides I'm  
also

curious about the answer


From: Griffith Feeney (gfee...@hawaii.edu)
Date: Fri 28 Jan 2000 - 07:48:45 EST   wrote (to R-help)
Constructing lists with

list(name1=name1, name2=name2, ...)

is tedious when there are many objects and names are long. Is there  
an R
function that takes a character vector of object names as an  
argument and

returns a list with each objected tagged by its name?


The idiom

lapply(ls(pat = "^name"), function(x) eval(as.name(x)))

makes the list, but (ironically)  doesn't assign the names to the  
components.




url:www.econ.uiuc.edu/~rogerRoger Koenker
emailrkoen...@uiuc.eduDepartment of Economics
vox: 217-333-4558University of Illinois
fax:   217-244-6678Urbana, IL 61801

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

Re: [R] Split plot analysis problems

[Mark Difford]

Hi Jean-Paul,

>> However, I've tried both solutions on my model, and I got different
>> residuals :...
>> What could be the difference between the two?

There is no difference. You have made a mistake.

##
tt <- data.frame(read.csv(file="tt.csv", sep=""))  ## imports your data set
T.aov <- aov(PH~Community*Mowing*Water + Error(Block), data=tt)

summary(T.aov)  ## This matches your output

   Df Sum Sq Mean Sq F valuePr(>F)
Mowing  1 40.007  40.007 21.1747 0.0002215 ***
Water   1 23.120  23.120 12.2370 0.0025673 ** 
Community:Mowing1 21.060  21.060 11.1467 0.0036554 ** 
Community:Water 1  6.901   6.901  3.6524 0.0720478 .  
Mowing:Water1  1.611   1.611  0.8527 0.3680090
Community:Mowing:Water  1  0.858   0.858  0.4542 0.5089331
Residuals  18 34.008   1.889  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

##
length(residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within))
[1] 24
length(proj(T.aov)[,"Residuals"])
[1] 24

## Details
residuals(T.aov <- aov(PH~Community*Mowing*Water + Error(Block),
data=tt)$Within)

   9   10   11   12   13  
14   15   16 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118 
  17   18   19   20   21  
22   23   24 
 1.241311954  1.498811954 -0.616188046  0.483811954 -0.477827381
-1.660327381  2.684672619  1.924672619 
  25   26   27   28   29  
30   31   32 
-0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358
-0.751015358  0.836484642  1.181484642 

proj(T.aov)[,"Residuals"]
   9   10   11   12   13  
14   15   16   17   18   19  
20   21 
-1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882 
0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954 -0.616188046 
0.483811954 -0.477827381 
  22   23   24   25   26  
27   28   29   30   31   32 
-1.660327381  2.684672619  1.924672619 -0.465198010 -1.735198010 
1.502301990  0.839801990 -0.428515358 -0.751015358  0.836484642  1.181484642 

Regards, Mark.


Jean-Paul Maalouf wrote:
> 
> Thanks Mark and Richard for your propositions on how to extract residuals.
> 
> However, I've tried both solutions on my model, and I got different  
> residuals :
> If we consider the within residuals :
> 
> Mark's solution gave me a vector of 24 residuals :
> 
>> as.vector(residuals(aov(PH~Community*Mowing*Water +
>> Error(Block))$Within))
>   [1] -1.107277057 -0.977277057 -0.007277057 -0.719777057 -1.021627882  
>   0.215872118 -1.621627882  0.508372118  1.241311954  1.498811954  
> -0.616188046
> [12]  0.483811954 -0.477827381 -1.660327381  2.684672619  1.924672619  
> -0.465198010 -1.735198010  1.502301990  0.839801990 -0.428515358  
> -0.751015358
> [23]  0.836484642  1.181484642
> 
> and Richard's solution gave me a vector 32 residuals :
> 
> do.call(data.frame,proj(aov(PH~Community*Water*Mowing +
> Error(Block->proj
> proj$Within.Residuals
>   [1] -0.216875  1.745625 -1.174375 -0.354375  0.800625  0.460625  
> -0.799375 -0.461875 -0.404375 -0.274375  0.695625 -0.016875 -0.541875   
> 0.695625 -1.141875
> [16]  0.988125  0.589375  0.846875 -1.268125 -0.168125 -1.095625  
> -2.278125  2.066875  1.306875 -0.500625 -1.770625  1.466875  0.804375  
> -0.638125 -0.960625
> [31]  0.626875  0.971875
> 
> What could be the difference between the two?
> 
> Regards
> 
> JPM
> 
> 
> Quoting "Richard M. Heiberger" :
> 
>> Jean-Paul Maalouf wrote:
>>> Do you have any idea on how can I verify preliminary assumptions in  
>>>  this model (normality of the residuals and variance homogeneity),   
>>> since R is not able to extract residuals?
>>
>> Of course, R extracts residuals.  Use the proj() function.  See ?proj
>> for the example
>> to get the projection of an aovlist object onto the terms of a linear
>> model
>>
>> proj(npk.aovE)
>>
>>
>> To get the projections into a simple data.frame, use
>> tmpE <- proj(npk.aovE)
>> tmpE.df <- do.call("data.frame", tmpE)
>> tmpE.df
>>
>> Mark Difford's solution effectively retrieved columns 3 and 10 from
>> tmpE.df
>>
>> Rich
> 
> 
> 
> -- 
> Jean-Paul Maalouf
> UMR 1202 BIOGECO
> Inra - Université Bordeaux 1
> Bâtiment B8, Avenue des Facultés
> 33405 Talence, France
> Tel : 05 40008772
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View th

Re: [R] setting up LMER for repeated measures and how do I get a p value for my fixed effect, group?

[Kingsford Jones]

On Thu, Jul 23, 2009 at 12:36 PM, Kingsford
Jones wrote:

>
> As you might have suspected from the lmer t-value close to 0, the
> associated p-value is about .5.

you can ignore the sentence above -- it's a two sided test and the
t-value is not that close to 0...

>
> hth,
>
> Kingsford Jones
>
>> My subjects are divided into two groups (variable GROUP), individual 
>> subjects are indicated by the variable SS, Value is the outcome measure, 
>> each subject has Value measured three times.
>>
>> I have used the following code:
>>
>> fit1<-lmer(Value~Group+(1|SS),data=smallDS)
>> print(fit1)
>>
>> Linear mixed model fit by REML
>> Formula: Value ~ Group + (1 | SS)
>>   Data: Dataset
>>   AIC   BIC logLik deviance REMLdev
>>  284.8 290.4 -138.4    291.4   276.8
>> Random effects:
>>  Groups   Name        Variance Std.Dev.
>>  SS       (Intercept) 1038.0   32.218
>>  Residual              552.7   23.510
>> Number of obs: 30, groups: SS, 10
>>
>> Fixed effects:
>>                        Estimate Std. Error t value
>> (Intercept)               152.87      15.63   9.778
>> Group[T.ContSIRNAnoEGF]   -15.87      22.11  -0.718
>>
>> Correlation of Fixed Effects:
>>            (Intr)
>> G[T.CSIRNAE -0.707
>>
>>
>> Questions:
>> (1) Have I set up lmer properly to analyze the repeated measures data?
>> (2) How can I get a p value for Group, my fixed effect?
>>
>>
>>
>>      Group SS Value
>>  [1,]     2  1   127
>>  [2,]     2  1   179
>>  [3,]     2  1   159
>>  [4,]     2  2   186
>>  [5,]     2  2   173
>>  [6,]     2  2   178
>>  [7,]     2  3   117
>>  [8,]     2  3   116
>>  [9,]     2  3    70
>> [10,]     2  4   176
>> [11,]     2  4   149
>> [12,]     2  4   138
>> [13,]     2  5   100
>> [14,]     2  5   105
>> [15,]     2  5    82
>> [16,]     1  6   142
>> [17,]     1  6   195
>> [18,]     1  6   222
>> [19,]     1  7   218
>> [20,]     1  7   178
>> [21,]     1  7   147
>> [22,]     1  8   135
>> [23,]     1  8   162
>> [24,]     1  8   177
>> [25,]     1  9   104
>> [26,]     1  9   102
>> [27,]     1  9   121
>> [28,]     1 10   135
>> [29,]     1 10   121
>> [30,]     1 10   134
>>
>>
>> John David Sorkin M.D., Ph.D.
>> Chief, Biostatistics and Informatics
>> University of Maryland School of Medicine Division of Gerontology
>> Baltimore VA Medical Center
>> 10 North Greene Street
>> GRECC (BT/18/GR)
>> Baltimore, MD 21201-1524
>> (Phone) 410-605-7119
>> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>>
>> Confidentiality Statement:
>> This email message, including any attachments, is for th...{{dropped:6}}
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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] Random # generator accuracy

[Jim Bouldin]

Thanks Greg, that most definitely was it.  So apparently the default is
sampling without replacement.  Fine, but this brings up a question I've had
for a bit now, which is, how do you know what the default settings are for
the arguments of any given function?  The HTML help files don't seem to
indicate in many (most) cases.  Thanks. 

> Try adding replace=TRUE to your call to sample, then you will get numbers
> closer to what you are expecting.
> 
> -- 
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> greg.s...@imail.org
> 801.408.8111
> 
> 
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> > project.org] On Behalf Of Jim Bouldin
> > Sent: Thursday, July 23, 2009 12:00 PM
> > To: r-help@r-project.org
> > Subject: [R] Random # generator accuracy
> > 
> > 
> > Dan Nordlund wrote:
> > 
> > "It would be necessary to see the code for your 'brief test' before
> > anyone
> > could meaningfully comment on your results.  But your results for a
> > single
> > test could have been a valid "random" result."
> > 
> > I've re-created what I did below.  The problem appears to be with the
> > weighting process: the unweighted sample came out much closer to the
> > actual
> > than the weighted sample (>1% error) did.  Comments?
> > Jim
> > 
> > > x
> >  [1]  1  2  3  4  5  6  7  8  9 10 11 12
> > > weights
> >  [1] 1 1 1 1 1 1 2 2 2 2 2 2
> > 
> > > a = mean(replicate(100,(sample(x, 3, prob = weights;a  # (1
> > million samples from x, of size 3, weighted by "weights"; the mean
> > should
> > be 7.50)
> > [1] 7.406977
> > > 7.406977/7.5
> > [1] 0.987597
> > 
> > > b = mean(replicate(100,(sample(x, 3;b  # (1 million samples
> > from
> > x, of size 3, not weighted this time; the mean should be 6.50)
> > [1] 6.501477
> > > 6.501477/6.5
> > [1] 1.000227
> > 
> > 
> > Jim Bouldin, PhD
> > Research Ecologist
> > Department of Plant Sciences, UC Davis
> > Davis CA, 95616
> > 530-554-1740
> > 
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/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 Bouldin, PhD
Research Ecologist
Department of Plant Sciences, UC Davis
Davis CA, 95616
530-554-1740

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

[Hongwei Dong]

Hi, Gabor, got it. Thanks a lot.
I have one more question about how the "predict" function works here,
especially for the lag(Y,-1) part.

In my model, I assume I know predictors x and z in the next two years, and
use them to predict Y. For each forecast step at time t, the lag(Y,-1) in
the model should be the estimated result from the last forecast step. For
example, the prediction of the Y in 2007 should be based on the estimated Y
for 2006, not the actual number already in my dataset. Before I report my
results, I want to make sure the "predict" function works in this way.
Thanks.

Harry



On Thu, Jul 23, 2009 at 10:28 AM, Gabor Grothendieck <
ggrothendi...@gmail.com> wrote:

> Please provide your code in a reproducible form (as requested previously).
>
> Here we fit with first 9 points and then add a point for prediction.  (Of
> course your model can only predict the current value of Y so you may
> have to rethink your model even aside from the implementation if you
> really want to predict future values of Y.)
>
> > library(dyn)
> > set.seed(123)
> > tt <- ts(cbind(Y = 1:10, x = rnorm(10), z = rnorm(10)))
> > L <- function(x, k = 1) lag(x, -k)
> > tt.zoo <- as.zoo(tt)
> > fit <- dyn$lm(Y ~ L(Y) + L(x, 0:2) + z, tt.zoo[-10, ])
> > predict(fit, tt.zoo)
>  1  2  3  4  5  6  7  8  9 10 11 12
> NA NA  3  4  5  6  7  8  9 10 NA NA
>
>
>
>
> On Thu, Jul 23, 2009 at 1:01 PM, Hongwei Dong wrote:
> > Hi, Gabor and Other R users,
> >   I'm re-posting my script and the results I got.
> >
> >   here is the dynamic model I used to estimate in-sample model
> (1996-2006)
> > and it works:
> >
> >   fit<-dyn$lm(Y~lag(Y,-1)+z+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4))
> >   Then I used this model to do out sample forecast with the following
> > scripts, which do not work:
> >
> >  z<-ts(Z[41:52],start=2006,frequency=4)
> >  x<-ts(X[41:52],start=2006,frequency=4)
> > newdata<-data.frame(cbind(z,x))
> > newdata<-ts(newdata,start=2006,frequency=4)
> > pred<-predict(fit,newdata)
> > Here is the results I got from R:
> >  Qtr1 Qtr2 Qtr3 Qtr4
> > 2006   NA   NA   NA   NA
> > 2007 3083.362   NA   NA   NA
> > 2008   NA   NA   NA   NA
> > 2009   NA   NA   NA   NA
> >
> > I got only one prediction for the first quarter in 2007. Intuitively,
> there
> > might be two problems: the definition of the newdata and how to define Y
> in
> > newdata. But I just can't figure this out. It will greatly appreciated if
> > someone can give me some help. Thanks.
> > Harry
> >
> >
> > On Thu, Jul 23, 2009 at 5:15 AM, Gabor Grothendieck
> >  wrote:
> >>
> >> You have to use consistent classes.  You can't start out using
> >> ts class and then suddenly switch to zoo class as if you had been
> >> using zoo all along.   Either use ts everywhere or zoo everywhere.
> >>
> >> Also in the future please post reproducible examples as requested
> >> at the bottom of every message to r-help.  That means include
> >> a minimal amount of data so we can get exactly what you
> >> are getting.
> >>
> >> On Thu, Jul 23, 2009 at 4:48 AM, Hongwei Dong wrote:
> >> > Thanks, Gabor. This is really helpful.
> >> > When the regressive part, lag(Y,-1), is not included, my sytax works
> >> > well.
> >> > However, when I include the lag(Y) in the model, the prediction part
> >> > does
> >> > not work. Here is my sytax for in-sample estimation and it works well:
> >> > fit<-dyn$lm(Y~lag(Y,-1)+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4)+z).
> >> > Then I use this model to do out of sample prediction:
> >> > x<-ts(X[41:52],start=2006,frequency=4)
> >> > z<-ts(Z[41:52],start=2006,frequency=4)
> >> > newdata<-data.frame(cbind(x,z))
> >> > newdata<-zooreg(newdata)
> >> > pred<-predict(fit,newdata)
> >> > With these, I got weird result, a prediction for each year from year 1
> >> > to
> >> > the first quarter of year 2007 (all "NA"). What I expect is a
> prediction
> >> > for
> >> > the 8 quarters from 2007 to 2008. Intuitively, I know there must be
> >> > something wrong with my newdata definition. But I can't figure it out.
> >> > I'll
> >> > appreciate it very much if you can give some suggestions to modify my
> >> > syntax. Thanks.
> >> >
> >> > Harry
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> >
> >> > On Wed, Jul 22, 2009 at 10:53 PM, Gabor Grothendieck
> >> >  wrote:
> >> >>
> >> >> Here is an example closer to yours.
> >> >>
> >> >> > library(dyn)
> >> >> > set.seed(123)
> >> >> > x <- zooreg(rnorm(10))
> >> >> > y <- zooreg(rnorm(10))
> >> >> > L <- function(x, k = 1) lag(x, k = -k)
> >> >> > mod <- dyn$lm(y ~ L(y) + L(x, 0:2))
> >> >> > mod
> >> >>
> >> >> Call:
> >> >> lm(formula = dyn(y ~ L(y) + L(x, 0:2)))
> >> >>
> >> >> Coefficients:
> >> >> (Intercept) L(y)   L(x, 0:2)1   L(x, 0:2)2   L(x, 0:2)3
> >> >>0.06355 -0.74540  0.63649  0.44957 -0.41360
> >> >>
> >> >> > newdata <- cbind(x = c(coredata(x), rnorm(1)), y = c(coredata

Re: [R] Elementary Symmetric Polynomials

[davidr]

How's this?

f2 <- function(v) {
  n <- length(v)
  s <- matrix(0,n,n)
  for(i in 1:n)
for(j in 1:i)
  s[j,i] <- ifelse(i>1, s[j,i-1]+v[i]*ifelse(j>1, s[j-1,i-1], 1),
v[i])
  c(1,s[,n])
}

> system.time(f2.result <- f2(1:20))
   user  system elapsed 
  0   0   0 
> system.time(f.result <- f(1:20))
   user  system elapsed 
  29.500.04   29.58 
> all.equal(f.result, f2.result)
[1] TRUE


-- David Reiner


-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Austin H. Jones
Sent: Tuesday, July 21, 2009 5:11 PM
To: r-help@r-project.org
Subject: [R] Elementary Symmetric Polynomials

We are interested in obtaining an efficient function that for a given 
vector of length t will output a vector of length t+1 that contains the 
associated values of the elementary symmetric polynomials in t 
variables. Below is what we have at the moment, but it is a little slow 
for our needs. Any suggestions?

Thanks ahead of time for any help you can offer,

Austin H. Jones
Department of Mathematics
Wake Forest University



f<-function(v)
{
prodsub<-function(v,q){prod(v[q])}
t<-length(v)
C<-vector("list",t)
for (i in 1:t)
{C[[i]]<-combn(1:t,i)}
e<-rep(0,t)
for (i in 1:(t))
{
e[i]<-sum(apply(C[[i]],2,prodsub,v=v))
}
e<-c(1,e)
e
}

Examples:


 > f(c(1,2,3))
[1]  1  6 11  6


 > f(c(4,6,3,1,9))
[1]1   23  193  729 1206  648

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


This e-mail and any materials attached hereto, including, without limitation, 
all content hereof and thereof (collectively, "Rho Content") are confidential 
and proprietary to Rho Trading Securities, LLC ("Rho") and/or its affiliates, 
and are protected by intellectual property laws.  Without the prior written 
consent of Rho, the Rho Content may not (i) be disclosed to any third party or 
(ii) be reproduced or otherwise used by anyone other than current employees of 
Rho or its affiliates, on behalf of Rho or its affiliates.

THE RHO CONTENT IS PROVIDED AS IS, WITHOUT REPRESENTATIONS OR WARRANTIES OF ANY 
KIND.  TO THE MAXIMUM EXTENT PERMISSIBLE UNDER APPLICABLE LAW, RHO HEREBY 
DISCLAIMS ANY AND ALL WARRANTIES, EXPRESS AND IMPLIED, RELATING TO THE RHO 
CONTENT, AND NEITHER RHO NOR ANY OF ITS AFFILIATES SHALL IN ANY EVENT BE LIABLE 
FOR ANY DAMAGES OF ANY NATURE WHATSOEVER, INCLUDING, BUT NOT LIMITED TO, 
DIRECT, INDIRECT, CONSEQUENTIAL, SPECIAL AND PUNITIVE DAMAGES, LOSS OF PROFITS 
AND TRADING LOSSES, RESULTING FROM ANY PERSON'S USE OR RELIANCE UPON, OR 
INABILITY TO USE, ANY RHO CONTENT, EVEN IF RHO IS ADVISED OF THE POSSIBILITY OF 
SUCH DAMAGES OR IF SUCH DAMAGES WERE FORESEEABLE.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 perform a calculation in each element of my list?

[Erik Iverson]

Mark, 

My example is essentially identical to Jorge's.  This is a good opportunity to 
compare two solutions to a problem, one using "for" loops, and one using the 
apply family of functions.  Compare this with Daniel's solution. 

## BEGIN EXAMPLE
## sample list of data.frames, different number of columns in each one
lst <- sapply(c(20, 25), function(x) as.data.frame(matrix(1:100, nrow = x)))

## define function to do what you want for one data.frame, needs at least 2 ## 
columns, no checks for that or for them being numeric...
mult2cols <- function(x) {
  nc <- ncol(x)
  x$product <- x[,nc] * x[,nc-1]
  x
}

## apply the function to your list of data.frames
lst <- lapply(lst, tmp)

## END EXAMPLE

--erik 

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Mark Na
Sent: Thursday, July 23, 2009 1:06 PM
To: r-help@r-project.org
Subject: [R] How to perform a calculation in each element of my list?

Hi R-helpers,

I have a list containing 10 elements, each of which is a dataframe. I wish
to add a new column to each list element (dataframe) containing the product
of the last two columns of each dataframe.

I'd appreciate any pointers, thanks!

Mark Na

[[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] setting up LMER for repeated measures and how do I get a p value for my fixed effect, group?

[Kingsford Jones]

On Thu, Jul 23, 2009 at 9:02 AM, John Sorkin wrote:
> R 2.8.1
> Windows XP
>
> I am trying to analyze repeated measures data (the data are listed at the end 
> of this Email message) and I need help to make sure that I have properly 
> specified my model, and would like to know why lmer does not return a p value 
> for Group, my fixed effect.
>

The model looks appropriate to me.  If you had more repeated measures
within subjects you might want to model the temporal correlation by
structuring the within-subject error covariance; the correlation
argument to nlme::lme accepts a variety of corClasses for this purpose
(e.g. AR1).

FAQ 7.35 provides a link to one explanation for why lmer omits p-values:

https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html

However in your balanced, normal case you can get (in this case
identical) p-values using eith aov or lme:

fit2 <- aov(Value ~ Group + Error(SS), data = smallDS)
#or
library(nlme)
fit3 <- lme(Value ~ Group, data = smallDS, random = ~1|SS)

As you might have suspected from the lmer t-value close to 0, the
associated p-value is about .5.

hth,

Kingsford Jones

> My subjects are divided into two groups (variable GROUP), individual subjects 
> are indicated by the variable SS, Value is the outcome measure, each subject 
> has Value measured three times.
>
> I have used the following code:
>
> fit1<-lmer(Value~Group+(1|SS),data=smallDS)
> print(fit1)
>
> Linear mixed model fit by REML
> Formula: Value ~ Group + (1 | SS)
>   Data: Dataset
>   AIC   BIC logLik deviance REMLdev
>  284.8 290.4 -138.4    291.4   276.8
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  SS       (Intercept) 1038.0   32.218
>  Residual              552.7   23.510
> Number of obs: 30, groups: SS, 10
>
> Fixed effects:
>                        Estimate Std. Error t value
> (Intercept)               152.87      15.63   9.778
> Group[T.ContSIRNAnoEGF]   -15.87      22.11  -0.718
>
> Correlation of Fixed Effects:
>            (Intr)
> G[T.CSIRNAE -0.707
>
>
> Questions:
> (1) Have I set up lmer properly to analyze the repeated measures data?
> (2) How can I get a p value for Group, my fixed effect?
>
>
>
>      Group SS Value
>  [1,]     2  1   127
>  [2,]     2  1   179
>  [3,]     2  1   159
>  [4,]     2  2   186
>  [5,]     2  2   173
>  [6,]     2  2   178
>  [7,]     2  3   117
>  [8,]     2  3   116
>  [9,]     2  3    70
> [10,]     2  4   176
> [11,]     2  4   149
> [12,]     2  4   138
> [13,]     2  5   100
> [14,]     2  5   105
> [15,]     2  5    82
> [16,]     1  6   142
> [17,]     1  6   195
> [18,]     1  6   222
> [19,]     1  7   218
> [20,]     1  7   178
> [21,]     1  7   147
> [22,]     1  8   135
> [23,]     1  8   162
> [24,]     1  8   177
> [25,]     1  9   104
> [26,]     1  9   102
> [27,]     1  9   121
> [28,]     1 10   135
> [29,]     1 10   121
> [30,]     1 10   134
>
>
> John David Sorkin M.D., Ph.D.
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>
> Confidentiality Statement:
> This email message, including any attachments, is for th...{{dropped:6}}
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Random # generator accuracy

[Greg Snow]

Try adding replace=TRUE to your call to sample, then you will get numbers 
closer to what you are expecting.

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Jim Bouldin
> Sent: Thursday, July 23, 2009 12:00 PM
> To: r-help@r-project.org
> Subject: [R] Random # generator accuracy
> 
> 
> Dan Nordlund wrote:
> 
> "It would be necessary to see the code for your 'brief test' before
> anyone
> could meaningfully comment on your results.  But your results for a
> single
> test could have been a valid "random" result."
> 
> I've re-created what I did below.  The problem appears to be with the
> weighting process: the unweighted sample came out much closer to the
> actual
> than the weighted sample (>1% error) did.  Comments?
> Jim
> 
> > x
>  [1]  1  2  3  4  5  6  7  8  9 10 11 12
> > weights
>  [1] 1 1 1 1 1 1 2 2 2 2 2 2
> 
> > a = mean(replicate(100,(sample(x, 3, prob = weights;a  # (1
> million samples from x, of size 3, weighted by "weights"; the mean
> should
> be 7.50)
> [1] 7.406977
> > 7.406977/7.5
> [1] 0.987597
> 
> > b = mean(replicate(100,(sample(x, 3;b  # (1 million samples
> from
> x, of size 3, not weighted this time; the mean should be 6.50)
> [1] 6.501477
> > 6.501477/6.5
> [1] 1.000227
> 
> 
> Jim Bouldin, PhD
> Research Ecologist
> Department of Plant Sciences, UC Davis
> Davis CA, 95616
> 530-554-1740
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] dimension trouble for a matrix

[Steve Lianoglou]

Hi Erin,

On Jul 23, 2009, at 2:12 PM, Erin Hodgess wrote:


Dear R People:

I'm having trouble with something that should be very simple.

I'm setting up a matrix outside of a loop and writing items into it
during the loop.

Here is the output:


glob3b("sites.info")

dim 27 3
[1] "/raid1/osg-app"
Error in xy[i, ] : incorrect number of dimensions


Here is the function:

glob3b

function(xx) {
 x.df <- read.table(xx,header=FALSE,as.is=TRUE,sep="\t")
n1 <- nrow(x.df)
xy <- matrix(rep("",n1*3),nrow=n1,ncol=3)
  cat("dim",dim(xy),"\n")
for(i in 1:n1) {
   xz <- paste("globus-job-run ",x.df[i,2]," /bin/sh -c 'echo  
$OSG_APP'",

   sep="")
xw <- system(xz,intern=TRUE)
if(length(xw)>0) {
xy[i,1:2] <- x.df[i,1:2]
print(xw)
cat(i,xy[i,],"\n")
cat(i,length(xw),dim(xw),"\n")
xy[i,3] <- xw
}
}
return(xy)
}


This is somehow hard for me to parse, but a few comments:

1. The object returned from your `system` call (xw) will be a  
character vector, calling dim(xw) will return NULL (doesn't look like  
that matters here, just sayin').


2. xw maybe of length > 1, in which case your final call `xy[,3] <-  
xw` is likely tripping you up (try: `xy[,3] <- xw[1]`.


That said, I expect some of the output from your calls to "cat" to be  
shown, but I'm not seeing it in the output you pasted.


Other randomness:

xy <- matrix(rep("",n1*3),nrow=n1,ncol=3)

can be:

xy <- matrix("", nrow=n1, ncol=3)


Also, are you sure you want to store this in a matrix? The fact that  
you want to save your `xw` var (which can be of length > 1) suggests  
to me that you might want to be stuffing whatever-it-is you're saving  
into a list of lists structure instead ... just a thought.


-steve

--
Steve Lianoglou
Graduate Student: Physiology, Biophysics and Systems Biology
Weill Medical College of Cornell University

Contact Info: http://cbio.mskcc.org/~lianos/contact

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

Re: [R] How to perform a calculation in each element of my list?

[Daniel Malter]

Hi Mark, study the following example. The two simulated dataframes are put
in a list called listdata. The loop iterates throught the elements in the
list and multiplies the last column "listdata[[i]][,length(listdata[[i]])]"
by the column before the last "listdata[[i]][,length(listdata[[i]])-1]".
This new columns is assigned to be the new column after the last of the
respective dataframe in the list "listdata[[i]][,length(listdata[[i]])+1]"

x1=rnorm(100,0,1)
x2=rnorm(100,0,1)
x3=rnorm(100,0,1)
x4=rnorm(100,0,1)

data1=data.frame(x1,x2)
data2=data.frame(x3,x4)

listdata=list(data1,data2)

for(i in 1:length(listdata)){
 
listdata[[i]][,length(listdata[[i]])+1]=listdata[[i]][,length(listdata[[i]])
]*listdata[[i]][,length(listdata[[i]])-1]
  }

HTH,
Daniel 


-
cuncta stricte discussurus
-

-Ursprüngliche Nachricht-
Von: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] Im
Auftrag von Mark Na
Gesendet: Thursday, July 23, 2009 2:06 PM
An: r-help@r-project.org
Betreff: [R] How to perform a calculation in each element of my list?

Hi R-helpers,

I have a list containing 10 elements, each of which is a dataframe. I wish
to add a new column to each list element (dataframe) containing the product
of the last two columns of each dataframe.

I'd appreciate any pointers, thanks!

Mark Na

[[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 perform a calculation in each element of my list?

[Jorge Ivan Velez]

Dear Mark,
Try this:

lapply(yourlistofdataframes, function(d){
   k <- ncol(d)
   d$product <- d[,k-1] * d[,k]
   d
}
   )

HTH,

Jorge


On Thu, Jul 23, 2009 at 2:05 PM, Mark Na  wrote:

> Hi R-helpers,
>
> I have a list containing 10 elements, each of which is a dataframe. I wish
> to add a new column to each list element (dataframe) containing the product
> of the last two columns of each dataframe.
>
> I'd appreciate any pointers, thanks!
>
> Mark Na
>
>[[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] dimension trouble for a matrix

[Erin Hodgess]

Dear R People:

I'm having trouble with something that should be very simple.

I'm setting up a matrix outside of a loop and writing items into it
during the loop.

Here is the output:

> glob3b("sites.info")
dim 27 3
[1] "/raid1/osg-app"
Error in xy[i, ] : incorrect number of dimensions


Here is the function:
> glob3b
function(xx) {
  x.df <- read.table(xx,header=FALSE,as.is=TRUE,sep="\t")
n1 <- nrow(x.df)
xy <- matrix(rep("",n1*3),nrow=n1,ncol=3)
   cat("dim",dim(xy),"\n")
for(i in 1:n1) {
xz <- paste("globus-job-run ",x.df[i,2]," /bin/sh -c 'echo $OSG_APP'",
   sep="")
xw <- system(xz,intern=TRUE)
if(length(xw)>0) {
xy[i,1:2] <- x.df[i,1:2]
print(xw)
cat(i,xy[i,],"\n")
cat(i,length(xw),dim(xw),"\n")
xy[i,3] <- xw   
}
}
return(xy)
}
>

I have no idea what's wrong.  This should run like clockwork.  Any
help is much appreciated.  By the way, this is NOT homework!

Thanks,
Sincerely,
Erin



-- 
Erin Hodgess
Associate Professor
Department of Computer and Mathematical Sciences
University of Houston - Downtown
mailto: erinm.hodg...@gmail.com

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

[R] How to perform a calculation in each element of my list?

[Mark Na]

Hi R-helpers,

I have a list containing 10 elements, each of which is a dataframe. I wish
to add a new column to each list element (dataframe) containing the product
of the last two columns of each dataframe.

I'd appreciate any pointers, thanks!

Mark Na

[[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] howto create a list row-by-row as input to function call?

[Mark Knecht]

Hi,
   I'm having trouble within my function CalcPos to get it to call
CalcHorz with values from each row. I *think* it's calling CalcHorz
with the final values of the inputs and not the values from each row.
How can I do this properly in R?

   The values aa,bb,cc,dd are inputs. CalcPos first calculates V1 and
V2 vertically, and then I attempt to call CalcHorz to handle H1, H2 &
H3 using a list L1 comprised of df[aa:V2]. I wanted that list to be
the values on that row and disposed of each pass, but I think all
calls are probably using the final values. Finally I use the returned
values H1:H3 for more vertical calculations in V3.

   How can I get list L1 to be the unique row values for each call to
CalcHorz? L1 would just be tossed when the call to CalcPos is
complete.

   PLEASE - do not combine any of the inputs or return values. I need
all of them aa through V3 and want to keep them separate. Thanks!

Cheers,
Mark


CalcHorz = function (L1) {
# Intended to calculate 3 values for each row
H1 <- with(L1, V1*cc)
H2 <- with(L1, V2*dd)
H3 <- with(L1, H1 - H2)
return(as.list(H1, H2, H3))
}

CalcPos = function (df) {

#calculate vertical initially
df$V1 <- with(df, aa*bb)
df$V2 <- with(df, aa*cc-dd)

#At this point do 3 horizontal calculations
#Goal is this is done for each row using
#the values aa,bb,cc,dd,V1,V2 for that row

L1 = as.list(df[1:dim(df)[2]])

R1 = as.list(with(df, CalcHorz(L1)))
df$H1 <- R1[1]
df$H2 <- R1[2]
df$H3 <- R1[3]

#Switch back to vertical

df$V3 = with (df, cumsum(H3))

return(df)
}


DF <- data.frame(aa=1:12, bb=12:15, cc=1:3, dd=1:6)
DF

DF <- CalcPos(DF)
DF

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

[Jim Bouldin]

Dan Nordlund wrote:

"It would be necessary to see the code for your 'brief test' before anyone
could meaningfully comment on your results.  But your results for a single
test could have been a valid "random" result."

I've re-created what I did below.  The problem appears to be with the
weighting process: the unweighted sample came out much closer to the actual
than the weighted sample (>1% error) did.  Comments?
Jim

> x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12
> weights
 [1] 1 1 1 1 1 1 2 2 2 2 2 2

> a = mean(replicate(100,(sample(x, 3, prob = weights;a  # (1
million samples from x, of size 3, weighted by "weights"; the mean should
be 7.50)
[1] 7.406977
> 7.406977/7.5
[1] 0.987597

> b = mean(replicate(100,(sample(x, 3;b  # (1 million samples from
x, of size 3, not weighted this time; the mean should be 6.50)
[1] 6.501477
> 6.501477/6.5
[1] 1.000227


Jim Bouldin, PhD
Research Ecologist
Department of Plant Sciences, UC Davis
Davis CA, 95616
530-554-1740

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

[Gabor Grothendieck]

Please provide your code in a reproducible form (as requested previously).

Here we fit with first 9 points and then add a point for prediction.  (Of
course your model can only predict the current value of Y so you may
have to rethink your model even aside from the implementation if you
really want to predict future values of Y.)

> library(dyn)
> set.seed(123)
> tt <- ts(cbind(Y = 1:10, x = rnorm(10), z = rnorm(10)))
> L <- function(x, k = 1) lag(x, -k)
> tt.zoo <- as.zoo(tt)
> fit <- dyn$lm(Y ~ L(Y) + L(x, 0:2) + z, tt.zoo[-10, ])
> predict(fit, tt.zoo)
 1  2  3  4  5  6  7  8  9 10 11 12
NA NA  3  4  5  6  7  8  9 10 NA NA




On Thu, Jul 23, 2009 at 1:01 PM, Hongwei Dong wrote:
> Hi, Gabor and Other R users,
>   I'm re-posting my script and the results I got.
>
>   here is the dynamic model I used to estimate in-sample model (1996-2006)
> and it works:
>
>   fit<-dyn$lm(Y~lag(Y,-1)+z+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4))
>   Then I used this model to do out sample forecast with the following
> scripts, which do not work:
>
>  z<-ts(Z[41:52],start=2006,frequency=4)
>  x<-ts(X[41:52],start=2006,frequency=4)
> newdata<-data.frame(cbind(z,x))
> newdata<-ts(newdata,start=2006,frequency=4)
> pred<-predict(fit,newdata)
> Here is the results I got from R:
>          Qtr1     Qtr2     Qtr3     Qtr4
> 2006       NA       NA       NA       NA
> 2007 3083.362       NA       NA       NA
> 2008       NA       NA       NA       NA
> 2009       NA       NA       NA       NA
>
> I got only one prediction for the first quarter in 2007. Intuitively, there
> might be two problems: the definition of the newdata and how to define Y in
> newdata. But I just can't figure this out. It will greatly appreciated if
> someone can give me some help. Thanks.
> Harry
>
>
> On Thu, Jul 23, 2009 at 5:15 AM, Gabor Grothendieck
>  wrote:
>>
>> You have to use consistent classes.  You can't start out using
>> ts class and then suddenly switch to zoo class as if you had been
>> using zoo all along.   Either use ts everywhere or zoo everywhere.
>>
>> Also in the future please post reproducible examples as requested
>> at the bottom of every message to r-help.  That means include
>> a minimal amount of data so we can get exactly what you
>> are getting.
>>
>> On Thu, Jul 23, 2009 at 4:48 AM, Hongwei Dong wrote:
>> > Thanks, Gabor. This is really helpful.
>> > When the regressive part, lag(Y,-1), is not included, my sytax works
>> > well.
>> > However, when I include the lag(Y) in the model, the prediction part
>> > does
>> > not work. Here is my sytax for in-sample estimation and it works well:
>> > fit<-dyn$lm(Y~lag(Y,-1)+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4)+z).
>> > Then I use this model to do out of sample prediction:
>> > x<-ts(X[41:52],start=2006,frequency=4)
>> > z<-ts(Z[41:52],start=2006,frequency=4)
>> > newdata<-data.frame(cbind(x,z))
>> > newdata<-zooreg(newdata)
>> > pred<-predict(fit,newdata)
>> > With these, I got weird result, a prediction for each year from year 1
>> > to
>> > the first quarter of year 2007 (all "NA"). What I expect is a prediction
>> > for
>> > the 8 quarters from 2007 to 2008. Intuitively, I know there must be
>> > something wrong with my newdata definition. But I can't figure it out.
>> > I'll
>> > appreciate it very much if you can give some suggestions to modify my
>> > syntax. Thanks.
>> >
>> > Harry
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > On Wed, Jul 22, 2009 at 10:53 PM, Gabor Grothendieck
>> >  wrote:
>> >>
>> >> Here is an example closer to yours.
>> >>
>> >> > library(dyn)
>> >> > set.seed(123)
>> >> > x <- zooreg(rnorm(10))
>> >> > y <- zooreg(rnorm(10))
>> >> > L <- function(x, k = 1) lag(x, k = -k)
>> >> > mod <- dyn$lm(y ~ L(y) + L(x, 0:2))
>> >> > mod
>> >>
>> >> Call:
>> >> lm(formula = dyn(y ~ L(y) + L(x, 0:2)))
>> >>
>> >> Coefficients:
>> >> (Intercept)         L(y)   L(x, 0:2)1   L(x, 0:2)2   L(x, 0:2)3
>> >>    0.06355     -0.74540      0.63649      0.44957     -0.41360
>> >>
>> >> > newdata <- cbind(x = c(coredata(x), rnorm(1)), y = c(coredata(y),
>> >> > rnorm(1)))
>> >> > newdata <- zooreg(newdata)
>> >> > predict(mod, newdata)
>> >>         1          2          3          4          5          6
>> >>  7
>> >>        NA         NA  0.9157808  0.6056333 -0.5496422  1.5984615
>> >> -0.2574875
>> >>         8          9         10         11         12         13
>> >> -1.6148859  0.3329285 -0.5284646 -0.1799693         NA         NA
>> >>
>> >>
>> >> On Thu, Jul 23, 2009 at 1:04 AM, Gabor
>> >> Grothendieck wrote:
>> >> > Use dyn.predict like this:
>> >> >
>> >> >> library(dyn)
>> >> >> x <- y <- zoo(1:5)
>> >> >> mod <- dyn$lm(y ~ lag(x, -1))
>> >> >> predict(mod, list(x = zoo(6:10, 6:10)))
>> >> >  7  8  9 10
>> >> >  7  8  9 10
>> >> >
>> >> >
>> >> > On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
>> >> > wrote:
>> >> >> I have a dynamic time series model like this:
>> >> >> dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
>> >> >>
>> >> >> I need 

[R] How to pass a character argument which contains expressions to arg.names in barplot?

[jcano]

Hi all

Can anybody help me with this? I am trying to include in an automatic way
the argument in arg.names in a barplot. I generate the labels I want to
appear below the bars with a for loop, and they contain subscripts, so I
need to use expression

  anch<-0.05
  esp<-4
  for (i in 1:dim(Ntot)[1])
{
naux<-Ntot[i,]
naux2<-naux[naux>0]
nind<-which(naux>0)
tit4<-character(0)
for (j in 1:length(nind))
  {
  tit4<-c(tit4,paste("expression(n[paste[",i,",",nind[j],"])",sep=""))
  }
windows()
barplot(naux2,xlab=eval(expression(substitute(n[i],list(i=i,
   
arg.names=tit4,col="gray",width=c(anch,anch),axes=TRUE,xlim=c(0,anch*(esp+2)*length(nind)),
space=esp,ylim=c(0,max(naux2)),main=paste("State
#",i),yaxp=c(0,M,Mint));lines(c(0,100),c(0,0))
}   # end of for
 
but I don't get what I expect. R plots literally the contents of tit4 below
each bar, i.e., (for the last value of i in the outer for loop)

"expression(n[paste[5,1])" 
"expression(n[paste[5,2])" 
"expression(n[paste[5,3])" 
"expression(n[paste[5,4])"

Of course, if I write directly
arg.names=c(expression(n[paste(5,1)]),expression(n[paste(5,2)]),expression(n[paste(5,3)]),expression(n[paste(5,4)]))
it works, but then I cannot do ii in an automatic way for all the graphics
coming from the for (i in 1:dim(Ntot)[1]).

I have also tried to store in another variable tit5 the whole expression as
> tit5
[1]
"c(expression(n[paste[5,1]),expression(n[paste[5,2]),expression(n[paste[5,3]),expression(n[paste[5,4]))"
but it doesn't work neither

I wonder if there is some way that barplot "tells" arg.names to "evaluate"
the contents of tit4 as expressions, not as characters

I hope I made myself clear enough

Thanks in advance to any response

Cheers

Javi
-- 
View this message in context: 
http://www.nabble.com/How-to-pass-a-character-argument-which-contains-expressions-to-arg.names-in-barplot--tp24630656p24630656.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] simple question about exporting data...

[John Kane]

I vote for it. 

--- On Thu, 7/23/09, Marc Schwartz  wrote:

> From: Marc Schwartz 
> On Jul 23, 2009, at 8:59 AM, Greg Snow wrote:
> 
> > Doing the computations in R then the graphs in Excel
> reminds me of the maxim:
> > 
> > Measure with a micrometer
> > Mark with chalk
> > Cut with an ax

> 
> Definitely a fortunes candidate...
> 
> Marc Schwartz
> 



  __
[[elided Yahoo spam]]

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

Re: [R] error message: .Random.seed is not an integer vector but

[Nordlund, Dan (DSHS/RDA)]

> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf Of Jim Bouldin
> Sent: Thursday, July 23, 2009 9:49 AM
> To: ted.hard...@manchester.ac.uk; ted.hard...@manchester.ac.uk; r-h...@r-
> project.org
> Subject: Re: [R] error message: .Random.seed is not an integer vector but
> 
> 
> Thanks much Ted.  I actually had just tried what you suggest here before
> you posted, and resolved the problem.  Thanks also for the other tips.  I
> wrote x = as.vector(c(1:12)) because I thought that the mode of x might be
> the problem, the error message pointing to .Random.seed notwithstanding.
> 
> On a related note, I did a brief test a couple weeks back where I ran a
> million random samples of 3 from the vector 1:12 and compared the mean
> against the known mean.  It was off by 1 percent, which indicated that the
> RNG was biased more than I'd have thought.  Comments?
> Jim
> 
<<>>

It would be necessary to see the code for your 'brief test' before anyone could 
meaningfully comment on your results.  But your results for a single test could 
have been a valid "random" result.

Dan

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

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

Re: [R] how to predict dynamic model in R

[Hongwei Dong]

Hi, Gabor and Other R users,
  I'm re-posting my script and the results I got.

  here is the dynamic model I used to estimate in-sample model (1996-2006)
and it works:

  fit<-dyn$lm(Y~lag(Y,-1)+z+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4))

  Then I used this model to do out sample forecast with the following
scripts, which do not work:

 z<-ts(Z[41:52],start=2006,frequency=4)
 x<-ts(X[41:52],start=2006,frequency=4)

newdata<-data.frame(cbind(z,x))
newdata<-ts(newdata,start=2006,frequency=4)
pred<-predict(fit,newdata)

Here is the results I got from R:

 Qtr1 Qtr2 Qtr3 Qtr4
2006   NA   NA   NA   NA
2007 3083.362   NA   NA   NA
2008   NA   NA   NA   NA
2009   NA   NA   NA   NA


I got only one prediction for the first quarter in 2007. Intuitively, there
might be two problems: the definition of the newdata and how to define Y in
newdata. But I just can't figure this out. It will greatly appreciated if
someone can give me some help. Thanks.

Harry



On Thu, Jul 23, 2009 at 5:15 AM, Gabor Grothendieck  wrote:

> You have to use consistent classes.  You can't start out using
> ts class and then suddenly switch to zoo class as if you had been
> using zoo all along.   Either use ts everywhere or zoo everywhere.
>
> Also in the future please post reproducible examples as requested
> at the bottom of every message to r-help.  That means include
> a minimal amount of data so we can get exactly what you
> are getting.
>
> On Thu, Jul 23, 2009 at 4:48 AM, Hongwei Dong wrote:
> > Thanks, Gabor. This is really helpful.
> > When the regressive part, lag(Y,-1), is not included, my sytax works
> well.
> > However, when I include the lag(Y) in the model, the prediction part does
> > not work. Here is my sytax for in-sample estimation and it works well:
> > fit<-dyn$lm(Y~lag(Y,-1)+x+lag(x,-1)+lag(x,-2)+lag(x,-3)+lag(x,-4)+z).
> > Then I use this model to do out of sample prediction:
> > x<-ts(X[41:52],start=2006,frequency=4)
> > z<-ts(Z[41:52],start=2006,frequency=4)
> > newdata<-data.frame(cbind(x,z))
> > newdata<-zooreg(newdata)
> > pred<-predict(fit,newdata)
> > With these, I got weird result, a prediction for each year from year 1 to
> > the first quarter of year 2007 (all "NA"). What I expect is a prediction
> for
> > the 8 quarters from 2007 to 2008. Intuitively, I know there must be
> > something wrong with my newdata definition. But I can't figure it out.
> I'll
> > appreciate it very much if you can give some suggestions to modify my
> > syntax. Thanks.
> >
> > Harry
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > On Wed, Jul 22, 2009 at 10:53 PM, Gabor Grothendieck
> >  wrote:
> >>
> >> Here is an example closer to yours.
> >>
> >> > library(dyn)
> >> > set.seed(123)
> >> > x <- zooreg(rnorm(10))
> >> > y <- zooreg(rnorm(10))
> >> > L <- function(x, k = 1) lag(x, k = -k)
> >> > mod <- dyn$lm(y ~ L(y) + L(x, 0:2))
> >> > mod
> >>
> >> Call:
> >> lm(formula = dyn(y ~ L(y) + L(x, 0:2)))
> >>
> >> Coefficients:
> >> (Intercept) L(y)   L(x, 0:2)1   L(x, 0:2)2   L(x, 0:2)3
> >>0.06355 -0.74540  0.63649  0.44957 -0.41360
> >>
> >> > newdata <- cbind(x = c(coredata(x), rnorm(1)), y = c(coredata(y),
> >> > rnorm(1)))
> >> > newdata <- zooreg(newdata)
> >> > predict(mod, newdata)
> >> 1  2  3  4  5  6
> >>  7
> >>NA NA  0.9157808  0.6056333 -0.5496422  1.5984615
> >> -0.2574875
> >> 8  9 10 11 12 13
> >> -1.6148859  0.3329285 -0.5284646 -0.1799693 NA NA
> >>
> >>
> >> On Thu, Jul 23, 2009 at 1:04 AM, Gabor
> >> Grothendieck wrote:
> >> > Use dyn.predict like this:
> >> >
> >> >> library(dyn)
> >> >> x <- y <- zoo(1:5)
> >> >> mod <- dyn$lm(y ~ lag(x, -1))
> >> >> predict(mod, list(x = zoo(6:10, 6:10)))
> >> >  7  8  9 10
> >> >  7  8  9 10
> >> >
> >> >
> >> > On Thu, Jul 23, 2009 at 12:54 AM, Hongwei Dong
> wrote:
> >> >> I have a dynamic time series model like this:
> >> >> dyn$lm( y ~ lag(y,-1) + x + lag(x,-1)+lag(x,-2) )
> >> >>
> >> >> I need to do an out of sample forecast with this model. Is there any
> >> >> way I
> >> >> can do this with R?
> >> >> It would be greatly appreciated if some one can give me an example.
> >> >> Thanks.
> >> >>
> >> >>
> >> >> Harry
> >> >>
> >> >>[[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-

[R] How to get w in SVR with e1071 package

[marlene marchena]

>
> Hi all,
>
> I need some help about how to calculate w in a SVR in package e1071.
>
> I have a regression y_i=f(x_i)+e
>
> where f(*x*)=(w,phi(x))+b
>
> then go on with the SVR calculation I know that w*=Sum_i=1^n [(á_i -
> á*_i)K(x,x_i) ] where á_i and á*_i are the lagrangian multipliers of the
> dual form.
>
> o.k but how I will get it in R?
>
> When I run my model I have information about the model but I don't know how
> to use it to get w. Where are the alphas? Is w a vector or a value?
>
> Someone can help me please.
>
> Regards,
>
> Marlene
>
> svm.m1 <- svm(q ~ ., data = train, cost = 100, gamma = 1e-05)
>
> svm.m1
>
> Call:
> svm(formula = q ~ ., data = train, cost = 100, gamma = 1e-05)
>
>
> Parameters:
>SVM-Type:  eps-regression
>  SVM-Kernel:  radial
>cost:  100
>   gamma:  1e-05
> epsilon:  0.1
>
>
> Number of Support Vectors:  187
>
>
>
>
>
>

[[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] error message: .Random.seed is not an integer vector but

[Jim Bouldin]

Thanks much Ted.  I actually had just tried what you suggest here before
you posted, and resolved the problem.  Thanks also for the other tips.  I
wrote x = as.vector(c(1:12)) because I thought that the mode of x might be
the problem, the error message pointing to .Random.seed notwithstanding.

On a related note, I did a brief test a couple weeks back where I ran a
million random samples of 3 from the vector 1:12 and compared the mean
against the known mean.  It was off by 1 percent, which indicated that the
RNG was biased more than I'd have thought.  Comments?
Jim
 
> Follow-up to my previous reply (just posted). Having read the other
> responses and your reactions, try the following:
> 
>   rm(.Random.seed)
>   set.seed(54321) ## (Or your favourite magic number) [*]
>   x = as.vector(c(1:12))  ## To reproduce your original code ... !
>   sample(x,3)
> 
> [*] When you did rm(.Random.seed) as suggested by Uwe, the variable
> .Random.seed was lost, so you have to create it again.
> 
> If, after the above, you still get the problem, then something is
> very seriously wrong.
> 
> Ted.
> 
> 
> E-Mail: (Ted Harding) 
> Fax-to-email: +44 (0)870 094 0861
> Date: 23-Jul-09   Time: 17:23:09
> -- XFMail --
> 

Jim Bouldin, PhD
Research Ecologist
Department of Plant Sciences, UC Davis
Davis CA, 95616
530-554-1740

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

[Greg Snow]

I am not sure that I fully understand what all you want to do (and I don't 
understand why you need the correlation and if a correlation based on 3 pairs 
is even meaningful), but here is a first stab at what you are trying to do:

tmp <- "Species Control_CR Damage_DR
A 10 2
A 9 3
A 7 4
A 9 2
A 8 3
B 5 4
B 6 3
B 4 2
B 5 4
B 6 3
C 8 1
C 6 2
C 4 4
C 7 1
C 6 2"

plants <- read.table( textConnection(tmp), header=TRUE )

plants2 <- rbind( data.frame( Species=plants$Species, growth=plants[[2]], 
group=1 ), 
data.frame(Species=plants$Species, growth=plants[[3]], group=2 ) )

tmpfun <- function() {
group2 <- ave( plants2$group, plants2$Species, FUN=sample )
cr <- tapply( plants2[ group2==1, 'growth'], plants2[group2==1, 
'Species'],
FUN=mean)
dr <- tapply( plants2[ group2==2, 'growth'], plants2[group2==2, 
'Species'],
FUN=mean)

ir <- cr - dr
cor( ir, cr )
}

out <- replicate( 1000, tmpfun() )

hist(out)


Is this what you wanted? Or at least helpful for going in the right direction?


-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Anne Kempel
> Sent: Thursday, July 23, 2009 9:02 AM
> To: r-help@r-project.org
> Subject: [R] help with randomisation test...
> 
> Dear R-people,
> I hope asking this is not too cheeky, but I do have a R Problem. I hope
> that some of you like to play around with R and can help me.
> 
> Its like this. I have several plant species (A,B,C) and 10 replicates
> per species. 5 plants per species are damaged, 5 not. I let a
> caterpillar feed on each plant and measured the growth of the
> caterpillars on control plants (CR) and on damaged plants (DR). The
> difference, DR-CR is an indicator for induced resistance IR.
> 
> The data could look like this...just as an example.
> 
> Species   Control_CR  Damage_DR
> A 10  2
> A 9   3
> A 7   4
> A 9   2
> A 8   3
> B 5   4
> B 6   3
> B 4   2
> B 5   4
> B 6   3
> C 8   1
> C 6   2
> C 4   4
> C 7   1
> C 6   2
> 
> 
> Now, i want to see if there is a trade-off (or a negative correlation)
> between CR and IR, but because IR is calculated out of CR there is a
> danger of spurious correlation. What I have to do to get rid of it is
> the following:
> 
> 1)Draw with replacement sets of 5 plants per treatment for each species
> 
> 2)Calculate of those sets the means (in the control and the damage
> treatment) and then abstract the mean of the damage minus the mean of
> the control to get a mean for IR.
> 
> 3)Randomize this IR (now called IRrandom)among all the species.
> 4)Go back to the data, calculate:  Each of the 5 CR values minus this
> IRrandom to get 5 new DR, then draw again with replacements sets of 5
> plants of each treatment for each species and calculate new means for
> the CR and the new DR. Now you have a mean CR and a new mean DR - you
> can calculate DR-CR to get a new IR (IRnew).
> 
> 5)Correlate this CR with the new IR. And store the result.
> 
> 6)Repeat step 1-5 many times. And make a histogram of the resulting
> corelation coefficients.
> 
> 
> You see, this is quite complicated. Does any of you have any idea how
> to
> start? I am very very grateful for any suggestions!
> The helper will get a parcel with swiss chocolate ;-)  !!
> 
> Best wishes from Bern,
> Anne
> 
> 
> --
> Anne Kempel
> Institute of Plant Sciences
> University of Bern
> Altenbergrain 21
> CH-3103 Bern
> kem...@ips.unibe.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.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] using k-means clustering in conjunction with heatmap.2 function

[Suzanne Matthews]

Hello,

I am trying to create a heatmap that clusters based on a k-means scheme
rather than a hierarchical clustering scheme.

Suppose I have the following input data, located in sample.table:
x1 x2 x3 x4
x1 17.198 16.306 16.806 16.374
x2 14.554 10.866 15.780 14.596
x3 14.374 14.118 14.569 17.352
x4 17.505 14.596 15.738 14.070

By using the heatmap.2 function as follows, I can create a heatmap using
hierarchical clustering:
library(gplots)
x=read.table('sample.table', header=TRUE)
mat=data.matrix(x)
heatmap.2(mat, col=redgreen(75), trace="none")

However, I am now interested in trying to show a heatmap in which the data
is clustered using k-means. By looking at some previous posts on the mailing
lists, and the documentation, I've attempted to cluster my data using kmeans
using the following modified script:

x=read.table('sample.table', header=TRUE)
mat=data.matrix(x)
km <-kmeans(t(scale(t(mat))), 3); km$cluster

heatmap.2(mat, col=redgreen(75), hclustfun=km, trace="none")

However, I'm not sure what I need to modify in the heatmap.2 list of
parameters to have the matrix reflect a k-means clustering. The above usage
will throw an error, since "km" is not a valid heirarchacal clustering
function. My question is, what should the heatmap.2() function look like in
order for it to show kmeans clustering?

I'm not even sure if I'm using the kmeans function in a way that's optimal
for producing heatmaps. I saw a post related to something like this in June
of 2006, but I was not able to get any further than I what I've shown you
above. I would greatly appreciate it if someone could point me in the right
direction.

One part of me also isn't even sure if k-means makes sense in this context.
I would really appreciate your opinions related to this.

Thanks in advance for your time!

-SM

[[alternative HTML version deleted]]

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

Re: [R] read.csv greater than

[stephen sefick]

here is the solution
x <- subset(read.csv("foo.csv"), VALUE>0)

On Thu, Jul 23, 2009 at 11:28 AM, stephen sefick wrote:
> I have a csv file that I am trying to read in and know that values <0
> are erroneous - is there a way to read only value grater than 0.
>
> thanks,
>
> --
> Stephen Sefick
>
> Let's not spend our time and resources thinking about things that are
> so little or so large that all they really do for us is puff us up and
> make us feel like gods.  We are mammals, and have not exhausted the
> annoying little problems of being mammals.
>
>                                                                -K. Mullis
>



-- 
Stephen Sefick

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

-K. Mullis

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

[R] read.csv greater than

[stephen sefick]

I have a csv file that I am trying to read in and know that values <0
are erroneous - is there a way to read only value grater than 0.

thanks,

-- 
Stephen Sefick

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

-K. Mullis

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

Re: [R] error message: .Random.seed is not an integer vector but

[Ted Harding]

On 23-Jul-09 16:08:17, Jim Bouldin wrote:
>> Jim Bouldin wrote:
>> > Thank you.  However, when I tried that, I got this message:
>> 
>> > Warning message:
>> > In rm(.Random.seed) : variable ".Random.seed" was not found
>> 
>> In that case, have you attached some package that has its own
>> .Random.seed?
>> Try to find where the current .random.seed comes from R complains
>> about.
>> 
>> Uwe Ligges
> 
> No, there are no attached packages, just the ones that load
> automatically. The R commander has some type of RNG but it is not
loaded.  Completely stumped.

Follow-up to my previous reply (just posted). Having read the other
responses and your reactions, try the following:

  rm(.Random.seed)
  set.seed(54321) ## (Or your favourite magic number) [*]
  x = as.vector(c(1:12))  ## To reproduce your original code ... !
  sample(x,3)

[*] When you did rm(.Random.seed) as suggested by Uwe, the variable
.Random.seed was lost, so you have to create it again.

If, after the above, you still get the problem, then something is
very seriously wrong.

Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 23-Jul-09   Time: 17:23:09
-- 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] error message: .Random.seed is not an integer vector but

[Ted Harding]

On 23-Jul-09 15:30:05, Jim Bouldin wrote:
> I'm trying to run this simple random sample procedure and keep
> getting the error message shown. I don't understand this; I've
> Thanks.
> 
>> x = as.vector(c(1:12));x
>  [1]  1  2  3  4  5  6  7  8  9 10 11 12
>> mode(x)
> [1] "numeric"
>> sample(x, 3)
> Error in sample(x, 3) : 
>   .Random.seed is not an integer vector but of type 'list'
> 
> Jim Bouldin, PhD

I think this error is not arising from the code you have given above.
I get:

  x=as.vector(c(1:12))
  sample(x,3)
  # [1]  8 12  2

which is fine. '.Random.seed' should be an integer vector to start with.
See '?set.seed'. Therefore it should not be a list. Try entering

  .Random.seed

and see what you get -- when I do it it is a vector of 626 assorted
integers. Also try:

  is.vector(.Random.seed)

and, if that is not TRUE, try

  str(.Random.seed)
  typeof(.Random.seed)

If .Random.seed turns out to be a list, then something has gone
seriously wrong somewhere. It may be that somwehere in your code
there is a command that meddles with .Random.seed (or perhaps it
was saved, in .Rdata, from a previous session where something
went wrong). Or, worst scenario, your R installation is broken ...

And, while I'm at it, your "x=as.vector(c(1:12))" is overkill!
Simply

  x <- (1:12) ; sample(x,3)

should be enough, or even

  sample((1:12),3)

since (1:12) is a vector (of integers).

Hoping this helps,
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 23-Jul-09   Time: 17:13:55
-- 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] error message: .Random.seed is not an integer vector but of type 'list'

[Jim Bouldin]

> 
> 
> Jim Bouldin wrote:
> > Thank you.  However, when I tried that, I got this message:
>  >
> > Warning message:
> > In rm(.Random.seed) : variable ".Random.seed" was not found
> 
> 
> In that case, have you attached some package that has its own
> .Random.seed?
> Try to find where the current .random.seed comes from R complains about.
> 
> Uwe Ligges

No, there are no attached packages, just the ones that load automatically.
The R commander has some type of RNG but it is not loaded.  Completely stumped.

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