Re: [R] Pseudo R2 for Tobit Regression

[Achim Zeileis]

On Sat, 16 Aug 2008, ldotero wrote:



Dear All:
I need some guidance in calculating a goodness-of-fit statistic for a Tobit
Regression model.
To develop the Tobit regression, I used the tobit() method from the AER
package, which is basically a simpler interface to the survreg() method.

I've read about pseudo R2 and C-index and was wondering if there is a
package that calculates this for me.  Also, is there a reason to select
pseudo R2 over C-Index, or viceversa?


McFadden's pseudo-R^2 can be computed via
  fm <- tobit(...)
  fm0 <- update(fm, . ~ 1)
  1 - as.vector(logLik(fm)/logLik(fm0))
using accessor functions or via computing on the internal structure via
  1 - fm$loglik[2]/fm$loglik[1]

hth,
Z


I will appreciate any suggestions.

Many thanks.
--
View this message in context: 
http://www.nabble.com/Pseudo-R2-for-Tobit-Regression-tp19014572p19014572.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] continuous coloring of a polygon

[Roger Leenders]

WinXP, R2.7.1

Thanks so much to all who have reponded to my inquery. The solutions are 
most helpful.

There is only one final thing that I can't seem to get right.
All of the proposed solutions yield a figure that is not quite round, 
but is wider than its height. This can be resolved by manually resizing 
the window in which it is plotted, but that becomes cumbersome when many 
of these figures need to be generated. Is there a way to enforce the 
figure to be perfectly round?
I have tried plotting directly to a windows device (using for example 
the code below) or to a pdf-device in which I specifically specify the 
square size of the plot, but the situation remains in which the figure 
is wider than it is high. Does anyone know how to solve this final piece 
of the puzzle?


Thanks, Roger


example code:

#pdf(file = ifelse(onefile, "Rplots.pdf", "Rplot%03d.pdf"), 
width=4,height=4,family="Palatino", title=c("De speedometer"), fonts = 
"AvantGarde",paper = "a4")


radius <- 3
x <- seq(-radius,radius,length=2000)
y <- sqrt(radius^2-x^2)
xx <- c(x,-x)
yy <- c(y,-y)
plot(xx,yy, xlim=c(-radius,radius),ylim=c(-radius,radius), type="l", 
ylab="", xlab="", axes=F)


tmp <- rainbow(1000, start=0/6, end=2/6)
theta <- seq(pi, 0, length=1000)

segments(2.7*cos(theta),2.7*sin(theta),
   2.0*cos(theta), 2.0*sin(theta), col=tmp, lwd=2)
dev.off()




Greg Snow schreef:

Here are a couple of other solutions, use whichever works best for you (after 
modifications, changing increments, etc.).

radius <- 3
x <- seq(-radius,radius,length=2000)
y <- sqrt(radius^2-x^2)
xx <- c(x,-x)
yy <- c(y,-y)
plot(xx,yy, xlim=c(-radius,radius),ylim=c(-radius,radius), type="l", ylab="", 
xlab="", axes=F)

radius <- 2.7
x1 <- seq(-radius,radius,length=2000)
y1 <- sqrt(radius^2-x1^2)
radius <- 2.0
x2 <- seq(radius,-radius,length=2000)
y2 <- sqrt(radius^2-x2^2)

#tmp <- rainbow(1000, start=2/6, end=0/6)
tmp <- rev(rainbow(1000, start=0/6, end=2/6))
theta <- seq(pi, 0, length=1000)

segments(2.7*cos(theta),2.7*sin(theta),
2.0*cos(theta), 2.0*sin(theta), col=tmp, lwd=2)

polygon(c(x1,x2),c(y1,y2))


library(TeachingDemos)



radius <- 3
x <- seq(-radius,radius,length=2000)
y <- sqrt(radius^2-x^2)
xx <- c(x,-x)
yy <- c(y,-y)
plot(xx,yy, xlim=c(-radius,radius),ylim=c(-radius,radius), type="l", ylab="", 
xlab="", axes=F)

radius <- 2.7
x1 <- seq(-radius,radius,length=2000)
y1 <- sqrt(radius^2-x1^2)
radius <- 2.0
x2 <- seq(radius,-radius,length=2000)
y2 <- sqrt(radius^2-x2^2)



tmpfun <- function(...){
   image( seq(-3,3,length=101), c(-3.24,3.24), matrix( 1:100, ncol=1 ),
col=rev(rainbow(100,start=0, end=1/3)), add=TRUE )
}

top <- approxfun( x1, y1 )
bottom <- approxfun( c(-3, x2, 3), c(min(y2), y2, min(y2) ) )

xx <- seq(-2.7,2.7, length.out=1000)
xxx <- embed(xx,2)[,2:1]
for(i in 1:999) {
clipplot(tmpfun(), xxx[i,], c( min(bottom(xxx[i,])), max(top(xxx[i,]))) 
)
}

polygon(c(x1,x2),c(y1,y2))





Hope this helps,

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



  

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Roger Leenders
Sent: Friday, August 15, 2008 6:01 AM
To: r-help@r-project.org
Subject: [R] continuous coloring of a polygon


R2.7.1, WinXP

Hi,

I have a polygon inside a circle as follows:

radius <- 3
x <- seq(-radius,radius,length=2000)
y <- sqrt(radius^2-x^2)
xx <- c(x,-x)
yy <- c(y,-y)
plot(xx,yy, xlim=c(-radius,radius),ylim=c(-radius,radius),
type="l", ylab="", xlab="", axes=F)

radius <- 2.7
x1 <- seq(-radius,radius,length=2000)
y1 <- sqrt(radius^2-x1^2)
radius <- 2.0
x2 <- seq(radius,-radius,length=2000)
y2 <- sqrt(radius^2-x2^2)

polygon(c(x1,x2),c(y1,y2))

(the graph much resembles a speed dial inside a car).
Now I want to fill the polygon with color, such that it
starts on the left with red and ends on the right with green,
following the coloring of the rainbow.
Preferably, the coloring should be "continuous", such that
colors naturally fade into each other.
I can draw the polygon as above, but I don't know how to do
the coloring. It is easy to give the polygon only one color
(e.g. through polygon(c(x1,x2),c(y1,y2), col="red")), but I
need a way in which to color the polygon such that the color
moves through the color spectrum from red (left) to green (right).
Can anyone help me to achieve this?

Thanks, Roger

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

[Barry Rowlingson]

2008/8/17 Roger Leenders <[EMAIL PROTECTED]>:
> WinXP, R2.7.1
>
> Thanks so much to all who have reponded to my inquery. The solutions are
> most helpful.
> There is only one final thing that I can't seem to get right.
> All of the proposed solutions yield a figure that is not quite round, but is
> wider than its height. This can be resolved by manually resizing the window
> in which it is plotted, but that becomes cumbersome when many of these
> figures need to be generated. Is there a way to enforce the figure to be
> perfectly round?
> I have tried plotting directly to a windows device (using for example the
> code below) or to a pdf-device in which I specifically specify the square
> size of the plot, but the situation remains in which the figure is wider
> than it is high. Does anyone know how to solve this final piece of the
> puzzle?
>

 add 'asp=1' to your plot() function to specify a 1:1 aspect ratio.
Should do the trick.

Barry

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

Re: [R] Lattice: problem using panel.superpose and panel.groups

[Dieter Menne]

Michael Braun  MIT.EDU> writes:

> 
> Dieter:
> 
> Thank you for your response.  As you requested, I created a self- 
> running example, pasted below.  It may be a little wordier than I  
> would like, but it runs.

.. Details removed
> 
> panel.ppc.plot <- function(...,group.number) {
> 
>if (group.number==1) {
>  panel.bwplot(...)
>} else {
> 
>  panel.lines(...)
>}
> }

Trellis graphics are a bit like hash functions: you can be close to the 
target, but get a far-off result. I admit that I do not know why 
panel.lines does not work in the above example, but  

panel.polygon(...)

works in your special case of ordered data. More generally, I would suggest 
to use

panel.xyplot(x,y,type="l")

or, if you want the ... notation, 

panel.xyplot(...)

but then you have to set type="l" in your bwplot calling function.

Dieter

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

[R] package building problem on windows

[Christophe Dutang]

Hi,

I'm trying to compile the package biglm, but when I build it with R  
CMD build biglm, it failed :


C:\LOCAL\c-dutang\code\R\biglm2>R CMD build biglm
* checking for file 'biglm/DESCRIPTION' ... OK
* preparing 'biglm':
* checking DESCRIPTION meta-information ...C:/DOCUME~1/c-dutang/Local:  
Can't op

n C:/DOCUME~1/c-dutang/Local: No such file or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
 OK
* cleaning src
* removing junk files
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
Error: cannot open file 'biglm/DESCRIPTION' for reading

But when I install it, the installation is performed but I still have  
'can't open' and 'no such file or directory'.


C:\LOCAL\c-dutang\code\R\biglm2>R CMD install biglm
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No  
such file

or directory
installing to ''


-- Making package biglm 
  adding build stamp to DESCRIPTION
  installing NAMESPACE file and metadata
  making DLL ...
gcc  -std=gnu99  -IC:/PROGRA~1/R/R-2.7.1/include -O3 -Wall  -c  
boundedQR.c -

o boundedQR.o
gfortran   -O3  -c boundedQRf.f -o boundedQRf.o
windres --preprocessor="gcc -E -xc -DRC_INVOKED" -I C:/PROGRA~1/R/ 
R-2.7.1/includ

e  -i biglm_res.rc -o biglm_res.o
gcc  -std=gnu99  -shared -s  -o biglm.dll biglm.def boundedQR.o  
boundedQRf.o big

lm_res.o  -LC:/PROGRA~1/R/R-2.7.1/bin   -lgfortran -lR
  ... DLL made
  installing DLL
  installing R files
  installing man source files
  installing indices
  installing help
 >>> Building/Updating help pages for package 'biglm'
 Formats: text html latex example chm
  adding MD5 sums

* DONE (biglm)

Actually, I have no directory local in 'document and settings' but a  
directory local settings exists. What's wrong with my installation? I  
install all tools as recommanded on Murdoch's page. Maybe R CMD build  
tries to use a temp file? in 'document and settings'?


I use R 2.7.1 with Windows xp pro.

Thanks in advance

Christophe

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

Re: [R] Wichmann-Hill Random Number Generator and the Birthday Problem

[Duncan Murdoch]

Shengqiao Li wrote:

Dear all,

Recently I am generating large random samples (10M) and any duplicated 
numbers are not desired.
We tried several RNGs in R and found Wichmann-Hill did not produce 
duplications.


The duplication problem is the interesting birthday problem. If there are 
M possible numbers, randomly draw N numbers from them,

the average number of dupilcations D = N(N-1)/2/M.

For Knuth-TAOCP and Knuth-TAOCP-2002, M=2^30, since this modulus is used. 
D = 46566.12 for N=10M samples.


For Marsaglia-Multicarry, Super-Duper and Mersene-Twister, M=2^32. D = 
11641.53 for N = 10M samples.


My testing results (see below) agree with above analysis. But for 
Wichmann-Hill, it wasn't. Wichmann-Hill's cycle is 6.9536e12 (refer to RNG 
help by ?RNG and Whichmann's correction in 1984). Thus M <= 6.9536e12. D 
  
= 7.19052 when N=1e7 and D>= 179.763 for N=5e7. 


But in my tests, duplications were surprisingly not observed.

It seems that Wichmann-Hill has a much larger cycle than the one 
documented!


Anybody can solve this puzzle?
  


As I told you, look at the code. 

The cycle length of Wichmann-Hill in R appears to be 2.8e13, and that 
also appears to be M in your notation.  Your birthday problem 
calculation does not apply here.
You could probably get a good approximation to it by rounding the 
Wichmann-Hill output (e.g. look at round(2^30 * runif()), or maybe more 
severe rounding).


M is larger than what was previously documented, but Brian Ripley has 
corrected the documentation.


Duncan Murdoch

Regards,

Shengqiao Li

Department of Statistics
West Virgina Unversity

==Testing===

  

RNGkind(kind="Knuth-TAOCP");
RNGkind();


[1] "Knuth-TAOCP" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 46216
  

RNGkind(kind="Knuth-TAOCP-2002");
RNGkind();


[1] "Knuth-TAOCP-2002" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 46373
  

  

RNGkind(kind="Marsaglia-Multicarry");
RNGkind();


[1] "Marsaglia-Multicarry" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 11671
  

RNGkind(kind="Super-Duper");
RNGkind();


[1] "Super-Duper" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 11704
  

RNGkind(kind="Mersenne-Twister");
RNGkind();


[1] "Mersenne-Twister" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 11508
  

RNGkind(kind="Wichmann-Hill");
RNGkind();


[1] "Wichmann-Hill" "Inversion"
  

sum(duplicated(runif(1e7)));


[1] 0

  

gc()


   used (Mb) gc trigger  (Mb)  max used   (Mb)
Ncells  199157  5.4 646480  17.3  10149018  271.1
Vcells 4497151 34.4  108714629 829.5 169585390 1293.9

  

sum(duplicated(runif(5e7)))


[1] 0


  

sessionInfo()


R version 2.7.1 (2008-06-23)
i386-pc-mingw32

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


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

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


==End of Testing===

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Making use of names of viewports (grid)

[Patrick Connolly]

The following code, though not brilliant, works on an A4 page.  It
might look odd on other devices of a very different size.

=X8--- cut here 
  require(grid)
  wide <- 15
  vps <- grid.layout(nrow = 3, ncol = 4,
 widths = unit(rep(1, 4), rep("null", 4)),
 heights = unit(c(99, 1, 99),
   c("mm", "null", "mm")))
  pushViewport(viewport(layout = vps))
  for(k in 1:4){# label 4 viewports in top row
cube.k <- ppaste("Cube", k)
pushViewport(viewport(layout = vps, name = cube.k,
  layout.pos.row = 1,
  layout.pos.col = k))
grid.rect(gp = gpar(lty = "dashed", lwd = .1))
grid.text(cube.k, y = .9, gp = gpar(cex = .8))
popViewport()
  }
## label first viewport in bottom row
  pushViewport(viewport(layout = vps, name = "Cube5",
layout.pos.row = 3,
layout.pos.col = 1))
  grid.rect(gp = gpar(lty = "dashed", lwd = .1))
  grid.text("Cube5", y = .95, gp = gpar(cex = .8))
## Set up a grid inside Cube5 viewport
  vpc <- grid.layout(nrow = 5, ncol = 4,
 widths = unit(rep(wide, 4), "mm"),
 heights = unit(rep(wide, 5), "mm"))
  pushViewport(viewport(layout = vpc))
  for(i in 1:5){
for(j in 1:4){
  pushViewport(viewport(layout = vpc,  clip = "on",
layout.pos.row = i,
layout.pos.col = j))
  grid.rect(gp = gpar(col = "red"))
  popViewport()
}
  }
  ## wipe out the outer edges 
  grid.rect(width = wide * 4, height = wide * 5, default.units = "mm",
gp = gpar(col = "white"))
  ## cover lines from centre squares
  grid.rect(width = wide * 2, height = wide * 3, default.units = "mm",
gp = gpar(col = "white", fill = "white"))

=X8--- cut here  

I've tried naming the viewports at the time they're pushed.  I had
hoped I could get back to them to do what I've done in 'Cube5' using
something like seekViewport(), but I see that won't work if the
viewport has been popped.

Reading through 'R Graphics' (the blue-toned book), I see there's such
a thing as vpList which I could make while I'm in that for loop,
but I don't see how I could do something like vpL <- c(vpL, cube.k)
sort of thing that I could do with a vector.  It's even less clear how
I could make use of vpStack or vpTree.

What is a smarter way to get back to those upper viewports?

-- 
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.   
   ___Patrick Connolly   
 {~._.~} Great minds discuss ideas
 _( Y )_Middle minds discuss events 
(:_~*~_:)Small minds discuss people  
 (_)-(_)   . Anon
  
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.

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

Re: [R] Axes in filled.contour plots

[Jim Lemon]

On Sat, 2008-08-16 at 14:28 -0700, hippie dream wrote:
> I am still struggling on how edit axes on a filled contour plot. I have
> managed to figure out how to place labels on the key of this graph but how
> to place the axes I want on this plot still eludes me. This command produces
> the plot I am looking for however as mentioned before these axes only go
> from 0 to 1.  
> 
> > filled.contour(contour, frame.plot=TRUE, color=terrain.colors)
> 
> If I try to set the xlim and ylim the plot itself no longer appears but
> simply the key. Moreover, when I set the limits and add the axes in this way
> the axes spill over into the key.
> 
> > axis(side=1, at=c(10, 30, 50, 70, 90, 110, 130, 150, 170, 190))
> 
> Is there any way to manually replace the given axis label with one of my
> choosing? That is, the 0 to 1 labels will correspond easily to the labels I
> would like to place on the graph. The tick marks are conveniently already in
> the right place. For example, can I make 0.6 appear as 100 instead? I don't
> really need to change the axis at all but rather just the labels.I have
> spent considerable time trying to figure this out but I am still having
> trouble. Sorry if this is ridiculously simple.
> 
Hi Sam,
Are you looking for something like axis.mult in the plotrix package?

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] continuous coloring of a polygon

[Roger Leenders]

Thanks Barry, that works beautifully!

Roger


Barry Rowlingson schreef:

2008/8/17 Roger Leenders <[EMAIL PROTECTED]>:
  

WinXP, R2.7.1

Thanks so much to all who have reponded to my inquery. The solutions are
most helpful.
There is only one final thing that I can't seem to get right.
All of the proposed solutions yield a figure that is not quite round, but is
wider than its height. This can be resolved by manually resizing the window
in which it is plotted, but that becomes cumbersome when many of these
figures need to be generated. Is there a way to enforce the figure to be
perfectly round?
I have tried plotting directly to a windows device (using for example the
code below) or to a pdf-device in which I specifically specify the square
size of the plot, but the situation remains in which the figure is wider
than it is high. Does anyone know how to solve this final piece of the
puzzle?




 add 'asp=1' to your plot() function to specify a 1:1 aspect ratio.
Should do the trick.

Barry



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

Re: [R] package building problem on windows

[Prof Brian Ripley]

Does your setting of TMPDIR, TEMP or TMP have a space in it?


From the R-admin manual:


  You may need to set TMPDIR to the absolute path to a suitable temporary
  directory: the default is 'c:/TEMP'.  (Use forward slashes and do not
  use a path including spaces.)

The documentation is that manual, not 'Murdoch's page', whatever that is.

On Sun, 17 Aug 2008, Christophe Dutang wrote:


Hi,

I'm trying to compile the package biglm, but when I build it with R CMD build 
biglm, it failed :


C:\LOCAL\c-dutang\code\R\biglm2>R CMD build biglm
* checking for file 'biglm/DESCRIPTION' ... OK
* preparing 'biglm':
* checking DESCRIPTION meta-information ...C:/DOCUME~1/c-dutang/Local: Can't 
op

n C:/DOCUME~1/c-dutang/Local: No such file or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
OK
* cleaning src
* removing junk files
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
Error: cannot open file 'biglm/DESCRIPTION' for reading

But when I install it, the installation is performed but I still have 'can't 
open' and 'no such file or directory'.


C:\LOCAL\c-dutang\code\R\biglm2>R CMD install biglm
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local: No such 
file

or directory
installing to ''


-- Making package biglm 
adding build stamp to DESCRIPTION
installing NAMESPACE file and metadata
making DLL ...
gcc  -std=gnu99  -IC:/PROGRA~1/R/R-2.7.1/include -O3 -Wall  -c 
boundedQR.c -

o boundedQR.o
gfortran   -O3  -c boundedQRf.f -o boundedQRf.o
windres --preprocessor="gcc -E -xc -DRC_INVOKED" -I 
C:/PROGRA~1/R/R-2.7.1/includ

e  -i biglm_res.rc -o biglm_res.o
gcc  -std=gnu99  -shared -s  -o biglm.dll biglm.def boundedQR.o boundedQRf.o 
big

lm_res.o  -LC:/PROGRA~1/R/R-2.7.1/bin   -lgfortran -lR
... DLL made
installing DLL
installing R files
installing man source files
installing indices
installing help

Building/Updating help pages for package 'biglm'

   Formats: text html latex example chm
adding MD5 sums

* DONE (biglm)

Actually, I have no directory local in 'document and settings' but a 
directory local settings exists. What's wrong with my installation? I install 
all tools as recommanded on Murdoch's page. Maybe R CMD build tries to use a 
temp file? in 'document and settings'?


I use R 2.7.1 with Windows xp pro.

Thanks in advance

Christophe

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


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

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

Re: [R] package building problem on windows

[Christophe Dutang]

Actually I was suspected that 'r cmd build' search for 'C:/DOCUME~1/c- 
dutang/Local' and not 'C:/DOCUME~1/c-dutang/Local Settings'.


I will check (tomorrow) if there are any spaces in environment  
variables TMP*.


Thanks for your suggestion.

Christophe


Le 17 août 08 à 14:08, Prof Brian Ripley a écrit :


Does your setting of TMPDIR, TEMP or TMP have a space in it?

From the R-admin manual:

 You may need to set TMPDIR to the absolute path to a suitable  
temporary

 directory: the default is 'c:/TEMP'.  (Use forward slashes and do not
 use a path including spaces.)

The documentation is that manual, not 'Murdoch's page', whatever  
that is.


On Sun, 17 Aug 2008, Christophe Dutang wrote:


Hi,

I'm trying to compile the package biglm, but when I build it with R  
CMD build biglm, it failed :


C:\LOCAL\c-dutang\code\R\biglm2>R CMD build biglm
* checking for file 'biglm/DESCRIPTION' ... OK
* preparing 'biglm':
* checking DESCRIPTION meta-information ...C:/DOCUME~1/c-dutang/ 
Local: Can't op

n C:/DOCUME~1/c-dutang/Local: No such file or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
OK
* cleaning src
* removing junk files
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
Error: cannot open file 'biglm/DESCRIPTION' for reading

But when I install it, the installation is performed but I still  
have 'can't open' and 'no such file or directory'.


C:\LOCAL\c-dutang\code\R\biglm2>R CMD install biglm
C:/DOCUME~1/c-dutang/Local: Can't open C:/DOCUME~1/c-dutang/Local:  
No such file

or directory
installing to ''


-- Making package biglm 
adding build stamp to DESCRIPTION
installing NAMESPACE file and metadata
making DLL ...
gcc  -std=gnu99  -IC:/PROGRA~1/R/R-2.7.1/include -O3 -Wall  -c  
boundedQR.c -

o boundedQR.o
gfortran   -O3  -c boundedQRf.f -o boundedQRf.o
windres --preprocessor="gcc -E -xc -DRC_INVOKED" -I C:/PROGRA~1/R/ 
R-2.7.1/includ

e  -i biglm_res.rc -o biglm_res.o
gcc  -std=gnu99  -shared -s  -o biglm.dll biglm.def boundedQR.o  
boundedQRf.o big

lm_res.o  -LC:/PROGRA~1/R/R-2.7.1/bin   -lgfortran -lR
... DLL made
installing DLL
installing R files
installing man source files
installing indices
installing help

Building/Updating help pages for package 'biglm'

  Formats: text html latex example chm
adding MD5 sums

* DONE (biglm)

Actually, I have no directory local in 'document and settings' but  
a directory local settings exists. What's wrong with my  
installation? I install all tools as recommanded on Murdoch's page.  
Maybe R CMD build tries to use a temp file? in 'document and  
settings'?


I use R 2.7.1 with Windows xp pro.

Thanks in advance

Christophe

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


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


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

[R] grangertest/lmtest ... what am I doing wrong ?

[tolga . i . uzuner]

Dear Achim, R Users,

What am I doing wrong in this example ?

a<-zoo(rnorm(100),order.by=1:100)
b<-lag(a)
regr<-na.exclude(merge(a,b))
plot(regr)
grangertest(regr[,1],regr[,2],3)

> a<-zoo(rnorm(100),order.by=1:100)
> b<-lag(a)
> regr<-na.exclude(merge(a,b))
> plot(regr)
> grangertest(regr[,1],regr[,2],3)
Error in solve(vc[ovar, ovar]) : subscript out of bounds
> 

This works fine in granger.test from the MSBVAR package:

> a<-zoo(rnorm(100),order.by=1:100)
> b<-lag(a)
> regr<-na.exclude(merge(a,b))
> plot(regr)
> granger.test(regr,3)
F-statistic   p-value
b -> a 1.259314e+34 0.000
a -> b 6.167865e-01 0.6059222
> 

Thanks in advance,
Tolga

Generally, this communication is for informational purposes only
and it is not intended as an offer or solicitation for the purchase
or sale of any financial instrument or as an official confirmation
of any transaction. In the event you are receiving the offering
materials attached below related to your interest in hedge funds or
private equity, this communication may be intended as an offer or
solicitation for the purchase or sale of such fund(s).  All market
prices, data and other information are not warranted as to
completeness or accuracy and are subject to change without notice.
Any comments or statements made herein do not necessarily reflect
those of JPMorgan Chase & Co., its subsidiaries and affiliates.

This transmission may contain information that is privileged,
confidential, legally privileged, and/or exempt from disclosure
under applicable law. If you are not the intended recipient, you
are hereby notified that any disclosure, copying, distribution, or
use of the information contained herein (including any reliance
thereon) is STRICTLY PROHIBITED. Although this transmission and any
attachments are believed to be free of any virus or other defect
that might affect any computer system into which it is received and
opened, it is the responsibility of the recipient to ensure that it
is virus free and no responsibility is accepted by JPMorgan Chase &
Co., its subsidiaries and affiliates, as applicable, for any loss
or damage arising in any way from its use. If you received this
transmission in error, please immediately contact the sender and
destroy the material in its entirety, whether in electronic or hard
copy format. Thank you.
Please refer to http://www.jpmorgan.com/pages/disclosures for
disclosures relating to UK legal entities.
[[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] grangertest/lmtest ... what am I doing wrong ?

[Gabor Grothendieck]

There is something wrong with the description below
since MSBVAR has granger.test, not grangertest,
and it takes 2 args, not 3.  This refers to 0.3-1 of
MSBVAR.

On Sun, Aug 17, 2008 at 9:07 AM,  <[EMAIL PROTECTED]> wrote:
> Dear Achim, R Users,
>
> What am I doing wrong in this example ?
>
> a<-zoo(rnorm(100),order.by=1:100)
> b<-lag(a)
> regr<-na.exclude(merge(a,b))
> plot(regr)
> grangertest(regr[,1],regr[,2],3)
>
>> a<-zoo(rnorm(100),order.by=1:100)
>> b<-lag(a)
>> regr<-na.exclude(merge(a,b))
>> plot(regr)
>> grangertest(regr[,1],regr[,2],3)
> Error in solve(vc[ovar, ovar]) : subscript out of bounds
>>
>
> This works fine in granger.test from the MSBVAR package:
>
>> a<-zoo(rnorm(100),order.by=1:100)
>> b<-lag(a)
>> regr<-na.exclude(merge(a,b))
>> plot(regr)
>> granger.test(regr,3)
>F-statistic   p-value
> b -> a 1.259314e+34 0.000
> a -> b 6.167865e-01 0.6059222
>>
>
> Thanks in advance,
> Tolga
>
> Generally, this communication is for informational purposes only
> and it is not intended as an offer or solicitation for the purchase
> or sale of any financial instrument or as an official confirmation
> of any transaction. In the event you are receiving the offering
> materials attached below related to your interest in hedge funds or
> private equity, this communication may be intended as an offer or
> solicitation for the purchase or sale of such fund(s).  All market
> prices, data and other information are not warranted as to
> completeness or accuracy and are subject to change without notice.
> Any comments or statements made herein do not necessarily reflect
> those of JPMorgan Chase & Co., its subsidiaries and affiliates.
>
> This transmission may contain information that is privileged,
> confidential, legally privileged, and/or exempt from disclosure
> under applicable law. If you are not the intended recipient, you
> are hereby notified that any disclosure, copying, distribution, or
> use of the information contained herein (including any reliance
> thereon) is STRICTLY PROHIBITED. Although this transmission and any
> attachments are believed to be free of any virus or other defect
> that might affect any computer system into which it is received and
> opened, it is the responsibility of the recipient to ensure that it
> is virus free and no responsibility is accepted by JPMorgan Chase &
> Co., its subsidiaries and affiliates, as applicable, for any loss
> or damage arising in any way from its use. If you received this
> transmission in error, please immediately contact the sender and
> destroy the material in its entirety, whether in electronic or hard
> copy format. Thank you.
> Please refer to http://www.jpmorgan.com/pages/disclosures for
> disclosures relating to UK legal entities.
>[[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] grangertest/lmtest ... what am I doing wrong ?

[tolga . i . uzuner]

Hi Gabor,

Perhaps, I should have been clearer. 

In the subject like, I put "grangertest/lmtest" to indicate the 
grangertest function in the lmtest package.

Later in the example, I am specifically talking about granger.test in the 
MSBVAR package (which is not failing, as in the example at the very 
bottom) to which you are referring, as a comparison to grangertest in the 
lmtest package (which is failing, in the example at the top). 

In the example at the very bottom of the email, I invoke granger.test from 
the MSBVAR package with two arguments, as you rightly put.  grangertest in 
the lmtest package is invoked with three arguments. This is the function 
which is failing.

Hope this clarifies my question.

Tolga




"Gabor Grothendieck" <[EMAIL PROTECTED]> 
17/08/2008 14:14

To
[EMAIL PROTECTED]
cc
"Achim Zeileis" <[EMAIL PROTECTED]>, R-HELP@r-project.org
Subject
Re: [R] grangertest/lmtest ... what am I doing wrong ?






There is something wrong with the description below
since MSBVAR has granger.test, not grangertest,
and it takes 2 args, not 3.  This refers to 0.3-1 of
MSBVAR.

On Sun, Aug 17, 2008 at 9:07 AM,  <[EMAIL PROTECTED]> wrote:
> Dear Achim, R Users,
>
> What am I doing wrong in this example ?
>
> a<-zoo(rnorm(100),order.by=1:100)
> b<-lag(a)
> regr<-na.exclude(merge(a,b))
> plot(regr)
> grangertest(regr[,1],regr[,2],3)
>
>> a<-zoo(rnorm(100),order.by=1:100)
>> b<-lag(a)
>> regr<-na.exclude(merge(a,b))
>> plot(regr)
>> grangertest(regr[,1],regr[,2],3)
> Error in solve(vc[ovar, ovar]) : subscript out of bounds
>>
>
> This works fine in granger.test from the MSBVAR package:
>
>> a<-zoo(rnorm(100),order.by=1:100)
>> b<-lag(a)
>> regr<-na.exclude(merge(a,b))
>> plot(regr)
>> granger.test(regr,3)
>F-statistic   p-value
> b -> a 1.259314e+34 0.000
> a -> b 6.167865e-01 0.6059222
>>
>
> Thanks in advance,
> Tolga
>
> Generally, this communication is for informational purposes only
> and it is not intended as an offer or solicitation for the purchase
> or sale of any financial instrument or as an official confirmation
> of any transaction. In the event you are receiving the offering
> materials attached below related to your interest in hedge funds or
> private equity, this communication may be intended as an offer or
> solicitation for the purchase or sale of such fund(s).  All market
> prices, data and other information are not warranted as to
> completeness or accuracy and are subject to change without notice.
> Any comments or statements made herein do not necessarily reflect
> those of JPMorgan Chase & Co., its subsidiaries and affiliates.
>
> This transmission may contain information that is privileged,
> confidential, legally privileged, and/or exempt from disclosure
> under applicable law. If you are not the intended recipient, you
> are hereby notified that any disclosure, copying, distribution, or
> use of the information contained herein (including any reliance
> thereon) is STRICTLY PROHIBITED. Although this transmission and any
> attachments are believed to be free of any virus or other defect
> that might affect any computer system into which it is received and
> opened, it is the responsibility of the recipient to ensure that it
> is virus free and no responsibility is accepted by JPMorgan Chase &
> Co., its subsidiaries and affiliates, as applicable, for any loss
> or damage arising in any way from its use. If you received this
> transmission in error, please immediately contact the sender and
> destroy the material in its entirety, whether in electronic or hard
> copy format. Thank you.
> Please refer to http://www.jpmorgan.com/pages/disclosures for
> disclosures relating to UK legal entities.
>[[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.
>



Generally, this communication is for informational purposes only
and it is not intended as an offer or solicitation for the purchase
or sale of any financial instrument or as an official confirmation
of any transaction. In the event you are receiving the offering
materials attached below related to your interest in hedge funds or
private equity, this communication may be intended as an offer or
solicitation for the purchase or sale of such fund(s).  All market
prices, data and other information are not warranted as to
completeness or accuracy and are subject to change without notice.
Any comments or statements made herein do not necessarily reflect
those of JPMorgan Chase & Co., its subsidiaries and affiliates.

This transmission may contain information that is privileged,
confidential, legally privileged, and/or exempt from disclosure
under applicable law. If you are not the intended recipient, you
are hereby notified that any disclosure, copying, 

[R] how to override/replace a function in a package namespace?

[Michael Friendly]

I'm trying to test an extension of mosaic() from the vcd package that 
requires a change to the
basic strucplot() function from that package.  I want to test my change 
by sourcing the

replacement function into my R session.

But when I do that,
source("c:/R/mosaics/strucplot-MF.R")
and run my extension, it is apparent that it is vcd:::strucplot that is 
eventually called, not
my replacement.  Thus, I'm looking for something with the semantics and 
effect of

source("c:/R/mosaics/strucplot-MF.R", namespace="vcd")
that would replace the package version with my testing one.

This question may have come up before, but searching the usual sources 
hasn't turned up anything

useful.

-Michael

--
Michael Friendly Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.

York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 override/replace a function in a package namespace?

[Gabor Grothendieck]

Check out:

?assignInNamespace


On Sun, Aug 17, 2008 at 12:42 PM, Michael Friendly <[EMAIL PROTECTED]> wrote:
> I'm trying to test an extension of mosaic() from the vcd package that
> requires a change to the
> basic strucplot() function from that package.  I want to test my change by
> sourcing the
> replacement function into my R session.
>
> But when I do that,
> source("c:/R/mosaics/strucplot-MF.R")
> and run my extension, it is apparent that it is vcd:::strucplot that is
> eventually called, not
> my replacement.  Thus, I'm looking for something with the semantics and
> effect of
> source("c:/R/mosaics/strucplot-MF.R", namespace="vcd")
> that would replace the package version with my testing one.
>
> This question may have come up before, but searching the usual sources
> hasn't turned up anything
> useful.
>
> -Michael
>
> --
> Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology
> Dept.
> York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
> 4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html
> Toronto, ONT  M3J 1P3 CANADA
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 override/replace a function in a package namespace?

[Ben Bolker]

Gabor Grothendieck  gmail.com> writes:

> 
> Check out:
> 
> ?assignInNamespace
> 
> On Sun, Aug 17, 2008 at 12:42 PM, Michael Friendly wrote:
> > I'm trying to test an extension of mosaic() from the vcd package that
> > requires a change to the
> > basic strucplot() function from that package.  I want to test my change by
> > sourcing the
> > replacement function into my R session.
> >

  Yes, but:

from the note in ?assignInNamespace

assignInNamespace and fixInNamespace change the copy in the namespace,
 but not any copies already exported from the namespace, in particular 
an object of that name in the package (if already attached) and any 
copies already imported into
other namespaces. They are really intended to be used only 
for objects which are
not exported from the namespace. They do attempt to alter 
a copy registered as
an S3 method if one is found. 

  So if strucplot is exported from the vcd namespace (which I guess
it is), this won't work.

  When I ran into a similar situation recently I couldn't find any
solution other than rebuilding the package with my changes.

  Ben Bolker

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

Re: [R] Axes in filled.contour plots

[hippie dream]

Jim,

Thanks so much for getting back to me. Axis.mult could definitely work so me
although there still appear to be a couple hiccups. I've included a .png
file of my plot rather than trying to explain it. I have for explanatory
purpose not called to plot with "axes=F". This illustrates that the
axis.mult call produces an axis with the proper labels but widens out the
axis for some reason. So the overlapping axes are here on purpose just to
show the problem.

Here is what I used:

filled.contour(contour, axes=FALSE, frame.plot=TRUE, color=terrain.colors)
axis.mult(side=1,mult=0.001,mult.label="Width (cm)")

So is there any way with axis.mult to produce the axes with the same scale
as the original? This is bringing me closer for sure. Thanks again.

Sam

Jim Lemon-2 wrote:
> 
> On Sat, 2008-08-16 at 14:28 -0700, hippie dream wrote:
>> I am still struggling on how edit axes on a filled contour plot. I have
>> managed to figure out how to place labels on the key of this graph but
>> how
>> to place the axes I want on this plot still eludes me. This command
>> produces
>> the plot I am looking for however as mentioned before these axes only go
>> from 0 to 1.  
>> 
>> > filled.contour(contour, frame.plot=TRUE, color=terrain.colors)
>> 
>> If I try to set the xlim and ylim the plot itself no longer appears but
>> simply the key. Moreover, when I set the limits and add the axes in this
>> way
>> the axes spill over into the key.
>> 
>> > axis(side=1, at=c(10, 30, 50, 70, 90, 110, 130, 150, 170, 190))
>> 
>> Is there any way to manually replace the given axis label with one of my
>> choosing? That is, the 0 to 1 labels will correspond easily to the labels
>> I
>> would like to place on the graph. The tick marks are conveniently already
>> in
>> the right place. For example, can I make 0.6 appear as 100 instead? I
>> don't
>> really need to change the axis at all but rather just the labels.I have
>> spent considerable time trying to figure this out but I am still having
>> trouble. Sorry if this is ridiculously simple.
>> 
> Hi Sam,
> Are you looking for something like axis.mult in the plotrix package?
> 
> 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.
> 
> 
http://www.nabble.com/file/p19021692/test.png test.png 
-- 
View this message in context: 
http://www.nabble.com/Axes-in-filled.contour-plots-tp18897760p19021692.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] how to override/replace a function in a package namespace?

[Michael Friendly]

Thanks, Gabor

For the record, what I was looking for is accomplished by
> source("c:/R/mosaics/strucplot-MF.R")
> assignInNamespace("strucplot",strucplot, ns="vcd")

-Michael

Gabor Grothendieck wrote:

Check out:

?assignInNamespace


On Sun, Aug 17, 2008 at 12:42 PM, Michael Friendly <[EMAIL PROTECTED]> wrote:
  

I'm trying to test an extension of mosaic() from the vcd package that
requires a change to the
basic strucplot() function from that package.  I want to test my change by
sourcing the
replacement function into my R session.

But when I do that,
source("c:/R/mosaics/strucplot-MF.R")
and run my extension, it is apparent that it is vcd:::strucplot that is
eventually called, not
my replacement.  Thus, I'm looking for something with the semantics and
effect of
source("c:/R/mosaics/strucplot-MF.R", namespace="vcd")
that would replace the package version with my testing one.

This question may have come up before, but searching the usual sources
hasn't turned up anything
useful.

-Michael

--
Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology
Dept.
York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA

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





--
Michael Friendly Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.

York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Streethttp://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 override/replace a function in a package namespace?

[Gabor Grothendieck]

In that case add this:

unlockBinding("strucplot", as.environment("package:vcd"))
assign("strucplot", my.strucplot, "package:vcd")



On Sun, Aug 17, 2008 at 1:22 PM, Ben Bolker <[EMAIL PROTECTED]> wrote:
> Gabor Grothendieck  gmail.com> writes:
>
>>
>> Check out:
>>
>> ?assignInNamespace
>>
>> On Sun, Aug 17, 2008 at 12:42 PM, Michael Friendly wrote:
>> > I'm trying to test an extension of mosaic() from the vcd package that
>> > requires a change to the
>> > basic strucplot() function from that package.  I want to test my change by
>> > sourcing the
>> > replacement function into my R session.
>> >
>
>  Yes, but:
>
> from the note in ?assignInNamespace
>
> assignInNamespace and fixInNamespace change the copy in the namespace,
>  but not any copies already exported from the namespace, in particular
> an object of that name in the package (if already attached) and any
> copies already imported into
> other namespaces. They are really intended to be used only
> for objects which are
> not exported from the namespace. They do attempt to alter
> a copy registered as
> an S3 method if one is found.
>
>  So if strucplot is exported from the vcd namespace (which I guess
> it is), this won't work.
>
>  When I ran into a similar situation recently I couldn't find any
> solution other than rebuilding the package with my changes.
>
>  Ben Bolker
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 override/replace a function in a package namespace?

[Henrik Bengtsson]

The few times I want to patch a function like this, I use:

  unlockBinding(name, env);
  assignInNamespace(name, value, ns=pkgName, envir=env);
  assign(name, value, envir=env);
  lockBinding(name, env);

/Henrik

On Sun, Aug 17, 2008 at 10:50 AM, Gabor Grothendieck
<[EMAIL PROTECTED]> wrote:
> In that case add this:
>
> unlockBinding("strucplot", as.environment("package:vcd"))
> assign("strucplot", my.strucplot, "package:vcd")
>
>
>
> On Sun, Aug 17, 2008 at 1:22 PM, Ben Bolker <[EMAIL PROTECTED]> wrote:
>> Gabor Grothendieck  gmail.com> writes:
>>
>>>
>>> Check out:
>>>
>>> ?assignInNamespace
>>>
>>> On Sun, Aug 17, 2008 at 12:42 PM, Michael Friendly wrote:
>>> > I'm trying to test an extension of mosaic() from the vcd package that
>>> > requires a change to the
>>> > basic strucplot() function from that package.  I want to test my change by
>>> > sourcing the
>>> > replacement function into my R session.
>>> >
>>
>>  Yes, but:
>>
>> from the note in ?assignInNamespace
>>
>> assignInNamespace and fixInNamespace change the copy in the namespace,
>>  but not any copies already exported from the namespace, in particular
>> an object of that name in the package (if already attached) and any
>> copies already imported into
>> other namespaces. They are really intended to be used only
>> for objects which are
>> not exported from the namespace. They do attempt to alter
>> a copy registered as
>> an S3 method if one is found.
>>
>>  So if strucplot is exported from the vcd namespace (which I guess
>> it is), this won't work.
>>
>>  When I ran into a similar situation recently I couldn't find any
>> solution other than rebuilding the package with my changes.
>>
>>  Ben Bolker
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] multiple t-tests

[B Nichols]

This is probably a basic data manipulation issue and an easy answer, but 
I can't seem to make it work. I am more used to working with data by 
participant, not by answer.


I have survey data, one answer per line (18000+ lines and counting), and 
each line includes variables for:


prepost (0 for pre, 1 for post)

question (the number of the question, from S1 to S85)

answer (a Likert range of 0-4)

I want to do t.tests (and Wilcox tests actually) on each question's pre 
and post means, but I can't figure out a single statement that will do 
this automatically, and I don't want to have to make 85 individual calls.


An individual call that seems to work with either t.test or wilcox.test 
looks like:


t.test(subset(answer, question=="S04") ~ subset(prepost, question=="S04"))

It would also be great to have a table including the difference between 
pre and posttest means. It seems these things should be doable, but I 
just can't find the right statement with the data in this format.


Thanks,

Bryan

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

[R] Error fitting overdispersed logistic regression: package dispmod

[Jason Curole]

Hi all,

First, a quick thank you for R; it's amazing.

I am trying to fit models for a count dataset following the overdispersed 
logisitic regression approach outlined in Baggerly et al. (BMC Bioinformatics, 
5:144; Annotated R code is given at the end of the paper) but R is returning an 
error with the data below.  Any help in understanding or overcoming this 
obstacle is appreciated.

>library(dispmod) # required for dispersion fitting

# Now the data

>counts<-matrix(c(2,1,3,1,2715597,3296062,2945864,2215143), ncol=2)
>temp<-c(25,20,25,20) # linear factor-Temperature
>trtmnt<-c(0,0,1,1) #categorical factor

#And the models

>fit1<-glm(counts~temp+trtmnt, family="binomial")
>fit2<-glm.binomial.disp(fit1)

Binomial overdispersed logit model fitting...
Iter.  1  phi: -3.615313e-07 
Error in glm(formula = counts ~ temp + trtmnt, family = "binomial", weights = 
disp.weights) : 
  negative weights not allowed

> disp.weights
[1]  54.865991  -5.218398 -15.379204   5.021180

So, clearly some dispersion weights are negative, which, according to my 
understanding of the model, would produce negative variances.  Is there a way 
around this?  

Thanks, Jason

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

Re: [R] Wichmann-Hill Random Number Generator and the Birthday Problem

[Shengqiao Li]

On Sun, 17 Aug 2008, Duncan Murdoch wrote:


Shengqiao Li wrote:

Dear all,

Recently I am generating large random samples (10M) and any duplicated 
numbers are not desired.
We tried several RNGs in R and found Wichmann-Hill did not produce 
duplications.


The duplication problem is the interesting birthday problem. If there are M 
possible numbers, randomly draw N numbers from them,

the average number of dupilcations D = N(N-1)/2/M.

For Knuth-TAOCP and Knuth-TAOCP-2002, M=2^30, since this modulus is used. D 
= 46566.12 for N=10M samples.


For Marsaglia-Multicarry, Super-Duper and Mersene-Twister, M=2^32. D = 
11641.53 for N = 10M samples.


My testing results (see below) agree with above analysis. But for 
Wichmann-Hill, it wasn't. Wichmann-Hill's cycle is 6.9536e12 (refer to RNG 
help by ?RNG and Whichmann's correction in 1984). Thus M <= 6.9536e12. D 
= 7.19052 when N=1e7 and D>= 179.763 for N=5e7. 

But in my tests, duplications were surprisingly not observed.

It seems that Wichmann-Hill has a much larger cycle than the one 
documented!


Anybody can solve this puzzle?



As I told you, look at the code. 
The cycle length of Wichmann-Hill in R appears to be 2.8e13, and that also 
appears to be M in your notation.  Your birthday problem calculation does not 
apply here.
You could probably get a good approximation to it by rounding the 
Wichmann-Hill output (e.g. look at round(2^30 * runif()), or maybe more 
severe rounding).


M is larger than what was previously documented, but Brian Ripley has 
corrected the documentation.


Wichmann claimed 2.78X10^13 in his 1982 original paper. He made correction 
in 1984, "The period of the generator was incorrectly given as 2.78 X 
10^13. In fact its period is only a quarter of that value: 6.95X10^12." 
In R version 2.7.1, Brian Ripley adopted Wichmann's 1984 correction in RNG 
document. That is the currently documented cycle is 6.95X10^12 not 
2.78X10^13 nor your approxmation 2.8e13.


So, is the cycle claimed in 1982 correct or the one claimed in 1984?

Even if the larger one 2.78e13 claimed in 1982 is correct,  around 45 
duplications were expected for 50M samples. But I got 0.  My testing 
method is based on the example code in the last third line in RNG help.


Shengqiao Li



Duncan Murdoch

Regards,

Shengqiao Li

Department of Statistics
West Virgina Unversity

==Testing===



RNGkind(kind="Knuth-TAOCP");
RNGkind();


[1] "Knuth-TAOCP" "Inversion"


sum(duplicated(runif(1e7)));


[1] 46216


RNGkind(kind="Knuth-TAOCP-2002");
RNGkind();


[1] "Knuth-TAOCP-2002" "Inversion"


sum(duplicated(runif(1e7)));


[1] 46373



RNGkind(kind="Marsaglia-Multicarry");
RNGkind();


[1] "Marsaglia-Multicarry" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11671


RNGkind(kind="Super-Duper");
RNGkind();


[1] "Super-Duper" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11704


RNGkind(kind="Mersenne-Twister");
RNGkind();


[1] "Mersenne-Twister" "Inversion"


sum(duplicated(runif(1e7)));


[1] 11508


RNGkind(kind="Wichmann-Hill");
RNGkind();


[1] "Wichmann-Hill" "Inversion"


sum(duplicated(runif(1e7)));


[1] 0



gc()


   used (Mb) gc trigger  (Mb)  max used   (Mb)
Ncells  199157  5.4 646480  17.3  10149018  271.1
Vcells 4497151 34.4  108714629 829.5 169585390 1293.9



sum(duplicated(runif(5e7)))


[1] 0




sessionInfo()


R version 2.7.1 (2008-06-23)
i386-pc-mingw32

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


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

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

==End of Testing===

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Making use of names of viewports (grid)

[Paul Murrell]

Hi


Patrick Connolly wrote:
> The following code, though not brilliant, works on an A4 page.  It
> might look odd on other devices of a very different size.
> 
> =X8--- cut here 
>   require(grid)
>   wide <- 15
>   vps <- grid.layout(nrow = 3, ncol = 4,
>  widths = unit(rep(1, 4), rep("null", 4)),
>  heights = unit(c(99, 1, 99),
>c("mm", "null", "mm")))
>   pushViewport(viewport(layout = vps))
>   for(k in 1:4){# label 4 viewports in top row
> cube.k <- ppaste("Cube", k)
> pushViewport(viewport(layout = vps, name = cube.k,
>   layout.pos.row = 1,
>   layout.pos.col = k))
> grid.rect(gp = gpar(lty = "dashed", lwd = .1))
> grid.text(cube.k, y = .9, gp = gpar(cex = .8))
> popViewport()
>   }
> ## label first viewport in bottom row
>   pushViewport(viewport(layout = vps, name = "Cube5",
> layout.pos.row = 3,
> layout.pos.col = 1))
>   grid.rect(gp = gpar(lty = "dashed", lwd = .1))
>   grid.text("Cube5", y = .95, gp = gpar(cex = .8))
> ## Set up a grid inside Cube5 viewport
>   vpc <- grid.layout(nrow = 5, ncol = 4,
>  widths = unit(rep(wide, 4), "mm"),
>  heights = unit(rep(wide, 5), "mm"))
>   pushViewport(viewport(layout = vpc))
>   for(i in 1:5){
> for(j in 1:4){
>   pushViewport(viewport(layout = vpc,  clip = "on",
> layout.pos.row = i,
> layout.pos.col = j))
>   grid.rect(gp = gpar(col = "red"))
>   popViewport()
> }
>   }
>   ## wipe out the outer edges 
>   grid.rect(width = wide * 4, height = wide * 5, default.units = "mm",
> gp = gpar(col = "white"))
>   ## cover lines from centre squares
>   grid.rect(width = wide * 2, height = wide * 3, default.units = "mm",
> gp = gpar(col = "white", fill = "white"))
> 
> =X8--- cut here  
> 
> I've tried naming the viewports at the time they're pushed.  I had
> hoped I could get back to them to do what I've done in 'Cube5' using
> something like seekViewport(), but I see that won't work if the
> viewport has been popped.
> 
> Reading through 'R Graphics' (the blue-toned book), I see there's such
> a thing as vpList which I could make while I'm in that for loop,
> but I don't see how I could do something like vpL <- c(vpL, cube.k)
> sort of thing that I could do with a vector.  It's even less clear how
> I could make use of vpStack or vpTree.
> 
> What is a smarter way to get back to those upper viewports?


Don't pop them ?  (use upViewport() instead)

Paul
-- 
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
[EMAIL PROTECTED]
http://www.stat.auckland.ac.nz/~paul/

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

[R] Optim stripping attributes from relistable objects

[rhelp . 20 . trevva]

Dear all,

The following code is inspired by the help file for the relist()
function (see?relist), which explicitly details how you can use a
relistable object in conjunction with optim to pass and reconstruct
complex parameter structures/groupings. The idea is that the optim()
function can only work with vectors, but in many cases you would like
to use a complex structure inside the objective function- relist is
one way to do that. The problem is that optim appears to be stripping
the attributes and therefore the example doesn't seem to run, giving
the error at the bottom.

> rb.banana <- function(params) {
+   #Params is initially a vector
+   cat("Params initially has the attributes:\n")
+   print(names(attributes(params)))
+   #Relisting it turns it into a list...
+   params <- relist(params)
+   cat("-\n")
+   #..which can then be called in the standard list manner
+   return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
+ }
>
> ipar<-  as.relistable(list(x=5,y=0))
> initial.params  <-  unlist(ipar)
>
> #Test to see if rb.banana works properly in the "normal" case
> rb.banana(initial.params)
Params initially has the attributes:
[1] "names""skeleton"
-
[1] 62516
>
> #OK, that's good. How about with optim though?
> optim(initial.params,rb.banana)
Params initially has the attributes:
[1] "names"
Error in relist(params) :
  The flesh argument does not contain a skeleton attribute.
Either ensure you unlist a relistable object, or specify the skeleton
separately.
>

What's going on here? Has this functionality been removed and the
documentation in relist() not updated? Or has the feature been broken?
Or have I misinterpreted something here (it wouldn't be the first
time!!)

Am running R version 2.7.1 (2008-06-23) under windows.

Cheers,

Mark

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

Re: [R] Axes in filled.contour plots [Solved]

[hippie dream]

I should let you guys know that I figured this whole axes on filled contour
plots issue out. This might not be the most elegant solution but it will
work me:

filled.contour(contour, axes=F, frame.plot=TRUE, color=terrain.colors, ylab=
"Length Along Flume (m)", key.title = title(main="Velocity\n(m/s)"),
key.axes = axis(4, seq(0, 0.6, by = 0.1)), plot.axes = {
axis.mult(side=1,mult=0.005,mult.label="Width (cm)")
axis(side=2, at=x, labels=colnames(contour)) })

Thanks for your help.

Sam

hippie dream wrote:
> 
> Jim,
> 
> Thanks so much for getting back to me. Axis.mult could definitely work so
> me although there still appear to be a couple hiccups. I've included a
> .png file of my plot rather than trying to explain it. I have for
> explanatory purpose not called to plot with "axes=F". This illustrates
> that the axis.mult call produces an axis with the proper labels but widens
> out the axis for some reason. So the overlapping axes are here on purpose
> just to show the problem.
> 
> Here is what I used:
> 
> filled.contour(contour, axes=FALSE, frame.plot=TRUE, color=terrain.colors)
> axis.mult(side=1,mult=0.001,mult.label="Width (cm)")
> 
> So is there any way with axis.mult to produce the axes with the same scale
> as the original? This is bringing me closer for sure. Thanks again.
> 
> Sam
> 
> Jim Lemon-2 wrote:
>> 
>> On Sat, 2008-08-16 at 14:28 -0700, hippie dream wrote:
>>> I am still struggling on how edit axes on a filled contour plot. I have
>>> managed to figure out how to place labels on the key of this graph but
>>> how
>>> to place the axes I want on this plot still eludes me. This command
>>> produces
>>> the plot I am looking for however as mentioned before these axes only go
>>> from 0 to 1.  
>>> 
>>> > filled.contour(contour, frame.plot=TRUE, color=terrain.colors)
>>> 
>>> If I try to set the xlim and ylim the plot itself no longer appears but
>>> simply the key. Moreover, when I set the limits and add the axes in this
>>> way
>>> the axes spill over into the key.
>>> 
>>> > axis(side=1, at=c(10, 30, 50, 70, 90, 110, 130, 150, 170, 190))
>>> 
>>> Is there any way to manually replace the given axis label with one of my
>>> choosing? That is, the 0 to 1 labels will correspond easily to the
>>> labels I
>>> would like to place on the graph. The tick marks are conveniently
>>> already in
>>> the right place. For example, can I make 0.6 appear as 100 instead? I
>>> don't
>>> really need to change the axis at all but rather just the labels.I have
>>> spent considerable time trying to figure this out but I am still having
>>> trouble. Sorry if this is ridiculously simple.
>>> 
>> Hi Sam,
>> Are you looking for something like axis.mult in the plotrix package?
>> 
>> 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.
>> 
>> 
>  http://www.nabble.com/file/p19021692/test.png test.png 
> 

-- 
View this message in context: 
http://www.nabble.com/Axes-in-filled.contour-plots-tp18897760p19023781.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] Optim stripping attributes from relistable objects

[Katharine Mullen]

>
> The following code is inspired by the help file for the relist()
> function (see?relist), which explicitly details how you can use a

The example is in the 'Details' section and, indeed, it looks like it no
longer works.

> relistable object in conjunction with optim to pass and reconstruct
> complex parameter structures/groupings. The idea is that the optim()
> function can only work with vectors, but in many cases you would like
> to use a complex structure inside the objective function- relist is
> one way to do that. The problem is that optim appears to be stripping
> the attributes and therefore the example doesn't seem to run, giving
> the error at the bottom.

You can get around this by specifying skeleton for relist:

rb.banana <- function(params) {
params <- relist(params, skeleton=list(x=NA,y=NA))
return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
}

ipar <-  as.relistable(list(x=5,y=0))

initial.params <- unlist(ipar)

xx <- optim(unlist(initial.params), rb.banana)

>
> > rb.banana <- function(params) {
> +   #Params is initially a vector
> +   cat("Params initially has the attributes:\n")
> +   print(names(attributes(params)))
> +   #Relisting it turns it into a list...
> +   params <- relist(params)
> +   cat("-\n")
> +   #..which can then be called in the standard list manner
> +   return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
> + }
> >
> > ipar<-  as.relistable(list(x=5,y=0))
> > initial.params  <-  unlist(ipar)
> >
> > #Test to see if rb.banana works properly in the "normal" case
> > rb.banana(initial.params)
> Params initially has the attributes:
> [1] "names""skeleton"
> -
> [1] 62516
> >
> > #OK, that's good. How about with optim though?
> > optim(initial.params,rb.banana)
> Params initially has the attributes:
> [1] "names"
> Error in relist(params) :
>   The flesh argument does not contain a skeleton attribute.
> Either ensure you unlist a relistable object, or specify the skeleton
> separately.
> >
>
> What's going on here? Has this functionality been removed and the
> documentation in relist() not updated? Or has the feature been broken?
> Or have I misinterpreted something here (it wouldn't be the first
> time!!)
>
> Am running R version 2.7.1 (2008-06-23) under windows.
>
> Cheers,
>
> Mark
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] grangertest/lmtest ... what am I doing wrong ?

[Achim Zeileis]

On Sun, 17 Aug 2008, [EMAIL PROTECTED] wrote:


Dear Achim, R Users,

What am I doing wrong in this example ?


The Granger test of x and y with lag k tests the significance of the x 
variables in

  y ~ Lags(y, 1:k) + Lags(x, 1:k)
If x is already a lag of y (as in your example), this is not really 
meaningful. In the "lmtest" implementation, it leads to collinearities. 
I'm not sure what MSBVAR does to resolve that.

Z


a<-zoo(rnorm(100),order.by=1:100)
b<-lag(a)
regr<-na.exclude(merge(a,b))
plot(regr)
grangertest(regr[,1],regr[,2],3)


a<-zoo(rnorm(100),order.by=1:100)
b<-lag(a)
regr<-na.exclude(merge(a,b))
plot(regr)
grangertest(regr[,1],regr[,2],3)

Error in solve(vc[ovar, ovar]) : subscript out of bounds




This works fine in granger.test from the MSBVAR package:


a<-zoo(rnorm(100),order.by=1:100)
b<-lag(a)
regr<-na.exclude(merge(a,b))
plot(regr)
granger.test(regr,3)

   F-statistic   p-value
b -> a 1.259314e+34 0.000
a -> b 6.167865e-01 0.6059222




Thanks in advance,
Tolga

Generally, this communication is for informational pur...{{dropped:30}}


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

[Nikolaos Lampadariou]

Hello everybody,

In am trying to analyse a BACI experiment and I really want to do it 
with R (which I find really exciting).  So, before moving on I though it 
would be a good idea to repeat some known experiments which are quite 
similar to my own. I tried to reproduce 2 published examples but without 
much success. The first one in particular is a published dataset 
analysed with SAS by McDonald et al. (2000). They also provide the 
original data as well as the SAS code. I don't know much about SAS and i 
really want to stick to R. So here follow 3 questions based on these 2 
papers.



Paper 1
(McDonald et al. 2000. Analysis of count data from before-after 
control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279)


Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples 
from 1984 are considered as before an impact and the remaining years as 
 after the impact. Each year, 96 transects were sampled (36 in the 
oiled area and 60 in the non-oiled area; "0" is for oiled and "1" for 
non-oiled). The authors compare 3 different ways of analysing the data 
including glmm. The data can be reproduced with the following commands 
(density is fake numbers but I can provide the original data since I've 
typed them in anyway):


>x<-c(rep("0",36),rep("1",60))
>oiled<-rep(x,6)
>year<-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), 
rep(1993,96), rep(1996,96))

>density<-runif(576, min=0, max=10)
>transect<-c(rep(1:96,6))
>oil<-data.frame("oiled"=oiled, "transect"=transect, "year"=year, 
"density"=density)


Question 1:
I can reproduce the results of the repeated measures anova with:
>oil.aov1<-aov(density~factor(year)*factor(oiled)+Error(factor(transect))

But why is the following command not working?
>oil.aov2<-aov(density~oiled*year + Error(oiled/transect), data=oil)

After reading the R-help archive, as well as Chambers and Hasties 
(Statistical Models in S) and Pinheiro's and Bates (Mixed effects models 
in S and S-plus) I would expect that the correct model is the oil.aov2. 
As you might see from the data, transect is nested within oiled and I 
still get the wrong results when I try +Error(transect) or 
+Error(oiled:transect) and many other variants.


Question 2 (on the same paper):
The authors conclude that it is better to analyse the data with a 
Generalized Linear (Mixed) Model Technique. I tried lme and after 
reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows:

>oil.lme<-lme(density~year*oiled, random=~1|oiled/transect)
or
>harvest.lmer<-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site))
but again no luck. I will get always some error messages or some 
interactions missing.



Paper 2
(Keough & Quinn, 2000. Legislative vs. practical protection of an 
intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881)


Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 
1997). the 1989-1991 are the before impact years and the other 4 years 
are the after the impact years. Eight sites were samples (2 protected 
and 6 impacted). In this dataset, site is nested within harvest (H) and 
year is nested within before-after (BA). Also, site is considered as 
random by the authors. The data (fake again) can be produced with the 
following commands:


>site<-c(rep(c("A1","A2", "RR1", "RR2", "WT1", "WT2", "WT3", "WT4"),8))
>H<-c(rep(c("exp", "exp", "prot", "pro", "exp", "exp", "exp", "exp"), 8))
>year<-c(rep(1989,8), rep(1990,8), rep(1991,8), rep(1993,8), 
rep(1994,8), rep(1995,8), rep(1996,8), rep(1997,8))

>BA<-c(rep("bef",24), rep("after",40))
>abund<-runif(64, min=0, max=10)
>harvest<-data.frame(abund, BA, H, site, year)

Question 3.
The authors use a complex anova design and here is part of their anova 
table which shows the design and the df.


source.of.var  df
Harvesting(H) 1, 6
before-after(BA)  1, 6
H x BA1, 6
Site(H)   6, 31
Year(BA)  6, 31
Site x BA 6, 31
Year x H  6, 31
Resid.31


My question again is: Why can't I reproduce the results? When I try a 
simple anova without any random factors:

>harvest.lm<-lm(abund~H*BA+H/site+BA/year+site:BA+year:H)
I get close enought but the results are different (at least I think they 
are different because the nomin. df are different).


So obviously I need to assign sites as a random factor somehow. So when 
I try to include site (which is nested in H) as a random factor and also 
year nested in BA (as the authors did) the best I can come up with is:

>harvest.lme<-lme(abund~H*BA+BA/year+BA/year:H, random=~1|H/site)
But I get a warning message (Warning message:In pf(q, df1, df2, 
lower.tail, log.p) : NaNs produced) and also I don't know where to put 
the site x BA term (whatever I try I get error messages).


I really apologise for the long post but after a week of studying and 
trying as many ideas and examples I could find and think of I felt that 
I really need advice from some more experienced users if I really want 
to

[R] trying to mimic page 70 of Software for Data Analysis

[markleeds]

I was trying to do what is on page 70 of John Chambers' new book namely 
using trace to invoke the browser by doing


trace(zapsmall, edit = TRUE)

but , typing above at an R prompt, I get

 trace(zapsmall, edit=TRUE)
sh: EMACS: command not found
Error in edit(name, file, title, editor) :
  problem with running editor EMACS

I looked in the archives but maybe it's a problem with the keywords I 
was using. If anyone could tell me what I need to do, it would be great. 
I'm thinking that I need to put something in my .Rprofile but I'm not 
sure what.  Thanks.


Oh, I use linux and my sessionInfo is below and I am using  GNU Emacs 
22.1.1. I usually use ESS and EMACS together but for this I figured it 
would be simpler to start from the R-prompt so that I was aligned with 
John's book as much as possible.



#

sessionInfo()
R version 2.7.0 (2008-04-22)
i686-redhat-linux-gnu

locale:

LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

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

other attached packages:
 [1] gsubfn_0.3-3   proto_0.3-8latticeExtra_0.4-1 
RColorBrewer_1.0-2 lattice_0.17-7 faraway_1.0.3  nnet_7.2-42 
car_1.2-8  reshape_0.8.0

[10] zoo_1.5-3  chron_2.3-22   MASS_7.2-42

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] i fixed my previous trace problem

[markleeds]

My previous problem with the emacs editor was that I had 
options(editor="EMACS",  when it should have been small letters as in


options(editor="emacs",papersize="letter", width=180, htmlhelp=FALSE, 
browser="/usr/bin/firefox")


Thanks for the help that I'm confident would have been there if I didn't 
figure it out myself.


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

Re: [R] trying to mimic page 70 of Software for Data Analysis

[Marc Schwartz]

on 08/17/2008 07:15 PM [EMAIL PROTECTED] wrote:
> I was trying to do what is on page 70 of John Chambers' new book namely
> using trace to invoke the browser by doing
> 
> trace(zapsmall, edit = TRUE)
> 
> but , typing above at an R prompt, I get
> 
>  trace(zapsmall, edit=TRUE)
> sh: EMACS: command not found
> Error in edit(name, file, title, editor) :
>   problem with running editor EMACS



Mark,

Take a look at:

  getOption("editor")

and see what it returns.  More than likely "EMACS" from what I see above.

It should be in all lower case. That is 'emacs'.

Take a look at your .Rprofile or perhaps at the 'EDITOR' or "VISUAL'
shell variables to see where you have this set and edit accordingly.
See ?edit for more information.

Enjoy the book. My copy came just before I left for useR! and I have not
had a chance to read it yet.

Regards,

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] Multiple Plotting help (lines don't always connect)

[stephen sefick]

d <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, NA, NA, NA, 14), .Dim = c(14L, 2L
), .Dimnames = list(NULL, c("a", "b")))
plot(d, type="b")

This is simplified, but Is there an option I am missing that will
force all of the points to be joined by a line?

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] Multiple Plotting help (lines don't always connect)

[jim holtman]

You have to get rid of the NAs since they indicate that there are no
values and will not draw a line:

> c.d <- d[complete.cases(d),]
> plot(c.d, type='b')
>


On Sun, Aug 17, 2008 at 6:55 PM, stephen sefick <[EMAIL PROTECTED]> wrote:
> d <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1,
> 2, 3, 4, 5, 6, 7, 8, 9, 10, NA, NA, NA, 14), .Dim = c(14L, 2L
> ), .Dimnames = list(NULL, c("a", "b")))
> plot(d, type="b")
>
> This is simplified, but Is there an option I am missing that will
> force all of the points to be joined by a line?
>
> 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.
>



-- 
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] Multiple Plotting help (lines don't always connect)

[markleeds]

I don't know if it's the best way but you can do

d<-d[complete.cases(d),]
plot(d, type="b")

plot doesn't connect the points when there are NAs between them.




On Sun, Aug 17, 2008 at  9:55 PM, stephen sefick wrote:


d <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, NA, NA, NA, 14), .Dim = c(14L, 2L
), .Dimnames = list(NULL, c("a", "b")))
plot(d, type="b")

This is simplified, but Is there an option I am missing that will
force all of the points to be joined by a line?

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

Re: [R] Multiple Plotting help (lines don't always connect)

[jim holtman]

You can also do:

plot(approx(d[,1], d[,2], xout=d[,1]), type='b')

This will interprete the values between the missing values, but it is
probably best to leave it as you first had it with the missing lines
so that you know there is something different in the data.

On Sun, Aug 17, 2008 at 6:55 PM, stephen sefick <[EMAIL PROTECTED]> wrote:
> d <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1,
> 2, 3, 4, 5, 6, 7, 8, 9, 10, NA, NA, NA, 14), .Dim = c(14L, 2L
> ), .Dimnames = list(NULL, c("a", "b")))
> plot(d, type="b")
>
> This is simplified, but Is there an option I am missing that will
> force all of the points to be joined by a line?
>
> 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.
>



-- 
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] Multiple Plotting help (lines don't always connect)

[Rolf Turner]

There was a thread on this recently.  Solutions were posted to allow you
to join interpolated points to ``real'' ones using a different line  
type.

See

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/127669.html

cheers,

Rolf Turner

On 18/08/2008, at 2:10 PM, jim holtman wrote:


You can also do:

plot(approx(d[,1], d[,2], xout=d[,1]), type='b')

This will interprete the values between the missing values, but it is
probably best to leave it as you first had it with the missing lines
so that you know there is something different in the data.

On Sun, Aug 17, 2008 at 6:55 PM, stephen sefick <[EMAIL PROTECTED]>  
wrote:

d <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, NA, NA, NA, 14), .Dim = c(14L, 2L
), .Dimnames = list(NULL, c("a", "b")))
plot(d, type="b")

This is simplified, but Is there an option I am missing that will
force all of the points to be joined by a line?

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.





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



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

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

[R] matrix row product and cumulative product

[Jeff Laake]

I spent a lot of time searching and came up empty handed on the 
following query. Is there an equivalent to rowSums that does product or 
cumulative product and avoids use of apply or looping? I found a rowProd 
in a package but it was a convenience function for apply. As part of a 
likelihood calculation called from optim, I’m computing products and 
cumulative products of rows of matrices with far more rows than columns. 
I started with apply and after some thought realized that a loop of 
columns might be faster and it was substantially faster (see below). 
Because the likelihood function is called many times I’d like to speed 
it up even more if possible.


Below is an example showing the cumprod.matrix and prod.matrix looping 
functions that I wrote and some timing comparisons to the use of apply 
for different column and row dimensions. At this point I’m better off 
with looping but I’d like to hear of any further suggestions.


Thanks –jeff

> prod.matrix=function(x)
+ {
+ y=x[,1]
+ for(i in 2:dim(x)[2])
+ y=y*x[,i]
+ return(y)
+ }

> cumprod.matrix=function(x)
+ {
+ y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
+ y[,1]=x[,1]
+ for (i in 2:dim(x)[2])
+ y[,i]=y[,i-1]*x[,i]
+ return(y)
+ }

> N=1000
> xmat=matrix(runif(N),ncol=10)
> system.time(cumprod.matrix(xmat))
user system elapsed
1.07 0.09 1.15
> system.time(t(apply(xmat,1,cumprod)))
user system elapsed
29.27 0.21 29.50
> system.time(prod.matrix(xmat))
user system elapsed
0.29 0.00 0.30
> system.time(apply(xmat,1,prod))
user system elapsed
30.69 0.00 30.72
> xmat=matrix(runif(N),ncol=100)
> system.time(cumprod.matrix(xmat))
user system elapsed
1.05 0.13 1.18
> system.time(t(apply(xmat,1,cumprod)))
user system elapsed
3.55 0.14 3.70
> system.time(prod.matrix(xmat))
user system elapsed
0.38 0.01 0.39
> system.time(apply(xmat,1,prod))
user system elapsed
2.87 0.00 2.89
> xmat=matrix(runif(N),ncol=1000)
> system.time(cumprod.matrix(xmat))
user system elapsed
1.30 0.18 1.46
> system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.77 0.27 2.05
> system.time(prod.matrix(xmat))
user system elapsed
0.46 0.00 0.47
> system.time(apply(xmat,1,prod))
user system elapsed
0.7 0.0 0.7
> xmat=matrix(runif(N),ncol=1)
> system.time(cumprod.matrix(xmat))
user system elapsed
1.28 0.00 1.29
> system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.19 0.08 1.26
> system.time(prod.matrix(xmat))
user system elapsed
0.40 0.00 0.41
> system.time(apply(xmat,1,prod))
user system elapsed
0.57 0.00 0.56
> xmat=matrix(runif(N),ncol=10)
> system.time(cumprod.matrix(xmat))
user system elapsed
3.18 0.00 3.19
> system.time(t(apply(xmat,1,cumprod)))
user system elapsed
1.42 0.21 1.64
> system.time(prod.matrix(xmat))
user system elapsed
1.05 0.00 1.05
> system.time(apply(xmat,1,prod))
user system elapsed
0.82 0.00 0.81
> R.Version()
$platform
[1] "i386-pc-mingw32"
.
.
.
$version.string
[1] "R version 2.7.1 (2008-06-23)"

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

Re: [R] matrix row product and cumulative product

[Charles C. Berry]

On Sun, 17 Aug 2008, Jeff Laake wrote:

I spent a lot of time searching and came up empty handed on the following 
query. Is there an equivalent to rowSums that does product or cumulative 
product and avoids use of apply or looping? I found a rowProd in a package 
but it was a convenience function for apply. As part of a likelihood 
calculation called from optim, I’m computing products and cumulative products 
of rows of matrices with far more rows than columns. I started with apply and 
after some thought realized that a loop of columns might be faster and it was 
substantially faster (see below). Because the likelihood function is called 
many times I’d like to speed it up even more if possible.



You might check out the 'inline' or 'jit' packages.

Otherwise, if you can as easily treat xmat as a list (or 
data.frame),


Reduce( "*", xmat.data.frame, accumulate=want.cumprod )

(where want.cumprod is FALSE for product, TRUE for cumulative product) 
will be a bit faster in many circumstances. However, this advantage is 
lost if you must retain xmat as a matrix since converting it to a 
data.frame seems to require more time than you save.


HTH,

Chuck



Below is an example showing the cumprod.matrix and prod.matrix looping 
functions that I wrote and some timing comparisons to the use of apply for 
different column and row dimensions. At this point I’m better off with 
looping but I’d like to hear of any further suggestions.


Thanks –jeff


 prod.matrix=function(x)

+ {
+ y=x[,1]
+ for(i in 2:dim(x)[2])
+ y=y*x[,i]
+ return(y)
+ }


 cumprod.matrix=function(x)

+ {
+ y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
+ y[,1]=x[,1]
+ for (i in 2:dim(x)[2])
+ y[,i]=y[,i-1]*x[,i]
+ return(y)
+ }


 N=1000
 xmat=matrix(runif(N),ncol=10)
 system.time(cumprod.matrix(xmat))

user system elapsed
1.07 0.09 1.15

 system.time(t(apply(xmat,1,cumprod)))

user system elapsed
29.27 0.21 29.50

 system.time(prod.matrix(xmat))

user system elapsed
0.29 0.00 0.30

 system.time(apply(xmat,1,prod))

user system elapsed
30.69 0.00 30.72

 xmat=matrix(runif(N),ncol=100)
 system.time(cumprod.matrix(xmat))

user system elapsed
1.05 0.13 1.18

 system.time(t(apply(xmat,1,cumprod)))

user system elapsed
3.55 0.14 3.70

 system.time(prod.matrix(xmat))

user system elapsed
0.38 0.01 0.39

 system.time(apply(xmat,1,prod))

user system elapsed
2.87 0.00 2.89

 xmat=matrix(runif(N),ncol=1000)
 system.time(cumprod.matrix(xmat))

user system elapsed
1.30 0.18 1.46

 system.time(t(apply(xmat,1,cumprod)))

user system elapsed
1.77 0.27 2.05

 system.time(prod.matrix(xmat))

user system elapsed
0.46 0.00 0.47

 system.time(apply(xmat,1,prod))

user system elapsed
0.7 0.0 0.7

 xmat=matrix(runif(N),ncol=1)
 system.time(cumprod.matrix(xmat))

user system elapsed
1.28 0.00 1.29

 system.time(t(apply(xmat,1,cumprod)))

user system elapsed
1.19 0.08 1.26

 system.time(prod.matrix(xmat))

user system elapsed
0.40 0.00 0.41

 system.time(apply(xmat,1,prod))

user system elapsed
0.57 0.00 0.56

 xmat=matrix(runif(N),ncol=10)
 system.time(cumprod.matrix(xmat))

user system elapsed
3.18 0.00 3.19

 system.time(t(apply(xmat,1,cumprod)))

user system elapsed
1.42 0.21 1.64

 system.time(prod.matrix(xmat))

user system elapsed
1.05 0.00 1.05

 system.time(apply(xmat,1,prod))

user system elapsed
0.82 0.00 0.81

 R.Version()

$platform
[1] "i386-pc-mingw32"
.
.
.
$version.string
[1] "R version 2.7.1 (2008-06-23)"

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



Charles C. Berry(858) 534-2098
Dept of Family/Preventive Medicine
E mailto:[EMAIL PROTECTED]  UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

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

[Applejus]

Hello,

I am trying to convert the following command from SPLUS to R:

graphsheet(pages = TRUE)

Does anyone have an idea what is the equivalent in R?

Thanks


-- 
View this message in context: 
http://www.nabble.com/graphsheet-tp19026010p19026010.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] Fucntion scope question. General non-linear solution help.

[rkevinburton]

I would like to solve the equation is is the sum from k = i to N of

choose(N,k) * MR ^ k * (1 - MR) ^ (N - k) - 0.50 = 0

I want to solve for MR. This seems like a non-linear equation to me. But I am 
having a hard time writing the function that implements the above. I could use 
'for(...) as a brute force appoarch but I would like a more "elegant" solution. 
The variables 'N' and 'i' are basically constant so the function has to take 
these from some kind of global space. So if I take t brute force apporach I 
came up with:

f <- function(MR)
{
k <- i:N
return sum(choose(N,k) * MR ^ k * (1 - MR) ^ (N - k)) - 0.5
}

Does this seem like a reasonable implemetation? How are 'N' and 'i' declare as 
"global"? For each equation N and I are constant but I want to be able to 
modify them. In other words solve the equantion after setting N to 6 and i to 5 
then again after setting i to 4.

The next question is regarding which 'R' function would be best suited to 
solving this equation? I looked at 'nls' but that seems to take data as an 
input. I want to solve the equation. What other options do I have? There must 
be an 'R' function to solve a non-linear equation. I did 
help.search("non-linear") and the closest match was nlm. But nlm minimizes the 
function rather than solving it.

Ideas?

Thank you.

Kevin

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

Re: [R] graphsheet

[Duncan Mackay]

At 14:53 18/08/2008, you wrote:


Hello,

I am trying to convert the following command from SPLUS to R:

graphsheet(pages = TRUE)

Does anyone have an idea what is the equivalent in R?

Thanks


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


Your intention in Splus is not clear; I assume you want to have tabbed 
window pages which is not possible in R.
I assume that you are using a MS OS and you did not specify your R version. 
There were changes for R 2.7

for windowing and using the graphics history.

See Prof Brian Ripley's email details:
Date: Mon, 21 Apr 2008 13:16:39 +0100 (BST)
From: Prof Brian Ripley <[EMAIL PROTECTED]>
To: Duncan Murdoch <[EMAIL PROTECTED]>
Subject: Re: [R] graphics history

Sorry I have not got a url for it.
also
see ?win.graph and ?Devices (postscript/pdf) if you want to combine plots

Regards

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
ARMIDALE NSW 2351
Email home: [EMAIL PROTECTED]

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

Re: [R] Optim stripping attributes from relistable objects

[Prof Brian Ripley]

On Mon, 18 Aug 2008, Katharine Mullen wrote:



The following code is inspired by the help file for the relist()
function (see?relist), which explicitly details how you can use a


The example is in the 'Details' section and, indeed, it looks like it no
longer works.


I don't think it ever did (after all there is a comment that dnorm does 
not have the parameters stated).


The subject line here illustrates the misunderstanding.  Optim *strips* 
nothing -- it creates a new object to pass to the function.  The help says


 Any names given to 'par' will be copied to the vectors passed to
 'fn' and 'gr'.

Note the use of 'copy' here.  You cannot assume that any other attributes 
will be copied, and this behaviour has not been changed (except to add 
names in March 2006, well before relist() was added), AFAICS.



relistable object in conjunction with optim to pass and reconstruct
complex parameter structures/groupings. The idea is that the optim()
function can only work with vectors, but in many cases you would like
to use a complex structure inside the objective function- relist is
one way to do that. The problem is that optim appears to be stripping
the attributes and therefore the example doesn't seem to run, giving
the error at the bottom.


You can get around this by specifying skeleton for relist:

rb.banana <- function(params) {
   params <- relist(params, skeleton=list(x=NA,y=NA))
   return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
}

ipar <-  as.relistable(list(x=5,y=0))

initial.params <- unlist(ipar)

xx <- optim(unlist(initial.params), rb.banana)




rb.banana <- function(params) {

+   #Params is initially a vector
+   cat("Params initially has the attributes:\n")
+   print(names(attributes(params)))
+   #Relisting it turns it into a list...
+   params <- relist(params)
+   cat("-\n")
+   #..which can then be called in the standard list manner
+   return( (1-params$x)^2 + 100*(params$y - params$x^2)^2)
+ }


ipar<-  as.relistable(list(x=5,y=0))
initial.params  <-  unlist(ipar)

#Test to see if rb.banana works properly in the "normal" case
rb.banana(initial.params)

Params initially has the attributes:
[1] "names""skeleton"
-
[1] 62516


#OK, that's good. How about with optim though?
optim(initial.params,rb.banana)

Params initially has the attributes:
[1] "names"
Error in relist(params) :
  The flesh argument does not contain a skeleton attribute.
Either ensure you unlist a relistable object, or specify the skeleton
separately.




What's going on here? Has this functionality been removed and the
documentation in relist() not updated? Or has the feature been broken?
Or have I misinterpreted something here (it wouldn't be the first
time!!)

Am running R version 2.7.1 (2008-06-23) under windows.

Cheers,

Mark

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



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



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

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

Re: [R] graphsheet

[Prof Brian Ripley]

You were probably referring to

https://stat.ethz.ch/pipermail/r-help/2008-April/160078.html

but a reference to ?windows would have sufficed.

On Mon, 18 Aug 2008, Duncan Mackay wrote:



At 14:53 18/08/2008, you wrote:


Hello,

I am trying to convert the following command from SPLUS to R:

graphsheet(pages = TRUE)

Does anyone have an idea what is the equivalent in R?

Thanks


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


Your intention in Splus is not clear; I assume you want to have tabbed window 
pages which is not possible in R.
I assume that you are using a MS OS and you did not specify your R version. 
There were changes for R 2.7

for windowing and using the graphics history.

See Prof Brian Ripley's email details:
Date: Mon, 21 Apr 2008 13:16:39 +0100 (BST)
From: Prof Brian Ripley <[EMAIL PROTECTED]>
To: Duncan Murdoch <[EMAIL PROTECTED]>
Subject: Re: [R] graphics history

Sorry I have not got a url for it.
also
see ?win.graph and ?Devices (postscript/pdf) if you want to combine plots

Regards

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
ARMIDALE NSW 2351
Email home: [EMAIL PROTECTED]

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



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

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

Re: [R] exporting adaBoost model

[Hans W. Borchers]

One way to port these kinds of models between applications is the 
Predictive Model Markup Language (PMML). The R package 'PMML' supports 
linear regression, rpart, SVM, and others, not adaBoost. On the other
side, not even the Python machine learning library Orange does have
an import function for PMML. 

Perhaps a more attractive possibility is to call R functions from Python
through the 'Rpy'  Python interface to R.
You could send data from Python to R for being handled by adaBoost and
get back the results. Of course, R and Python need be installed on your
machine and you cannot generate a single executable.,

//  Hans Werner Borchers


Bob Flagg wrote:
> 
> Dear all,
> 
> I'm using adaBoost from the ada package to build a classification
> model.  After training the model in R I'd like to use it in a
> Python application.  Is it possible to export the model in some
> way to make translating into python easier?  
> 
> Any help would be greatly appreciated. Thanks.
> Bob
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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/exporting-adaBoost-model-tp19005415p19026563.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] Fucntion scope question. General non-linear solution help.

[Prof Brian Ripley]

Looks like you want to solve pbinom(k-1, N, p) = 0.5 for p.  That is easy: 
use uniroot.


testit <- function(k, N)
{
fn <- function(p, k, N) pbinom(k-1, N, p) - 0.5
uniroot(fn, c(0,1), k=k, N=N)
}
testit(6, 10)

On Sun, 17 Aug 2008, [EMAIL PROTECTED] wrote:


I would like to solve the equation is is the sum from k = i to N of

choose(N,k) * MR ^ k * (1 - MR) ^ (N - k) - 0.50 = 0


That's not what you have below: I presume that you want the sum to be 0.5.
The sum is 1 - pbinom(k-1, N, MR).


I want to solve for MR. This seems like a non-linear equation to me. But I am having a 
hard time writing the function that implements the above. I could use 'for(...) as a 
brute force appoarch but I would like a more "elegant" solution. The variables 
'N' and 'i' are basically constant so the function has to take these from some kind of 
global space. So if I take t brute force apporach I came up with:

f <- function(MR)
{
   k <- i:N
   return sum(choose(N,k) * MR ^ k * (1 - MR) ^ (N - k)) - 0.5
}

Does this seem like a reasonable implemetation? How are 'N' and 'i' declare as 
"global"? For each equation N and I are constant but I want to be able to 
modify them. In other words solve the equantion after setting N to 6 and i to 5 then 
again after setting i to 4.

The next question is regarding which 'R' function would be best suited to solving this 
equation? I looked at 'nls' but that seems to take data as an input. I want to solve the 
equation. What other options do I have? There must be an 'R' function to solve a 
non-linear equation. I did help.search("non-linear") and the closest match was 
nlm. But nlm minimizes the function rather than solving it.

Ideas?

Thank you.

Kevin

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



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

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

Re: [R] matrix row product and cumulative product

[Moshe Olshansky]

Hi Jeff,

If I understand correctly, the overhead of a loop is that at each iteration the 
command must be interpreted, and this time is independent of the number of rows 
N. So if N is small this overhead may be very significant but when N is large 
this should be very small compared to the time needed to multiply N pairs of 
numbers. 


--- On Mon, 18/8/08, Jeff Laake <[EMAIL PROTECTED]> wrote:

> From: Jeff Laake <[EMAIL PROTECTED]>
> Subject: [R] matrix row product and cumulative product
> To: r-help@r-project.org
> Received: Monday, 18 August, 2008, 12:49 PM
> I spent a lot of time searching and came up empty handed on
> the 
> following query. Is there an equivalent to rowSums that
> does product or 
> cumulative product and avoids use of apply or looping? I
> found a rowProd 
> in a package but it was a convenience function for apply.
> As part of a 
> likelihood calculation called from optim, I’m computing
> products and 
> cumulative products of rows of matrices with far more rows
> than columns. 
> I started with apply and after some thought realized that a
> loop of 
> columns might be faster and it was substantially faster
> (see below). 
> Because the likelihood function is called many times I’d
> like to speed 
> it up even more if possible.
> 
> Below is an example showing the cumprod.matrix and
> prod.matrix looping 
> functions that I wrote and some timing comparisons to the
> use of apply 
> for different column and row dimensions. At this point
> I’m better off 
> with looping but I’d like to hear of any further
> suggestions.
> 
> Thanks –jeff
> 
>  > prod.matrix=function(x)
> + {
> + y=x[,1]
> + for(i in 2:dim(x)[2])
> + y=y*x[,i]
> + return(y)
> + }
> 
>  > cumprod.matrix=function(x)
> + {
> + y=matrix(1,nrow=dim(x)[1],ncol=dim(x)[2])
> + y[,1]=x[,1]
> + for (i in 2:dim(x)[2])
> + y[,i]=y[,i-1]*x[,i]
> + return(y)
> + }
> 
>  > N=1000
>  > xmat=matrix(runif(N),ncol=10)
>  > system.time(cumprod.matrix(xmat))
> user system elapsed
> 1.07 0.09 1.15
>  > system.time(t(apply(xmat,1,cumprod)))
> user system elapsed
> 29.27 0.21 29.50
>  > system.time(prod.matrix(xmat))
> user system elapsed
> 0.29 0.00 0.30
>  > system.time(apply(xmat,1,prod))
> user system elapsed
> 30.69 0.00 30.72
>  > xmat=matrix(runif(N),ncol=100)
>  > system.time(cumprod.matrix(xmat))
> user system elapsed
> 1.05 0.13 1.18
>  > system.time(t(apply(xmat,1,cumprod)))
> user system elapsed
> 3.55 0.14 3.70
>  > system.time(prod.matrix(xmat))
> user system elapsed
> 0.38 0.01 0.39
>  > system.time(apply(xmat,1,prod))
> user system elapsed
> 2.87 0.00 2.89
>  > xmat=matrix(runif(N),ncol=1000)
>  > system.time(cumprod.matrix(xmat))
> user system elapsed
> 1.30 0.18 1.46
>  > system.time(t(apply(xmat,1,cumprod)))
> user system elapsed
> 1.77 0.27 2.05
>  > system.time(prod.matrix(xmat))
> user system elapsed
> 0.46 0.00 0.47
>  > system.time(apply(xmat,1,prod))
> user system elapsed
> 0.7 0.0 0.7
>  > xmat=matrix(runif(N),ncol=1)
>  > system.time(cumprod.matrix(xmat))
> user system elapsed
> 1.28 0.00 1.29
>  > system.time(t(apply(xmat,1,cumprod)))
> user system elapsed
> 1.19 0.08 1.26
>  > system.time(prod.matrix(xmat))
> user system elapsed
> 0.40 0.00 0.41
>  > system.time(apply(xmat,1,prod))
> user system elapsed
> 0.57 0.00 0.56
>  > xmat=matrix(runif(N),ncol=10)
>  > system.time(cumprod.matrix(xmat))
> user system elapsed
> 3.18 0.00 3.19
>  > system.time(t(apply(xmat,1,cumprod)))
> user system elapsed
> 1.42 0.21 1.64
>  > system.time(prod.matrix(xmat))
> user system elapsed
> 1.05 0.00 1.05
>  > system.time(apply(xmat,1,prod))
> user system elapsed
> 0.82 0.00 0.81
>  > R.Version()
> $platform
> [1] "i386-pc-mingw32"
> .
> .
> .
> $version.string
> [1] "R version 2.7.1 (2008-06-23)"
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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.