[R] nls() and numerical integration (e.g. integrate()) working together?

2007-08-22 Thread Michael Lindner
Dear List-Members,

since 3 weeks I have been heavily working on reproducing the results of an
economic paper. The method there uses the numerical solution of an integral
within nonlinear least squares. Within the integrand there is also some
parameter to estimate. Is that in the end possible to implement in R
[Originally it was done in GAUSS]? I'm nearly into giving up.

I constucted an example to showing the problems I face.

I have three questions - related to three errors shown below: 
1) How to make clear that in the integrand z is the integration variable and
b1 is a parameter while x1 is a data variable
2) and 3) How to set up a correct estimation of the integral?

library(stats)
y <- c(2,15,24,21,5,6,)
x1 <- c(2.21,5,3,5,2,1)
x2 <- c(4.51,6,2,11,0.4,3)
f <- function(z) {z + b1*x1}
vf <- Vectorize(f)
g <- function(z) {z + x1}
vg <- Vectorize(f)

Error 1: 
> nls(y ~ integrate(vf,0,1)+b2*x2,start=list(b1=0.5,b2=2))
Error in function (z)  : object "b1" not found

Error 2:
> nls(y ~ integrate(vg,0,1)+b2*x2,start=list(b1=0.5,b2=2))
Error in integrate(vg, 0, 1) : REAL() can only be applied to a 'numeric',
not a 'list'

Error 3:
> nls(y ~ integrate(g,0,1)+b2*x2,start=list(b1=0.5,b2=2))
Error in integrate(g, 0, 1) + b2 * x2 : non-numeric argument to binary
operator
In addition: Warning messages:
1: longer object length
is not a multiple of shorter object length in: z + x1

With a lot of thanks in advance,

Michael

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


[R] two labels on x-axis (year and month)

2007-08-22 Thread mgilgen
hej

i'm plotting time-series and label the x-axis as follows:

r <- as.POSIXct(round(range(p1$time), "month"))

to define the time range for labeling the xaxis

plot(p1$time,p1$ p1, type="l", xaxt="n")

plots p1 against time

axis.POSIXct(1, at=seq(r[1], r[2], by="month"), format="%m")

labels the axis in mothly steps.

what I want do do now, is a second label for the x-axis, that stands lower and 
indicates the years.

like 05 06 07 08 09 10 11 12 01 02 03 04 05 etc...
 2004  2005
I don't know how to proceed with all the possibilities in axis, par and plot

thanks for your help

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


Re: [R] gWidgets (tcltk): problem extracting values f rom widgets in glayout grid

2007-08-22 Thread j verzani
Dear Evan, 
 epamail.epa.gov> writes:
> 
> 
> Hello,
> 
> I haven't been able to find an example for the second case below -- or
> perhaps I didn't recognize it when I saw it.
> Is there a value for x such that svalue(x) will return "bbb", either by
> itself or as part of an array? Or do I need to do something else
> entirely?
> (R2.5.1; Windows XP)
> 

You will need to store the widget. See below.

> > ### Case 2: this makes an identical window & widget, with "b"
> replacing "a"
> > ### but does't return the value
> > wtesta = gwindow("wtestb",visible=TRUE)
> > wtestb = glayout(visible=TRUE,container=wtesta)
> > wtestb[1,1]=gedit("bbb",container=wtestb)
> > print(svalue(wtestb[1,1]))


the glayout container doesn't work that way. You don't get a [
method, just a [<- one. Try something like this instead to store the widget:

wtestb[1,1]=(bbb <- gedit("bbb",container=wtestb) )
print(svalue(bbb))

--John


> Evan Englund
> U.S. EPA
> 702-798-2248
> 
> __
> R-help  stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] degrees of freedom question

2007-08-22 Thread Greg Tarpinian
R2.3, WinXP


Dear all,

I am using the following functions:

 f1 = Phi1+(Phi2-Phi1)/(1+exp((log(Phi3)-log(x))/exp(log(Phi4)))
 f2 = Phi1+(Phi2-Phi1)/(1+exp((log(Phi3)-log(r)-log(x))/exp(log(Phi4)))

subject to the residual weighting

 Var(e[i]) = sigma^2 * abs( E(y) )^(2*Delta)


Here is my question, in steps:

 1.  Function f1 is separately fitted to two different datasets
 corresponding to two different dose response curves.  These
 fits are unweighted.
 2.  Function f2 is fitted to the pooled data such that the two
 dose response curves are assumed to differ _only_ in log(r).
 This fit is also unweighted.
 3.  The residuals from #2 are used to estimate an appropriate
 sigma^2 and Delta to use in weighting.
 4.  The functions described in #1 and #2 are refitted, but this
 time weighted using the information gathered in #3.
 5.  How many degrees of freedom should be allocated to the 
 weighted residual sums of squares?  (There are three such
 SSE's, one for each individual model, and one for the overall
 joint model)


Much thanks in advance,

   Greg

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


Re: [R] Need a variant of rbind for datasets with different numbers of columns

2007-08-22 Thread jim holtman
Where is the data coming from since it has a variable number of
columns in each row?  Is it coming from a text file?  If so, you can
use the "fill=TRUE" option when reading to fill out empty columns.
You need to provide at least a subset of the data so we can see what
you are working with.

On 8/22/07, Kirsten Beyer <[EMAIL PROTECTED]> wrote:
> Hello.  I am looking for a function that will allow me to paste rows
> together without regard for the numbers of columns in the datasets to
> be joined.  The only columns where it matters if they are aligned
> correctly are at the beginning - the rest of the columns represent
> differing numbers of ICD9 (disease) codes reported by each
> person(record) at a health visit.  They are in no particular order.
>
> For example, a result would look like this:
>
> patient  ICD91  ICD92  ICD93
> patient A   12345  67891543
> patient B3469   9090
> patient C   1234
>
> I am trying to accomplish this inside a loop which first identifies
> the codes associated with the person and then joins them to the
> person.  I have the code working so that it can create a row for each
> person, but I can't figure out how to join these rows together!  FYI,
> my dataset has 200,000+ people.
>
> Thanks
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


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

What is the problem you are trying to solve?

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] gWidgets (tcltk): problem extracting values from widgets in glayout grid

2007-08-22 Thread englund . evan

Hello,

I haven't been able to find an example for the second case below -- or
perhaps I didn't recognize it when I saw it.
Is there a value for x such that svalue(x) will return "bbb", either by
itself or as part of an array? Or do I need to do something else
entirely?
(R2.5.1; Windows XP)

>   gWidgets test
> options("guiToolkit"="tcltk")
> require(gWidgets)
[1] TRUE
>
> ## Case 1: this works 
> wtesta = gwindow("wtesta",visible=TRUE)
> wtesta1 = gedit("aaa",container=wtesta)
> print(svalue(wtesta1))
[1] "aaa"
>
> ### Case 2: this makes an identical window & widget, with "b"
replacing "a"
> ### but does't return the value
> wtesta = gwindow("wtestb",visible=TRUE)
> wtestb = glayout(visible=TRUE,container=wtesta)
> wtestb[1,1]=gedit("bbb",container=wtestb)
> print(svalue(wtestb[1,1]))
Error in function (classes, fdef, mtable)  :
unable to find an inherited method for function ".leftBracket",
for signature "gLayouttcltk", "guiWidgetsToolkittcltk"
Error in svalue(wtestb[1, 1]) : error in evaluating the argument 'obj'
in selecting a method for function 'svalue'
> print(svalue(wtestb))
Error in function (classes, fdef, mtable)  :
unable to find an inherited method for function ".svalue", for
signature "gLayouttcltk", "NULL", "NULL", "guiWidgetsToolkittcltk"

Evan Englund
U.S. EPA
702-798-2248

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


Re: [R] integrate

2007-08-22 Thread Ravi Varadhan
You can divide your domain of integration into smaller intervals and then
add up the individual contributions.  This could improve the speed of
adaptive Gauss-Kronrod quadrature used in integrate().

Ravi.


---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: [EMAIL PROTECTED]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 





-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Santanu Pramanik
Sent: Wednesday, August 22, 2007 6:41 PM
To: [EMAIL PROTECTED]
Cc: r-help@stat.math.ethz.ch; [EMAIL PROTECTED]
Subject: Re: [R] integrate

Hi Andy,

Thank your very much for your input. I also tried something like that
which gives a value close to 20, basically using the same trapezoidal
rule.

> sum(apply(as.matrix(seq(-10,10,by=0.1)),1,my.fcn))*0.1
[1] 20.17385

Actually my function is much more complicated than the posted example
and it is taking a long time...

Anyway, thanks again,
Santanu


JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916



>>> <[EMAIL PROTECTED]> 8/22/2007 5:20:05 PM >>>
As Duncan Murdoch mentioned in his reply, the problem is with the fact
that
your function is not really a properly defined function in the sense
of
"assigning a unique y to each x".  The integrate function uses an
adaptive
quadrature routine which probably makes multiple calls to the function
being integrated and expects to get the same y's for the same x's
every
time.

If you want to get a number close to 20 (for your example) you need an
integration routine which will use single evaluation of your "function"
at
each value of x.  A simple method like rectangular approximation on a
grid
or the so-called trapezoidal rule will do just that.

Here is a very crude prototype of such an integrator:

integrate1 <- function(f, lower, upper){
   f <- match.fun(f)
   xx <- seq(lower, upper, length=100)
   del <- xx[2] - xx[1]
   yy <- f(xx[-100])
return(del*sum(yy))

Now when you run integrate1(my.fun, -10, 10) you will get a number
close to
20 but, of course, every time you do it you will get a different
value.

Hope this helps,

Andy

__
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-
E-mail: [EMAIL PROTECTED] 
Tel:  (651) 733-6092
Fax:  (651) 736-3122


   
   
 "Santanu  
   
 Pramanik" 
   
 <[EMAIL PROTECTED] 
To 
 .umd.edu>   
   
 Sent by:  
cc 
 [EMAIL PROTECTED] 
   
 at.math.ethz.ch  
Subject 
   [R] integrate   
   
   
   
 08/22/2007 02:56  
   
 PM
   
   
   
   
   
   
   




Hi,
I am trying to integrate a function which is approximately constant
over the range of the integration. The function is as follows:

> my.fcn = function(mu){
+ m = 1000
+ z = 0
+ z.mse = 0
+ for(i in 1:m){
+ z[i] = rnorm(1, mu, 1)
+ z.mse = z.mse + (z[i] - mu)^2
+ }
+ return(z.mse/m)
+ }

> my.fcn(-10)
[1] 1.021711
> my.fcn(10)
[1] 0.9995235
> my.fcn(-5)
[1] 1.012727
> my.fcn(5)
[1] 1.033595
> my.fcn(0)
[1] 1.106282
>
The function takes the value (approx) 1 over the range of mu. So the
integration result should be close to 20 if we integrate over (-10,
10).
But R gives:

> integrate(my.fcn, -10, 10)
685.4941 with absolute error < 7.6e-12

> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate

I have seen in the "?integrate" page it is clearly written:

If the function is approximately constant (in particular, zero) over
nearly all its range it is possible that the result and error estimate
may be seriously wrong.

But this doesn't help in solving the problem.
Thanks,
Santanu



JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916



 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list

Re: [R] integrate

2007-08-22 Thread Santanu Pramanik
Hi Andy,

Thank your very much for your input. I also tried something like that
which gives a value close to 20, basically using the same trapezoidal
rule.

> sum(apply(as.matrix(seq(-10,10,by=0.1)),1,my.fcn))*0.1
[1] 20.17385

Actually my function is much more complicated than the posted example
and it is taking a long time...

Anyway, thanks again,
Santanu


JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916



>>> <[EMAIL PROTECTED]> 8/22/2007 5:20:05 PM >>>
As Duncan Murdoch mentioned in his reply, the problem is with the fact
that
your function is not really a properly defined function in the sense
of
"assigning a unique y to each x".  The integrate function uses an
adaptive
quadrature routine which probably makes multiple calls to the function
being integrated and expects to get the same y's for the same x's
every
time.

If you want to get a number close to 20 (for your example) you need an
integration routine which will use single evaluation of your "function"
at
each value of x.  A simple method like rectangular approximation on a
grid
or the so-called trapezoidal rule will do just that.

Here is a very crude prototype of such an integrator:

integrate1 <- function(f, lower, upper){
   f <- match.fun(f)
   xx <- seq(lower, upper, length=100)
   del <- xx[2] - xx[1]
   yy <- f(xx[-100])
return(del*sum(yy))

Now when you run integrate1(my.fun, -10, 10) you will get a number
close to
20 but, of course, every time you do it you will get a different
value.

Hope this helps,

Andy

__
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-
E-mail: [EMAIL PROTECTED] 
Tel:  (651) 733-6092
Fax:  (651) 736-3122


   
   
 "Santanu  
   
 Pramanik" 
   
 <[EMAIL PROTECTED] 
To 
 .umd.edu>   
   
 Sent by:  
cc 
 [EMAIL PROTECTED] 
   
 at.math.ethz.ch  
Subject 
   [R] integrate   
   
   
   
 08/22/2007 02:56  
   
 PM
   
   
   
   
   
   
   




Hi,
I am trying to integrate a function which is approximately constant
over the range of the integration. The function is as follows:

> my.fcn = function(mu){
+ m = 1000
+ z = 0
+ z.mse = 0
+ for(i in 1:m){
+ z[i] = rnorm(1, mu, 1)
+ z.mse = z.mse + (z[i] - mu)^2
+ }
+ return(z.mse/m)
+ }

> my.fcn(-10)
[1] 1.021711
> my.fcn(10)
[1] 0.9995235
> my.fcn(-5)
[1] 1.012727
> my.fcn(5)
[1] 1.033595
> my.fcn(0)
[1] 1.106282
>
The function takes the value (approx) 1 over the range of mu. So the
integration result should be close to 20 if we integrate over (-10,
10).
But R gives:

> integrate(my.fcn, -10, 10)
685.4941 with absolute error < 7.6e-12

> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate

I have seen in the "?integrate" page it is clearly written:

If the function is approximately constant (in particular, zero) over
nearly all its range it is possible that the result and error estimate
may be seriously wrong.

But this doesn't help in solving the problem.
Thanks,
Santanu



JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916



 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] triviality solved: Unwanted but added additional first column when using write.csv2()

2007-08-22 Thread Delcour Libertus
Hello!

When I export a dataframe to a csv-file there is always an unwanted
column added as the following demonstrates:


> system("cat example.csv")
ID;Name;Value1;Value2
100;Tom;11,90;32,90
101;Martha;15,90;49,00
102;Ute;20,00;300,00

> write.csv2(read.csv2 ("example.csv"), "examplecopy.csv")
> system("cat examplecopy.csv")
"";"ID";"Name";"Value1";"Value2"
"1";100;"Tom";11,9;32,9
"2";101;"Martha";15,9;49
"3";102;"Ute";20;300
>

I never want the column

""
"1"
"2"
"3"

to be written into a csv-file.

I could find any hint neither in the manual nor in the archive, but then
I found the solution: row.names = FALSE

> write.csv2(read.csv2 ("example.csv"), "examplecopy.csv", row.names =
FALSE)
> system("cat examplecopy.csv")
"ID";"Name";"Value1";"Value2"
100;"Tom";11,9;32,9
101;"Martha";15,9;49
102;"Ute";20;300

But (perhaps because I am not a native in English) it was more a
solution by trial and chance than by reading and understanding the
manual. I did not expect this unwanted column to be called "row.names"
(the values are numbers which at last becomes after a subset()
completely unusuable in most cases). So I did never search for row names
but for "unwanted column write."

I did not dare to ask before because of the triviality of the question
(even I stumble over this problem for weeks) but I want the next
beginner who searches this list archive for "unwanted column" to get a
result.

Greetings

Delcour

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


Re: [R] integrate

2007-08-22 Thread apjaworski
As Duncan Murdoch mentioned in his reply, the problem is with the fact that
your function is not really a properly defined function in the sense of
"assigning a unique y to each x".  The integrate function uses an adaptive
quadrature routine which probably makes multiple calls to the function
being integrated and expects to get the same y's for the same x's every
time.

If you want to get a number close to 20 (for your example) you need an
integration routine which will use single evaluation of your "function" at
each value of x.  A simple method like rectangular approximation on a grid
or the so-called trapezoidal rule will do just that.

Here is a very crude prototype of such an integrator:

integrate1 <- function(f, lower, upper){
   f <- match.fun(f)
   xx <- seq(lower, upper, length=100)
   del <- xx[2] - xx[1]
   yy <- f(xx[-100])
return(del*sum(yy))

Now when you run integrate1(my.fun, -10, 10) you will get a number close to
20 but, of course, every time you do it you will get a different value.

Hope this helps,

Andy

__
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-
E-mail: [EMAIL PROTECTED]
Tel:  (651) 733-6092
Fax:  (651) 736-3122


   
 "Santanu  
 Pramanik" 
 <[EMAIL PROTECTED]  To 
 .umd.edu>   
 Sent by:   cc 
 [EMAIL PROTECTED] 
 at.math.ethz.ch   Subject 
   [R] integrate   
   
 08/22/2007 02:56  
 PM
   
   
   




Hi,
I am trying to integrate a function which is approximately constant
over the range of the integration. The function is as follows:

> my.fcn = function(mu){
+ m = 1000
+ z = 0
+ z.mse = 0
+ for(i in 1:m){
+ z[i] = rnorm(1, mu, 1)
+ z.mse = z.mse + (z[i] - mu)^2
+ }
+ return(z.mse/m)
+ }

> my.fcn(-10)
[1] 1.021711
> my.fcn(10)
[1] 0.9995235
> my.fcn(-5)
[1] 1.012727
> my.fcn(5)
[1] 1.033595
> my.fcn(0)
[1] 1.106282
>
The function takes the value (approx) 1 over the range of mu. So the
integration result should be close to 20 if we integrate over (-10, 10).
But R gives:

> integrate(my.fcn, -10, 10)
685.4941 with absolute error < 7.6e-12

> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate

I have seen in the "?integrate" page it is clearly written:

If the function is approximately constant (in particular, zero) over
nearly all its range it is possible that the result and error estimate
may be seriously wrong.

But this doesn't help in solving the problem.
Thanks,
Santanu



JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916



 [[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Need a variant of rbind for datasets with different numbers of columns

2007-08-22 Thread Kirsten Beyer
Hello.  I am looking for a function that will allow me to paste rows
together without regard for the numbers of columns in the datasets to
be joined.  The only columns where it matters if they are aligned
correctly are at the beginning - the rest of the columns represent
differing numbers of ICD9 (disease) codes reported by each
person(record) at a health visit.  They are in no particular order.

For example, a result would look like this:

patient  ICD91  ICD92  ICD93
patient A   12345  67891543
patient B3469   9090
patient C   1234

I am trying to accomplish this inside a loop which first identifies
the codes associated with the person and then joins them to the
person.  I have the code working so that it can create a row for each
person, but I can't figure out how to join these rows together!  FYI,
my dataset has 200,000+ people.

Thanks

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


Re: [R] plotting lda results

2007-08-22 Thread Silvia Lomascolo

Thanks for trying to answer even if my question isn't very clear. Please
allow me to try again with an example:
I want to predict group membership by the three variables below, and then
plot the results in one figure that compares the frequency of scores of both
groups on the LD obtained:
Group<- c(1,1,1,1,1,1,1,2,2,2,2,2,2,2,2)
var1<- c(2,3,1,2,3,2,3,3,5,6,7,6,8,5,5)
var2<- c(9,9,9,8,7,8,9,3,2,2,1,1,2,3,3)
var3<- c(6,7,6,6,5,6,7,1,2,1,2,3,1,1,2)
data.df <- as.data.frame (cbind(Group, var1, var2, var3)) 
data.lda <- lda(Group~., data=data.df)
predict(data.lda)
plot(data.lda)

As you said, I get two separate figures, one for each group, but I need them
both in one figure.  I tried ldahist by reading the help because it has an
argument sep, but I'm not managing to get it going.

I hope this code is reproducible but I'll certainly read the book suggested
to make my questions easier to answer.

Thanks! Silvia.


Prof Brian Ripley wrote:
> 
> Read ?plot.lda, which tells you the ... arguments are (for dimen=1, the 
> only option for two groups) passed to ldahist, so then read its help page.
> 
> I tried what is suggested plot.lda but 
> 
> I don't know what you want (and your example is not reproducible): I would 
> expect you to get a single plot with two panels (figures), but there are 
> options to have a single panel.  (Reading 'An Introduction to R' may help 
> you to use standard terminology that others will be able to follow.)
> 
> On Wed, 22 Aug 2007, Silvia Lomascolo wrote:
> 
>>
>> Hi all,
>> I am trying to plot the results of a discriminant analysis done with
>> lda(MASS) but my groups appear in two different plots (in the same
>> graphics
>> device) and I want to combine them in one plot. My code looks like:
>>
>> BirdTrain.lda <- lda(Bdisperser~., data=BirdTrain.mx)
>> predict(BirdTrain.lda)
>> plot(BirdTrain.lda)
>>
>> I have two types of Bdisperser, so I only get one linear discriminant
>> function. Can anyone please tell me how to combine the data in one plot?
>>
>> I work with R 2.4.1 using Windows.
> 
> But the version of MASS is what is relevant, and it would have been in 
> the sessionInfo() output the R posting guide asked you for.
> 
> -- 
> 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@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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/plotting-lda-results-tf4312870.html#a12282028
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] distance by vegan

2007-08-22 Thread Victor Landeiro
Hi,
In vegan the function to calculate distance is vegdist in which the default
distance is Sorensen (Bray-Curtis).
So, try:
vegdist(x, method="bray")
Also see ?vegdist

victor

On 8/22/07, [EMAIL PROTECTED] <[EMAIL PROTECTED]> wrote:
>
> How to calculate sorensen (bray-curtis) distance by dist  function
> within the vegan package?
>
> Cheers
> Duccio
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Victor Lemes Landeiro
Instituto Nacional de Pesquisas da Amazônia - INPA
Manaus, Amazonas, Brasil
Skypename: landeiro

[[alternative HTML version deleted]]

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


Re: [R] integrate

2007-08-22 Thread Duncan Murdoch
On 8/22/2007 3:54 PM, Santanu Pramanik wrote:
> Hi,
> I am trying to integrate a function which is approximately constant
> over the range of the integration. The function is as follows:

That's not a function of the input mu.  It includes a random component:

 > my.fcn(10)
[1] 0.9786558
 > my.fcn(10)
[1] 1.022467

You can't expect integrate() to return a sensible answer if you don't 
give it a function that returns consistent results.

Duncan Murdoch

>  
>> my.fcn = function(mu){
> + m = 1000
> + z = 0
> + z.mse = 0
> + for(i in 1:m){
> + z[i] = rnorm(1, mu, 1)
> + z.mse = z.mse + (z[i] - mu)^2
> + }
> + return(z.mse/m)
> + }
> 
>> my.fcn(-10)
> [1] 1.021711
>> my.fcn(10)
> [1] 0.9995235
>> my.fcn(-5)
> [1] 1.012727
>> my.fcn(5)
> [1] 1.033595
>> my.fcn(0)
> [1] 1.106282
>> 
> The function takes the value (approx) 1 over the range of mu. So the
> integration result should be close to 20 if we integrate over (-10, 10).
> But R gives:
> 
>> integrate(my.fcn, -10, 10)
> 685.4941 with absolute error < 7.6e-12
> 
>> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate
>  
> I have seen in the "?integrate" page it is clearly written:
>  
> If the function is approximately constant (in particular, zero) over
> nearly all its range it is possible that the result and error estimate
> may be seriously wrong. 
>  
> But this doesn't help in solving the problem.
> Thanks,
> Santanu
> 
> 
>  
> JPSM, 1218J Lefrak Hall
> University of Maryland, College Park
> Phone: 301-314-9916
>  
>  
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] integrate

2007-08-22 Thread Santanu Pramanik
Hi,
I am trying to integrate a function which is approximately constant
over the range of the integration. The function is as follows:
 
> my.fcn = function(mu){
+ m = 1000
+ z = 0
+ z.mse = 0
+ for(i in 1:m){
+ z[i] = rnorm(1, mu, 1)
+ z.mse = z.mse + (z[i] - mu)^2
+ }
+ return(z.mse/m)
+ }

> my.fcn(-10)
[1] 1.021711
> my.fcn(10)
[1] 0.9995235
> my.fcn(-5)
[1] 1.012727
> my.fcn(5)
[1] 1.033595
> my.fcn(0)
[1] 1.106282
> 
The function takes the value (approx) 1 over the range of mu. So the
integration result should be close to 20 if we integrate over (-10, 10).
But R gives:

> integrate(my.fcn, -10, 10)
685.4941 with absolute error < 7.6e-12

> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate
 
I have seen in the "?integrate" page it is clearly written:
 
If the function is approximately constant (in particular, zero) over
nearly all its range it is possible that the result and error estimate
may be seriously wrong. 
 
But this doesn't help in solving the problem.
Thanks,
Santanu


 
JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916
 
 

[[alternative HTML version deleted]]

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

2007-08-22 Thread Santanu Pramanik
Hi,
I am trying to integrate a function which is approximately constant
over the range of the integration. The function is as follows:
 
> my.fcn = function(mu){
+ m = 1000
+ z = 0
+ z.mse = 0
+ for(i in 1:m){
+ z[i] = rnorm(1, mu, 1)
+ z.mse = z.mse + (z[i] - mu)^2
+ }
+ return(z.mse/m)
+ }

> my.fcn(-10)
[1] 1.021711
> my.fcn(10)
[1] 0.9995235
> my.fcn(-5)
[1] 1.012727
> my.fcn(5)
[1] 1.033595
> my.fcn(0)
[1] 1.106282
> 
The function takes the value (approx) 1 over the range of mu. So the
integration result should be close to 20 if we integrate over (-10, 10).
But R gives:

> integrate(my.fcn, -10, 10)
685.4941 with absolute error < 7.6e-12

> integrate(Vectorize(my.fcn), -10, 10)  # this code never terminate
 
I have seen in the "?integrate" page it is clearly written:
 
If the function is approximately constant (in particular, zero) over
nearly all its range it is possible that the result and error estimate
may be seriously wrong. 
 
But this doesn't help in solving the problem.
Thanks,
Santanu


 
JPSM, 1218J Lefrak Hall
University of Maryland, College Park
Phone: 301-314-9916
 
 

[[alternative HTML version deleted]]

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


Re: [R] All subsets regression

2007-08-22 Thread Thomas Lumley
On Wed, 22 Aug 2007, Alan Harrison wrote:

> Hey folks,
>
> I'm trying to do all subsets on a zero-inflated poisson regression. 
> I'm aware of the leaps and regsubsets functions but I don't know if they 
> work for ZIP regressions or how the syntax fits in for them.

They don't.

-thomas

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

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


Re: [R] tackle memory insufficiency for large dataset using save() & load()??

2007-08-22 Thread Charles C. Berry
On Tue, 21 Aug 2007, Jessica Z wrote:


[snip]


I did not notice a comment on this bit in the other replies:

>>
>> newdata <- load ("compact_d.Rdata")
>>
>> summary(newdata)
>   Length Class  Mode
>1 character character

newdata is a string whose value is 'd'

try print( newdata )

ls() should tell you there are two objects - 'd' and 'newdata'

So just continue using 'd', e.g.

summary( d )

HTH,

Chuck

[snip]

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


Re: [R] Output from while and for loop

2007-08-22 Thread Charles C. Berry
On Tue, 21 Aug 2007, Ryan Briscoe Runquist wrote:

>
> Hello,
>
> I am new and am having a hard time getting the proper syntax for output
> from loops.  I am working on a simulation to generate a null expectation of
> bee behavior.  Pieces of it work.  The part that I am having specific
> difficulty is in output of a vector from within the while loop that I am
> using.


Hmmm. I think I can guess at your problem.

You need to assign the result to a vector or list element along these 
lines:

my.list <- list()
k <- 0
while( keep.going ){
   k <- k+1
   ## do lots of stuff
   ## details omitted
   my.list[[ k ]] <- result.of.doing.lots.of.stuff
}

unlist( my.list ) returns a vector if result.of.doing.lots.of.stuff is 
atomic.

There are lots of variations on this. If the result ('my.list') will be 
very long and each element is atomic, it helps speed to set it up a vector 
ahead of time and then fill one element per pass thru the loop.


Basically the simulation works as such:  I have a starting point
> and a neighbor matrix and a certain threshold distance for travel.  In the
> while loop the "bee" moves to a randomly chosen neighbor location.  I want
> to be able to record the elevations of these points (including the starting
> point) so that I can look at variance in elevation and mean elevation.  The
> loop itself works as does the calculation of the final elevation list,
> change in elevation list, and true total distance traveled.  I have looked
> in all of the email archives but have not come across a correct way of
> doing it.  Code below:
>
> start.elev.list<-list()
> final.mean.elev.list<-list()
> final.elev.list<-list()
> final.distance.list<-list()
> final.delta.elev.list<-list()
> final.var.elev<-list()
>
>
> b<-length(Bees.Day.1$bee)
> for (bee in 1:b){
>   #this is for number of bees that are trackable in the day with starting
> points and threshold distances
>   elev.current.vector<-vector(mode="numeric", length=0)
>   count<-1
>   ElevSS<-0
>   d.traveled<-0
>   thresh<-Bees.Day.1$cum.dist[bee]
>   n<-Bees.Day.1$grid.pt[bee]
>   #I'm making this up for the threshold, want to be bee specific
>   #current.point<-round(runif(1,1,n)) #random starting point
>   current.point<-Day.1.neighbor.matrix[1,n]
>   #I want to specify the first point in the matrix
>   Elev.Sum<-Day.1.elev.vector[current.point]
>
>
>   while(d.traveled   #which of the four options will be selection
>   transition<-round(runif(1,1,4))
>
>   #so, what's the new point?
>   new.point<- Day.1.neighbor.matrix[transition,n]
>
>   #what is the variance in elevation changed
>   Elev.current<-Day.1.elev.vector[current.point]
>   elev.current.vector[i]<-Elev.current
>   Elev.new<-Day.1.elev.vector[new.point]
>   Elev.Sum<-(Elev.Sum+Elev.new)
>
>   #how far will bee travelled
>   current.travel<- Day.1.distance.matrix[current.point, new.point]
>   d.traveled<- current.travel + d.traveled
>   current.point<- new.point
>
>   #Number of iterations until we reach the threshold
>   count<-count+1
>
>   }
>
>   print(count)
>   print(elev.current.vector)
>   mean.elev<-Elev.Sum/count
>   print(paste("Final mean elev for bee", bee, "is", mean.elev, sep=" "))
>   final.mean.elev.list[bee]<-list(mean.elev)
>
>   #What was the start elevation?
>   start.elev<-Day.1.elev.vector[n]
>   print(paste("Start elev for bee",bee,"is",start.elev, sep=" "))
>   start.elev.list[bee]<-list(start.elev)
>
>   #what is the final elevation?
>   final.elev<-Day.1.elev.vector[current.point]
>   print(paste("Final elev for bee",bee,"is", final.elev,sep=" "))
>   final.elev.list[bee]<-list(final.elev)
>
>   print(paste("Final travel distance for bee", bee,"is", d.traveled, sep=" "))
>   final.distance.list[bee]<-list(d.traveled)
>
>   net.delta.elev<-(final.elev-Day.1.elev.vector[n])
>   print(paste("Final net change in elevation for bee",bee,"is",
> net.delta.elev,sep=" "))
>   final.delta.elev.list[bee]<-list(net.delta.elev)
>
> }
> ~~
> Ryan D. Briscoe Runquist
> Population Biology Graduate Group
> University of California, Davis
> [EMAIL PROTECTED]
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] plotting lda results

2007-08-22 Thread Prof Brian Ripley
Read ?plot.lda, which tells you the ... arguments are (for dimen=1, the 
only option for two groups) passed to ldahist, so then read its help page.

I don't know what you want (and your example is not reproducible): I would 
expect you to get a single plot with two panels (figures), but there are 
options to have a single panel.  (Reading 'An Introduction to R' may help 
you to use standard terminology that others will be able to follow.)

On Wed, 22 Aug 2007, Silvia Lomascolo wrote:

>
> Hi all,
> I am trying to plot the results of a discriminant analysis done with
> lda(MASS) but my groups appear in two different plots (in the same graphics
> device) and I want to combine them in one plot. My code looks like:
>
> BirdTrain.lda <- lda(Bdisperser~., data=BirdTrain.mx)
> predict(BirdTrain.lda)
> plot(BirdTrain.lda)
>
> I have two types of Bdisperser, so I only get one linear discriminant
> function. Can anyone please tell me how to combine the data in one plot?
>
> I work with R 2.4.1 using Windows.

But the version of MASS is what is relevant, and it would have been in 
the sessionInfo() output the R posting guide asked you for.

-- 
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] All subsets regression

2007-08-22 Thread Alan Harrison
Hey folks,

I'm trying to do all subsets on a zero-inflated poisson regression.  I'm aware 
of the leaps and regsubsets functions but I don't know if they work for ZIP 
regressions or how the syntax fits in for them.  Can anyone help?
The model syntax is:
> zip.zc <- zicounts(resp=Scars~.,x=~Location + Lar + Mass + Lac + Lar:Mass + 
> Location:Mass, data = alan2, distr="ZIP")

Many thanks

Alan Harrison

Quercus
Queen's University Belfast
MBC, 97 Lisburn Road
Belfast

BT9 7BL

T: 02890 972219
M: 07798615682


[[alternative HTML version deleted]]

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


[R] plotting lda results

2007-08-22 Thread Silvia Lomascolo

Hi all,
I am trying to plot the results of a discriminant analysis done with
lda(MASS) but my groups appear in two different plots (in the same graphics
device) and I want to combine them in one plot. My code looks like:

BirdTrain.lda <- lda(Bdisperser~., data=BirdTrain.mx)
predict(BirdTrain.lda)
plot(BirdTrain.lda)

I have two types of Bdisperser, so I only get one linear discriminant
function. Can anyone please tell me how to combine the data in one plot?

I work with R 2.4.1 using Windows.
-- 
View this message in context: 
http://www.nabble.com/plotting-lda-results-tf4312870.html#a12279005
Sent from the R help mailing list archive at Nabble.com.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] distance by vegan

2007-08-22 Thread [EMAIL PROTECTED]
How to calculate sorensen (bray-curtis) distance by dist  function 
within the vegan package?

Cheers
Duccio

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Rdonlp2 control parameter question

2007-08-22 Thread Talbot Katz
Hi.

I have been trying to work with the Rdonlp2 package.  There is a control 
parameter called epsx that has the following description:

epsx(1e-5) - successful temination is indicated if the Kuhn-Tucker criteria 
are satisfied within the value

I have tried to use the donlp2.control() function to adjust the value of 
epsx, e.g.:

>cntlx1_3<-donlp2.control(epsx=1e-03)
>cntlx1_3$epsx
[1] 0.001
>

I have adjusted it larger (as above) and smaller, but it has not had any 
effect on the solution when I run donlp2(...,control=cntlx1_3).  And the 
console output from running the donlp2() function reports a value of epsx, 
which is always "epsx= 1.000e-05".

Am I doing something wrong?  Is there a way to adjust the value of epsx?  
I've been comparing donlp2() with quadprog() on a quadratic programming 
problem; the objective function values are pretty close, but Id like to see 
if I can get them any closer by adjusting epsx.  Any help is appreciated.  
Thanks!

--  TMK  --
212-460-5430home
917-656-5351cell

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] 2nd label category on x axis

2007-08-22 Thread mgilgen
i'm plotting time-series and label the x-axis as follows:

r <- as.POSIXct(round(range(p1$time), "month"))

plot(p1$time,p1$ p1, type="l", xaxt="n")

axis.POSIXct(1, at=seq(r[1], r[2], by="month"), format="%Y.%m")

what I want do do now, is a second label for the x-axis, that stands lower and 
indicates the years.

I don't know how to proceed with all the possibilities in axis, par and plot

thanks for your help

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


Re: [R] Help with vector gymnastics

2007-08-22 Thread Erik Iverson
Philip -

I don't know if this is the "best" way, but it gives you the output you 
want.

Using your tf,

vals <- rle(ifelse(tf, 5*which(tf), 0))
vals$values[vals$values == 0] <- vals$values[which(vals$values==0) - 1]
inverse.rle(vals)
[1]  5  5  5  5 25 30 30

Gladwin, Philip wrote:
> Hello,
> 
> What is the best way of solving this problem?
> 
> answer <- ifelse(tf=TRUE, i * 5, previous answer)
> where as an initial condition 
> tf[1] <- TRUE
> 
> 
> For example if,
> tf <- c(T,F,F,F,T,T,F)
> over i = 1 to 7
> then the output of the function will be
> answer = 5 5 5 5 25 30 30 
> 
> Thank you.
> 
> Phil,
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Help with vector gymnastics

2007-08-22 Thread Felix Andrews
library(zoo)
tf <- c(T,F,F,F,T,T,F)
i <- seq(7)
answer <- ifelse(tf, i*5, NA)
answer <- na.locf(answer)


On 8/23/07, Gladwin, Philip <[EMAIL PROTECTED]> wrote:
> Hello,
>
> What is the best way of solving this problem?
>
> answer <- ifelse(tf=TRUE, i * 5, previous answer)
> where as an initial condition
> tf[1] <- TRUE
>
>
> For example if,
> tf <- c(T,F,F,F,T,T,F)
> over i = 1 to 7
> then the output of the function will be
> answer = 5 5 5 5 25 30 30
>
> Thank you.
>
> Phil,
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Felix Andrews / 安福立
PhD candidate
Integrated Catchment Assessment and Management Centre
The Fenner School of Environment and Society
The Australian National University (Building 48A), ACT 0200
Beijing Bag, Locked Bag 40, Kingston ACT 2604
http://www.neurofractal.org/felix/
voice:+86_1051404394 (in China)
mobile:+86_13522529265 (in China)
mobile:+61_410400963 (in Australia)
xmpp:[EMAIL PROTECTED]
3358 543D AAC6 22C2 D336  80D9 360B 72DD 3E4C F5D8

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


Re: [R] within-subject factors in lme

2007-08-22 Thread lorenz.gygax
> > If I understand correctly, you want to include the interactions  
> > between the random and fixed terms?
> 
> Yes that is exactly I wanted to model.
> 
> > This is done like:
> >
> > model.lme <- lme(Beta ~ Trust*Sex*Freq,
> >  random = ~Trust*Sex*Freq|Subj, Model)
> >
> > But this needs a lot of observations as quite a few 
> > parameters need to be estimated!
> 
> Well, I tried this as well, but it seems R kept hanging there and  
> never finished the modeling. It is very likely due to some  
> singularity as you suspected about the large number of parameters  
> needed to estimate. But this is not a problem with aov. So does
> it mean that I can't run a similar model to that in aov with lme?

It depends what you mean by 'similar'. You could still include some of the 
interactions, e.g. by random = ~(Trust+Sex+Freq)^2|Subj, or even further 
reduced such as ~Trust+Sex+Freq|Subj. I am not very familiar with aov, but I 
would suspect that the model you calcualted in aov is not really the same than 
the one with all possible interactions in lme. In any case, I would personally 
trust lme much more than aov.

> but I feel this is not good enough to account for cross-subject  
> variations for those interactions. Why wouldn't those patterned  
> variance-covariance matrix specifications work as I mentioned in
> my previous mail? Any more thoughts and suggestions?

Sorry, I have never really worked with those.

Lorenz
- 
Lorenz Gygax
Centre for proper housing of ruminants and pigs
Agroscope Reckenholz-Tänikon Research Station ART

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


[R] Help with vector gymnastics

2007-08-22 Thread Gladwin, Philip
Hello,

What is the best way of solving this problem?

answer <- ifelse(tf=TRUE, i * 5, previous answer)
where as an initial condition 
tf[1] <- TRUE


For example if,
tf <- c(T,F,F,F,T,T,F)
over i = 1 to 7
then the output of the function will be
answer = 5 5 5 5 25 30 30 

Thank you.

Phil,

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 do i use the get function to obtain an element from alist...

2007-08-22 Thread Juan Manuel Barreneche
To Mark Leeds and the others, thank you for solving my problem, and so quickly,

JM

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


Re: [R] within-subject factors in lme

2007-08-22 Thread Gang Chen
Hi Lorenz,

I really appreciate your comments.

> If I understand correctly, you want to include the interactions  
> between the random and fixed terms?

Yes that is exactly I wanted to model.

> This is done like:
>
> model.lme <- lme(Beta ~ Trust*Sex*Freq,
>  random = ~Trust*Sex*Freq|Subj, Model)
>
> But this needs a lot of observations as quite a few parameters need  
> to be estimated!

Well, I tried this as well, but it seems R kept hanging there and  
never finished the modeling. It is very likely due to some  
singularity as you suspected about the large number of parameters  
needed to estimate. But this is not a problem with aov. So does it  
mean that I can't run a similar model to that in aov with lme?

Sure I can simply run the following model

fit.lme <- lme(Beta ~ Trust*Sex*Freq, random = ~1|Subj, Model)

but I feel this is not good enough to account for cross-subject  
variations for those interactions. Why wouldn't those patterned  
variance-covariance matrix specifications work as I mentioned in my  
previous mail? Any more thoughts and suggestions?

> Possibly, you can not include the variable Sex in this, because I  
> assume that Subj is nested within Sex.

No, Sex is NOT a subject classifying factor. Instead it is a task- 
related within-subject factor.

Again thanks a lot for your help,
Gang

On Aug 22, 2007, at 6:52 AM, <[EMAIL PROTECTED]>  
<[EMAIL PROTECTED]> wrote:

> I don't think, this has been answered:
>
>> I'm trying to run a 3-way within-subject anova in lme with 3
>> fixed factors (Trust, Sex, and Freq), but get stuck with handling
>> the  random effects. As I want to include all the possible random
>> effects in the model, it would be something more or less
>> equivalent to using aov
>>
>>> fit.aov <- aov(Beta ~
>> Trust*Sex*Freq+Error(Subj/(Trust*Sex*Freq)),
>> Model)
>>
>> However I'm not so sure what I should do in lme. Sure
>>
>>> lme(Beta ~ Trust*Sex*Freq, random = ~1|Subj, Model)
>>
>> works fine, but it only models the random effect of the
>> intercept. I tried the following 4 possibilities:
>
> If I understand correctly, you want to include the interactions  
> between the random and fixed terms? This is done like:
>
> model.lme <- lme(Beta ~ Trust*Sex*Freq,
>  random = ~Trust*Sex*Freq|Subj, Model)
>
> But this needs a lot of observations as quite a few parameters need  
> to be estimated! Possibly, you can not include the variable Sex in  
> this, because I assume that Subj is nested within Sex. If you just  
> refer to within and between subject effects and their corresponding  
> degrees of freedom: you should see this being handled automatically  
> and correctly by lme e.g. in the output of anova (model.lme)
>
> Lorenz
> -
> Lorenz Gygax
> Centre for proper housing of ruminants and pigs
> Agroscope Reckenholz-Tänikon Research Station ART
> Tänikon, CH-8356 Ettenhausen / Switzerland

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


Re: [R] Evaluating f(x(2,3)) on a function f<- function(a,b){a+b}

2007-08-22 Thread Paul Smith
On 8/22/07, Søren Højsgaard <[EMAIL PROTECTED]> wrote:
> I have a function and a vector, say
> f <- function(a,b){a+b}
> x <- c(2,3)
> I want to "evaluate f on x" in the sense of computing f(x[1],x[2]). I would 
> like it to be so that I can write f(x). (I know I can write a wrapper 
> function g <- function(x){f(x[1],x[2])}, but this is not really what I am 
> looking for). Is there a general way doing this (programmatically)? (E.g. by 
> "unpacking" the elements of x and putting them in the "right places" when 
> calling f...)
> I've looked under formals, alist etc. but so far without luck.

I hope that the following helps:

> f <- function(x) {sum(x)}
> f(c(2,3))
[1] 5
> f(c(2,3,5))
[1] 10
>

Paul

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


Re: [R] Evaluating f(x(2,3)) on a function f<- function(a,b){a+b}

2007-08-22 Thread Gabor Grothendieck
Try this:

do.call(f, as.list(x))

On 8/22/07, Søren Højsgaard <[EMAIL PROTECTED]> wrote:
> Dear list
> I have a function and a vector, say
>f <- function(a,b){a+b}
>x <- c(2,3)
> I want to "evaluate f on x" in the sense of computing f(x[1],x[2]). I would 
> like it to be so that I can write f(x). (I know I can write a wrapper 
> function g <- function(x){f(x[1],x[2])}, but this is not really what I am 
> looking for). Is there a general way doing this (programmatically)? (E.g. by 
> "unpacking" the elements of x and putting them in the "right places" when 
> calling f...)
> I've looked under formals, alist etc. but so far without luck.
>
> Regards
> Søren
>
>
>[[alternative HTML version deleted]]
>
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Evaluating f(x(2,3)) on a function f<- function(a,b){a+b}

2007-08-22 Thread Søren Højsgaard
Dear list
I have a function and a vector, say
f <- function(a,b){a+b}
x <- c(2,3)
I want to "evaluate f on x" in the sense of computing f(x[1],x[2]). I would 
like it to be so that I can write f(x). (I know I can write a wrapper function 
g <- function(x){f(x[1],x[2])}, but this is not really what I am looking for). 
Is there a general way doing this (programmatically)? (E.g. by "unpacking" the 
elements of x and putting them in the "right places" when calling f...)
I've looked under formals, alist etc. but so far without luck.
 
Regards
Søren
 

[[alternative HTML version deleted]]

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


Re: [R] Formatting Sweave in R-News

2007-08-22 Thread Gregor Gorjanc
Hi Duncan,

Duncan Murdoch wrote:
>> Slightly off Arjuns topic is a problem with Schunk, Sinput and Soutput
>> environments. They just use to much inline space. I have tweade this
>> problem with sed (see bellow for Makefile content), but wonder if there
>> is a better solution.
> 
> You can change the spacing by redefining those environments.  I use this:
> 
> % This removes the extra spacing after code and output chunks in Sweave,
> % but keeps the spacing around the whole block.
> 
> \fvset{listparameters={\setlength{\topsep}{0pt}}}
> \renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}

Nice trick. Thank you.

I am CCing to a friend that had the same issue.

-- 
Lep pozdrav / With regards,
 Gregor Gorjanc
--
University of Ljubljana PhD student
Biotechnical Facultywww: http://www.bfro.uni-lj.si/MR/ggorjan
Zootechnical Department blog: http://ggorjan.blogspot.com
Groblje 3   mail: gregor.gorjanc  bfro.uni-lj.si
SI-1230 Domzale fax: +386 (0)1 72 17 888
Slovenia, Europetel: +386 (0)1 72 17 861

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


Re: [R] Synchronzing workspaces

2007-08-22 Thread Barry Rowlingson
Eric Turkheimer wrote:
> How do people go about synchronizing multiple workspaces on different
> workstations?  I tend to wind up with projects spread around the various
> machines I work on.  I find that placing the directories on a server and
> reading them remotely tends to slow things down.

  If R were to store all its workspace data objects in individual files 
instead of one big .RData file, then you could use a revision control 
system like SVN.  Check out the data, work on it, check it in, then on 
another machine just update to get the changes.

  However SVN doesn't work too well for binary files - conflicts being 
hard to resolve without someone backing down - so maybe its not such a 
good idea anyway...

  On unix boxes and derivatives, you can keep things in sync efficiently 
with the 'rsync' command.  I think there are GUI addons for it, and 
Windows ports.

Barry

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


Re: [R] Synchronzing workspaces

2007-08-22 Thread Duncan Murdoch
On 8/22/2007 9:20 AM, Eric Turkheimer wrote:
> How do people go about synchronizing multiple workspaces on different
> workstations?  I tend to wind up with projects spread around the various
> machines I work on.  I find that placing the directories on a server and
> reading them remotely tends to slow things down.

I use Subversion.  It works best with text, but can handle binary files too.

So the idea would be to write scripts that either generate or load the 
data I want to work on, and source those at the beginning of a session. 
  I always try to start with a clean workspace.

Separate the code that does complex calculations into functions (and 
consider organizing it into a package, if it's big enough); put those 
functions in separate files from the scripts that run them on particular 
examples.  This way you can run source("fns.R") to update the function 
definitions without re-doing all the simulations in your project.

Duncan Murdoch

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


Re: [R] Synchronzing workspaces

2007-08-22 Thread James W. MacDonald
You could use a versioning system like Subversion or Git on the server 
but work with the local repository. I believe Git is the easier of the 
two to set up.

Best,

Jim



Eric Turkheimer wrote:
> How do people go about synchronizing multiple workspaces on different
> workstations?  I tend to wind up with projects spread around the various
> machines I work on.  I find that placing the directories on a server and
> reading them remotely tends to slow things down.
> 
> thanks,
> Eric
> 

-- 
James W. MacDonald, M.S.
Biostatistician
Affymetrix and cDNA Microarray Core
University of Michigan Cancer Center
1500 E. Medical Center Drive
7410 CCGC
Ann Arbor MI 48109
734-647-5623

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Synchronzing workspaces

2007-08-22 Thread Eric Turkheimer
How do people go about synchronizing multiple workspaces on different
workstations?  I tend to wind up with projects spread around the various
machines I work on.  I find that placing the directories on a server and
reading them remotely tends to slow things down.

thanks,
Eric

-- 
Eric Turkheimer, PhD
Department of Psychology
University of Virginia
PO Box 400400
Charlottesville, VA  22904-4400

434-982-4732
434-982-4766 (FAX)

[[alternative HTML version deleted]]

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


Re: [R] Processing Sweave documents without Sweave

2007-08-22 Thread Duncan Murdoch
On 8/22/2007 6:40 AM, Renaud Lancelot wrote:
> Sweave produces a TeX file which can be exchanged like any other such
> file, providing that you give the possible graphic files (and bib,
> etc.) called from the TeX code.

Yes, and since Sweave doesn't touch the formatting of the TeX code, it's 
not hard to recognize and merge changes to the .tex file back into the 
Rnw file.

Duncan Murdoch
> 
> Best,
> 
> Renaud
> 
> 2007/8/22, Werner Wernersen <[EMAIL PROTECTED]>:
>> Hi,
>>
>> I am very intrigued by the idea of integrating
>> statistical analysis directly with a paper as Sweave
>> does it. But as I am collaborating with several people
>> and they won't have set up R and also they're very
>> unlikely interested in learning it, I am concerned
>> about how such an Latex/Sweave document will work on
>> their side.
>>
>> Can they Latex-compile the document without having R /
>> Sweave installed and still see the numbers and figures
>> I produced by an earlier run of Latex / R / Sweave? Is
>> there any option or so planned for this case? Or how
>> are you coping with this problem otherwise?
>>
>> Many thanks,
>>   Werner
>>
>> __
>> R-help@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Formatting Sweave in R-News

2007-08-22 Thread Duncan Murdoch
On 8/22/2007 3:20 AM, Gregor Gorjanc wrote:
> Arjun Ravi Narayan  arjunnarayan.com> writes:
>> I am editing a document for submission to the R-news newsletter, and
>> in my article my Sweave code inserts a dynamically generated PDF
>> report that my R program generates.
>> 
> 
> Slightly off Arjuns topic is a problem with Schunk, Sinput and Soutput 
> environments. They just use to much inline space. I have tweade this
> problem with sed (see bellow for Makefile content), but wonder if there 
> is a better solution.

You can change the spacing by redefining those environments.  I use this:

% This removes the extra spacing after code and output chunks in Sweave,
% but keeps the spacing around the whole block.

\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}


> Thanks, Gregor
> 
> default: Sweave fixTex
>   texi2pdf --clean wrapper.tex && evince wrapper.pdf
> 
> Sweave: # Sweave myPaper.Rnw
>   R CMD Sweave myPaper.Rnw
> 
> fixTex: Sweave # Change all S* environments to smallverbatim
>   @cat myPaper.tex | sed -e '/\\begin{Sinput}/d' \
>  -e '/\\end{Sinput}/d' \
>  -e '/\\begin{Soutput}/d' \
>  -e '/\\end{Soutput}/d' \
>  -e 's/\\begin{Schunk}/\\begin{smallverbatim}/g' \
>  -e 's/\\end{Schunk}/\\end{smallverbatim}/g' \
>   > myPaper2.tex
>   @mv -f myPaper2.tex myPaper.tex
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] C code generators

2007-08-22 Thread Duncan Murdoch
On 8/21/2007 10:20 PM, [EMAIL PROTECTED] wrote:
> Dear R-helpers
> 
> Are there any established R packages that include a C code generator -- 
> that generates new C language files and compiles them?

Oleg Sklyar's inline package does that.  It takes input in the form of a 
fragment of C (or C++, Fortran, etc.), completes the source code, 
compiles it, loads it, and creates an R function to call it.  See the 
?cfunction man page for examples.

I wouldn't recommend this for extensive use:  you're better off putting 
your own code in a package instead.  But it's good for quick tests, or 
writing the first version of the code that goes in the package.

> To be precise what I'm looking for is a process that takes text input in
> some format (it might be pseudocode, fragments of C code, etc) and creates
> a valid C language source file that can be compiled by R CMD COMPILE.
> Ideally the procedure should also cause the C code to be compiled and
> dynamically loaded.
> 
> To give a trivial example, suppose I want to be able to perform image
> filtering operations on matrices. All filters have the same structure:
> each position [i,j] in the matrix is visited in a double loop; a
> calculation (depending on the filter) is performed using the value at
> [i,j] and also the values at neighbouring positions [i+1,j] , [i+1,j+1]
> etc; the result is written to the output matrix at the position [i,j]. The
> C code for the loop is always the same; only a few lines of code that
> perform the filter calculation will change. I would like a procedure that
> accepts a text file containing just a few lines of C code or pseudocode,
> and inserts these lines into the appropriate place in the loop, producing
> a valid C language routine which can then be compiled by R CMD COMPILE and
> dynamically loaded.

As a matter of fact, one of Oleg's examples is an image filter.

Duncan Murdoch

> 
> (Of course, it can get more complicated than just inserting a single text
> fragment!)
> 
> I once implemented such a feature in an image processing package, so I
> know it's not hard. Before dusting off this ancient code I would like to
> learn whether there's an R package or other open source program that
> already does it.
> 
> Any pointers would be welcome
> 
> thanks
> Adrian Baddeley
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] IRT model nonparametric from Ramsay in R

2007-08-22 Thread Xavier Giovanni Ordóñez Camacho
Hi:
I need to apply the IRT model from Ramsay (1991), He apply the Smoothing Kernel 
to multiple choice test. Is possible in R?.

Thank you,

Xavier G. Ordonez
Doctoral Student
Universidad Complutense de Madrid

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


Re: [R] Optimization problem

2007-08-22 Thread Gabor Grothendieck
Try this.

1. following Ben remove the Randalstown point and reset the levels of the
Location factor.

2. then replace solve with ginv so it uses the generalized inverse to calculate
the hessian:

alan2 <- subset(alan, subset = Location != "Randalstown")
alan2$Location <- factor(as.character(alan2$Location))

library(MASS)
solve <- ginv

zinb.zc <- zicounts(resp=Scars~.,x =~Location + Lar + Mass + Lar:Mass
+ Location:Mass,z =~Location + Lar + Mass + Lar:Mass + Location:Mass,
data = alan2)

rm(solve)

On 8/21/07, Ben Bolker <[EMAIL PROTECTED]> wrote:
>
>  (Hope this gets threaded properly.  Sorry if it doesn't.)
>
>   Gabor: Lac and Lacfac being the same is irrelevant, wouldn't
> produce NAs (but would produce something like a singular Hessian
> and maybe other problems) -- but they're not even specified in this
> model.
>
>  The bottom line is that you have a location with a single
> observation, so the GLM that zicounts runs to get the initial
> parameter values has an unestimable location:mass interaction
> for one location, so it gives an NA, so optim complains.
>
>  In gruesome detail:
>
> ## set up  data
> scardat = read.table("scars.dat",header=TRUE)
> library(zicounts)
> ## try to run model
> zinb.zc <- zicounts(resp=Scars~.,
>x =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>z =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>data=scardat)
> ## tried to debug this by dumping zicounts.R to a file, modifying
> ## it to put a "trace" argument in that would print out the parameters
> ## and log-likelihood for every call to the log-likelihood function.
> dump("zicounts",file="zicounts.R")
> source("zicounts.R")
> zinb.zc <- zicounts(resp=Scars~.,
>x =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>z =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>data=scardat,trace=TRUE)
> ## this actually didn't do any good because the negative log-likelihood
> ## function never gets called -- as it turns out optim() barfs when it
> ## gets its initial values, before it ever gets to evaluating the
> log-likelihood
>
> ## check the glm -- this is the equivalent of what zicounts does to
> ## get the initial values of the x parameters
> p1 <- glm(Scars~Location + Lar + Mass + Lar:Mass + Location:Mass,
>  data=scardat,family="poisson")
> which(is.na(coef(p1)))
>
> ## find out what the deal is
> table(scardat$Location)
>
> scar2 = subset(scardat,Location!="Randalstown")
> ## first step to removing the bad point from the data set -- but ...
> table(scar2$Location)
> ## it leaves the Location factor with the same levels, so
> ##  now we have ZERO counts for one location:
> ## redefine the factor to drop unused levels
> scar2$Location <- factor(scar2$Location)
> ## OK, looks fine now
> table(scar2$Location)
>
> zinb.zc <- zicounts(resp=Scars~.,
>x =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>z =~Location + Lar + Mass + Lar:Mass + Location:Mass,
>data=scar2)
> ## now we get another error ("system is computationally singular" when
> ## trying to compute Hessian -- overparameterized?)   Not in any
> ## trivial way that I can see.  It would be nice to get into the guts
> ## of zicounts and stop it from trying to invert the Hessian, which is
> ## I think where this happens.
>
>  In the meanwhile, I have some other  ideas about this analysis (sorry,
> but you started it ...)
>
>  Looking at the data in a few different ways:
>
> library(lattice)
> xyplot(Scars~Mass,groups=Location,data=scar2,jitter=TRUE,
>   auto.key=list(columns=3))
> xyplot(Scars~Mass|Location,data=scar2,jitter=TRUE)
>
> xyplot(Scars~Lar,groups=Location,data=scar2,
>   auto.key=list(columns=3))
> xyplot(Scars~Mass|Lar,data=scar2)
> xyplot(Scars~Lar|Location,data=scar2)
>
>   Some thoughts: (1) I'm not at all sure that
> zero-inflation is necessary (see Warton 2005, Environmentrics).
> This is a fairly small, noisy data set without huge numbers
> of zeros -- a plain old negative binomial might be fine.
>
>   I don't actually see a lot of signal here, period (although there may
> be some) ...
> there's not a huge range in Lar (whatever it is -- the rest of the
> covariates I
> think I can interpret).  It would be tempting to try to fit location as
> a random
> effect, because fitting all those extra degrees of freedom is going to
> kill you.
> On the other hand, GLMMs are a bit hairy.
>
>   cheers
>   Ben
>
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Ans on: How do i print a "main title" on a win.graph with several plots?

2007-08-22 Thread Tom Willems
Thanks Mr.Pr. Ripley,
 
For pointing me in the good direction,
i use win.graph() so i get an overview of the different graphs.
Until now, i v had no problems with it, hope it stays that way.


for those whom have the same porblem, this is what i do:

 win.graph();   op <- par(mfrow = c(1,2)  , oma= c(1,1,1,1))
 boxplot(lg_value~labo, main="Test 1 at day 1",ylab="log(x) ",
 xlab="different lab's", data=dataset,ylim=c(-0.05,5))
 mtext("At day 1",line=-1,outer=TRUE)
 boxplot(lg_value~labo, main="Test 2 at day 1",ylab="log(x)",
 xlab="different lab's", data=dataset,ylim=c(-0.05,5))
 par(op);

Kind regards,
Tom W.


Disclaimer: click here
[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] ordering of variables in plots

2007-08-22 Thread elin sagulin
I'm making plots where the x-axis variables are called l hl s (for
"long" "half-long" and "short").  I'd like them to appear in that
order, but they turn up in alphabetical order: hl l s.
Except for fiddling with the names of the variables, how can I fix
this? That is, how can I specify what order I want?

Thanks,

Elin

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


Re: [R] within-subject factors in lme

2007-08-22 Thread lorenz.gygax
I don't think, this has been answered:

> I'm trying to run a 3-way within-subject anova in lme with 3
> fixed factors (Trust, Sex, and Freq), but get stuck with handling
> the  random effects. As I want to include all the possible random
> effects in the model, it would be something more or less
> equivalent to using aov
> 
>  > fit.aov <- aov(Beta ~ 
> Trust*Sex*Freq+Error(Subj/(Trust*Sex*Freq)),  
> Model)
> 
> However I'm not so sure what I should do in lme. Sure
> 
>  > lme(Beta ~ Trust*Sex*Freq, random = ~1|Subj, Model)
> 
> works fine, but it only models the random effect of the
> intercept. I tried the following 4 possibilities:

If I understand correctly, you want to include the interactions between the 
random and fixed terms? This is done like:

model.lme <- lme(Beta ~ Trust*Sex*Freq,
 random = ~Trust*Sex*Freq|Subj, Model)

But this needs a lot of observations as quite a few parameters need to be 
estimated! Possibly, you can not include the variable Sex in this, because I 
assume that Subj is nested within Sex. If you just refer to within and between 
subject effects and their corresponding degrees of freedom: you should see this 
being handled automatically and correctly by lme e.g. in the output of anova 
(model.lme)

Lorenz
- 
Lorenz Gygax
Centre for proper housing of ruminants and pigs
Agroscope Reckenholz-Tänikon Research Station ART
Tänikon, CH-8356 Ettenhausen / Switzerland

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


Re: [R] Processing Sweave documents without Sweave

2007-08-22 Thread Werner Wernersen
That sound's perfect, that's all I needed to know.

Thanks for your kind help,
  Werner

--- Renaud Lancelot <[EMAIL PROTECTED]>
schrieb:

> Sweave produces a TeX file which can be exchanged
> like any other such
> file, providing that you give the possible graphic
> files (and bib,
> etc.) called from the TeX code.
> 
> Best,
> 
> Renaud
> 
> 2007/8/22, Werner Wernersen
> <[EMAIL PROTECTED]>:
> > Hi,
> >
> > I am very intrigued by the idea of integrating
> > statistical analysis directly with a paper as
> Sweave
> > does it. But as I am collaborating with several
> people
> > and they won't have set up R and also they're very
> > unlikely interested in learning it, I am concerned
> > about how such an Latex/Sweave document will work
> on
> > their side.
> >
> > Can they Latex-compile the document without having
> R /
> > Sweave installed and still see the numbers and
> figures
> > I produced by an earlier run of Latex / R /
> Sweave? Is
> > there any option or so planned for this case? Or
> how
> > are you coping with this problem otherwise?
> >
> > Many thanks,
> >   Werner
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained,
> reproducible code.
> >
> 
> 
> -- 
> Renaud LANCELOT
> Département Systèmes Biologiques du CIRAD
> CIRAD, Biological Systems Department
> 
> Campus International de Baillarguet
> TA 30 / B
> F34398 Montpellier
> Tel   +33 (0)4 67 59 37 17
> Secr. +33 (0)4 67 59 37 37
> Fax   +33 (0)4 67 59 37 95
> 



  

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


Re: [R] Processing Sweave documents without Sweave

2007-08-22 Thread Renaud Lancelot
Sweave produces a TeX file which can be exchanged like any other such
file, providing that you give the possible graphic files (and bib,
etc.) called from the TeX code.

Best,

Renaud

2007/8/22, Werner Wernersen <[EMAIL PROTECTED]>:
> Hi,
>
> I am very intrigued by the idea of integrating
> statistical analysis directly with a paper as Sweave
> does it. But as I am collaborating with several people
> and they won't have set up R and also they're very
> unlikely interested in learning it, I am concerned
> about how such an Latex/Sweave document will work on
> their side.
>
> Can they Latex-compile the document without having R /
> Sweave installed and still see the numbers and figures
> I produced by an earlier run of Latex / R / Sweave? Is
> there any option or so planned for this case? Or how
> are you coping with this problem otherwise?
>
> Many thanks,
>   Werner
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Renaud LANCELOT
Département Systèmes Biologiques du CIRAD
CIRAD, Biological Systems Department

Campus International de Baillarguet
TA 30 / B
F34398 Montpellier
Tel   +33 (0)4 67 59 37 17
Secr. +33 (0)4 67 59 37 37
Fax   +33 (0)4 67 59 37 95

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 do i print a "main title" on a win.graph with several plots?

2007-08-22 Thread Prof Brian Ripley
?title, look at the 'outer' argument.

You can see further discussion of the outer margins in 'An Introduction to 
R'.

I don't know why you are using win.graph(): it is a deprecated form of 
windows() with many of the arguments taking unchangable defaults.

On Wed, 22 Aug 2007, Tom Willems wrote:

> Good Mornig All,
>
> How R you today? ;o)
>
> I have lots of questions, but i l start with the simplest one,
> to wich i am shy to say, i did not find the answer.
>
> It is the following:
>
> When i make a summary plot like for example plot( summary(glm)),
> i get one window, one main title, and 4 graph's in that window.
>
> Now i do know how to get several graphs in one window,
> buit i don't manage geting one main title , in the top middel of the
> window.
> I can give every plot a different main and subtitle, but i can't put no
> title in the "win.graph()" box.
>
> this is what i do:
>
> win.graph();   op <- par(mfrow = c(1,2))
> boxplot(lg_value~labo, main="Test 1 at day 1",ylab="log(x) ",
> xlab="different lab's", data=dataset,ylim=c(-0.05,5))
> boxplot(lg_value~labo, main="Test 2 at day 1",ylab="log(x)",
> xlab="different lab's", data=dataset,ylim=c(-0.05,5))
> par(op);
>
> Then i get one graph window, two plots each with their onw main title.
> What i'd like to have is, one main title saying " At day 1", and then two
> plots with the different tests.
> How can i do this, pls?
>
> Kind regards,
> Tom W.
>
>
> Disclaimer: click here
>   [[alternative HTML version deleted]]
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Processing Sweave documents without Sweave

2007-08-22 Thread Werner Wernersen
Hi,

I am very intrigued by the idea of integrating
statistical analysis directly with a paper as Sweave
does it. But as I am collaborating with several people
and they won't have set up R and also they're very
unlikely interested in learning it, I am concerned
about how such an Latex/Sweave document will work on
their side. 

Can they Latex-compile the document without having R /
Sweave installed and still see the numbers and figures
I produced by an earlier run of Latex / R / Sweave? Is
there any option or so planned for this case? Or how
are you coping with this problem otherwise?

Many thanks,
  Werner

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


Re: [R] Recurrence plots for CGH

2007-08-22 Thread Ido M. Tamir
>CGH

Hi,

I guess you should start looking at the bioconductor project.

http://www.bioconductor.org/packages/release/Software.html

best wishes,

ido

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


Re: [R] divided scatter plots

2007-08-22 Thread Jim Lemon
Caroline Nganga wrote:
> I have a data set which contains two columns. The first column is a
> list of countries, and the second column contains their political risk
> ratings. I would like to create  one large plot that contains 5
> different sections, each with a scatter plot. To clarify, I have
> divided the countries into 5 groups. For each group (continent), I
> would like to have the name of the continent on the x-axis, and points
> representing countries and  their risk rating on the y-axis. However,
> I want all 5 scatter plots to be in one large plot. What function
> should I use to do this? Also, is it possible to label each point?
> thanks for any help!
> 
Hi Caroline,
If I understand your request, you might be able to use the axis.break 
function in the plotrix package. That is, you make one big plot with the 
points in five columns and then put gap style axis breaks between the 
columns. Here's a toy example:

library(plotrix)
prr.df<-data.frame(country=c("us","mx","ca","br","ar","pe",
  "ch","mn","in","nl","fr","es","na","mz","rw"),
  continent=c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5),
  prr=rnorm(15)+4)
plot(prr.df$continent,prr.df$prr,main="Political risk ratings",
  xlim=c(0.7,5.3),xlab="Continent",ylab="Risk rating",type="n")
text(prr.df$continent,prr.df$prr,prr.df$country)
axis.break(1,1.5,style="gap")
axis.break(1,2.5,style="gap")
axis.break(1,3.5,style="gap")
axis.break(1,4.5,style="gap")

Jim

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


Re: [R] plot in for loop with "i" indexed variables does not work

2007-08-22 Thread Dimitris Rizopoulos
you need something like this:

par(mfrow = c(5, 5))
for (i in 1:10) {
nam <- paste("Var.", i, sep = "")
plot(x = time, y = mat[, nam], xlab = "Time",
ylab = nam)
}

where `mat' is the matrix containing Var.1, Var.2, Var.3, etc.


I hope it helps.

Best,
Dimitris


Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
 http://www.student.kuleuven.be/~m0390867/dimitris.htm


- Original Message - 
From: <[EMAIL PROTECTED]>
To: 
Sent: Wednesday, August 22, 2007 10:56 AM
Subject: [R] plot in for loop with "i" indexed variables does not work


> Hi there
>
> Does anyone know, why this is not working?:
>
>
> par(mfrow = c(5, 5))
>
> for (i in 1:10){
> plot(x=time, y=Var.i)
> }
>
> (x-variable is time
> y-variable is Var out of matrix with columns Var.1 Var.2 Var.3 
> etc...)
>
> even
>
> for (i in 1:10){
> plot(x=time, y=paste("Var", i, sep=".")
> }
>
> does not work
>
> Thanks
>
> marc
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] plot in for loop with "i" indexed variables does not work

2007-08-22 Thread mgilgen
Hi there

Does anyone know, why this is not working?:


par(mfrow = c(5, 5))

for (i in 1:10){
plot(x=time, y=Var.i)
}

(x-variable is time
y-variable is Var out of matrix with columns Var.1 Var.2 Var.3 etc...)

even

for (i in 1:10){
plot(x=time, y=paste("Var", i, sep=".")
}

does not work

Thanks 

marc

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


Re: [R] library(fCalendar) timeDate("12.03.2005", format="%d.%m.%Y")

2007-08-22 Thread Martin Becker
Martin Maechler wrote:
>> "OL" == Ola Lindqvist <[EMAIL PROTECTED]>
>> on Tue, 21 Aug 2007 14:32:19 +0200 writes:
>> 
>
> OL> Thanks!
> OL> Seems to work fine now!
>
> Well, for your example.
>
> But sorry to say, the patch breaks other cases.
>
> I'm investigating further
> (and will hopefully contribute to a new CRAN release of fCalendar
>  once Diethelm Wuertz is back from wherever; I've already made
>  more changes)
>
>   

Good to see experts investigating :-)

Sorry, I should have mentioned, at least, that I am not an Rmetrics 
developer, that I only took a quick glance at the code (focused on the 
reported problem), and that the patched version is neither official nor 
tested. Now I see that I completely messed up the case for "Date" input 
and sometimes .midnightStandard() was not working correct any more.

One further issue, which was brought to my attention by another R-user 
off list: I think, that all occurences of

  if (sum(lt$sec+lt$min+lt$hour) == 0) isoFormat = "%Y-%m-%d"

in R/3A-TimeDateClass.R could safely be replaced by

  if (sum(lt$sec+lt$min+lt$hour,na.rm=TRUE) == 0) isoFormat = "%Y-%m-%d"

in order to make fCalendar more NA-friendly. Currently, at least for 
some cases, a single NA entry in the input vector suffices to let 
timeDate() break down. Again, I am not sure, if this change breaks some 
special cases, or if breaking down in case of NAs is intended.

Nevertheless, I have built a new, unofficial and untested patched 
(windows-)version of fCalendar (same weblink), which makes the "Date" 
input case and .midnightStandard() working again (I hope so, at least), 
and incorporates  the (improved?) NA-handling. But of course: only use 
this version with care (and only if the official version doesn't work 
for your particular problem). Hopefully, we see an official update on 
CRAN soon :-)

Kind regards,

  Martin Becker


> Martin Maechler,
> ETH Zurich  [but different department than D.Wuertz]
>
>
> OL> Martin Becker wrote:
> >> Dear Ola,
> >> 
> >> I think you spotted a small bug in *package* fCalendar.
> >> Explicit specification should prevent "autodetection" of the date 
> >> format, which is not the case for fCalendar v251.70, instead 
> >> autodetection is done at least once (twice, if actually appropriate).
> >> With the following patch, things should work ok:
> >> 
> >> diff --recursive fCalendar.orig/R/3A-TimeDateClass.R 
> >> fCalendar/R/3A-TimeDateClass.R
> >> 433c433
> >> < charvec = format(strptime(charvec, .whichFormat(charvec)), 
> >> isoFormat)
> >> ---
> >> > charvec = format(strptime(charvec, format), isoFormat)
> >> 
> >> You did not provide the output of sessionInfo() (which you are asked 
> >> for in the posting guide). If you are using Windows and don't know how 
> >> to apply the patch, you can download a patched binary version here:
> >> http://www.saar-gate.net/download/fCalendar_251.70.zip
> >> 
> >> Regards,
> >> 
> >> Martin
> >> 
> >> PS: Maybe r-sig-finance is more appropriate for questions concerning 
> >> Rmetrics.
> >> 
> >> 
> >> Ola Lindqvist wrote:
> >>> Dear R users,
> >>> I have problem with the library fCalendar.
> >>> 
> >>> I am not using the US standard format notations. It seems like it is 
> >>> not possible to have different format than the US standards.
> >>> Anyone how knows a way to go around this problem?
> >>> 
> >>> Here is the code I enter:
> >>> myDate = "12.03.2005"
> >>> timeDate(myDate, format = "%d.%m.%Y")
> >>> 
> >>> And I get following error message:
> >>> Error in if (sum(lt$sec + lt$min + lt$hour) == 0) isoFormat = 
> >>> "%Y-%m-%d" :
> >>> missing value where TRUE/FALSE needed
> >>> 
> >>> Thanks,
> >>> Ola
> >>> 
> >>> __
> >>> R-help@stat.math.ethz.ch mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide 
> >>> http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible code.
> >>> 
> >> 
>
> OL> __
> OL> R-help@stat.math.ethz.ch mailing list
> OL> https://stat.ethz.ch/mailman/listinfo/r-help
> OL> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> OL> and provide commented, minimal, self-contained, reproducible code.
>

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


[R] How do i print a "main title" on a win.graph with several plots?

2007-08-22 Thread Tom Willems
Good Mornig All,

How R you today? ;o)

I have lots of questions, but i l start with the simplest one,
to wich i am shy to say, i did not find the answer.

It is the following: 

When i make a summary plot like for example plot( summary(glm)),
i get one window, one main title, and 4 graph's in that window.

Now i do know how to get several graphs in one window,
buit i don't manage geting one main title , in the top middel of the 
window.
I can give every plot a different main and subtitle, but i can't put no 
title in the "win.graph()" box.

this is what i do:

win.graph();   op <- par(mfrow = c(1,2))
boxplot(lg_value~labo, main="Test 1 at day 1",ylab="log(x) ", 
xlab="different lab's", data=dataset,ylim=c(-0.05,5))
boxplot(lg_value~labo, main="Test 2 at day 1",ylab="log(x)", 
xlab="different lab's", data=dataset,ylim=c(-0.05,5))
 par(op);

Then i get one graph window, two plots each with their onw main title.
What i'd like to have is, one main title saying " At day 1", and then two 
plots with the different tests.
How can i do this, pls?

Kind regards,
Tom W.


Disclaimer: click here
[[alternative HTML version deleted]]

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


Re: [R] Random Sampling from a Matrix

2007-08-22 Thread Uwe Ligges


Anup Nandialath wrote:
> Dear Friends,
> 
> I have a matrix of size 5000 X 20. The first two columns are indicator 
> variables taking the value of either 0 or 1. Let us call the first two 
> columns Y1 and Y2. 
> 
> I need to randomly sample 1000 rows with all the associated columns, in other 
> words my new matrix should be of size 1000 X 20. I realize that using this 
> command
> 
>  newmat <- mainmat[sample(1000,replace=F),]
> 
> achieves this. However, I would like to make sure that both Y1 and Y2 have 
> more or less an equal amount of 0's and 1's. At present when I sample, I get 
> cases where sometimes all my Y2's are 0. Is there any way to accomodate this 
> problem.


Yes, for example drawing samples stratified by the values of Y1 or Y2.

Uwe Ligges



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


Re: [R] rectify a program of seasonal dummies matrix

2007-08-22 Thread Uwe Ligges


Friedrich Schuster wrote:
> Hello, 
>  
> the main problem seems to be the "if else", should be "else if". 
> 
> Your code is hard to read, maybe you should consider using more () {}: 
> 
> T <- 100;
> br <- matrix(0,T,4);

Thanks for the contribution. Please note:
a) It is a bad idea to have a variable called T. Some people still use 
it as a logical value even if they should not.
b) R does not need any ";" at the end of a line.

Uwe Ligges


> for (i in 1:T) {
>for (j in 1:4) {
>   if (i==j) {
>  br[i,j] <- 1;
>}
>else if ((abs(i-j)%%4)==0) {
>   br[i,j] <- 1;
> }
> else {
>   br[i,j] <- 0;
> }
> }
>  }
> 
> A simpler approach is creating a diagonal matrix and multply it : 
> 
> # create small diagonal matrix
> mat = diag(x=1, nrow=4, ncol=4);
> mat
> # multiply diagonal matrix and re-dimension it   to 4 cols
> br <- rep(mat, 25);
> dim(br) <- c(100, 4);
> br;
> 
> Hope this helps, 
> FS
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] Odp: displaying the means on a scatter diagram

2007-08-22 Thread Petr PIKAL
Hi

[EMAIL PROTECTED] napsal dne 22.08.2007 09:10:26:

> 
> Hello,
> 
> I created a simple scatter diagram:
> 
> x= c(53, 52, 55, 65, 71, 69, 75, 78, 88, 70)
> y= c(162, 165, 165, 171, 173, 175, 179, 181, 184, 176)
> plot(x, y)
> 
> Now I would like to display the mean on that diagram.

do you think

points(mean(x), mean(y))

or e.g.

abline(h=mean(y))
abline(v=mean(x))

Regards
Petr


> 
> Can anyone tell me what possibilities I have to do that?
> Thanks in advance
> 
> Tobias
> 
> -- 
> View this message in context: 
http://www.nabble.com/displaying-the-means-on-a-
> scatter-diagram-tf4309832.html#a12269270
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Formatting Sweave in R-News

2007-08-22 Thread Gregor Gorjanc
Arjun Ravi Narayan  arjunnarayan.com> writes:
> I am editing a document for submission to the R-news newsletter, and
> in my article my Sweave code inserts a dynamically generated PDF
> report that my R program generates.
> 

Slightly off Arjuns topic is a problem with Schunk, Sinput and Soutput 
environments. They just use to much inline space. I have tweade this
problem with sed (see bellow for Makefile content), but wonder if there 
is a better solution.

Thanks, Gregor

default: Sweave fixTex
texi2pdf --clean wrapper.tex && evince wrapper.pdf

Sweave: # Sweave myPaper.Rnw
R CMD Sweave myPaper.Rnw

fixTex: Sweave # Change all S* environments to smallverbatim
@cat myPaper.tex | sed -e '/\\begin{Sinput}/d' \
   -e '/\\end{Sinput}/d' \
   -e '/\\begin{Soutput}/d' \
   -e '/\\end{Soutput}/d' \
   -e 's/\\begin{Schunk}/\\begin{smallverbatim}/g' \
   -e 's/\\end{Schunk}/\\end{smallverbatim}/g' \
> myPaper2.tex
@mv -f myPaper2.tex myPaper.tex

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/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] displaying the means on a scatter diagram

2007-08-22 Thread squall44

Hello,

I created a simple scatter diagram:

x= c(53, 52, 55, 65, 71, 69, 75, 78, 88, 70)
y= c(162, 165, 165, 171, 173, 175, 179, 181, 184, 176)
plot(x, y)

Now I would like to display the mean on that diagram.

Can anyone tell me what possibilities I have to do that?
Thanks in advance

Tobias

-- 
View this message in context: 
http://www.nabble.com/displaying-the-means-on-a-scatter-diagram-tf4309832.html#a12269270
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Time conversion problems

2007-08-22 Thread Peter Dalgaard
[EMAIL PROTECTED] wrote:
> Hi there
>
> I have precipitation data from 2004 to 2006 in varying resolutions (10 to 
> 20min intervals) with time in seconds from beginnig of the year (summation) 
> and a second variable as year.
>
> I applied follwing code to convert the time into a date:
>
> times<-strptime("2004-01-01", "%Y-%m-%d", tz="GMT") + precipitation$time1
>
> everytihng went well, except that every year, the seconds-counter starts by 
> zero, therefore I have now three 2004 series instead of going further from 04 
> to 05 etc.
>
> I tried to sum the last seconds-values of 2004 to the first of 2005 with an 
> if command like:
>
> if (year=2005) time2=time1+632489 ;(seconds)
>
> but it doesn't work.
>
> thanks for a solution
>   
Can't you just do strptime(paste(year, "01-01", sep="-"), ? (or use 
ISOdatetime(year,1,1,0,0,0,tz="GMT"))

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


[R] [R-pkgs] brew 1.0-1

2007-08-22 Thread Jeffrey Horner
brew implements a templating framework for mixing text and R code for 
report generation. brew template syntax is similar to PHP, Ruby's erb 
module, Java Server Pages, and Python's psp module.

brew is written in R with no package dependencies, and it's not just for 
the web. It can be used as an alternative to Sweave in a limited 
context. See the brew-test-1.brew file in the distribution for some 
salient differences between the two. brew can also complement Sweave 
since it can be written to do conditional inclusion of or loop over 
Sweave code chunks.

The 1.0-1 version should show up on the CRAN mirrors shortly, but in the 
mean time it can be got from:

http://www.rforge.net/brew/

Best,

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

___
R-packages mailing list
[EMAIL PROTECTED]
https://stat.ethz.ch/mailman/listinfo/r-packages

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