Re: [R] searching for data.frame rows / processing of rows

2006-10-05 Thread Gabor Grothendieck
Grouping the data frame by the first two columns, apply colMeans
and then rbind the resulting by-structure together:

do.call(rbind, by(DF, DF[2:1], colMeans, na.rm = TRUE))


On 10/5/06, Greg Tarpinian [EMAIL PROTECTED] wrote:
 R 2.3.1, WinXP:

 I have a puzzling problem that I suspect may be solved using
 grep or a regular expression but am at a loss how to actually do it...
 My data.frame looks like

  Location   TimeXY
        ---  ---
  1  0  1.6  9.3
  1  3  4.2  10.4
  1  6  2.7  16.3
  2  0  0.5  2.1
  2  3  NA   3.6
  2  3  5.0  0.06
  2  6  3.4  14.0

 and so forth.  I would like to search for duplicate Time values
 within a Location and take the numerical average (where possible)
 of the elements in X and Y.  These numerical averages should
 then be used to create a single row where multiple rows once
 existed.  So, I would like to obtain

  2  3  5.0  1.83

 for the two Time = 3 rows for Location = 2 and use it to replace
 these two rows.  Ideally, avoiding for(i in 1:blah) loops would be
 nice because the data.frame has about 10,000 rows that need to
 be searched and processed.  My intent is to do some comparing of
 SAS to R -- the DATA step processing in SAS is quite fast and
 using the RETAIN statement along with the LAG( ) function allows
 this sort of thing to be done rapidly.


 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.


__
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 script

2006-10-05 Thread Gabor Grothendieck
On 10/4/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
 Also see package caTools.

 On 10/4/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
  See:
 
  http://tolstoy.newcastle.edu.au/R/help/04/10/5161.html
 
  On 10/4/06, JOHN VOIKLIS [EMAIL PROTECTED] wrote:
   Hello,
  
   I wrote the function, below, in the hope of _quickly_ generating a
   sliding-window time series of the mean, sd, median, and mad values
   from a matrix of data. The script below is anything but quick; it has
   been working away on a 2500 x 2500 matrix with a sliding window of 100
   x 100 for five days, with no end in sight. Obviously, I am not the
   best of programmers. Can anyone offer suggestions as to how I might
   optimize this script.
  
   Thank you,
  
   John
  
   common.ground-function(inMatrix,window){
  
  cleanMatrix-as.matrix(inMatrix)
  diag(cleanMatrix)-NA
  outMatrix-matrix(0,dim(inMatrix)[1]-window,4)
  colnames(outMatrix)-c(mean,SD, median,mad)
  
  for(i in 1:dim(outMatrix)[1]){
  for(j in window:dim(cleanMatrix)[2]){
  
   outMatrix[i,1]-mean(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE)
  
   outMatrix[i,2]-sd(c(cleanMatrix[c(i:j),c(i:j)]),na.rm=TRUE)
  
   outMatrix[i,3]-median(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE)
  
   outMatrix[i,4]-mad(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE)
  }
  }
  return(outMatrix)
   }

Also you could look at rollmax, rollmean and rollmedian in the zoo package.
rollapply in zoo with FUN = sd could be used for the sd.

__
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] A statement over multiple lines (i.e. the ... feature in Matlab)

2006-10-05 Thread Wee-Jin Goh
Hello again list,

I thought I'd start a new thread, seeing as it's completely different  
from my previous question. Some functions I have written require many  
parameters, and so do not fit nicely into an 80 column width display.  
This is usually avoided, by spreading that particular statement over  
a few lines. This is something that I do in Matlab with the following:

myFunc( parameter1, ...
parameter2, ...
parameter3, ...
parameter4)

The  ... operator in Matlab allows me to spread a statement over  
several lines. The ... operator in R seems to be more like the ...  
operator in C, which allows for a variable argument list.

How do I accomplish this task in R? It's not a show stopper, but it  
would make reading my code much much easier.

Cheers,
Wee-Jin

__
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] Variables in RODBC environment

2006-10-05 Thread Thorsten Muehge

Hello Experts,
how can I use variables in the RODBC environment.

Example which does not work:

Thanks for your help.

Thorsten

pn -  '39R5238';

library(RODBC);
odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx);
sql - select
u.unitid,
from test
where part in ('pn')
;
parameter - sqlQuery(odbcobj,sql);
odbcClose(odbcobj);

__
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] A statement over multiple lines (i.e. the ... feature in Matlab)

2006-10-05 Thread Robin Hankin
Hi  Wee-Jin

you can block out bits of R code with

if(FALSE){
  code not executed
}

For the line breaking, R deals with incomplete lines by not executing  
the statement
until you finish it.  In the function case, it waits for you to close  
a bracket.


If you type:

myFunc(a=3,
b=5,
c=6,
d=7
)

myFunc() will only execute
when you close the bracket


HTH

rksh


On 5 Oct 2006, at 08:28, Wee-Jin Goh wrote:

 Hello again list,

 I thought I'd start a new thread, seeing as it's completely different
 from my previous question. Some functions I have written require many
 parameters, and so do not fit nicely into an 80 column width display.
 This is usually avoided, by spreading that particular statement over
 a few lines. This is something that I do in Matlab with the following:

 myFunc( parameter1, ...
 parameter2, ...
 parameter3, ...
 parameter4)

 The  ... operator in Matlab allows me to spread a statement over
 several lines. The ... operator in R seems to be more like the ...
 operator in C, which allows for a variable argument list.

 How do I accomplish this task in R? It's not a show stopper, but it
 would make reading my code much much easier.

 Cheers,
 Wee-Jin

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

--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

__
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] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX

2006-10-05 Thread Stefan Grosse

 pdf(sgr6100.pdf, horizontal=FALSE, onefile=FALSE, 
 height=3, width=3, pointsize=6)

 Reducing point-size below 6 does not seem to make any difference to 
 the size of text and symbols. Any suggestions to get smaller font sizes?
   

I have never used the pointsize option. Increasing the height and the
width decreases the relative size of the font so if you scale down the
figure in LaTeX you would get smaller fonts.

 I am using WinEdt with MikTeX set-up. 
   

as well do I ...



 Latex complains about the [scale...] part for any scale in 

 \includegraphics[scale=1]{} 
   

A scale of one does not make any sense since it is 100% or with other
words you could leave that out.

If you want to scale the figure to the text width you could write
\begin{figure}
\includegraphic[width=\textwidth]{sgr6100.pdf }
\end{figure}

btw. I use the postscript device to produce an eps file and dvipdfm to
convert the dvi to pdf

postscript(sgr6100.eps, width = 6, height = 6, horizontal = FALSE,
onefile =T, paper=special)

and then in LaTeX:

\begin{figure}
\includegraphic[width=\textwidth]{sgr6100}
\end{figure}



 and pauses. Suggestions will be appreciated about what is the best way to 
 scale
 R graphics for inclusion in LaTeX.

   
see above.

Have a look at a good LaTeX guide like lshort
http://tobi.oetiker.ch/lshort/lshort.pdf
(e.g. p.73/73)

Stefan Grosse

__
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] A statement over multiple lines (i.e. the ... feature in Matlab)

2006-10-05 Thread Wee-Jin Goh

On 5 Oct 2006, at 08:39, Robin Hankin wrote:

 Hi  Wee-Jin

 you can block out bits of R code with

 if(FALSE){
  code not executed
 }

 For the line breaking, R deals with incomplete lines by not  
 executing the statement
 until you finish it.  In the function case, it waits for you to  
 close a bracket.


 If you type:

 myFunc(a=3,
 b=5,
 c=6,
 d=7
 )

 myFunc() will only execute
 when you close the bracket


 HTH

 rksh


Great stuff. That's exactly what I was looking for.

Thanks a bunch :)
Wee-Jin

__
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] A statement over multiple lines (i.e. the ... feature in

2006-10-05 Thread Ted Harding
On 05-Oct-06 Wee-Jin Goh wrote:
 Hello again list,
 
 I thought I'd start a new thread, seeing as it's completely different  
 from my previous question. Some functions I have written require many  
 parameters, and so do not fit nicely into an 80 column width display.  
 This is usually avoided, by spreading that particular statement over  
 a few lines. This is something that I do in Matlab with the following:
 
 myFunc( parameter1, ...
 parameter2, ...
 parameter3, ...
 parameter4)
 
 The  ... operator in Matlab allows me to spread a statement over  
 several lines. The ... operator in R seems to be more like the ...  
 operator in C, which allows for a variable argument list.
 
 How do I accomplish this task in R? It's not a show stopper, but it  
 would make reading my code much much easier.
 
 Cheers,
 Wee-Jin

Just spread it over several lines, without the ..., which in R
you do not need for this purpose:

  myFunc(parameter1,
  parameter2,
  parameter3,
  parameter4)

The point is that R uses its syntactic rules to determine when
a component of an expression is complete. Therefore (in this instance)
the opening ( must at some stage be matched by a ), so, whether
or not you interpolate newlines, so R will continue to expect you
to input more items until ) is encountered.

Since you can (at any stage of input) have incomplete expressions
within incomplete expressions, the same principle applies until
every thing is closed up, to the outermost level.

Since R *does* use ... for a variable argument list, you
can of course use this as well, for that purpose.

The important thing to remember, though, is that an expression
which is incomplete from the point of view of your intentions
may look complete to R. So, if you want to break it in the
middle, make sure you do so at a point where it is incomplete
from R's point of view.

For example, suppose you intend

  x - 3*4*5*6*7*8

If, in R, you put

  x - 2*3*4*
  5*6*7*8

in your input file, it will work:

   x - 2*3*4*
  + 5*6*7*8
   x
  [1] 40320
   

because, syntactically, the * at the end of the first line
is expected to be followed by something, so R waits for that.

But if you put it in as

  x - 2*3*4
  *5*6*7*8

it will not, because the first line is already complete when
the newline is encountered, and then the second line will
generate a syntax error because the initial * is required
to be preceded by something. Having come to the end of a
complete valid expression, R is now waiting for the start of
a new valid expression; and an expression which starts with
a * is not valid.

Hence:

   x - 2*3*4
   *5*6*7*8
  Error: syntax error
   x
  [1] 24

Best wishes,
Ted.


E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 05-Oct-06   Time: 09:18:51
-- XFMail --

__
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] A statement over multiple lines (i.e. the ... feature in Matlab)

2006-10-05 Thread Bjørn-Helge Mevik
Robin Hankin wrote:

 For the line breaking, R deals with incomplete lines by not
 executing the statement until you finish it.

Beware, however, that syntactically valid lines do get executed
immediately (at least at the prompt).  So

1 + 2
- 3

will be interpreted as two commands (returning 3 and -3,
respectively), while

1 + 2 -
3

will be interpreted as a single command (returnig 0).

-- 
Bjørn-Helge Mevik

__
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] A statement over multiple lines (i.e. the ... feature in Matlab)

2006-10-05 Thread Robin Hankin

On 5 Oct 2006, at 09:34, Bjørn-Helge Mevik wrote:

 Robin Hankin wrote:

 For the line breaking, R deals with incomplete lines by not
 executing the statement until you finish it.

 Beware, however, that syntactically valid lines do get executed
 immediately (at least at the prompt).


yes, you're quite right.  Your observation
leads to my number 1 gotcha:

if(12)  print(12)
else
print(12)

which (properly) returns an error because the first line is  
syntactically
complete, so the else statement is processed ab initio.

So now I always always always always type


if(12){
   print(12)
} else {
   print(12)
}


if by some chance I ignore emacs/ESS's indentation telling me I've  
forgotten
a brace.


best wishes


rksh




 So

 1 + 2
 - 3

 will be interpreted as two commands (returning 3 and -3,
 respectively), while

 1 + 2 -
 3

 will be interpreted as a single command (returnig 0).

 -- 
 Bjørn-Helge Mevik

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

--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

__
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] Block comments in R?

2006-10-05 Thread Uwe Ligges


Wee-Jin Goh wrote:
 Hello list,
 
 Is there any way to perform a block comment in R? In C++, anything in  
 between a /* and */ is considered a comment, and it allows  
 programmers to comment out chunks of code for testing and debugging.  
 Is there such a feature in R?


This has frequently been asked on the list. Try to google for it. You 
will find answers like

if(FALSE){
 code block
 commented
}


or use a good editor that supports commenting and uncommenting blocks.


Uwe Ligges


 Cheers,
 Wee-Jin
 
 __
 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] compiling rgdal package on windows / macos

2006-10-05 Thread Roger Bivand
On Wed, 4 Oct 2006, Dylan Beaudette wrote:

 Greetings:
 
 As I am not a windows user, I cannot try this: is it possible to install
 rgdal on windows without having to compile it from source ?

Andy Jaworski already replied that rgdal Windows binaries are available 
from CRAN mirrors, thanks to Uwe Ligges.

 
 Compilation on MacOS is within my abilities, however each time i try and 
 install the rgdal package it dies complaining that it cannot find 
 gdal-config --- which was recently installed with GRASS. I have updated my 
 PATH environment variable, logged out, but R still cannot find the 
 gdal-config program.
 

For Mac OSX, please consult the detailed instructions on:

http://www.r-project.org/Rgeo - Maps (in the left navigation bar).

Hint: see the configure.args= argument of install.packages() and use 
--with-gdal-config= to get the address right.

No Mac OSX binaries will be provided from CRAN, there are simply too many 
Mac OSX varieties out there to be sure that the installed external 
software harmonises with the rgdal binary build. However, one way of 
getting there is referenced in Rgeo, using William Kyngesburye's 
frameworks (rgdal binaries are provided in harmony with PROJ.4 and GDAL).

By the way:

RSiteSearch(rgdal macosx)

first hit takes you to a recent reply on this topic on the R-sig-geo list.

 any tips on getting the rgdal package up and running on MacOS or Windows
 would be greatly appreciated.
 
 Cheers,
 
 
 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: [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.


Re: [R] Block comments in R?

2006-10-05 Thread Robin Hankin

On 5 Oct 2006, at 10:05, Uwe Ligges wrote:



 Wee-Jin Goh wrote:
 Hello list,

 Is there any way to perform a block comment in R? In C++, anything in
 between a /* and */ is considered a comment, and it allows
 programmers to comment out chunks of code for testing and debugging.
 Is there such a feature in R?


 This has frequently been asked on the list. Try to google for it. You
 will find answers like

 if(FALSE){
  code block
  commented
 }




That method doesn't work for me:

if(FALSE){

if(12)
   print(12)
else
   print(12)

}


returns an error.  How would I comment out that block of (incorrect)  
code?







 or use a good editor that supports commenting and uncommenting  
 blocks.


 Uwe Ligges



--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

__
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] Block comments in R?

2006-10-05 Thread Uwe Ligges


Robin Hankin wrote:
 
 On 5 Oct 2006, at 10:05, Uwe Ligges wrote:
 


 Wee-Jin Goh wrote:
 Hello list,

 Is there any way to perform a block comment in R? In C++, anything in
 between a /* and */ is considered a comment, and it allows
 programmers to comment out chunks of code for testing and debugging.
 Is there such a feature in R?


 This has frequently been asked on the list. Try to google for it. You
 will find answers like

 if(FALSE){
  code block
  commented
 }


 
 
 That method doesn't work for me:
 
 if(FALSE){
 
 if(12)
   print(12)
 else
   print(12)
 
 }


Use an editor that comments out a whole block which is what I do all the 
time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them.

Or you can go ahead and use

if(FALSE){

if(12)
   print(12)
else
   print(12)

}


or

if(FALSE){'

if(12)
   print(12)
else
   print(12)

'}


Uwe


 
 returns an error.  How would I comment out that block of (incorrect) code?
 
 
 
 
 
 
 
 or use a good editor that supports commenting and uncommenting blocks.


 Uwe Ligges


 
 -- 
 Robin Hankin
 Uncertainty Analyst
 National Oceanography Centre, Southampton
 European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743
 
 
 


__
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] searching for data.frame rows / processing of rows

2006-10-05 Thread chao gai
I tought that aggregate was the way to go, but only for large dataframes it is 
faster.

   df - read.table(stdin(),header=TRUE)
0: Location   TimeXY
1:   1  0  1.6  9.3
2:   1  3  4.2  10.4
3:   1  6  2.7  16.3
4:   2  0  0.5  2.1
5:   2  3  NA   3.6
6:   2  3  5.0  0.06
7:   2  6  3.4  14.0
8:
 aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE)
  Location Time   X Y
110 1.6  9.30
220 0.5  2.10
313 4.2 10.40
423 5.0  1.83
516 2.7 16.30
626 3.4 14.00
  system.time(  aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE)  )
[1] 0.008 0.000 0.008 0.000 0.000

  system.time(  do.call(rbind, by(df, df[2:1], colMeans, na.rm = TRUE)))
[1] 0.005 0.000 0.005 0.000 0.000


 df - data.frame(Location=rep(1:50,50),
+  Time=sample(rep(1:50,each=10),2500,replace=TRUE),
+  X=runif(2500),Y=runif(2500))

  system.time(  aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE)  )
[1] 0.162 0.000 0.163 0.000 0.000

  system.time(  do.call(rbind, by(df, df[2:1], colMeans, na.rm = TRUE)))
[1] 2.179 0.006 2.216 0.000 0.000


Kees

Gabor Grothendieck wrote:
 Grouping the data frame by the first two columns, apply colMeans
 and then rbind the resulting by-structure together:
 
 do.call(rbind, by(DF, DF[2:1], colMeans, na.rm = TRUE))
 
 
 On 10/5/06, Greg Tarpinian [EMAIL PROTECTED] wrote:
 R 2.3.1, WinXP:

 I have a puzzling problem that I suspect may be solved using
 grep or a regular expression but am at a loss how to actually do it...
 My data.frame looks like

  Location   TimeXY
        ---  ---
  1  0  1.6  9.3
  1  3  4.2  10.4
  1  6  2.7  16.3
  2  0  0.5  2.1
  2  3  NA   3.6
  2  3  5.0  0.06
  2  6  3.4  14.0

 and so forth.  I would like to search for duplicate Time values
 within a Location and take the numerical average (where possible)
 of the elements in X and Y.  These numerical averages should
 then be used to create a single row where multiple rows once
 existed.  So, I would like to obtain

  2  3  5.0  1.83

 for the two Time = 3 rows for Location = 2 and use it to replace
 these two rows.  Ideally, avoiding for(i in 1:blah) loops would be
 nice because the data.frame has about 10,000 rows that need to
 be searched and processed.  My intent is to do some comparing of
 SAS to R -- the DATA step processing in SAS is quite fast and
 using the RETAIN statement along with the LAG( ) function allows
 this sort of thing to be done rapidly.


 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.

 
 __
 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] Variables in RODBC environment

2006-10-05 Thread Søren Merser
hi, use paste:
# a string/number/date
pn -  '39R5238';
sql - paste(select u.unitid from test where part =,  pn)
^^

# an array of strings/numbers/dates
pn=c(1,2,3,4)
sql - paste(select u.unitid from test where part in (,  paste(pn, 
collapse=','), ))
^^
regards soren
obs the use of '=' and in

- Original Message - 
From: Thorsten Muehge [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Thursday, October 05, 2006 9:35 AM
Subject: [R] Variables in RODBC environment



 Hello Experts,
 how can I use variables in the RODBC environment.

 Example which does not work:

 Thanks for your help.

 Thorsten

 pn -  '39R5238';

 library(RODBC);
 odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx);
 sql - select
 u.unitid,
 from test
 where part in ('pn')
 ;
 parameter - sqlQuery(odbcobj,sql);
 odbcClose(odbcobj);

 __
 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] VGAM Package ?

2006-10-05 Thread KOITA Lassana - STAC/ACE




Hi!  R users
I would like to ask you where could we find the VGAM Package. I don't find
it in the list of packages.

Thak you for your help


Lassana KOITA
Etudes de Sécurité et d'Exploitation aéroportuaires / Safety Study 
Statistical analysis
Service Technique de l'Aviation Civile (STAC) / Civil Aviation Technical
Department
Direction Générale de l'Aviation Civile (DGAC) / French Civil Aviation
Authority
Tel: 01 49 56 80 60
Fax: 01 49 56 82 14
E-mail: [EMAIL PROTECTED]
http://www.stac.aviation-civile.gouv.fr/

__
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] VGAM Package ?

2006-10-05 Thread David Scott


google VGAM


On Thu, 5 Oct 2006, KOITA Lassana - STAC/ACE wrote:






Hi!  R users
I would like to ask you where could we find the VGAM Package. I don't find
it in the list of packages.

Thak you for your help


Lassana KOITA
Etudes de Sécurité et d'Exploitation aéroportuaires / Safety Study 
Statistical analysis
Service Technique de l'Aviation Civile (STAC) / Civil Aviation Technical
Department
Direction Générale de l'Aviation Civile (DGAC) / French Civil Aviation
Authority
Tel: 01 49 56 80 60
Fax: 01 49 56 82 14
E-mail: [EMAIL PROTECTED]
http://www.stac.aviation-civile.gouv.fr/

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



_
David Scott Visiting (Until January 07)
Department of Probability and Statistics
The University of Sheffield
The Hicks Building
Hounsfield Road
Sheffield S3 7RH
United Kingdom
Phone:  +44 114 222 3908
Email:  [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.


Re: [R] (not so real) function (more like #include)

2006-10-05 Thread Meinhard Ploner
Thanks a lot! longfun2  longfun3 work perfect for my very  
untypical problem. As I have many local variables, usual functions  
with parameters are very uncomfortable, but the code you gave me is  
great!
Meinhard


On Oct 4, 2006, at 5:25 PM, Gabor Grothendieck wrote:

 longfun could just pass a, b and d to each of the individual
 functions and each of the individual functions could pass
 out back as a return value.

 f1 - f2 - function(a, b, d) a+b+d

 longfun1 - function() {
   a - b - d - 1
   out - f1(a, b, d)
   out - f2(a, b, d) + out
   out
 }

 longfun1() # 6

 If the problem is that a, b and d actually represent 100 variables,  
 say,
 and its a nuisance to pass them all explicitly you could do this:

 f1 - f2 - function() with(parent.frame(), a + b + d)

 longfun2 - function() {
   a - b - d - 1
   out - f1()
   out - f2() + out
   out
 }

 longfun2() # 6


 or you could do it this way if you would prefer to have longfun
 control the scoping rather than having it done within each
 individual function:


 f1 - f2 - function() a+b+d

 longfun3 - function() {
   a - b - d - 1
   environment(f1) - environment(f2) - environment()
   out - f1()
   f2() + out
 }

 longfun3() # 6


 On 10/4/06, Meinhard Ploner [EMAIL PROTECTED] wrote:
 Hello R users,

 I have a very long function with parts of that looking like a list of
 jobs. The first part of the function defines a lot of local
 variables, which are used by the jobs. Now I extracted the job part
 und putted them into an external file, used by source(, local=T),
 which is more comfortable to me. Apart from that it gives me the
 opportunity that more users can work on the project. Is it possible
 to define a function for that code and passing to it the environment
 of f() without save(list=ls()) and load()  ?

 Thanks in advance
 Meinhard Ploner

 longfun - f() {

## lot of local variables:
a - ...
b - ...
d - ...
...

out - ...

source(job1.txt, local=TRUE)   #changes out, uses a, b,  
 d, ...

source(job2.txt, local=TRUE)  # changes out, uses a, b,  
 d, ...
...
 }



 version
_
 platform   i386-apple-darwin8.6.1
 arch   i386
 os darwin8.6.1
 system i386, darwin8.6.1
 status
 major  2
 minor  3.1
 year   2006
 month  06
 day01
 svn rev38247
 language   R
 version.string Version 2.3.1 (2006-06-01)

 __
 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] Variables in RODBC environment

2006-10-05 Thread Bonfigli Sandro
 how can I use variables in the RODBC environment.
 
 Example which does not work:
 
 Thanks for your help.
 
 Thorsten
 
 pn -  '39R5238';
 
 library(RODBC);
 odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx);
 sql - select
 u.unitid,
 from test
 where part in ('pn')
 ;
 parameter - sqlQuery(odbcobj,sql);
 odbcClose(odbcobj);
 

You can compose your query simply using paste.

In your example it would be

pn -  '39R5238';

sql - paste(select
u.unitid,
from test
where part in (', 
pn, 
'));

or, to avoid problems with the newline character:

sql - paste(select,
u.unitid,,
from test,
where part in (', 
pn, 
'));

In this case you'd have:
 sql
[1] select u.unitid, from test where part in (' 39R5238 ')

It's clear that, if the spaces in (' 39R5238 ') are a problem for you ,
you can use 

sql - paste(select ,
u.unitid, ,
from test ,
where part in (', 
pn, 
'), sep = );

which lends to 

 sql
[1] select u.unitid, from test where part in ('39R5238')

HTH

  Sandro Bonfigli

__
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] VGAM Package ?

2006-10-05 Thread Gavin Simpson
On Thu, 2006-10-05 at 11:39 +0200, KOITA Lassana - STAC/ACE wrote:
 
 
 
 Hi!  R users
 I would like to ask you where could we find the VGAM Package. I don't find
 it in the list of packages.
 
 Thak you for your help
 
 
 Lassana KOITA

That's because it isn't on CRAN yet. You can find a link to VGAM on the
Environmetrics Task View:

http://www.stats.bris.ac.uk/R/src/contrib/Views/Environmetrics.html

HTH

G

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

__
R-help@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] format of data for use in clim.pact

2006-10-05 Thread isidora k
Goodmorning everyone!
Just a quick question: 
I want to apply eof analysis and downscaling using the
clim.pact package.
What form should my data have? Only netCDF?
Right now my data are of the form of ascii files for
each month with 3 columns corresponding to latitude,
longitude and value of my observations.
Any idea on how to handle them would be welcome.
Thanks a lot
Isidora

__
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] glm with nesting

2006-10-05 Thread Jeffrey Stratford
I just had a manuscript returned with the biggest problem being the
analysis.  Instead of using principal components in a regression I've
been asked to analyze a few variables separately. So that's what I'm
doing.

I pulled a feather from young birds and we quantified certain aspects of
the color of those feathers.  Since I often have more than one sample
from a nest, I thought I should use a nested design. 

Here's the code I've been using and I'd appreciate if someone could look
it over and see if it was correct. 

bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale,
family=gaussian, na.action=na.omit)

where rtot = total reflectance, box = nest box (i.e., birdhouse), 
julian = day of the year and purbank = the proportion of urban cover in
a 1 km buffer around the nest box.  I'm not interested in the box effect
and I've seperated males and female chicks.   

 I've asked about nestedness before and I was given code that included
| to indicate nestedness but this indicates a grouping does it not?  I
suspect that there is something wrong.  In the summary I get 

Coefficients:
  Estimate Std. Error t value Pr(|t|)
(Intercept)  2.880e-01  3.224e-03  89.322   2e-16 ***
box -3.219e-05  6.792e-05  -0.4740.636
box:julian   7.093e-08  3.971e-07   0.1790.859
box:purbank -1.735e-05  1.502e-04  -0.1150.908   

The other question I have is how do I test a null hypothesis - no
explanatory variables?  [rtot ~ NULL?]

Many thanks,

Jeff 




Jeffrey A. Stratford, Ph.D.
Postdoctoral Associate
331 Funchess Hall
Department of Biological Sciences
Auburn University
Auburn, AL 36849
334-329-9198
FAX 334-844-9234
http://www.auburn.edu/~stratja

__
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] glm with nesting

2006-10-05 Thread Peter Dalgaard
Jeffrey Stratford [EMAIL PROTECTED] writes:

 I just had a manuscript returned with the biggest problem being the
 analysis.  Instead of using principal components in a regression I've
 been asked to analyze a few variables separately. So that's what I'm
 doing.
 
 I pulled a feather from young birds and we quantified certain aspects of
 the color of those feathers.  

 Since I often have more than one sample
 from a nest, I thought I should use a nested design. 

Notwithstanding comments below, that quote could be aiming for the
fortunes package... 

 
 Here's the code I've been using and I'd appreciate if someone could look
 it over and see if it was correct. 
 
 bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale,
 family=gaussian, na.action=na.omit)
 
 where rtot = total reflectance, box = nest box (i.e., birdhouse), 
 julian = day of the year and purbank = the proportion of urban cover in
 a 1 km buffer around the nest box.  I'm not interested in the box effect
 and I've seperated males and female chicks.   
 
  I've asked about nestedness before and I was given code that included
 | to indicate nestedness but this indicates a grouping does it not?  I
 suspect that there is something wrong.  In the summary I get 
 
 Coefficients:
   Estimate Std. Error t value Pr(|t|)
 (Intercept)  2.880e-01  3.224e-03  89.322   2e-16 ***
 box -3.219e-05  6.792e-05  -0.4740.636
 box:julian   7.093e-08  3.971e-07   0.1790.859
 box:purbank -1.735e-05  1.502e-04  -0.1150.908   

Several things look wrong here. 

Most importantly, you appear to have single-degree of freedom effects
(t tests) of things that appear not to be linear effects: Certainly,
you have more than two nest boxes, but also day of year as a linear
term looks suspicious to me. Unless there is something I have missed
completely, box should be a factor variable, and you might also need
trigonometric terms for the julian effect (depending on what sort of
time spans we are talking about.)

Secondly, notation like box/julian suggests that julian only makes
sense within a nest box i.e. 1st of March in one box is completely
different from 1st of March in another box (the notation is more
commonly used to describe bird number within nests and the like). And
with purbank presumably constant for measurements from the same box,
the box:purbank term looks strange indeed.

If you want to take account of a between-box variation in the effect
of covariates, you probably need to add them as variance components,
but this requires non-glm software, either lme() or lmer(). However,
instructing you on those is outside the scope of this mailing list,
and you may need to find a local consultant.

 The other question I have is how do I test a null hypothesis - no
 explanatory variables?  [rtot ~ NULL?]
 
 Many thanks,
 
 Jeff 
 
 
 
 
 Jeffrey A. Stratford, Ph.D.
 Postdoctoral Associate
 331 Funchess Hall
 Department of Biological Sciences
 Auburn University
 Auburn, AL 36849
 334-329-9198
 FAX 334-844-9234
 http://www.auburn.edu/~stratja
 
 __
 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.
 

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

__
R-help@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] [Fwd: Re: Block comments in R?]

2006-10-05 Thread Philippe Grosjean
Ooops! Sorry, I send it only to Uwe Ligges the first time.
Best,

Philippe Grosjean

This is perhaps another solution, more elegant in the way the block
comment is written... but it requires to redefine `!` and slows it a
little bit because it tests first its arguments before calling
.Primitive(!):

It takes advantage of `!` being not defined for character arguments:

 !2
[1] FALSE
 !some text
Error in !some text : invalid argument type

So, now, we will define it for character arguments (invisibly returns
the argument)

 `!`- function(x)
+   if (inherits(x, character) == FALSE)
+   .Primitive(!)(x) else invisible(x)


Now, `!` can be used to construct a block comment:

 A R script with block comments =

1+1 # This is a line comment

!'
This is a block comment spread
on several lines...
'

ls()

!!'
This is another block comment, possibly of higher
importance than the previous one
'
search()

!!!'
For color syntax highlighting and to better detect
the end of block comments, one may also decide to
use the same code for opening and closing the comments
like it is the present case
!!!'

!'
Note that the only constraint is to escape single quotes
in block comments (or not to use single quotes)
Of course, one could also decide to use double quotes
instead of single quotes
!'

!'
Now, it would be nice to have a little patch of .Primitive(!)
that simply displays no error message in case the argument of
`!`is a character sting. So, the hock would not be required
any more
!'

 And of the R script =

Best,

Philippe Grosjean


Uwe Ligges wrote:
 
 Robin Hankin wrote:
 On 5 Oct 2006, at 10:05, Uwe Ligges wrote:


 Wee-Jin Goh wrote:
 Hello list,

 Is there any way to perform a block comment in R? In C++, anything in
 between a /* and */ is considered a comment, and it allows
 programmers to comment out chunks of code for testing and debugging.
 Is there such a feature in R?

 This has frequently been asked on the list. Try to google for it. You
 will find answers like

 if(FALSE){
  code block
  commented
 }



 That method doesn't work for me:

 if(FALSE){

 if(12)
   print(12)
 else
   print(12)

 }
 
 
 Use an editor that comments out a whole block which is what I do all the 
 time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them.
 
 Or you can go ahead and use
 
 if(FALSE){
 
 if(12)
print(12)
 else
print(12)
 
 }
 
 
 or
 
 if(FALSE){'
 
 if(12)
print(12)
 else
print(12)
 
 '}
 
 
 Uwe
 
 
 returns an error.  How would I comment out that block of (incorrect) code?







 or use a good editor that supports commenting and uncommenting blocks.


 Uwe Ligges


 -- 
 Robin Hankin
 Uncertainty Analyst
 National Oceanography Centre, Southampton
 European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743




 
 __
 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] Partition into quantiles

2006-10-05 Thread Alberto Monteiro
Is there any function that divides a sample into N quantiles?

For example, for N = 2, this would be the solution:

x - rnorm(100)
m - median(x)
q - ifelse(x = median, 1, 2)

Alberto Monteiro

__
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] Partition into quantiles

2006-10-05 Thread Peter Dalgaard
Alberto Monteiro [EMAIL PROTECTED] writes:

 Is there any function that divides a sample into N quantiles?
 
 For example, for N = 2, this would be the solution:
 
 x - rnorm(100)
 m - median(x)
 q - ifelse(x = median, 1, 2)

Have a look at

 N - 2
 table(cut(x,quantile(x,seq(0,1,1/N)), include.lowest=TRUE))

[-2.78,0.205]  (0.205,2.22]
   5050

 table(cut(x,c(-Inf,quantile(x,(1:(N-1))/N),Inf)))

(-Inf,0.205] (0.205, Inf]
  50   50

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

__
R-help@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] Partition into quantiles

2006-10-05 Thread Richard M. Heiberger
cut.quantile - function(x, N, use.ppoints=TRUE) {
  qq - if (use.ppoints) c(0,ppoints(N-1),1)
  else seq(0, 1, 1/N)
  breaks - quantile(x, qq)
  breaks[1] - breaks[1] - 1 
  breaks[N+1] - breaks[N+1]+1
  cut(x, breaks)
}

tmpT - cut.quantile(rnorm(100), 10, TRUE)
table(tmpT)
tmpF - cut.quantile(rnorm(100), 10, FALSE)
table(tmpF)

__
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] Partition into quantiles

2006-10-05 Thread David Barron
You might also want to look at the function quantcut in the gtools
package (part of the gregmisc bundle).



On 05/10/06, Alberto Monteiro [EMAIL PROTECTED] wrote:
 Is there any function that divides a sample into N quantiles?

 For example, for N = 2, this would be the solution:

 x - rnorm(100)
 m - median(x)
 q - ifelse(x = median, 1, 2)

 Alberto Monteiro

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



-- 
=
David Barron
Said Business School
University of Oxford
Park End Street
Oxford OX1 1HP

__
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] [Fwd: Re: Block comments in R?]

2006-10-05 Thread Barry Rowlingson
Philippe Grosjean wrote:

 
 It takes advantage of `!` being not defined for character arguments:
 

  *gasp*

  how much R code is destined to feature on www.thedailywtf.com in the 
future?

  whats the chances of block commenting being included in a future 
version? and more generally, is there a feature request system, or 
discussion of whats going in new releases? the nearest I can find is the 
'ideas' file:

http://developer.r-project.org/ideas.txt

  which has a bit of a 1990's feel to it (okay, it is in the 'Older 
Material' section).


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] glm with nesting

2006-10-05 Thread Jeffrey Stratford
Peter and list,

Thanks for the response.  A did add box as a factor (box -
factor(box)).  Julian should be linear - bluebird chicks are bluer as
the season progresses from March to August.  

I did try the following

rtot.lme - lmer(rtot ~ sex +(purban|box:chick) + (purban|box), data=bb,
na.action=na.omit) # from H. Doran

and

rtot.lme2 - lme(fixed=rtot ~ sex + purban + sexv:purban, data = bb,
random = ~1 |box) # from K. Jones [EMAIL PROTECTED]

but these did not work (months ago and I don't remember exactly why) and
I have since seperated males and females and added day of the year
(julian).  But | does indicate grouping not nested, correct?

Could someone suggest some coding that might work?

Thanks again,

Jeff


 Peter Dalgaard [EMAIL PROTECTED] 10/05/06 7:14 AM 
Jeffrey Stratford [EMAIL PROTECTED] writes:

 I just had a manuscript returned with the biggest problem being the
 analysis.  Instead of using principal components in a regression I've
 been asked to analyze a few variables separately. So that's what I'm
 doing.
 
 I pulled a feather from young birds and we quantified certain aspects of
 the color of those feathers.  

 Since I often have more than one sample
 from a nest, I thought I should use a nested design. 

Notwithstanding comments below, that quote could be aiming for the
fortunes package... 

 
 Here's the code I've been using and I'd appreciate if someone could look
 it over and see if it was correct. 
 
 bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale,
 family=gaussian, na.action=na.omit)
 
 where rtot = total reflectance, box = nest box (i.e., birdhouse), 
 julian = day of the year and purbank = the proportion of urban cover in
 a 1 km buffer around the nest box.  I'm not interested in the box effect
 and I've seperated males and female chicks.   
 
  I've asked about nestedness before and I was given code that included
 | to indicate nestedness but this indicates a grouping does it not?  I
 suspect that there is something wrong.  In the summary I get 
 
 Coefficients:
   Estimate Std. Error t value Pr(|t|)
 (Intercept)  2.880e-01  3.224e-03  89.322   2e-16 ***
 box -3.219e-05  6.792e-05  -0.4740.636
 box:julian   7.093e-08  3.971e-07   0.1790.859
 box:purbank -1.735e-05  1.502e-04  -0.1150.908   

Several things look wrong here. 

Most importantly, you appear to have single-degree of freedom effects
(t tests) of things that appear not to be linear effects: Certainly,
you have more than two nest boxes, but also day of year as a linear
term looks suspicious to me. Unless there is something I have missed
completely, box should be a factor variable, and you might also need
trigonometric terms for the julian effect (depending on what sort of
time spans we are talking about.)

Secondly, notation like box/julian suggests that julian only makes
sense within a nest box i.e. 1st of March in one box is completely
different from 1st of March in another box (the notation is more
commonly used to describe bird number within nests and the like). And
with purbank presumably constant for measurements from the same box,
the box:purbank term looks strange indeed.

If you want to take account of a between-box variation in the effect
of covariates, you probably need to add them as variance components,
but this requires non-glm software, either lme() or lmer(). However,
instructing you on those is outside the scope of this mailing list,
and you may need to find a local consultant.

 The other question I have is how do I test a null hypothesis - no
 explanatory variables?  [rtot ~ NULL?]
 
 Many thanks,
 
 Jeff 
 
 
 
 
 Jeffrey A. Stratford, Ph.D.
 Postdoctoral Associate
 331 Funchess Hall
 Department of Biological Sciences
 Auburn University
 Auburn, AL 36849
 334-329-9198
 FAX 334-844-9234
 http://www.auburn.edu/~stratja
 
mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

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

__
R-help@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] [Fwd: Re: Block comments in R?]

2006-10-05 Thread Duncan Murdoch
On 2006-10-5 9:20, Barry Rowlingson wrote:
 Philippe Grosjean wrote:
 
 
 It takes advantage of `!` being not defined for character arguments:
 
 
   *gasp*
 
   how much R code is destined to feature on www.thedailywtf.com in the 
 future?
 
   whats the chances of block commenting being included in a future 
 version? and more generally, is there a feature request system, or 
 discussion of whats going in new releases? the nearest I can find is the 
 'ideas' file:

There's a feature request system (one of the categories of bug report is 
  wishlist).

I think a specific proposal with some work behind it would probably get 
action in a case like this.  That is, block comments would be nice 
likely won't get acted on, but block comments in the form ... which is 
used by other language ..., which has the fantastic properties ... and 
no downside, and which could be implemented by the following patch to 
the parser,, likely would.

Duncan Murdoch


 
 http://developer.r-project.org/ideas.txt
 
   which has a bit of a 1990's feel to it (okay, it is in the 'Older 
 Material' section).
 
 
 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.

__
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] glm with nesting

2006-10-05 Thread Doran, Harold
It's not really possible to help without knowing what errors you received and 
maybe some reproducible code. I think I remember this, though. From what I 
recall, there was no distinction between box and chick, so you cannot estimate 
both variance components. 

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 Jeffrey Stratford
 Sent: Thursday, October 05, 2006 9:27 AM
 To: [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] glm with nesting
 
 Peter and list,
 
 Thanks for the response.  A did add box as a factor (box - 
 factor(box)).  Julian should be linear - bluebird chicks are 
 bluer as the season progresses from March to August.  
 
 I did try the following
 
 rtot.lme - lmer(rtot ~ sex +(purban|box:chick) + 
 (purban|box), data=bb,
 na.action=na.omit) # from H. Doran
 
 and
 
 rtot.lme2 - lme(fixed=rtot ~ sex + purban + sexv:purban, 
 data = bb, random = ~1 |box) # from K. Jones [EMAIL PROTECTED]
 
 but these did not work (months ago and I don't remember 
 exactly why) and I have since seperated males and females and 
 added day of the year (julian).  But | does indicate 
 grouping not nested, correct?
 
 Could someone suggest some coding that might work?
 
 Thanks again,
 
 Jeff
 
 
  Peter Dalgaard [EMAIL PROTECTED] 10/05/06 7:14 AM 
 Jeffrey Stratford [EMAIL PROTECTED] writes:
 
  I just had a manuscript returned with the biggest problem being the 
  analysis.  Instead of using principal components in a 
 regression I've 
  been asked to analyze a few variables separately. So that's 
 what I'm 
  doing.
  
  I pulled a feather from young birds and we quantified 
 certain aspects 
  of the color of those feathers.
 
  Since I often have more than one sample from a nest, I thought I 
  should use a nested design.
 
 Notwithstanding comments below, that quote could be aiming 
 for the fortunes package... 
 
  
  Here's the code I've been using and I'd appreciate if someone could 
  look it over and see if it was correct.
  
  bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale, 
  family=gaussian, na.action=na.omit)
  
  where rtot = total reflectance, box = nest box (i.e., birdhouse), 
  julian = day of the year and purbank = the proportion of 
 urban cover 
  in a 1 km buffer around the nest box.  I'm not interested 
 in the box effect
  and I've seperated males and female chicks.   
  
   I've asked about nestedness before and I was given code 
 that included 
  | to indicate nestedness but this indicates a grouping 
 does it not?  
  I suspect that there is something wrong.  In the summary I get
  
  Coefficients:
Estimate Std. Error t value Pr(|t|)
  (Intercept)  2.880e-01  3.224e-03  89.322   2e-16 ***
  box -3.219e-05  6.792e-05  -0.4740.636
  box:julian   7.093e-08  3.971e-07   0.1790.859
  box:purbank -1.735e-05  1.502e-04  -0.1150.908   
 
 Several things look wrong here. 
 
 Most importantly, you appear to have single-degree of freedom 
 effects (t tests) of things that appear not to be linear 
 effects: Certainly, you have more than two nest boxes, but 
 also day of year as a linear term looks suspicious to me. 
 Unless there is something I have missed completely, box 
 should be a factor variable, and you might also need 
 trigonometric terms for the julian effect (depending on what 
 sort of time spans we are talking about.)
 
 Secondly, notation like box/julian suggests that julian only 
 makes sense within a nest box i.e. 1st of March in one box is 
 completely different from 1st of March in another box (the 
 notation is more commonly used to describe bird number within 
 nests and the like). And with purbank presumably constant for 
 measurements from the same box, the box:purbank term looks 
 strange indeed.
 
 If you want to take account of a between-box variation in the 
 effect of covariates, you probably need to add them as 
 variance components, but this requires non-glm software, 
 either lme() or lmer(). However, instructing you on those is 
 outside the scope of this mailing list, and you may need to 
 find a local consultant.
 
  The other question I have is how do I test a null hypothesis - no 
  explanatory variables?  [rtot ~ NULL?]
  
  Many thanks,
  
  Jeff
  
  
  
  
  Jeffrey A. Stratford, Ph.D.
  Postdoctoral Associate
  331 Funchess Hall
  Department of Biological Sciences
  Auburn University
  Auburn, AL 36849
  334-329-9198
  FAX 334-844-9234
  http://www.auburn.edu/~stratja
  
 mailman/listinfo/r-help
  PLEASE do read the posting guide 
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
  
 
 -- 
O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
   c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
  (*) \(*) -- University of Copenhagen   Denmark  Ph:  
 (+45) 35327918
 ~~ - ([EMAIL PROTECTED])  

[R] Changes in console preferences won't take

2006-10-05 Thread David Kaplan
Greetings,

I've attempted to change the font size directly in Rconsole as well as 
through the GUI preferences.  Neither seems to take even when I save the 
preference changes.  When I make the change through the gui preference, 
and click apply the change is made, but when I save it, close, and 
reopen, there is no change.  Any thouughts?  I'm running R.2.40.

David
-- 
===
David Kaplan, Ph.D.
Professor
Department of Educational Psychology
University of Wisconsin - Madison
Educational Sciences, Room, 1061
1025 W. Johnson Street
Madison, WI 53706

email: [EMAIL PROTECTED]
homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm
Phone: 608-262-0836

__
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] Block comments in R?

2006-10-05 Thread Seth Falcon
Uwe Ligges [EMAIL PROTECTED] writes:
 Use an editor that comments out a whole block which is what I do all the 
 time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of
 them.

This, of course, works.  The if(FALSE) approach does not because it
requires the comment to be syntactically correct.  The multiline
string trick is _almost_ there, the problem is that both  and ' are
fairly common in text and having to escape that is a pain.  

My wtf feature request is to add a multiline string delimiter ala
Python like .

8-P

+ seth

__
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] Partition into quantiles

2006-10-05 Thread Frank E Harrell Jr
David Barron wrote:
 You might also want to look at the function quantcut in the gtools
 package (part of the gregmisc bundle).

Also, cut2 in Hmisc will do this and will label the intervals compactly.

Frank

 
 
 
 On 05/10/06, Alberto Monteiro [EMAIL PROTECTED] wrote:
 Is there any function that divides a sample into N quantiles?

 For example, for N = 2, this would be the solution:

 x - rnorm(100)
 m - median(x)
 q - ifelse(x = median, 1, 2)

 Alberto Monteiro

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

 
 


-- 
Frank E Harrell Jr   Professor and Chair   School of Medicine
  Department of Biostatistics   Vanderbilt University

__
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] research

2006-10-05 Thread aziz tomi
  dear sir
   iam algerian student iam 30 years old  i finished my postgraduate in applied 
economic and statistics .i prepare my doctorat about effect of kyoto protocol 
on our economic.
  i need the document and the articals and the adress of doctor in this domain 
in
  your university 
   i shall be looking forward and hoping that give my request forward 
consideration 
   rachid toumache 
  adress: 11 rue doudou mothtar ben aknoun algeria  
 institut national de la planification et de la 
statistique inps 



-

[[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] Partition into quantiles

2006-10-05 Thread Alberto Monteiro
Sorry, folks. Thanks for the help, but none of the suggestions 
worked quite as I wanted :-/

So I wrote the code. It seems to work:

gera_particao - function(x, n) {
  y - sort(x)
  icut - as.integer(seq(1, length(x)+1, length = n + 1))
  icut - icut[c(-(n+1))]
  ycut - y[icut]
  for (i in 1:length(x))
xpart[i] - sum(x[i] = ycut)
  return(xpart)
}

First, x is sorted to y. Then, I select the minima y
of each of the n segments. Then, an ugly and slow loop,
counts for each x how many ycut lie below them.

 x - runif(12)
 x
 [1] 0.2971455266 0.6112766485 0.5571073645 0.5886481798 0.7499023860
 [6] 0.1681732289 0.6319822536 0.0005354732 0.8055324992 0.8841625380
[11] 0.0726578285 0.6250309648
 gera_particao(x, 2)
 [1] 1 2 1 1 2 1 2 1 2 2 1 2
 gera_particao(x, 3)
 [1] 1 2 2 2 3 1 3 1 3 3 1 2
 gera_particao(x, 4)
 [1] 2 3 2 2 4 1 3 1 4 4 1 3
 gera_particao(x, 12)
 [1]  4  7  5  6 10  3  9  1 11 12  2  8

Alberto Monteiro

__
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] warning message in nlm

2006-10-05 Thread Spencer Graves
  Without 'msc39', I can't say for sure, but I doubt if it's 
anything to worry about.  It looks like 'nlm' tests values for theta and 
len that produce either NA or Inf in computing 'loglikcs'.  Since you 
got an answer, it does not look to me like it's anything worth worrying 
about;  just make sure you are minimizing (-log(likelihood)), not 
maximizing it.  For more detail, you can run 'nlm' with print.level = 2 
rather than the default of 0 -- and make contour plots of loglikcs in an 
appropriate region of the optimum, as suggested in Venables and Ripley 
(2002) Modern Applied Statistics with S (Springer), discussing 
'expand.grid' and 'contour'. 

  Hope this helps. 
  Spencer Graves

singyee ling wrote:
 Dear R-users,

 I am trying to find the MLEs for a loglikelihood function (loglikcs39) and
 tried using both optim and nlm.

 fredcs39-function(b1,b2,x){return(exp(b1+b2*x))}
 loglikcs39-function(theta,len){
 sum(mcs39[1:len]*fredcs39(theta[1],theta[2],c(8:(7+len))) - pcs39[1:len] *
 log(fredcs39(theta[1],theta[2],c(8:(7+len)
 }
 theta.start-c(0.1,0.1)


 1. The output from using optim is as follow
 --

 optcs39-optim(theta.start,loglikcs39,len=120,method=BFGS)
   
 optcs39
 
 $par
 [1] -1.27795226 -0.03626846

 $value
 [1] 7470.551

 $counts
 function gradient
  133   23

 $convergence
 [1] 0

 $message
 NULL

 2. The output from using nlm is as follow
 ---

   
 outcs39-nlm(loglikcs39,theta.start,len=120)
 
 Warning messages:
 1: NA/Inf replaced by maximum positive value
 2: NA/Inf replaced by maximum positive value
 3: NA/Inf replaced by maximum positive value
 4: NA/Inf replaced by maximum positive value
 5: NA/Inf replaced by maximum positive value
 6: NA/Inf replaced by maximum positive value
 7: NA/Inf replaced by maximum positive value
   
 outcs39
 
 $minimum
 [1] 7470.551

 $estimate
 [1] -1.27817854 -0.03626027

 $gradient
 [1] -8.933577e-06 -1.460512e-04

 $code
 [1] 1

 $iterations
 [1] 40


 As you can see, the values obtained from using both functions are very
 similar.  But, what puzzled is the warning message that i got from using
 nlm. Could anyone please shed some light on how this warning message come
 about and whether it is a cause for concern?


 Many thanks in advance for any advice!

 singyee

   [[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] Block comments in R?

2006-10-05 Thread Peter Dalgaard
Seth Falcon [EMAIL PROTECTED] writes:

 Uwe Ligges [EMAIL PROTECTED] writes:
  Use an editor that comments out a whole block which is what I do all the 
  time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of
  them.
 
 This, of course, works.  The if(FALSE) approach does not because it
 requires the comment to be syntactically correct.  The multiline
 string trick is _almost_ there, the problem is that both  and ' are
 fairly common in text and having to escape that is a pain.  
 
 My wtf feature request is to add a multiline string delimiter ala
 Python like .

Still gives you problems with nested comments. *IF* we want to go down
that route, we should have directional symbols like


dsaldfysdfk


(which would likely drive version control users up the wall...)

or 

/# kdlsfk 
   dfkjsdfl #/

(but those character combinations have a marginal risk of being used
in existing code).

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

__
R-help@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] matrix multiplication

2006-10-05 Thread Ya-Hsiu Chuang
Dear all,

I have 2 matrices, one is a 8x3 matrix, called X; the other matrix is a 3x3 
indicator matrix with the diagonal element as 0 or 1. when a variable is 
included in the model, the corresponding diagonal element is 1, otherwise, it 
is 0.  Let A be a set of matrices that contain the possible indicator matrix. 
e.g.,
A= [A1, A2, A3], 
where A1= [1,0,0;0,0,0,0,0,0],
A2 =[1,0,0;0,1,0,0,0,0], 
A3 =[1,0,0;0,0,0,0,1,0]

In order to derive the new X matries depending on the indicator matrices, I use 
for loops to multiply both X and A as following

p-3
for ( i in 1:p){
  XX- X%*%A[i]
}

However, it only shows the result when i=p. How can I derive the results which 
include all possible value of XX[i], i=1,..,p, rather than just i=p

thanks for any help

Ya-Hsiu



[[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] Changes in console preferences won't take

2006-10-05 Thread Duncan Murdoch
On 10/5/2006 9:41 AM, David Kaplan wrote:
 Greetings,
 
 I've attempted to change the font size directly in Rconsole as well as 
 through the GUI preferences.  Neither seems to take even when I save the 
 preference changes.  When I make the change through the gui preference, 
 and click apply the change is made, but when I save it, close, and 
 reopen, there is no change.  Any thouughts?  I'm running R.2.40.

I would guess you've got conflicting settings in multiple Rconsole 
files, or you're saving the change to the wrong place. See the ?Rconsole 
man page for a description of how R looks for the file (and the example 
section of it to see what would normally be loaded on your own system).

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] glm with nesting

2006-10-05 Thread Jeffrey Stratford
Harold and list,

I've changed a few things since the last time so I'm really starting
from scratch.  

I start with

bbmale - read.csv(c:\\eabl\\2004\\feathers\\male_feathers2.csv,
header=TRUE)
box -factor(box)
chick - factor(chick)

Here's a sample of the data

box,chick,julian,cltchsz,mrtot,cuv,cblue,purbank,purban2,purban1,pgrassk,pgrass2,pgrass1,grassdist,grasspatchk
1,2,141,2,21.72290152,0.305723811,0.327178868,0.003813435,0.02684564,0.06896552,0.3282487,0.6845638,0.7586207,0,3.73
4,1,164,4,18.87699007,0.281863299,0.310935559,0.06072162,0.2080537,0.06896552,0.01936052,0,0,323.1099,0.2284615
4,2,164,4,19.64359348,0.294117388,0.316049817,0.06072162,0.2080537,0.06896552,0.01936052,0,0,323.1099,0.2284615
7,1,118,4,13.48699876,0.303649408,0.31765218,0.3807568,0.4362416,0.6896552,0.06864183,0.03355705,0,94.86833,0.468
12,1,180,4,21.42196378,0.289731361,0.317562294,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032
12,2,180,4,18.79487905,0.286052077,0.316367349,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032
12,3,180,4,12.83127682,0.260197475,0.292636914,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032
15,1,138,4,20.07161467,0.287632782,0.318671887,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
15,2,138,4,17.61146256,0.305581768,0.315848051,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
15,3,138,4,20.36397134,0.271795667,0.30539683,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
15,4,138,4,20.81940158,0.269468041,0.304160648,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818

As you can see I have multiple boxes ( 70).  Sometimes I have multiple
chicks per box each having their own response  to mrtot, cuv, and cblue
but the same landscape variables for that box.  Chick number is randomly
assigned and is not an effect I'm interested in.  I'm really not
interested in the box effect either.  I would like to know if landscape
affects the color of chicks (which may be tied into chick
health/physiology).  We also know that chicks get bluer as the season
progresses and that clutch size (cltchsz) has an effect so I'm including
that as covariates.  

Hopefully, this clears things up a bit. 

I do have the MASS and MEMS (Pineiro's) texts in hand. 

Many thanks,

Jeff

__
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 text with lattice

2006-10-05 Thread Ritwik Sinha
Hi,

On a related note, if I wanted to add different texts to different
panels, should I stick to using trellis.focus() for each text in each
panel? I cannot figure out a way to do it using a panel function.

Ritwik.

On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote:
 On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
  Here are two possibilities.  The first use trellis.focus/trellis/unfocus to
  add text subsequent to drawing the xyplot and the second uses a
  custom panel:
 
  xyplot(x ~ x, data = data.frame(x = 1:10))
  trellis.focus(panel, 1, 1)
  panel.text(x=2, y=4, labels=Text)
  trellis.unfocus()
 
  xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) {
  panel.xyplot(...)
  panel.text(x=2, y=4, labels=Text)
  })

 Right, on a more general note, this is necessary because it is not clear what

 ltext(x=2, y=4, labels=Text)

 should do for a multi-panel plot. On an even more general note, you
 will keep getting in trouble when using lattice if you (1) try to
 follow the incremental addition approach of standard graphics or (2)
 use the par() system in any way.

 Deepayan

  On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote:
   Hello,
   I've decided to take the leap and try my hand at the lattice package,
   though I am getting stuck at what one might consider a trivial problem,
   plotting text at a point in a graph. Apologies in advance if (that) I'm
   missing something extremely basic.
  
   Consider in base graphics:
plot(1:10)
text(2, 4, Text)
   In the above you will see text centered at the point (2, 4) on the
   graph.
  
   Now I would like to try to do the same thing using the lattice package:
  
xyplot(x ~ x, data = data.frame(x = 1:10))
ltext(x=2, y=4, labels=Text)
panel.text(x=2, y=4, labels=Text)
grid.text(label=Text, x=2, y=4)
grid.text(label=Text, x=unit(2, native), y=unit(4, native))
  
   None of the above four commands puts the Text at the (2, 4) point on
   the graph. Any help with this would be appreciated! Also, if I have more
   than one panel and would like to place text at different points on
   different panels how would I do this?
  
   Also, note that I'm hoping to use text to label interesting points in a
   levelplot, but am using the above xyplot as an example.
  
   Thanks,
   Robert
  
  
   Robert McGehee
   Quantitative Analyst
   Geode Capital Management, LLC
   53 State Street, 5th Floor | Boston, MA | 02109
   Tel: 617/392-8396Fax:617/476-6389
   mailto:[EMAIL PROTECTED]
  
  
  
   This e-mail, and any attachments hereto, are intended for us...{{dropped}}
  
   __
   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.
 


 --
 http://www.stat.wisc.edu/~deepayan/

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



-- 
Ritwik Sinha
Graduate Student
Epidemiology and Biostatistics
Case Western Reserve University
[EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha

__
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] Block comments in R?

2006-10-05 Thread hadley wickham
 Still gives you problems with nested comments. *IF* we want to go down
 that route, we should have directional symbols like

 
 dsaldfysdfk
 

What about heredocs? (http://en.wikipedia.org/wiki/Heredoc)

Hadley

__
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] Block comments in R?

2006-10-05 Thread Martin Maechler
 PD == Peter Dalgaard [EMAIL PROTECTED]
 on 05 Oct 2006 16:12:50 +0200 writes:

PD Seth Falcon [EMAIL PROTECTED] writes:
 Uwe Ligges [EMAIL PROTECTED] writes:
  Use an editor that comments out a whole block which is what I do all 
the 
  time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of
  them.
 
 This, of course, works.  The if(FALSE) approach does not because it
 requires the comment to be syntactically correct.  The multiline
 string trick is _almost_ there, the problem is that both  and ' are
 fairly common in text and having to escape that is a pain.  
 
 My wtf feature request is to add a multiline string delimiter ala
 Python like .

PD Still gives you problems with nested comments. *IF* we want to go down
PD that route, we should have directional symbols like

PD 
PD dsaldfysdfk
 

PD (which would likely drive version control users up the wall...)

i.e. all of R-core ...

PD or 

PD /# kdlsfk 
PD dfkjsdfl #/

PD (but those character combinations have a marginal risk of being used
PD in existing code).

or

##  
  ...
 ##

but I'm not yet fond of the general idea.
Martin


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

__
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] treatment effect at specific time point within mixed effects model

2006-10-05 Thread Afshartous, David
Hi Spencer,

Thanks for your reply.
I don't think this answers my question.

If I understand correctly, your model simply removes the intercept and
thus the intercept in fm1 is the same as the first time factor in 
fm1a ... but am I confused as to why the other coefficient estimates are
now different for the time factor if this is just a re-naming.  
The coefficient estimates for the interactions are the same for fm1 and
fm1a, as expected.

But my question relates to the signifcance of drug at a specific time
point,
e.g., time = 3.  The coeffecieint for say factor(time)3:drugP measures
the 
interaction of the effect of drug=P and time=3, which is not testing
what 
I want to test.  Based on the info below, I want to compare 3) versus
4).

1) time=1, Drug=I : Intercept
2) time=1, Drug=P : Intercept + DrugP
3) time=3, Drug=I : Intercept + factor(time)3
4) time=3, Drug=P : Intercept + factor(time)3 + DrugP +
factor(time)3:drugP

I'm surprised this isn't simple or maybe I'm missing something
competely.

thanks
dave





-Original Message-
From: Spencer Graves [mailto:[EMAIL PROTECTED] 
Sent: Wednesday, October 04, 2006 7:11 PM
To: Afshartous, David
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] treatment effect at specific time point within mixed
effects model

  Consider the following modification of your example: 

fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random =
list(Patient = ~ 1) )

summary(fm1a)
snip
 Value Std.Error DFt-value p-value
factor(Time)1   -0.6238472 0.7170161 10 -0.8700602  0.4047
factor(Time)2   -1.0155283 0.7170161 10 -1.4163256  0.1871
factor(Time)30.1446512 0.7170161 10  0.2017405  0.8442
factor(Time)40.7751736 0.7170161 10  1.0811105  0.3050
factor(Time)50.1566588 0.7170161 10  0.2184871  0.8314
factor(Time)60.0616839 0.7170161 10  0.0860286  0.9331
drugP1.2781723 1.0140139  3  1.2605077  0.2966
factor(Time)2:drugP  0.4034690 1.4340322 10  0.2813528  0.7842
factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104  0.6477
factor(Time)4:drugP -1.8149720 1.4340322 10 -1.2656424  0.2343
factor(Time)5:drugP -0.6416580 1.4340322 10 -0.4474502  0.6641
factor(Time)6:drugP -2.1396105 1.4340322 10 -1.4920240  0.1666

  Does this answer your question? 
  Hope this helps. 
  Spencer Graves

Afshartous, David wrote:
  
 All,

 The code below is for a pseudo dataset of repeated measures on 
 patients where there is also a treatment factor called drug.  Time 
 is treated as categorical.

 What code is necessary to test for a treatment effect at a single time

 point,
 e.g., time = 3?   Does the answer matter if the design is a crossover
 design,
 i.e, each patient received drug and placebo?

 Finally, what would be a good response to someone that suggests to do 
 a simple t-test (paired in crossover case) instead of the test above 
 within a mixed model?

 thanks!
 dave



 z = rnorm(24, mean=0, sd=1)
 time = rep(1:6, 4)
 Patient = rep(1:4, each = 6)
 drug = factor(rep(c(I, P), each = 6, times = 2)) ## P = placebo, I

 = Ibuprofen dat.new = data.frame(time, drug, z, Patient) data.grp = 
 groupedData(z ~ time | Patient, data = dat.new)
 fm1 = lme(z ~ factor(time) + drug + factor(time):drug, data = 
 data.grp, random = list(Patient = ~ 1) )

 __
 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] Block comments in R?

2006-10-05 Thread Duncan Murdoch
On 10/5/2006 10:57 AM, Martin Maechler wrote:

 
 ##  
   ...
  ##
 
 but I'm not yet fond of the general idea.
 Martin

Is that just because you don't need it, or that you see something 
objectionable in it?

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] Plotting text with lattice

2006-10-05 Thread Gabor Grothendieck
This requires R 2.4.0.  Its the same as my earlier example except
the labels= arg has been changed to labels=which.packet().


xyplot(x ~ x | g, data = data.frame(x = 1:12), panel = function(...) {
   panel.xyplot(...)
   panel.text(x=2, y=4, labels=which.packet())

})


On 10/5/06, Ritwik Sinha [EMAIL PROTECTED] wrote:
 Hi,

 On a related note, if I wanted to add different texts to different
 panels, should I stick to using trellis.focus() for each text in each
 panel? I cannot figure out a way to do it using a panel function.

 Ritwik.

 On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote:
  On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
   Here are two possibilities.  The first use trellis.focus/trellis/unfocus 
   to
   add text subsequent to drawing the xyplot and the second uses a
   custom panel:
  
   xyplot(x ~ x, data = data.frame(x = 1:10))
   trellis.focus(panel, 1, 1)
   panel.text(x=2, y=4, labels=Text)
   trellis.unfocus()
  
   xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) {
   panel.xyplot(...)
   panel.text(x=2, y=4, labels=Text)
   })
 
  Right, on a more general note, this is necessary because it is not clear 
  what
 
  ltext(x=2, y=4, labels=Text)
 
  should do for a multi-panel plot. On an even more general note, you
  will keep getting in trouble when using lattice if you (1) try to
  follow the incremental addition approach of standard graphics or (2)
  use the par() system in any way.
 
  Deepayan
 
   On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote:
Hello,
I've decided to take the leap and try my hand at the lattice package,
though I am getting stuck at what one might consider a trivial problem,
plotting text at a point in a graph. Apologies in advance if (that) I'm
missing something extremely basic.
   
Consider in base graphics:
 plot(1:10)
 text(2, 4, Text)
In the above you will see text centered at the point (2, 4) on the
graph.
   
Now I would like to try to do the same thing using the lattice package:
   
 xyplot(x ~ x, data = data.frame(x = 1:10))
 ltext(x=2, y=4, labels=Text)
 panel.text(x=2, y=4, labels=Text)
 grid.text(label=Text, x=2, y=4)
 grid.text(label=Text, x=unit(2, native), y=unit(4, native))
   
None of the above four commands puts the Text at the (2, 4) point on
the graph. Any help with this would be appreciated! Also, if I have more
than one panel and would like to place text at different points on
different panels how would I do this?
   
Also, note that I'm hoping to use text to label interesting points in a
levelplot, but am using the above xyplot as an example.
   
Thanks,
Robert
   
   
Robert McGehee
Quantitative Analyst
Geode Capital Management, LLC
53 State Street, 5th Floor | Boston, MA | 02109
Tel: 617/392-8396Fax:617/476-6389
mailto:[EMAIL PROTECTED]
   
   
   
This e-mail, and any attachments hereto, are intended for 
us...{{dropped}}
   
__
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.
  
 
 
  --
  http://www.stat.wisc.edu/~deepayan/
 
  __
  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.
 


 --
 Ritwik Sinha
 Graduate Student
 Epidemiology and Biostatistics
 Case Western Reserve University
 [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha

 __
 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] Block comments in R?

2006-10-05 Thread Barry Rowlingson
Seth Falcon wrote:

 
 My wtf feature request is to add a multiline string delimiter ala
 Python like .
 

  Anything that makes R more like Python syntax gets a +1 from me. How 
about getting rid of curly brackets for blocks and using indenting in R?

  Time to attack gram.y with a sledgehammer...

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.


[R] lmer BIC changes between output and anova

2006-10-05 Thread Darren M. Ward
list,

i am using lmer to fit multilevel models and trying to use anova to compare the 
models.  however, whenever i run the anova, the AIC, BIC and loglik are 
different from the original model output- as below.  can someone help me out 
with why this is happening?  (i'm hoping the output assocaited with the anova 
is right!).

thank you,

darren

 unconditional-lmer(log50 ~ 1 + (1 | Stream:Site) + (1|Stream), data)
 summary(unconditional)
Linear mixed-effects model fit by REML 
Formula: log50 ~ 1 + (1 | Stream:Site) + (1 | Stream) 
   Data: data 
AICBIC logLik MLdeviance REMLdeviance
 -138.8 -132.8  72.42 -150.4   -144.8


 nosection-lmer(log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) + 
 (1|Stream), data)
 summary(nosection)
Linear mixed-effects model fit by REML 
Formula: log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) +  (1 
| Stream) 
   Data: data 
AICBIC logLik MLdeviance REMLdeviance
 -140.8 -130.7   75.4 -168.9   -150.8


 anova(unconditional, nosection)
Data: data
Models:
unconditional: log50 ~ 1 + (1 | Stream:Site) + (1 | Stream)
nosection: log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) + 
unconditional: (1 | Stream)
  Df  AIC  BIC   logLik  Chisq Chi Df Pr(Chisq)
unconditional  3 -144.370 -138.294   75.185 
nosection  5 -158.861 -148.734   84.430 18.491  2  9.657e-05 ***

__
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] glm with nesting

2006-10-05 Thread Berton Gunter
Jeffrey:

Please... May I repeat what Peter Dalgaard already said: consult a local
statistician. The structure of your study is sufficiently complicated that
your stat 101 training is inadequate. Get professional help, which this list
is not set up to provide (though it often does, through the good offices and
patience of many wise contributors).

Bert Gunter
Genentech Nonclinical Statistics
South San Francisco, CA 94404

 

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 Jeffrey Stratford
 Sent: Thursday, October 05, 2006 7:46 AM
 To: [EMAIL PROTECTED]; r-help@stat.math.ethz.ch
 Subject: Re: [R] glm with nesting
 
 Harold and list,
 
 I've changed a few things since the last time so I'm really starting
 from scratch.  
 
 I start with
 
 bbmale - read.csv(c:\\eabl\\2004\\feathers\\male_feathers2.csv,
 header=TRUE)
 box -factor(box)
 chick - factor(chick)
 
 Here's a sample of the data
 
 box,chick,julian,cltchsz,mrtot,cuv,cblue,purbank,purban2,purba
 n1,pgrassk,pgrass2,pgrass1,grassdist,grasspatchk
 1,2,141,2,21.72290152,0.305723811,0.327178868,0.003813435,0.02
 684564,0.06896552,0.3282487,0.6845638,0.7586207,0,3.73
 4,1,164,4,18.87699007,0.281863299,0.310935559,0.06072162,0.208
 0537,0.06896552,0.01936052,0,0,323.1099,0.2284615
 4,2,164,4,19.64359348,0.294117388,0.316049817,0.06072162,0.208
 0537,0.06896552,0.01936052,0,0,323.1099,0.2284615
 7,1,118,4,13.48699876,0.303649408,0.31765218,0.3807568,0.43624
 16,0.6896552,0.06864183,0.03355705,0,94.86833,0.468
 12,1,180,4,21.42196378,0.289731361,0.317562294,0.09238011,0.13
 42282,0,0.2430127,0.8322148,1,0,1.199032
 12,2,180,4,18.79487905,0.286052077,0.316367349,0.09238011,0.13
 42282,0,0.2430127,0.8322148,1,0,1.199032
 12,3,180,4,12.83127682,0.260197475,0.292636914,0.09238011,0.13
 42282,0,0.2430127,0.8322148,1,0,1.199032
 15,1,138,4,20.07161467,0.287632782,0.318671887,0.07046477,0.03
 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
 15,2,138,4,17.61146256,0.305581768,0.315848051,0.07046477,0.03
 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
 15,3,138,4,20.36397134,0.271795667,0.30539683,0.07046477,0.033
 55705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
 15,4,138,4,20.81940158,0.269468041,0.304160648,0.07046477,0.03
 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818
 
 As you can see I have multiple boxes ( 70).  Sometimes I 
 have multiple
 chicks per box each having their own response  to mrtot, cuv, 
 and cblue
 but the same landscape variables for that box.  Chick number 
 is randomly
 assigned and is not an effect I'm interested in.  I'm really not
 interested in the box effect either.  I would like to know if 
 landscape
 affects the color of chicks (which may be tied into chick
 health/physiology).  We also know that chicks get bluer as the season
 progresses and that clutch size (cltchsz) has an effect so 
 I'm including
 that as covariates.  
 
 Hopefully, this clears things up a bit. 
 
 I do have the MASS and MEMS (Pineiro's) texts in hand. 
 
 Many thanks,
 
 Jeff
 
 __
 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] treatment effect at specific time point within mixedeffects model

2006-10-05 Thread Doran, Harold
Hi David:

In looking at your original post it is a bit difficult to ascertain
exactly what your null hypothesis was. That is, you want to assess
whether there is a treatment effect at time 3, but compared to what. I
think your second post clears this up. You should refer to pages 224-
225 of Pinhiero and Bates for your answer. This shows how to specify
contrasts.

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 Afshartous, David
 Sent: Thursday, October 05, 2006 11:08 AM
 To: Spencer Graves
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point 
 within mixedeffects model
 
 Hi Spencer,
 
 Thanks for your reply.
 I don't think this answers my question.
 
 If I understand correctly, your model simply removes the 
 intercept and thus the intercept in fm1 is the same as the 
 first time factor in fm1a ... but am I confused as to why the 
 other coefficient estimates are now different for the time 
 factor if this is just a re-naming.  
 The coefficient estimates for the interactions are the same 
 for fm1 and fm1a, as expected.
 
 But my question relates to the signifcance of drug at a 
 specific time point, e.g., time = 3.  The coeffecieint for 
 say factor(time)3:drugP measures the interaction of the 
 effect of drug=P and time=3, which is not testing what I want 
 to test.  Based on the info below, I want to compare 3) versus 4).
 
 1) time=1, Drug=I : Intercept
 2) time=1, Drug=P : Intercept + DrugP
 3) time=3, Drug=I : Intercept + factor(time)3
 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + 
 factor(time)3:drugP
 
 I'm surprised this isn't simple or maybe I'm missing 
 something competely.
 
 thanks
 dave
 
 
 
 
 
 -Original Message-
 From: Spencer Graves [mailto:[EMAIL PROTECTED]
 Sent: Wednesday, October 04, 2006 7:11 PM
 To: Afshartous, David
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point 
 within mixed effects model
 
   Consider the following modification of your example: 
 
 fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random 
 = list(Patient = ~ 1) )
 
 summary(fm1a)
 snip
  Value Std.Error DFt-value p-value
 factor(Time)1   -0.6238472 0.7170161 10 -0.8700602  0.4047
 factor(Time)2   -1.0155283 0.7170161 10 -1.4163256  0.1871
 factor(Time)30.1446512 0.7170161 10  0.2017405  0.8442
 factor(Time)40.7751736 0.7170161 10  1.0811105  0.3050
 factor(Time)50.1566588 0.7170161 10  0.2184871  0.8314
 factor(Time)60.0616839 0.7170161 10  0.0860286  0.9331
 drugP1.2781723 1.0140139  3  1.2605077  0.2966
 factor(Time)2:drugP  0.4034690 1.4340322 10  0.2813528  
 0.7842 factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104 
  0.6477 factor(Time)4:drugP -1.8149720 1.4340322 10 
 -1.2656424  0.2343 factor(Time)5:drugP -0.6416580 1.4340322 
 10 -0.4474502  0.6641 factor(Time)6:drugP -2.1396105 
 1.4340322 10 -1.4920240  0.1666
 
   Does this answer your question? 
   Hope this helps. 
   Spencer Graves
 
 Afshartous, David wrote:
   
  All,
 
  The code below is for a pseudo dataset of repeated measures on 
  patients where there is also a treatment factor called 
 drug.  Time 
  is treated as categorical.
 
  What code is necessary to test for a treatment effect at a 
 single time
 
  point,
  e.g., time = 3?   Does the answer matter if the design is a 
 crossover
  design,
  i.e, each patient received drug and placebo?
 
  Finally, what would be a good response to someone that 
 suggests to do 
  a simple t-test (paired in crossover case) instead of the 
 test above 
  within a mixed model?
 
  thanks!
  dave
 
 
 
  z = rnorm(24, mean=0, sd=1)
  time = rep(1:6, 4)
  Patient = rep(1:4, each = 6)
  drug = factor(rep(c(I, P), each = 6, times = 2)) ## P = 
 placebo, I
 
  = Ibuprofen dat.new = data.frame(time, drug, z, Patient) data.grp = 
  groupedData(z ~ time | Patient, data = dat.new)
  fm1 = lme(z ~ factor(time) + drug + factor(time):drug, data = 
  data.grp, random = list(Patient = ~ 1) )
 
  __
  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-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] unexpected behavior of boxplot(x, notch=TRUE, log=y)

2006-10-05 Thread bogdan romocea
A function I've been using for a while returned a surprising [to me,
given the data] error recently:
   Error in plot.window(xlim, ylim, log, asp, ...) :
   Logarithmic axis must have positive limits

After some digging I realized what was going on:
x - c(10460.97, 10808.67, 29499.98, 1, 35818.62, 48535.59, 1, 1,
   42512.1, 1627.39, 1, 7571.06, 21479.69, 25, 1, 16143.85, 12736.96,
   1, 7603.63, 1, 33155.24, 1, 1, 50, 3361.78, 1, 37781.84, 1, 1,
   1, 46492.05, 22334.88, 1, 1)
summary(x)
boxplot(x,notch=TRUE,log=y)  #unexpected
boxplot(x)  #ok
boxplot(x,log=y)  #ok
boxplot(x,notch=TRUE)  #aha

I can get around this, but thought that maybe boxplot() should be
adjusted to deal with something like this on its own.

Thank you,
b.

platform   i386-pc-mingw32
arch   i386
os mingw32
system i386, mingw32
status
major  2
minor  4.0
year   2006
month  10
day03
svn rev39566
language   R
version.string R version 2.4.0 (2006-10-03)

__
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] Fw: Bivariate Weibull distribution -- Copula

2006-10-05 Thread Jenny Stadt
I repost this in order for more responses. Thanks!




From: Jenny Stadt
Sent:  2006-10-04 16:49:33
To:  r-help@stat.math.ethz.ch
CC:  
Subject:  [R] Bivariate Weibull distribution -- Copula

Hi All,

I am struggling in a bivariate Weibull distribution although I searched 
R-Site-Help and found suggestion with Copula. Seems the maximum likelihood 
estimate is beyond what I can understand.

My case is: given two known marginal distribution (both are Weibull), and the 
correlation between them. How can I get the bivariate distribution based on 
these two? I wish if there is anybody has the experience in bivariate 
distribution could give me some hints. Thanks!

Jen.

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

[[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] Block comments in R?

2006-10-05 Thread Martin Maechler
 Duncan == Duncan Murdoch [EMAIL PROTECTED]
 on Thu, 05 Oct 2006 11:13:03 -0400 writes:

Duncan On 10/5/2006 10:57 AM, Martin Maechler wrote:
 
 ##  
 ...
 ##
 
 but I'm not yet fond of the general idea.
 Martin

Duncan Is that just because you don't need it, or that you see something 
Duncan objectionable in it?

The combination: It complicates the S language (slightly) which
I think is too expensive given that I don't see the need for
it (smart editors ...;  if(FALSE) { .. }).

{But my feelings are not very strong here probably mostly
 because I haven't thought very long about it.}

Martin

__
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] mixed models: correlation between fixed and random effects??

2006-10-05 Thread Bruno L. Giordano
Hello,
I built 4 mixed models using different data sets and standardized variables 
as predictors.

In all the models each of the fixed effects has an associated random effect
(same predictor).

What I find is that fixed effects with larger (absolute) standardized
parameter estimates have also a higher estimate of the related random
effect. In other words, the higher the average of the absolute value of the
BLUPs for a given standardized parameter, the higher its variance.

Is this a common situation or I am missing some additional normalization
factor necessary to compare the different random effects?

Thanks a lot,
Bruno



~~
Bruno L. Giordano, Ph.D.
CIRMMT
Schulich School of Music, McGill University
555 Sherbrooke Street West
Montréal, QC H3A 1E3
Canada
http://www.music.mcgill.ca/~bruno/

__
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] Block comments in R?

2006-10-05 Thread Gabor Grothendieck
There are two places that I find the current way it works to be
less than ideal although its not that bad either:

1 .when one wants to define strings that have quotes then
one must be careful to use double quoted strings to contain
single quotes and single quoted strings to contain double quotes
or if one needs both then one is into the complexities of multiple
backslashes.  Enhancing sQuote and dQuote might help somewhat
in this regards to allow them to be used in more situations but
having an extended syntax might be the best.

2. Sometime I like to define a file inline so that it is
defined like this:

Lines - 1 2
3 4
read.table(textConnection(Lines))

I haven't hit the size limit on the character string yet (I usually do this
with small config style files) but suspect that I may if I start using this
for data.  What is the limit anyways?

On 10/5/06, Martin Maechler [EMAIL PROTECTED] wrote:
  Duncan == Duncan Murdoch [EMAIL PROTECTED]
  on Thu, 05 Oct 2006 11:13:03 -0400 writes:

Duncan On 10/5/2006 10:57 AM, Martin Maechler wrote:

 ##  
 ...
 ##

 but I'm not yet fond of the general idea.
 Martin

Duncan Is that just because you don't need it, or that you see something
Duncan objectionable in it?

 The combination: It complicates the S language (slightly) which
 I think is too expensive given that I don't see the need for
 it (smart editors ...;  if(FALSE) { .. }).

 {But my feelings are not very strong here probably mostly
  because I haven't thought very long about it.}

 Martin

 __
 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] treatment effect at specific time point within mixedeffects model

2006-10-05 Thread Afshartous, David

Hi Harold,

Thanks for your response. 
I'll check out p.224 in PB, thanks.

The null hypothesis is that there is no difference between say
A=[time=3, drug=I] 
and B=[time=3, drug=P], or mu_A = mu_B.  If the study is a crossover
design, i.e., 
each patient receives drug=I and drug=P, I assume that a simple paried
t-test could
also be employed at time=3.  

However, I'd like to test this within a mixed effects model; With
respect to 3) and 4) below, 
it seems somewhat difficult to express this specific hypothesis in terms

of the model paramaters.   Ways in which this null are violated under 
the mixed effects models could be:

1) there is no interaction between time and Drug, i.e., there is a drug
effect but 
it is the same at all time points. (the specific interaction in 3) below
represents 
the shift of the effect of drug=P from time=1 to time=3 ... so the lack
of significance of
the paramater factor(time)3:drugP doesn't capture what I want)
2) there is neither interaction nor drug effect (variable Drug not
significant).

But both these violations are more general than my null; 
I think testing fixed effects 3) versus 4) below is what I want, but
this 
also seems strange since possibly the drug effect and drug:time
interaction as defined in the model
are signicant (with time=1 as the reference baseline).

Regardless, I assume I would need to employ coef() and vcov() to obtain
the needed 
info ... but I notice that coef() produces 4 values for the intercept of
fm1
below, does anyone know why this occurs?

I apologize if my explanation above is confusing, I've tried to make it
as clear as possible.

thanks again,
dave



-Original Message-
From: Doran, Harold [mailto:[EMAIL PROTECTED] 
Sent: Thursday, October 05, 2006 11:40 AM
To: Afshartous, David; Spencer Graves
Cc: r-help@stat.math.ethz.ch
Subject: RE: [R] treatment effect at specific time point within
mixedeffects model

Hi David:

In looking at your original post it is a bit difficult to ascertain
exactly what your null hypothesis was. That is, you want to assess
whether there is a treatment effect at time 3, but compared to what. I
think your second post clears this up. You should refer to pages 224-
225 of Pinhiero and Bates for your answer. This shows how to specify
contrasts.

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of Afshartous, 
 David
 Sent: Thursday, October 05, 2006 11:08 AM
 To: Spencer Graves
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point within 
 mixedeffects model
 
 Hi Spencer,
 
 Thanks for your reply.
 I don't think this answers my question.
 
 If I understand correctly, your model simply removes the intercept and

 thus the intercept in fm1 is the same as the first time factor in fm1a

 ... but am I confused as to why the other coefficient estimates are 
 now different for the time factor if this is just a re-naming.
 The coefficient estimates for the interactions are the same for fm1 
 and fm1a, as expected.
 
 But my question relates to the signifcance of drug at a specific time 
 point, e.g., time = 3.  The coeffecieint for say factor(time)3:drugP

 measures the interaction of the effect of drug=P and time=3, which is 
 not testing what I want to test.  Based on the info below, I want to 
 compare 3) versus 4).
 
 1) time=1, Drug=I : Intercept
 2) time=1, Drug=P : Intercept + DrugP
 3) time=3, Drug=I : Intercept + factor(time)3
 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + 
 factor(time)3:drugP
 
 I'm surprised this isn't simple or maybe I'm missing something 
 competely.
 
 thanks
 dave
 
 
 
 
 
 -Original Message-
 From: Spencer Graves [mailto:[EMAIL PROTECTED]
 Sent: Wednesday, October 04, 2006 7:11 PM
 To: Afshartous, David
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point within mixed 
 effects model
 
   Consider the following modification of your example: 
 
 fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = 
 list(Patient = ~ 1) )
 
 summary(fm1a)
 snip
  Value Std.Error DFt-value p-value
 factor(Time)1   -0.6238472 0.7170161 10 -0.8700602  0.4047
 factor(Time)2   -1.0155283 0.7170161 10 -1.4163256  0.1871
 factor(Time)30.1446512 0.7170161 10  0.2017405  0.8442
 factor(Time)40.7751736 0.7170161 10  1.0811105  0.3050
 factor(Time)50.1566588 0.7170161 10  0.2184871  0.8314
 factor(Time)60.0616839 0.7170161 10  0.0860286  0.9331
 drugP1.2781723 1.0140139  3  1.2605077  0.2966
 factor(Time)2:drugP  0.4034690 1.4340322 10  0.2813528
 0.7842 factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104
  0.6477 factor(Time)4:drugP -1.8149720 1.4340322 10
 -1.2656424  0.2343 factor(Time)5:drugP -0.6416580 1.4340322 10 
 -0.4474502  0.6641 factor(Time)6:drugP -2.1396105
 1.4340322 10 -1.4920240  0.1666
 
   Does this answer your question? 
   Hope this helps. 
   

Re: [R] (not so real) function (more like #include)

2006-10-05 Thread Gabor Grothendieck
I thought of a few more possibilities.

If options 1, 2 and 3 correspond to longfun, longfun2 and longfun3
then:

4. We could arrange it so that the variables are stored outside of
longfun and then longfun, f1 and f2 reference them symmetrically.

a - b - d - 1
longfun4 - function() {
   out - f1()
   # since a is not in longfun4 use - to set it
   a - a + 1
   f2() + out
}
f1 - f2 - function() a + b + d
longfun4()

Given the problem as stated I don't think this approach is really
the best but if the problem gets extended to one where there is
not just one longfun but several longfuns then it will not be feasible
to store all the variables in longfun since the other longfuns won't
be able to access them and this approach becomes desirable.

5. We could also use the proto package for this purpose placing
a, b, d, longfun, f1 and f2 into an object p.  If we wanted to be
able to override f1, f2, a, b and d in delegated subobjects and
be able to refer to f1 and f2 outside of longfun, i.e. outside of
the context of object p, then we would have
to refer to them as .$f1, .$f2, .$a, .$b and .$d and also pass
the target object as arg 1 of f1 and f2 but we have not bothered
with that here since it seems not to be a requirement:

library(proto)
p - proto(a = 1, b = 1, d = 1)
p$longfun - function(.) {
   out - f1()
   # since a is not in longfun4 use - to set it
   a - a + 1
   f2() + out
}
p$f1 - p$f2 - function() a + b + d
p$longfun() # 7

This is probably overkill here but it might be that if one thinks about
the design some more that there are some opportunities to define
subobjects of p in which case the ability to do inheritance would come
into play.

6. I think the R.oo package would work too although I am not sure
its as suitable here since it would require the definition of classes
which you don't really need.  proto lets you define objects directly
even without classes.


On 10/5/06, Meinhard Ploner [EMAIL PROTECTED] wrote:
 Thanks a lot! longfun2  longfun3 work perfect for my very
 untypical problem. As I have many local variables, usual functions
 with parameters are very uncomfortable, but the code you gave me is
 great!
 Meinhard


 On Oct 4, 2006, at 5:25 PM, Gabor Grothendieck wrote:

  longfun could just pass a, b and d to each of the individual
  functions and each of the individual functions could pass
  out back as a return value.
 
  f1 - f2 - function(a, b, d) a+b+d
 
  longfun1 - function() {
a - b - d - 1
out - f1(a, b, d)
out - f2(a, b, d) + out
out
  }
 
  longfun1() # 6
 
  If the problem is that a, b and d actually represent 100 variables,
  say,
  and its a nuisance to pass them all explicitly you could do this:
 
  f1 - f2 - function() with(parent.frame(), a + b + d)
 
  longfun2 - function() {
a - b - d - 1
out - f1()
out - f2() + out
out
  }
 
  longfun2() # 6
 
 
  or you could do it this way if you would prefer to have longfun
  control the scoping rather than having it done within each
  individual function:
 
 
  f1 - f2 - function() a+b+d
 
  longfun3 - function() {
a - b - d - 1
environment(f1) - environment(f2) - environment()
out - f1()
f2() + out
  }
 
  longfun3() # 6
 
 
  On 10/4/06, Meinhard Ploner [EMAIL PROTECTED] wrote:
  Hello R users,
 
  I have a very long function with parts of that looking like a list of
  jobs. The first part of the function defines a lot of local
  variables, which are used by the jobs. Now I extracted the job part
  und putted them into an external file, used by source(, local=T),
  which is more comfortable to me. Apart from that it gives me the
  opportunity that more users can work on the project. Is it possible
  to define a function for that code and passing to it the environment
  of f() without save(list=ls()) and load()  ?
 
  Thanks in advance
  Meinhard Ploner
 
  longfun - f() {
 
 ## lot of local variables:
 a - ...
 b - ...
 d - ...
 ...
 
 out - ...
 
 source(job1.txt, local=TRUE)   #changes out, uses a, b,
  d, ...
 
 source(job2.txt, local=TRUE)  # changes out, uses a, b,
  d, ...
 ...
  }
 
 
 
  version
 _
  platform   i386-apple-darwin8.6.1
  arch   i386
  os darwin8.6.1
  system i386, darwin8.6.1
  status
  major  2
  minor  3.1
  year   2006
  month  06
  day01
  svn rev38247
  language   R
  version.string Version 2.3.1 (2006-06-01)
 
  __
  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 

[R] any package can visualize original affymetrix data?

2006-10-05 Thread Baoqiang Cao
Dear All,

I'm thinking of visualizing the original binary data from 
affymetrix.  Would you recommend any package? Not necessarily limited 
to R packages. Thanks!

Best,
  -Cao

__
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] convert day of week from number to character and include in lm

2006-10-05 Thread Spencer Jones
All,

I am trying to include a day of week variable (1-7) in in a regression
model. I would like to have the day of week treated as a categorical
variable rather than a number

the code looks like

lm( dep ~ WKDY)

I know this is a basic question, but help would be appreciated

thanks

spencer

[[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] treatment effect at specific time point within mixedeffects model

2006-10-05 Thread Chuck Cleland
Afshartous, David wrote:
 Hi Harold,
 
 Thanks for your response. 
 I'll check out p.224 in PB, thanks.
 
 The null hypothesis is that there is no difference between say
 A=[time=3, drug=I] 
 and B=[time=3, drug=P], or mu_A = mu_B.  If the study is a crossover
 design, i.e., 
 each patient receives drug=I and drug=P, I assume that a simple paried
 t-test could
 also be employed at time=3.  
 
 However, I'd like to test this within a mixed effects model; With
 respect to 3) and 4) below, 
 it seems somewhat difficult to express this specific hypothesis in terms
 
 of the model paramaters.   Ways in which this null are violated under 
 the mixed effects models could be:
 
 1) there is no interaction between time and Drug, i.e., there is a drug
 effect but 
 it is the same at all time points. (the specific interaction in 3) below
 represents 
 the shift of the effect of drug=P from time=1 to time=3 ... so the lack
 of significance of
 the paramater factor(time)3:drugP doesn't capture what I want)

  If there is no evidence for an interaction between drug and time, do
you still want to ask about the drug effect at time=3, or would you then
want to ask about the time-averaged drug effect?

 2) there is neither interaction nor drug effect (variable Drug not
 significant).

 But both these violations are more general than my null; 
 I think testing fixed effects 3) versus 4) below is what I want, but
 this 
 also seems strange since possibly the drug effect and drug:time
 interaction as defined in the model
 are signicant (with time=1 as the reference baseline).

  If you fit a model with an intercept, main effects for drug and time,
and an interaction, would'nt the coefficient for the drug main effect
test the drug effect at a particular time?  Perhaps you only need to
change the contrasts for time so that time=3 is the reference category?

 Regardless, I assume I would need to employ coef() and vcov() to obtain
 the needed 
 info ... but I notice that coef() produces 4 values for the intercept of
 fm1
 below, does anyone know why this occurs?

  I think Harold was getting at the fact that you could get an estimate
of the drug effect at time=3 simply by setting the contrasts for time in
the right way.

 I apologize if my explanation above is confusing, I've tried to make it
 as clear as possible.
 
 thanks again,
 dave
 
 
 
 -Original Message-
 From: Doran, Harold [mailto:[EMAIL PROTECTED] 
 Sent: Thursday, October 05, 2006 11:40 AM
 To: Afshartous, David; Spencer Graves
 Cc: r-help@stat.math.ethz.ch
 Subject: RE: [R] treatment effect at specific time point within
 mixedeffects model
 
 Hi David:
 
 In looking at your original post it is a bit difficult to ascertain
 exactly what your null hypothesis was. That is, you want to assess
 whether there is a treatment effect at time 3, but compared to what. I
 think your second post clears this up. You should refer to pages 224-
 225 of Pinhiero and Bates for your answer. This shows how to specify
 contrasts.
 
 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of Afshartous, 
 David
 Sent: Thursday, October 05, 2006 11:08 AM
 To: Spencer Graves
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point within 
 mixedeffects model

 Hi Spencer,

 Thanks for your reply.
 I don't think this answers my question.

 If I understand correctly, your model simply removes the intercept and
 
 thus the intercept in fm1 is the same as the first time factor in fm1a
 
 ... but am I confused as to why the other coefficient estimates are 
 now different for the time factor if this is just a re-naming.
 The coefficient estimates for the interactions are the same for fm1 
 and fm1a, as expected.

 But my question relates to the signifcance of drug at a specific time 
 point, e.g., time = 3.  The coeffecieint for say factor(time)3:drugP
 
 measures the interaction of the effect of drug=P and time=3, which is 
 not testing what I want to test.  Based on the info below, I want to 
 compare 3) versus 4).

 1) time=1, Drug=I : Intercept
 2) time=1, Drug=P : Intercept + DrugP
 3) time=3, Drug=I : Intercept + factor(time)3
 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + 
 factor(time)3:drugP

 I'm surprised this isn't simple or maybe I'm missing something 
 competely.

 thanks
 dave





 -Original Message-
 From: Spencer Graves [mailto:[EMAIL PROTECTED]
 Sent: Wednesday, October 04, 2006 7:11 PM
 To: Afshartous, David
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] treatment effect at specific time point within mixed 
 effects model

   Consider the following modification of your example: 

 fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = 
 list(Patient = ~ 1) )

 summary(fm1a)
 snip
  Value Std.Error DFt-value p-value
 factor(Time)1   -0.6238472 0.7170161 10 -0.8700602  0.4047
 factor(Time)2   -1.0155283 0.7170161 10 -1.4163256  0.1871
 

Re: [R] Block comments in R?

2006-10-05 Thread BBands
+1 for Python style comments.


your comments here
and some more
...
and some related info


jab
-- 
John Bollinger, CFA, CMT
www.BollingerBands.com

If you advance far enough, you arrive at the beginning.

__
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] The W statistic in wilcox.exact

2006-10-05 Thread Jue.Wang2
Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as 
indicated below.

12 is the rank sum of group 0 of x, which is the linear statistic computed by 
wilcox_test.

y-c(1,2,3,4,5)
x-c(1,1,0,0,0)


(a) wilcox.exact

wilcox.exact(y~x)
Exact Wilcoxon rank sum test
data:  y by x 
W = 6, p-value = 0.2
alternative hypothesis: true mu is not equal to 0 


(b) wilcox_test

tt-wilcox_test(y~factor(x),distribution=exact)
statistic(tt,linear)
   
0 12



Jue Wang, Biostatistician
Contracted Position for Preclinical  Research Biostatistics
PrO Unlimited
(908) 231-3022

__
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] Multivariate AR - prediction

2006-10-05 Thread Fleischer Alexander 0969 SPI
Hi,

does anybody know how to predict a multivariate AR within R?
If I just estimate a multi AR-object and plug it into predict I get an
error from the aperm - just works for univariates.

thx 
alex

__
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 text with lattice

2006-10-05 Thread Gabor Grothendieck
I seem to have omitted g

library(lattice)
xyplot(x ~ x | g, data = data.frame(x = 1:12, g = gl(3,4)), panel =
function(...) {
   panel.xyplot(...)
   panel.text(x=2, y=4, labels=which.packet())
})


On 10/5/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
 This requires R 2.4.0.  Its the same as my earlier example except
 the labels= arg has been changed to labels=which.packet().


 xyplot(x ~ x | g, data = data.frame(x = 1:12), panel = function(...) {
   panel.xyplot(...)
   panel.text(x=2, y=4, labels=which.packet())

 })


 On 10/5/06, Ritwik Sinha [EMAIL PROTECTED] wrote:
  Hi,
 
  On a related note, if I wanted to add different texts to different
  panels, should I stick to using trellis.focus() for each text in each
  panel? I cannot figure out a way to do it using a panel function.
 
  Ritwik.
 
  On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote:
   On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
Here are two possibilities.  The first use 
trellis.focus/trellis/unfocus to
add text subsequent to drawing the xyplot and the second uses a
custom panel:
   
xyplot(x ~ x, data = data.frame(x = 1:10))
trellis.focus(panel, 1, 1)
panel.text(x=2, y=4, labels=Text)
trellis.unfocus()
   
xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) {
panel.xyplot(...)
panel.text(x=2, y=4, labels=Text)
})
  
   Right, on a more general note, this is necessary because it is not clear 
   what
  
   ltext(x=2, y=4, labels=Text)
  
   should do for a multi-panel plot. On an even more general note, you
   will keep getting in trouble when using lattice if you (1) try to
   follow the incremental addition approach of standard graphics or (2)
   use the par() system in any way.
  
   Deepayan
  
On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote:
 Hello,
 I've decided to take the leap and try my hand at the lattice package,
 though I am getting stuck at what one might consider a trivial 
 problem,
 plotting text at a point in a graph. Apologies in advance if (that) 
 I'm
 missing something extremely basic.

 Consider in base graphics:
  plot(1:10)
  text(2, 4, Text)
 In the above you will see text centered at the point (2, 4) on the
 graph.

 Now I would like to try to do the same thing using the lattice 
 package:

  xyplot(x ~ x, data = data.frame(x = 1:10))
  ltext(x=2, y=4, labels=Text)
  panel.text(x=2, y=4, labels=Text)
  grid.text(label=Text, x=2, y=4)
  grid.text(label=Text, x=unit(2, native), y=unit(4, native))

 None of the above four commands puts the Text at the (2, 4) point on
 the graph. Any help with this would be appreciated! Also, if I have 
 more
 than one panel and would like to place text at different points on
 different panels how would I do this?

 Also, note that I'm hoping to use text to label interesting points in 
 a
 levelplot, but am using the above xyplot as an example.

 Thanks,
 Robert


 Robert McGehee
 Quantitative Analyst
 Geode Capital Management, LLC
 53 State Street, 5th Floor | Boston, MA | 02109
 Tel: 617/392-8396Fax:617/476-6389
 mailto:[EMAIL PROTECTED]



 This e-mail, and any attachments hereto, are intended for 
 us...{{dropped}}

 __
 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.
   
  
  
   --
   http://www.stat.wisc.edu/~deepayan/
  
   __
   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.
  
 
 
  --
  Ritwik Sinha
  Graduate Student
  Epidemiology and Biostatistics
  Case Western Reserve University
  [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
 
  __
  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 

[R] [R-pkgs] new package proptest

2006-10-05 Thread David Kraus
Dear R users,

The package proptest (http://www.davidkraus.net/proptest/) has been
recently released to CRAN.

The package provides functions for testing the proportional hazards
assumption in the Cox model for right censored survival data. Two types of
tests for identifying nonproportional covariates are implemented:

* data-driven Neyman type smooth tests (using both nested and all-subsets
variant of the BIC),

* tests based on the score process (KS, CM, and AD type, with the LWY
simulation).

Bug reports, comments or suggestions are welcome.

Regards,

David.


-- 

David Kraus

Institute of Information Theory and Automation
Pod Vodarenskou vezi 4
18208 Praha 8
Czechia



Charles University in Prague
Department of Statistics
Sokolovska 83
18675 Praha 8
Czechia

kraus #at# karlin.mff.cuni.cz
http://www.davidkraus.net/

___
R-packages mailing list
R-packages@stat.math.ethz.ch
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.


Re: [R] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX

2006-10-05 Thread Jens Scheidtmann
Anupam Tyagi [EMAIL PROTECTED] writes:

[...]
 What is the best format to save R graphics for inclusion into a
 LaTeX documents?

When using pdflatex use pdf for graphics as reference format.  Using
ps2pdf or some such may have some problems when it comes to alpha
channels and transparency.

[...]
 I will be grateful if someone can share their experience.

Indispensable is having a good editing cycle.  I.e. compile the latex,
jump to the position where you are editing the file in the preview,
double click somewhere in the preview, move to that position in your
editor (or at least near that position).  As I am not aware of any
tools doing this with pdf, you may want to render dvi instead.  When
rendering DVI, you probably want to generate ps.

My 2c,

Jens

__
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] littler release 0.0.6

2006-10-05 Thread Jeffrey Horner
For those of you who are littler fans, you'll be pleased to know that 
we've already decommissioned version 0.0.6 in light of the new and 
improved 0.0.7 (actually 0.0.6 was broken because of a missing file).

You can get it from the links below.

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

Dirk Eddelbuettel wrote:
 What ?
 --
 
We are pleased to announce version 0.0.6 of littler 
 
 
 What's new ?
 
 
This version includes a bug fix or two as well as a number of small
enhancements to the documentation.
 
For OS X and the r/R confusion, our recommended suggestion is to call
configure using either the --program-suffix=X  or  --program-prefix=Y
option to have the binary and manual page renamed when 'make install'
runs. In this example, r would be renamed to rX and Yr, resptively.
Suitable choices of X and Y are left to the user.

 
 Where ?
 ---
 
You can find the source tarball at both
 
   http://dirk.eddelbuettel.com/code/littler/
 
and 
 
   http://biostat.mc.vanderbilt.edu/LittleR
 
 
 
 Comments are welcome, as are are suggestions, bug fixes, or patches.
 
   - Jeffrey Horner [EMAIL PROTECTED]
   - Dirk Eddelbuettel [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.


[R] Help on plot multiple plot(hexbin) in the same page

2006-10-05 Thread Li, Sue
Dear R-users,

I try to plot multiple plots of hexbin in the same page. However,  par
does not work when hexbin is used. Neither do xlim and ylim. Any ideas?


Thanks,

Sue

[[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] Writing Text into Plots

2006-10-05 Thread Lothar Botelho-Machado
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hello,


I have a simple question on the text() method in plots, e.g.:

text(3,4,adj=1,cex=1.0, expression(alpha == beta))

I know there exists a lot more like frac(), etc which could be used for
expression. But a help(frac) doesn't return any results - where do I
have to look for a documentation of possible text commands?

More in detail I'm searching for how to set an index like mu_{1} and
mu_{2} ? And how to get something like a bar over a character like \bar x ?


TIA,
 Lothar
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.5 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFFJWX0HRf7N9c+X7sRAtyiAJ4yssSZtt/DGzWVTfVI2qvzoO9FigCfe7Fm
NTOUZN7jL/KdthkNEzM8S0M=
=tBvl
-END PGP SIGNATURE-

__
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] matrix multiplication

2006-10-05 Thread Ferdinand Alimadhi
Is this what you want ?

  A1-matrix(c(1,0,0,0,0,0,0,0,0),3)
  A2-matrix(c(1,0,0,0,1,0,0,0,0),3)
  A3-matrix(c(1,0,0,0,1,0,0,1,0),3)
  X - matrix(1:24,8)
  XX-list()
  for(i in 1:3){
+ XX[[i]]-X%*%A[[i]]
+ }
  XX
[[1]]
 [,1] [,2] [,3]
[1,]100
[2,]200
[3,]300
[4,]400
[5,]500
[6,]600
[7,]700
[8,]800

[[2]]
 [,1] [,2] [,3]
[1,]190
[2,]2   100
[3,]3   110
[4,]4   120
[5,]5   130
[6,]6   140
[7,]7   150
[8,]8   160

[[3]]
 [,1] [,2] [,3]
[1,]199
[2,]2   10   10
[3,]3   11   11
[4,]4   12   12
[5,]5   13   13
[6,]6   14   14
[7,]7   15   15
[8,]8   16   16

 



Ya-Hsiu Chuang wrote:

Dear all,

I have 2 matrices, one is a 8x3 matrix, called X; the other matrix is a 3x3 
indicator matrix with the diagonal element as 0 or 1. when a variable is 
included in the model, the corresponding diagonal element is 1, otherwise, it 
is 0.  Let A be a set of matrices that contain the possible indicator matrix. 
e.g.,
A= [A1, A2, A3], 
where A1= [1,0,0;0,0,0,0,0,0],
A2 =[1,0,0;0,1,0,0,0,0], 
A3 =[1,0,0;0,0,0,0,1,0]

In order to derive the new X matries depending on the indicator matrices, I 
use for loops to multiply both X and A as following

p-3
for ( i in 1:p){
  XX- X%*%A[i]
}

However, it only shows the result when i=p. How can I derive the results which 
include all possible value of XX[i], i=1,..,p, rather than just i=p

thanks for any help

Ya-Hsiu



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

  



-- 
Ferdinand Alimadhi
Programmer / Analyst
Harvard University
The Institute for Quantitative Social Science
(617) 496-0187
[EMAIL PROTECTED]
www.iq.harvard.edu


[[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] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX

2006-10-05 Thread Mike Prager
Jens Scheidtmann [EMAIL PROTECTED] wrote:
 
 Indispensable is having a good editing cycle.  I.e. compile the latex,
 jump to the position where you are editing the file in the preview,
 double click somewhere in the preview, move to that position in your
 editor (or at least near that position).  As I am not aware of any
 tools doing this with pdf, you may want to render dvi instead.  

Agreed.  As a point of interest, the newer Windows versions of
vTeX (by MicroPress) will do this (in both directions) with PDF.
This is a commercial TeX distribution with value added.  They
have a free distribution suitable for Unix or Linux, but I doubt
that it does that particular trick.

-- 
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.

__
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 to get the function names

2006-10-05 Thread Søren Højsgaard
I've defined the function
 
getFunNames - function(FUN){
  if (!is.list(FUN)) 
fun.names - paste(deparse(substitute(FUN)), collapse =  )
  else
fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a)))
  fun.names
}

which gives what I want :
 getFunNames(mean)
[1] mean
 getFunNames(ff)
[1] ff
 getFunNames(c(mean,ff))
[1] mean ff  
 
If I call this within a function, things go wrong:
1] FUN
 foo(ff)
[1] FUN
 foo(c(mean,ff))
Error in substitute(FUN)[-1] : object is not subsettable

Obviously there are some things (quite a few things) which I have not 
understood. Can anyone help?
Thanks
Søren

__
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] The W statistic in wilcox.exact

2006-10-05 Thread Christos Hatzis
Probably because of the offset:

U = W - n*(n+1)/2 

In your example, W=12 (=3+4+5) as reported by wilcox.test.  The offset is 6
(=3*4/2) and therefore U=6.
I am not certain as I haven't installed the exactRankTests package, but it
seems that wilcox.exact reports U instead of W.

-Christos Hatzis

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
[EMAIL PROTECTED]
Sent: Thursday, October 05, 2006 1:35 PM
To: r-help@stat.math.ethz.ch
Subject: [R] The W statistic in wilcox.exact

Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as
indicated below.

12 is the rank sum of group 0 of x, which is the linear statistic computed
by wilcox_test.

y-c(1,2,3,4,5)
x-c(1,1,0,0,0)


(a) wilcox.exact

wilcox.exact(y~x)
Exact Wilcoxon rank sum test
data:  y by x
W = 6, p-value = 0.2
alternative hypothesis: true mu is not equal to 0 


(b) wilcox_test

tt-wilcox_test(y~factor(x),distribution=exact)
statistic(tt,linear)
   
0 12



Jue Wang, Biostatistician
Contracted Position for Preclinical  Research Biostatistics
PrO Unlimited
(908) 231-3022

__
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] xyplot

2006-10-05 Thread Osman Al-Radi
Hi,

for the data below:
time-c(rep(1:10,5))
y-time+rnorm(50,5,2)
subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10))
group-c(rep('A',30),rep('B',20))
df-data.frame(subject,group,time,y)

I'd like to produce a plot with a single pannel with two loess curves one
for each group. the code below does that but on two different pannels, I'd
like to have both curves on the same pannel with different colors.

xyplot(y~time|group, groups=subject,
panel=function(...){
panel.loess(...)
}
,data=df)


-- 
Osman O. Al-Radi, MD, MSc, FRCSC
Fellow, Cardiovascular Surgery
The Hospital for Sick Children
University of Toronto, Canada

[[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] Writing Text into Plots

2006-10-05 Thread Duncan Murdoch
On 10/5/2006 4:10 PM, Lothar Botelho-Machado wrote:
 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1
 
 Hello,
 
 
 I have a simple question on the text() method in plots, e.g.:
 
 text(3,4,adj=1,cex=1.0, expression(alpha == beta))
 
 I know there exists a lot more like frac(), etc which could be used for
 expression. But a help(frac) doesn't return any results - where do I
 have to look for a documentation of possible text commands?
 
 More in detail I'm searching for how to set an index like mu_{1} and
 mu_{2} ? And how to get something like a bar over a character like \bar x ?

See ?plotmath, and demo(plotmath).

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] Writing Text into Plots

2006-10-05 Thread Marc Schwartz
On Thu, 2006-10-05 at 22:10 +0200, Lothar Botelho-Machado wrote:
 Hello,
 
 
 I have a simple question on the text() method in plots, e.g.:
 
 text(3,4,adj=1,cex=1.0, expression(alpha == beta))
 
 I know there exists a lot more like frac(), etc which could be used for
 expression. But a help(frac) doesn't return any results - where do I
 have to look for a documentation of possible text commands?
 
 More in detail I'm searching for how to set an index like mu_{1} and
 mu_{2} ? And how to get something like a bar over a character like \bar x ?

See ?plotmath, which BTW, is listed in the See Also section of ?text

HTH,

Marc Schwartz

__
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] which command

2006-10-05 Thread AgusSusanto
I obtained error messages when I run these commands in UNIX, but I 
obtained correct result when I run these command in WINDOWS.  Can 
somebody point out the problem and give the solution. Thanks.

  dt-read.table(file=Fall.dat)
  dim(dt)
[1] 19415
  table(dt$V2)
  0   1   2   3
220 989 639  93
 
  Favg-as.matrix(c(1:max(dt$V2)))
 
  for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])}
  Favg
 [,1]
[1,]   NA
[2,]   NA
[3,]   NA
Warning messages:
argument is not numeric or logical: returning NA (etc)

__
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] Matrix input problem

2006-10-05 Thread Bill Wyatt

Hi,

 Included is an R.script that came from a much large date set being 
read in to R from a .txt file. Why is it that the first line codes with 
an error, but the second line works fine? Is there some line length 
limit in R? This happens at random places all through my data. Any help 
would be appreciated.



Bill Wyatt
Associate Instructor
Ergonomics Graduate Student
Department of Kinesiology
School of Health, Physical Education, and Recreation
Indiana University
O:(812)856-5924
[EMAIL PROTECTED]





//orginal line

b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-67.3
 5,-67.35,-61.8,-61.8,-55.65,-55.65))

//copy of line with end )) moved one data point forward

b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-67.3
 5,-67.35,-61.8,-61.8,-55.65)) ,-55.65
__
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] xyplot

2006-10-05 Thread Sundar Dorai-Raj


Osman Al-Radi said the following on 10/5/2006 3:43 PM:
 Hi,
 
 for the data below:
 time-c(rep(1:10,5))
 y-time+rnorm(50,5,2)
 subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10))
 group-c(rep('A',30),rep('B',20))
 df-data.frame(subject,group,time,y)
 
 I'd like to produce a plot with a single pannel with two loess curves one
 for each group. the code below does that but on two different pannels, I'd
 like to have both curves on the same pannel with different colors.
 
 xyplot(y~time|group, groups=subject,
 panel=function(...){
 panel.loess(...)
 }
 ,data=df)
 
 

Simple:

xyplot(y ~ time | group, data = df, groups = subject,
panel = function(...) {
  panel.superpose(...)
  panel.superpose(panel.groups = panel.loess, ...)
})

HTH,

--sundar

__
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] which command

2006-10-05 Thread AgusSusanto
I obtained error messages when I run these commands in UNIX, but I 
obtained correct result when I run these command in WINDOWS.  Can 
somebody point out the problem and give the solution. Thanks.

  dt-read.table(file=Fall.dat)
  dim(dt)
[1] 19415
  table(dt$V2)
  0   1   2   3
220 989 639  93
 
  Favg-as.matrix(c(1:max(dt$V2)))
 
  for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])}
  Favg
 [,1]
[1,]   NA
[2,]   NA
[3,]   NA
Warning messages:
argument is not numeric or logical: returning NA (etc)


Cheers,
@

__
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 to get the function names

2006-10-05 Thread Duncan Murdoch
On 10/5/2006 4:41 PM, Søren Højsgaard wrote:
 I've defined the function
  
 getFunNames - function(FUN){
   if (!is.list(FUN)) 
 fun.names - paste(deparse(substitute(FUN)), collapse =  )
   else
 fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a)))
   fun.names
 }
 
 which gives what I want :
 getFunNames(mean)
 [1] mean
 getFunNames(ff)
 [1] ff
 getFunNames(c(mean,ff))
 [1] mean ff  
  
 If I call this within a function, things go wrong:
 1] FUN
 foo(ff)
 [1] FUN
 foo(c(mean,ff))
 Error in substitute(FUN)[-1] : object is not subsettable
 
 Obviously there are some things (quite a few things) which I have not 
 understood. Can anyone help?

I don't think you'll be able to do what you want.  The problem is that R 
objects don't know their own name.  substitute just gives you the 
unevaluated expression that was passed as an arg; so getFunNames is 
reporting on the expression in foo() that was used when it was called.

You'll need to call substitute directly on the args to get the 
expressions, rather than passing them to another function to do that.

In C code you have a bit more flexibility for non-standard evaluation, 
but not in R code.

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] xyplot

2006-10-05 Thread Deepayan Sarkar
On 10/5/06, Osman Al-Radi [EMAIL PROTECTED] wrote:
 Hi,

 for the data below:
 time-c(rep(1:10,5))
 y-time+rnorm(50,5,2)
 subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10))
 group-c(rep('A',30),rep('B',20))
 df-data.frame(subject,group,time,y)

 I'd like to produce a plot with a single pannel with two loess curves one
 for each group. the code below does that but on two different pannels, I'd
 like to have both curves on the same pannel with different colors.

 xyplot(y~time|group, groups=subject,
 panel=function(...){
 panel.loess(...)
 }
 ,data=df)

Well, the fix for your code is

xyplot(y ~ time | group, data = df, groups = subject,

   panel = panel.superpose,
   panel.groups = function(...) panel.loess(...)

   )

Which, incidentally is effectively equivalent to

xyplot(y ~ time | group, data = df, groups = subject, type = smooth)

However, this is not what you describe as your goal. If you want a
single smooth for each group combining all subjects within that group,
use something like

xyplot(y ~ time, data = df, groups = group. ...

The panel part does not change.

-Deepayan

__
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] Matrix input problem

2006-10-05 Thread Duncan Murdoch
On 10/5/2006 4:52 PM, Bill Wyatt wrote:
 Hi,
 
   Included is an R.script that came from a much large date set being 
 read in to R from a .txt file. Why is it that the first line codes with 
 an error, but the second line works fine? 

I get syntax errors on both, due to the missing comma at the end of the 
very long line.

Is there some line length
 limit in R? 

Yes, there's a 1024 character length limit (which is documented in the 
Intro to R manual in the current release; not sure how long it's been 
there).  Why not put in some line breaks?  R source is supposed to be 
human readable, not just machine readable.

Duncan Murdoch

This happens at random places all through my data. Any help
 would be appreciated.
 
 
 Bill Wyatt
 Associate Instructor
 Ergonomics Graduate Student
 Department of Kinesiology
 School of Health, Physical Education, and Recreation
 Indiana University
 O:(812)856-5924
 [EMAIL PROTECTED]
 
 
 
 
 
 
 
 
 
 //orginal line
 
 b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-6
7.3
  5,-67.35,-61.8,-61.8,-55.65,-55.65))
 
 //copy of line with end )) moved one data point forward
 
 b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-6
7.3
  5,-67.35,-61.8,-61.8,-55.65)) ,-55.65
 
 
 
 
 __
 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] How to get the function names

2006-10-05 Thread Jerome Asselin
On Thu, 2006-10-05 at 22:41 +0200, Søren Højsgaard wrote:
 I've defined the function
  
 getFunNames - function(FUN){
   if (!is.list(FUN)) 
 fun.names - paste(deparse(substitute(FUN)), collapse =  )
   else
 fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a)))
   fun.names
 }

Hi,

Try this:
 getFunNames - function(x)
+   sapply(as.list(sys.call()[[2]][-1]),as.character)
 getFunNames(c(mean,ff))
[1] mean ff
 foo - function() getFunNames(c(mean,ff))
 foo()
[1] mean ff

HTH,
Jerome

-- 
Jerome Asselin, M.Sc., Agent de recherche, RHCE
CHUM -- Centre de recherche
3875 rue St-Urbain, 3e etage // Montreal QC  H2W 1V1
Tel.: 514-890-8000 Poste 15914; Fax: 514-412-7106

__
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] Writing Text into Plots

2006-10-05 Thread Lothar Botelho-Machado
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Duncan Murdoch wrote:
 On 10/5/2006 4:10 PM, Lothar Botelho-Machado wrote:
 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1

 Hello,


 I have a simple question on the text() method in plots, e.g.:

 text(3,4,adj=1,cex=1.0, expression(alpha == beta))

 I know there exists a lot more like frac(), etc which could be used for
 expression. But a help(frac) doesn't return any results - where do I
 have to look for a documentation of possible text commands?

 More in detail I'm searching for how to set an index like mu_{1} and
 mu_{2} ? And how to get something like a bar over a character like
 \bar x ?
 
 See ?plotmath, and demo(plotmath).
 
 Duncan Murdoch
 

Thanx a lot!!
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.5 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFFJXsfHRf7N9c+X7sRAkfGAKCUMsskOqH5pDA+ToiQRDAEaxpw9gCghDjd
X5wETK17s21+u255hQYcqWU=
=1JXU
-END PGP SIGNATURE-

__
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] The W statistic in wilcox.exact

2006-10-05 Thread Christos Hatzis
Jue,

On a second look, it appears that wilcox.test does report the
offset-adjusted statistic U, as also mentioned in the help page.

wilcox.test returns W=6 (instead of 12 as your example showed, unless
wilcox_test is a different function).

 wilcox.test( 1:5 ~ c(1,1,0,0,0) )$statistic  # or wilcox.test( 1:5 ~
factor(c(1,1,0,0,0)) )$statistic
W 
6  

So there does not appear to be a difference between the two methods.  Did I
miss something?

-Christos

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Christos Hatzis
Sent: Thursday, October 05, 2006 4:41 PM
To: [EMAIL PROTECTED]; r-help@stat.math.ethz.ch
Subject: Re: [R] The W statistic in wilcox.exact

Probably because of the offset:

U = W - n*(n+1)/2 

In your example, W=12 (=3+4+5) as reported by wilcox.test.  The offset is 6
(=3*4/2) and therefore U=6.
I am not certain as I haven't installed the exactRankTests package, but it
seems that wilcox.exact reports U instead of W.

-Christos Hatzis

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
[EMAIL PROTECTED]
Sent: Thursday, October 05, 2006 1:35 PM
To: r-help@stat.math.ethz.ch
Subject: [R] The W statistic in wilcox.exact

Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as
indicated below.

12 is the rank sum of group 0 of x, which is the linear statistic computed
by wilcox_test.

y-c(1,2,3,4,5)
x-c(1,1,0,0,0)


(a) wilcox.exact

wilcox.exact(y~x)
Exact Wilcoxon rank sum test
data:  y by x
W = 6, p-value = 0.2
alternative hypothesis: true mu is not equal to 0 


(b) wilcox_test

tt-wilcox_test(y~factor(x),distribution=exact)
statistic(tt,linear)
   
0 12



Jue Wang, Biostatistician
Contracted Position for Preclinical  Research Biostatistics PrO Unlimited
(908) 231-3022

__
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-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] which command

2006-10-05 Thread Don MacQueen
Are you sure Fall.dat is identical on both platforms, and that dt is also?
Did you intend to ignore the cases where dt$V2 equals zero?

Try, for example,
   table(which(dt[,2] == 0))   ## also ==1, ==2, ==3
   unique(  dt[which(dt[,2] == 0),5]## also ==1, ==2, ==3
   unique( dt$V5 )
  adding na.rm=TRUE to mean()


Also, try this:
   tmp - lapply( split( dt$V5, dt$V2) , mean )
   Favg - unlist(tmp)

-Don

At 4:57 PM -0400 10/5/06, AgusSusanto wrote:
I obtained error messages when I run these commands in UNIX, but I
obtained correct result when I run these command in WINDOWS.  Can
somebody point out the problem and give the solution. Thanks.

   dt-read.table(file=Fall.dat)
   dim(dt)
[1] 19415
   table(dt$V2)
   0   1   2   3
220 989 639  93
  
   Favg-as.matrix(c(1:max(dt$V2)))
  
   for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])}
   Favg
  [,1]
[1,]   NA
[2,]   NA
[3,]   NA
Warning messages:
argument is not numeric or logical: returning NA (etc)

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


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

__
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] emacs-ess tab indentation

2006-10-05 Thread Ferdinand Alimadhi
Hi,

I have a 2 tab-identation in my emacs-ess and I would like to make it 8.  Can 
somebody help me with this ?

thanks 
Johan

__
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] which command

2006-10-05 Thread Gavin Simpson
On Thu, 2006-10-05 at 16:57 -0400, AgusSusanto wrote:
 I obtained error messages when I run these commands in UNIX, but I 
 obtained correct result when I run these command in WINDOWS.  Can 
 somebody point out the problem and give the solution. Thanks.

Not without Fall.dat I suspect. please read the posting guide and act
accordingly. (And why post twice?)

I suspect some of your data are missing in Fall.dat (or end up missing
because you didn't check what had been read in by R) and that this has
propagated warnings in say mean, when called within na.rm = TRUE. But
I'm guessing as you fail to provide a reproducible example.

Why this is different between Unix/windows versions of R, I don't know,
but are they running the same versions of R?

G

 
   dt-read.table(file=Fall.dat)
   dim(dt)
 [1] 19415
   table(dt$V2)
   0   1   2   3
 220 989 639  93
  
   Favg-as.matrix(c(1:max(dt$V2)))
  
   for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])}
   Favg
  [,1]
 [1,]   NA
 [2,]   NA
 [3,]   NA
 Warning messages:
 argument is not numeric or logical: returning NA (etc)
 
 
 Cheers,
 @
 
 __
 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.
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 *Note new Address and Fax and Telephone numbers from 10th April 2006*
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Gavin Simpson [t] +44 (0)20 7679 0522
ECRC  [f] +44 (0)20 7679 0565
UCL Department of Geography
Pearson Building  [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street
London, UK[w] http://www.ucl.ac.uk/~ucfagls/cv/
WC1E 6BT  [w] http://www.ucl.ac.uk/~ucfagls/
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
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] [OT] testing for synchronicity

2006-10-05 Thread Andrew Robinson
Greetings, friends in the R community,

this is an OT question about statistics. Given four time series of
events, what possibilities do I have to test for synchronicity?

e.g.

times - data.frame(year=   c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
event.1=c(1, 0, 0, 1, 2, 4, 1, 0, 0, 0),
event.2=c(0, 0, 0, 1, 0, 2, 1, 0, 0, 0),
event.3=c(1, 0, 0, 0, 1, 2, 4, 1, 0, 0),
event.4=c(0, 1, 0, 1, 0, 0, 1, 1, 1, 0))

I have about 100 years of each, and my null hypothesis is that the
events are not synchronous. 

Right now my thinking is to focus on pairwise comparisons:

1) ignore the magnitude and convert the series to binary  
 
2) the sum of the product of the events of any two years is then the
   number of overlapping occurences.

3) I think that, in the absence of temporal autocorrelation, I could
   assume that this sum is hypergeometrically distributed.

4) I can test for statistical significance of this number by a moving
   blocks monte-carlo simulation.  I will do this by taking blocks of
   contiguous years with a random start and reordering them
   randomly.  This conditions on the number of event occurrences,
   which I would rather do than have it be random, and partially
   preserves the temporal autocorrelation.
   

If anyone has any thoughts, or pointers, they'd be very welcome.

Cheers

Andrew
-- 
Andrew Robinson  
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au

__
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] Aggregate Values for All Levels of a Factor

2006-10-05 Thread Kaom Te
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hello,

I'm a novice user trying to figure out how to retain NA aggregate
values. For example, given a data frame with data for 3 of the 4
possible factor colors(orange is omitted from the data frame), I want
to calculate the average height by color, but I'd like to retain the
knowledge that orange is a possible factor, its just missing. Here is
the example code:

 data - data.frame(color = factor(c(blue,red,red,green,blue),
levels = c(blue,red,green,orange)),
height = c(2,8,4,4,5))
 aggregate(data$height, list(color = data$color), mean)
  color   x
1  blue 3.5
2   red 6.0
3 green 4.0


Instead I would like to get

   color   x
1   blue 3.5
2red 6.0
3  green 4.0
4 orange  NA

Is this possible. I've read as much documentation as I can find, but am
unable to find the solution. It seems like something people would need
to do. So I would assume it must be built in somewhere or do I need to
write my own version of aggregate?

Thanks in advance,
Kaom
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.5 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFFJYrLaaZgZdCbWv4RApNoAJ9jqKXne3IlQnd+PprS+7Kz1l4oRACfeu5I
Nv/xYWVsSGJD5+fdCP+02jk=
=b5TI
-END PGP SIGNATURE-

__
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] solaris 64 build?

2006-10-05 Thread roger koenker
We have a solaris/sparc machine that has been running an old version
of R-devel:  Version 2.2.0 Under development (unstable) (2005-06-04  
r34577)
which was built as m64 from sources.  Attempting to upgrade to 2.4.0  
the configure step
goes ok, but I'm getting early on from make:

 gcc -m64  -L/opt/sfw/lib/sparcv9  -L/usr/lib/sparcv9
 -L/usr/openwin/lib/sparcv9  -L/usr/local/lib -o R.bin Rmain.o
 CConverters.o CommandLineArgs.o  Rdynload.o Renviron.o RNG.o  apply.o
 arithmetic.o apse.o array.o attrib.o  base.o bind.o builtin.o
 character.o coerce.o colors.o complex.o connections.o context.o  cov.o
 cum.o  dcf.o datetime.o debug.o deparse.o deriv.o  dotcode.o dounzip.o
 dstruct.o duplicate.o  engine.o envir.o errors.o eval.o  format.o
 fourier.o  gevents.o gram.o gram-ex.o graphics.o  identical.o  
 internet.o
 iosupport.o  lapack.o list.o localecharset.o logic.o  main.o mapply.o
 match.o memory.o model.o  names.o  objects.o optim.o optimize.o
 options.o  par.o paste.o pcre.o platform.o  plot.o plot3d.o plotmath.o
 print.o printarray.o printvector.o printutils.o qsort.o  random.o
 regex.o registration.o relop.o rlocale.o  saveload.o scan.o seq.o
 serialize.o size.o sort.o source.o split.o  sprintf.o startup.o
 subassign.o subscript.o subset.o summary.o sysutils.o  unique.o util.o
 version.o vfonts.o xxxpr.o  mkdtemp.o ../unix/libunix.a
 ../appl/libappl.a ../nmath/libnmath.a -L../../lib -lRblas
 -L/usr/local/encap/gf7764-3.4.3+2/lib/gcc/sparc64-sun-solaris2.9/3.4.3
 -L/usr/ccs/bin/sparcv9 -L/usr/ccs/bin -L/usr/ccs/lib
 -L/usr/local/encap/gf7764-3.4.3+2/lib/sparcv9
 -L/usr/local/encap/gf7764-3.4.3+2/lib -lg2c -lm -lgcc_s
 ../extra/zlib/libz.a  ../extra/bzip2/libbz2.a ../extra/pcre/libpcre.a
 ../extra/intl/libintl.a  -lreadline -ltermcap -lnsl -lsocket -ldl -lm


 Undefined   first referenced
 symbol in file
 __builtin_isnan arithmetic.o
 ld: fatal: Symbol referencing errors. No output written to R.bin
 collect2: ld returned 1 exit status

I've tried to look at the difference in outcomes in the old R-devel
version --  if I touch arithmetic.c  there and then type make I get a  
something
almost the same as above except for the following  bits that are new  
to 2.4.0
(this diff is after replacing spaces with linebreaks obviously.)

ysidro.econ.uiuc.edu% diff t0 t1
54a55
  localecharset.o
81a83
  rlocale.o
101a104
  mkdtemp.o
104a108,109
  -L../../lib
  -lRblas


Has there been some change in the way that Rblas is used, or in
isnan?  It didn't seem so from a look at arithmetic.c, but this is well
beyond me.

I hope that someone sees something suspicious, or could point me
toward a better diagnostic.  Thanks,

Roger


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

__
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] xyplot

2006-10-05 Thread David Barron
If you include a conditioning factor (here it is the variable
groups) you will get more than one panel. To get loess curves for
the two groups, the easiest way is to use the type argument:

 xyplot(y ~ time, groups=group,data=df,type=c(p,smooth))


On 05/10/06, Osman Al-Radi [EMAIL PROTECTED] wrote:
 Hello,

 Thanks for your answer, however, the resulting graph has two panels with a
 loess curve per subject. I need a single panel with a loess curve per
 group..

 Osman

 On 10/5/06, Sundar Dorai-Raj [EMAIL PROTECTED] wrote:
 
 
 
  Osman Al-Radi said the following on 10/5/2006 3:43 PM:
   Hi,
  
   for the data below:
   time-c(rep(1:10,5))
   y-time+rnorm(50,5,2)
   subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10))
   group-c(rep('A',30),rep('B',20))
   df-data.frame(subject,group,time,y)
  
   I'd like to produce a plot with a single pannel with two loess curves
  one
   for each group. the code below does that but on two different pannels,
  I'd
   like to have both curves on the same pannel with different colors.
  
   xyplot(y~time|group, groups=subject,
   panel=function(...){
   panel.loess(...)
   }
   ,data=df)
  
  
 
  Simple:
 
  xyplot(y ~ time | group, data = df, groups = subject,
  panel = function(...) {
panel.superpose(...)
panel.superpose(panel.groups = panel.loess, ...)
  })
 
  HTH,
 
  --sundar
 



 --
 Osman O. Al-Radi, MD, MSc, FRCSC
 Fellow, Cardiovascular Surgery
 The Hospital for Sick Children
 University of Toronto, Canada

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



-- 
=
David Barron
Said Business School
University of Oxford
Park End Street
Oxford OX1 1HP

__
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] Block comments in R?

2006-10-05 Thread Richard A. O'Keefe
Please let's not waste any time trying to figure out how to add
block comments to R.  In any guise they are highly error prone.
Although its precursors (PL/I and Pascal) had block comments and did not
have end-of-line comments, the programming language Ada was explicitly
designed to have end-of-line comments and not block comments, as block
comments were thought too easy to write incorrectly and too easy to misread.

As just one example, the version of vim I have does syntax highlighting,
BUT Vim doesn't parse the whole file, so if you have code in the middle
of a block comment it may be highlighed incorrectly.

Commenting code out and providing documentation comments are easily
done with a good editor, although R documentation comments really belong
in files where help() can find them.

In the text editor I use, comment-region (using end-of-line comments)
and uncomment-region (ditto) are just two keystrokes each; it's
LESS effort than writing block comment delimiters would be.
Also, when writing a big comment, the newline key automatically
copies white space*(end of line comment mark white space*)?
to the beginning of a new line, so again it is less typing to make
a block comment using end-of-line comments than /* ... */ or anything
like that, and of course fill-paragraph and justify-paragraph (two
keystrokes each) can cope with the comment marks (a change from lines
with comment marks to ones without, or vice versa, constitutes a
paragraph break).

__
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 to get the function names

2006-10-05 Thread Gabor Grothendieck
Probably the best you can hope for is to cover
most cases.  This one uses match.call and handles
a number of cases and perhaps if you spend more time
on it might be able to add some cases where it fails
such as the second L below:

f - function(x) {
if (!is.list(x)) x - list(x)
if (is.null(names(x))) names(x) - 
names(x)[names(x) == ] - NA
mc - match.call()[-1][[1]]
if (length(mc)  1) mc - mc[-1]
ifelse(is.na(names(x)), as.character(mc), names(x))
}

f(c(a = mean))
f(list(a = mean, b = sd))
f(c(f = function(x)x*x))
f(list(f = function(x)x*x, function(x)1-x))
L - list(a = mean, b = sd)
f(L)
L - list(a = mean, function(x)x)
f(L)

f(mean)
f(list(a = mean, sd))
f(list(mean, sd))
f(function(x)x*x)
f(list(function(x)x*x, function(y)y-1))


On 10/5/06, Søren Højsgaard [EMAIL PROTECTED] wrote:
 I've defined the function

 getFunNames - function(FUN){
  if (!is.list(FUN))
fun.names - paste(deparse(substitute(FUN)), collapse =  )
  else
fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a)))
  fun.names
 }

 which gives what I want :
  getFunNames(mean)
 [1] mean
  getFunNames(ff)
 [1] ff
  getFunNames(c(mean,ff))
 [1] mean ff

 If I call this within a function, things go wrong:
 1] FUN
  foo(ff)
 [1] FUN
  foo(c(mean,ff))
 Error in substitute(FUN)[-1] : object is not subsettable

 Obviously there are some things (quite a few things) which I have not 
 understood. Can anyone help?
 Thanks
 Søren

 __
 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] Aggregate Values for All Levels of a Factor

2006-10-05 Thread Marc Schwartz
On Thu, 2006-10-05 at 15:44 -0700, Kaom Te wrote:
 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1
 
 Hello,
 
 I'm a novice user trying to figure out how to retain NA aggregate
 values. For example, given a data frame with data for 3 of the 4
 possible factor colors(orange is omitted from the data frame), I want
 to calculate the average height by color, but I'd like to retain the
 knowledge that orange is a possible factor, its just missing. Here is
 the example code:
 
  data - data.frame(color = factor(c(blue,red,red,green,blue),
 levels = c(blue,red,green,orange)),
   height = c(2,8,4,4,5))
  aggregate(data$height, list(color = data$color), mean)
   color   x
 1  blue 3.5
 2   red 6.0
 3 green 4.0
 
 
 Instead I would like to get
 
color   x
 1   blue 3.5
 2red 6.0
 3  green 4.0
 4 orange  NA
 
 Is this possible. I've read as much documentation as I can find, but am
 unable to find the solution. It seems like something people would need
 to do. So I would assume it must be built in somewhere or do I need to
 write my own version of aggregate?
 
 Thanks in advance,
 Kaom

If you review the Details section of ?aggregate, you will note:

  Empty subsets are removed, ...

Thus, one approach is:

tmp - tapply(data$height, data$color, mean, na.rm = TRUE)

 tmp
  bluered  green orange
   3.56.04.0 NA

DF - data.frame(color = names(tmp), mean.height = tmp, 
 row.names = seq(along = tmp))

 DF
   color mean.height
1   blue 3.5
2red 6.0
3  green 4.0
4 orange  NA


HTH,

Marc Schwartz

__
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 to get the function names

2006-10-05 Thread Gabor Grothendieck
I should have mentioned is that the way it works is that
it uses the name of the list component, if any, otherwise
it uses the name of the function if its given as a number
and otherwise it uses the function itself or possibly the
name of the list.

On 10/5/06, Gabor Grothendieck [EMAIL PROTECTED] wrote:
 Probably the best you can hope for is to cover
 most cases.  This one uses match.call and handles
 a number of cases and perhaps if you spend more time
 on it might be able to add some cases where it fails
 such as the second L below:

 f - function(x) {
if (!is.list(x)) x - list(x)
if (is.null(names(x))) names(x) - 
names(x)[names(x) == ] - NA
mc - match.call()[-1][[1]]
if (length(mc)  1) mc - mc[-1]
ifelse(is.na(names(x)), as.character(mc), names(x))
 }

 f(c(a = mean))
 f(list(a = mean, b = sd))
 f(c(f = function(x)x*x))
 f(list(f = function(x)x*x, function(x)1-x))
 L - list(a = mean, b = sd)
 f(L)
 L - list(a = mean, function(x)x)
 f(L)

 f(mean)
 f(list(a = mean, sd))
 f(list(mean, sd))
 f(function(x)x*x)
 f(list(function(x)x*x, function(y)y-1))


 On 10/5/06, Søren Højsgaard [EMAIL PROTECTED] wrote:
  I've defined the function
 
  getFunNames - function(FUN){
   if (!is.list(FUN))
 fun.names - paste(deparse(substitute(FUN)), collapse =  )
   else
 fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a)))
   fun.names
  }
 
  which gives what I want :
   getFunNames(mean)
  [1] mean
   getFunNames(ff)
  [1] ff
   getFunNames(c(mean,ff))
  [1] mean ff
 
  If I call this within a function, things go wrong:
  1] FUN
   foo(ff)
  [1] FUN
   foo(c(mean,ff))
  Error in substitute(FUN)[-1] : object is not subsettable
 
  Obviously there are some things (quite a few things) which I have not 
  understood. Can anyone help?
  Thanks
  Søren
 
  __
  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.


  1   2   >