Re: [R] Generation of Normal Random Numbers

2011-01-11 Thread Michael Bedward
I realized after posting that you can also just do this...

x <- matrix(rnorm(25 * length(m), m, s), ncol=4, byrow=TRUE)

Michael


On 11 January 2011 19:13, Amelia Vettori  wrote:

> Dear Mr Bedward and Mr Newmiller,
>
> Thanks a lot for the guidance. Really appreciated.
>
> Sagga K
>
> --- On *Tue, 11/1/11, Michael Bedward * wrote:
>
>
> From: Michael Bedward 
> Subject: Re: [R] Generation of Normal Random Numbers
> To: "saggak" 
> Cc: r-help@r-project.org
> Received: Tuesday, 11 January, 2011, 7:24 AM
>
> m <- c(1004.1, 1028.3, 1044.3, 861.4)
> s <- c(194.5899, 158.7052, 123.3000, 285.8695)
> x <- mapply(function(mi, si) rnorm(25, mi, si), m, s)
>
> Hope this helps,
>
> Michael
>
>
> On 11 January 2011 17:44, saggak  wrote:
> > Dear R helpers
> >
> > I have a data frame as given below
> >
> > df = data.frame(A = c(776,827,836,995,855,1026,1203,1363,965,1195),
> >B =
> c(806,953,1049,1056,1243,764,1148,1162,948,1154),
> >C =
> c(959,1155,1193,1163,863,1070,1087,877,1132,944),
> >D = c(906,760,978,1170,1009,883,1007,960,828,113))
> >
> > # Actually the real data has number of vectors and not only A, B, C and
> D.
> >
> > m = as.numeric(lapply(df, mean))
> > s = as.numeric(lapply(df, sd))
> >
> > gives
> >
> >> m
> > [1] 1004.1 1028.3 1044.3  861.4
> >> s
> > [1] 194.5899 158.7052 123.3000 285.8695
> >
> > I need to generate 25 (normal) random numbers for each of these mean and
> corresponding standard deviation combination (m[i], s[i]). i.e. I need to
> have a table (dim 25 X  4) giving me random numbers.
> >
> >
> > rnorm(25, m[1], s[1])   rnorm(25, m[2], s[2])   rnorm(25, m[3], s[3])
> rnorm(25, m[4], s[4])
> >   .
> ... .
> >   .
> ... .
> >   .
> ... .
> >   .
> ... .
> >
> >
> >   .
> ... .
> >
> >
> > Kindly guide
> >
> >
> > Thanking in advance
> >
> > Sagga K
> >
> >
> >
> >
> >
> >
> >
> >[[alternative HTML version deleted]]
> >
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
>

[[alternative HTML version deleted]]

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


Re: [R] Help with Data Transformation

2011-01-11 Thread ONKELINX, Thierry
Dear Guy,

Have a look at cast() from the reshape package. You'll need something
like

cast(fldsampleid ~ Analysis, value = "Result", data = your.data.frame)

Best regards,

Thierry



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey
  

> -Oorspronkelijk bericht-
> Van: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] Namens Guy Jett
> Verzonden: maandag 10 januari 2011 22:00
> Aan: r-help@r-project.org
> Onderwerp: [R] Help with Data Transformation
> 
> Greetings,
> I am new to R and am having trouble with parsing a file with 
> the following characteristics:
> 
> * Individual results for a single sample are written 
> to multiple lines.
> 
> * First 16 columns are constant from sample to sample.
> 
> * Remaining 10 need to be matched up (cross-tabbed?)
> 
> o   (the exact contents for the remaining 10 vary from sample 
> to sample, as indicated in the extract below)
> 
> * Ultimate goal is to run various comparisons between 
> the variable columns, compare samples from separate 
> populations, and graph samples from the separate populations.
> 
> * (An extract is provided below)
> 
> The data is initially extracted from an SQL database into 
> Excel, then saved as a tab-delimited text file for use in R.
> I have been successful in using subset() to extract specific 
> sample types, but have not yet been able to transform the 
> data so that all the data needed is on a single line.  I have 
> looked at several R manuals, read through 'R in a Nutshell', 
> prowled the help resources (R Site Search and the Google 
> link), tried stack(), subset(), reshape(), and several other 
> functions, to no avail.
> 
> Thank you very much for your help.  This seems like a 
> wonderful community,
> Guy Jett, R.G.
> Project Geologist
> gj...@itsi.com
> 
> Example Data Input (subset):
> fldsampidCLP_ID  sacode  matrix   
> etc...  prccodeLab EXMCODE
>AnalysisPARLABELPARVQ Result
> 2268   LHR020GW-01E2   N  
>WGINOBRLS  NONE
> E300   CL   = 23590.9
> 2269   LHR020GW-01E2   N  
>WGINOBRLS  NONE
> E300   PO4ND  50
> 2270   LHR020GW-01E2   N  
>WGINOBRLS  NONE
> E300   SO4= 22460
> 2272   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1631HG  = 0.00171
> 2273   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638AG  = 2.57
> 2274   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638AL   = 122
> 2275   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638AS   = 317
> 2276   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638B = 9970
> 2289   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638V = 131
> 2290   LHR020GW-01E2   N  
>WGMET   BRLS  FLDFLT  
> E1638Zn   = 1.76
> 2291   LHR020GW-01E2   N  
>WGMET   BRLS  METHOD   
>E1638PB   ND0.008
> 2292   LHR020GW-01E2   N  
>WGMI  BRLS  NONE   
>  A2320ALK= 807000
> 2293   LHR020GW-01E2   N   

[R] estimates in subpopulations (package survey)

2011-01-11 Thread andrija djurovic
Dear all,

I have a example from library survey about estimating in subpopulations and
I need to reproduce it step by step.  This is the example (from paper:
Estimates in subpopulation of Thomas Lumley and just instead of svymean I
used svytotal) that I did in R:

library(survey)

data(fpc)
dfpc<- svydesign(id=~psuid, strata=~stratid, weight=~weight,
data=fpc, nest=TRUE)

dsub<- subset(dfpc, x>4)

svytotal(~x, design=dsub)
  total SE
x 123.9 33.004
svyby(~x, ~I(x>4), design = dfpc, svytotal)
I(x > 4) x se.x
FALSEFALSE  23.2 17.01764
TRUE  TRUE 123.9 33.00409

and what I can't figure out is which formula or approximation is used for
standard error calculation. I tried with some famouse literature written by
Cohran, Sarndal, Lohr  but I still can't figure out this.
 I would really appreciate if anyone can explain me in few steps formula
which is used to obtain value of 33.00409 for standard error.

Best regards,
Andrija

[[alternative HTML version deleted]]

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


[R] [R-pkgs] Package Rd2roxygen v0.1-5

2011-01-11 Thread Yihui Xie
Hi,

A new version (0.1-5) of Rd2roxygen is on CRAN now.

Rd2roxygen is a package to help developers transit from the raw Rd
files to roxygen-style documentations, and it also has utilities to
help build the package (with options to reformat the usage and
examples sections).

Updates:

 o a new way to reformat the usage and examples sections generated by
roxygen; it is more reliable than earlier versions, e.g. the
\subsection{}{} macro will no longer be destroyed when reformatting
the Rd and building the package.

See the package vignette for more details:

library(Rd2roxygen)
vignette('Rd2roxygen')
## or http://cran.r-project.org/web/packages/Rd2roxygen/vignettes/Rd2roxygen.pdf

Regards,
Yihui
--
Yihui Xie 
Phone: 515-294-2465 Web: http://yihui.name
Department of Statistics, Iowa State University
2215 Snedecor Hall, Ames, IA

___
R-packages mailing list
r-packa...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-packages

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


Re: [R] Assumptions for ANOVA: the right way to check the normality

2011-01-11 Thread DrorD


Frodo Jedi wrote:
> 
> 
>>What is the question you are really trying to find the answer for? 
Knowing that 
>>may help us give more meaningful answers.
> 
> Concerning your question I thought to have been clear. 
> 
> I want to understand which analysis I have to use in order to understand
> if the differences I am having are statistically significant or not. 
> 
> 

Dear Frodo, 

I would like to suggest that the question required is a concrete
specification of an hypothesis. 

Something like: 
"I hypothesize that the responses to condition-A would be different in
magnitude from the responses to condition-AH, across all stimuli."

Perhaps after having a detailed formulation of your hypothesis the required
analysis will be clearer for yourself, or at least make it easier for
experts to guide you. 

Best, 
dror
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Assumptions-for-ANOVA-the-right-way-to-check-the-normality-tp3176073p3208596.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] how to coerce part of each column of a matrix to a vector and merge them

2011-01-11 Thread Michael Bedward
Hello,

The answer to this one drops out of the answer to your previous question...

m <- matrix(1:16, nrow=4)
end <- c(2,3,1,3)
ii <- cbind(sequence(end), rep(1:length(end), end))
x <- m[ ii ]

Hope this helps,
Michael


2011/1/11 zhaoxing731 :
> Hello
>
> Suppose I have a matrix mat=(1:16,2)
> [,1] [,2] [,3] [,4]
> [1,]159   13
> [2,]26   10   14
> [3,]37   11   15
> [4,]48   12   16
>
> there is a vector end=c(2,3,1,3)
>
> #coerce the 1st 2 numbers of the 1st column to a vector  [1]  1  2
> #coerce the 1st 3 numbers of the 2nd column and append it to the previous 
> vector  [1]  1  2  5  6  7
> #coerce the 1st number of the 3rd column and append it to the previous vector 
>  [1]  1  2  5  6  7  9
> #coerce the 1st 3 numbers of the 4th column and append it to the previous 
> vector  [1]  1  2  5  6  7  9 13 14 15
> #they are specified by vector end
> a for loop
> mat<-matrix(1:16,4)
> end<-c(2,3,1,3)
> vec<-numeric()
> for (i in 1:4)
>{vec<-c(vec,mat[1:end[i],i])
>}
> #the result
>  > vec
> [1]  1  2  5  6  7  9 13 14 15
>
>  but when I want to do it to a large dataset, it's inefficiency becomes a 
> problem, so need vectorization
>
>  Thank you in advance
>  Yours sincerely
>
> ZhaoXing
> Department of Health Statistics
> West China School of Public Health
> Sichuan University
> No.17 Section 3, South Renmin Road
> Chengdu, Sichuan 610041
> P.R.China
>
> __
> 赶快注册雅虎超大容量免费邮箱?
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] Generation of Normal Random Numbers

2011-01-11 Thread Amelia Vettori
Dear Mr Bedward and Mr Newmiller,

Thanks a lot for the guidance. Really appreciated.

Sagga K

--- On Tue, 11/1/11, Michael Bedward  wrote:

From: Michael Bedward 
Subject: Re: [R] Generation of Normal Random Numbers
To: "saggak" 
Cc: r-help@r-project.org
Received: Tuesday, 11 January, 2011, 7:24 AM

m <- c(1004.1, 1028.3, 1044.3, 861.4)
s <- c(194.5899, 158.7052, 123.3000, 285.8695)
x <- mapply(function(mi, si) rnorm(25, mi, si), m, s)

Hope this helps,

Michael


On 11 January 2011 17:44, saggak  wrote:
> Dear R helpers
>
> I have a data frame as given below
>
> df = data.frame(A = c(776,827,836,995,855,1026,1203,1363,965,1195),
>    B = c(806,953,1049,1056,1243,764,1148,1162,948,1154),
>    C = c(959,1155,1193,1163,863,1070,1087,877,1132,944),
>    D = c(906,760,978,1170,1009,883,1007,960,828,113))
>
> # Actually the real data has number of vectors and not only A, B, C and
 D.
>
> m = as.numeric(lapply(df, mean))
> s = as.numeric(lapply(df, sd))
>
> gives
>
>> m
> [1] 1004.1 1028.3 1044.3  861.4
>> s
> [1] 194.5899 158.7052 123.3000 285.8695
>
> I need to generate 25 (normal) random numbers for each of these mean and 
> corresponding standard deviation combination (m[i], s[i]). i.e. I need to 
> have a table (dim 25 X  4) giving me random numbers.
>
>
> rnorm(25, m[1], s[1])   rnorm(25, m[2], s[2])   rnorm(25, m[3], s[3])  
> rnorm(25, m[4], s[4])
>       .   
> ...                 .
>   
    .   
...                 .
>       .   
> ...                 .
>       .   
> ...                 .
>
>
>      
 .   ...                 
.
>
>
> Kindly guide
>
>
> Thanking in advance
>
> Sagga K
>
>
>
>
>
>
>
>        [[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible
 code.
>
>

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



  
[[alternative HTML version deleted]]

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


[R] using Smolyak to genarate intervals

2011-01-11 Thread Lara Aleluia
 Hello

I am running a program from value intervals for a set of variables.
I thought I could use Smolyak to get my intervals adding points as I need
them , so instead of running the program for a whole interval I would run
only for the most important points.
I am new to Smolyak algorithm though, and I can not understand how can i use
smolyak.quad(d,k), to solve my problem.

Anybody can give me some lights?

[[alternative HTML version deleted]]

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


Re: [R] Problems producing quantreg-Tables

2011-01-11 Thread Jan Henckens


I think I found out about the existing problem with generation of tables 
done by latex():


Checking the source code of 'table.R' I noticed any optional arguments 
specified seem to be discarded when calling the function latex.table() 
after preprocessing the summary.rqs-output


So unlike stated in the reference manual (p. 25)

> latex.summary.rqs Make a latex table from a table of rq results
[...]
> ...   optional arguments for latex.table

it is not possible to pass any preferences regarding the layout to 
latex.table().



The solution is extremely simple:

by changing line 36 in the code 'table.R':

latex.table(table, caption = caption, rowlabel = rowlabel, file = file)

into:

latex.table(table, caption = caption, rowlabel = rowlabel, file = file, ...)

the funcionality described in the reference manual is implemented.


Maybe this helps someone who is confronted with the same problem...


Jan


--
jan.henckens | jöllenbecker str. 58 | 33613 bielefeld
tel 0521-5251970

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


Re: [R] problem with packages

2011-01-11 Thread Duncan Murdoch

On 11-01-10 5:43 PM, R Heberto Ghezzo, Dr wrote:

Hello, I am on a laptop with Win7, running R-2.12.1
if I click on Packages/InstallPackages I get :


utils:::menuInstallPkgs()

Warning: unable to access index for repository 
http://cran.skazkaforyou.com/bin/windows/contrib/2.12
Warning: unable to access index for repository 
http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.12
Error in install.packages(NULL, .libPaths()[1L], dependencies = NA, type = 
type) :
   no packages were specified



if I try update packages :

update.packages(ask='graphics')

Warning: unable to access index for repository 
http://cran.skazkaforyou.com/bin/windows/contrib/2.12
Warning: unable to access index for repository 
http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.12



My profile.site file is :
# Things you might want to change
# options(papersize="a4")
# options(editor="notepad")
# options(pager="internal")
# set the default help type
# options(help_type="text")
# options(help_type="html")
# set a CRAN mirror
  local({r<- getOption("repos")
r["CRAN"]<- "http://cran.skazkaforyou.com";
options(repos=r)})
# set library paths
  aaa<- strsplit(R.home(),"R-")[[1]][1]
.libPaths(c(paste(aaa,"cran",sep=""),paste(aaa,"Others",sep="")))
#  set for 'jags' , 'WinBUGS' , 'Sweave' , 'EBImage' , 'JGR'
  aa<- strsplit(R.home(),"R")[[1]][1]
local({
  Sys.setenv("JAGS_HOME"=paste(aa,"R/JAGS-2.2.0",sep=""))
  Sys.setenv("GTK_BASEPATH"=paste(aa,"R/GTK",sep=""))
  Sys.setenv("R_HOME"=paste(aa,"R/R-2.12.1",sep=""))
  Sys.setenv("R_LIBS"=paste(aa,"R/Cran",sep=""))
  Sys.setenv("DYLD_LIBRARY_PATH"=paste(aa,"R/R-2.12.1/bin/i386",sep=""))
  PathNew<- 
paste(aa,"R/GTK/bin;",aa,"R/ImageMagick;",Sys.getenv("PATH"),";",aa,"PortableUSB/PortableApps/MikeTex/miktex/bin",sep="")
  Sys.setenv("DirWinBugs"=paste(aa,"PortableUSB/PortableApps/WinBUGS14",sep=""))
  Sys.setenv("PATH"=PathNew)
#
  })


Any ideas what can be wrong?


Your system can't get through to those two web sites, probably because 
you have a firewall blocking something, or need to use a proxy.  You 
might try running


setInternet2()

first, and R will use the Internet Explorer proxy settings, assuming IE 
can get through to those sites.  If it can't, then you've got to fix 
your Internet connection.


Duncan Murdoch


Thanks

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


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


Re: [R] problem with packages

2011-01-11 Thread Uwe Ligges



On 10.01.2011 23:43, R Heberto Ghezzo, Dr wrote:

Hello, I am on a laptop with Win7, running R-2.12.1
if I click on Packages/InstallPackages I get :


utils:::menuInstallPkgs()

Warning: unable to access index for repository 
http://cran.skazkaforyou.com/bin/windows/contrib/2.12
Warning: unable to access index for repository 
http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.12
Error in install.packages(NULL, .libPaths()[1L], dependencies = NA, type = 
type) :
   no packages were specified



if I try update packages :

update.packages(ask='graphics')

Warning: unable to access index for repository 
http://cran.skazkaforyou.com/bin/windows/contrib/2.12
Warning: unable to access index for repository 
http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.12



No internet connection? Behind a proxy?

Uwe Ligges



My profile.site file is :
# Things you might want to change
# options(papersize="a4")
# options(editor="notepad")
# options(pager="internal")
# set the default help type
# options(help_type="text")
# options(help_type="html")
# set a CRAN mirror
  local({r<- getOption("repos")
r["CRAN"]<- "http://cran.skazkaforyou.com";
options(repos=r)})
# set library paths
  aaa<- strsplit(R.home(),"R-")[[1]][1]
.libPaths(c(paste(aaa,"cran",sep=""),paste(aaa,"Others",sep="")))
#  set for 'jags' , 'WinBUGS' , 'Sweave' , 'EBImage' , 'JGR'
  aa<- strsplit(R.home(),"R")[[1]][1]
local({
  Sys.setenv("JAGS_HOME"=paste(aa,"R/JAGS-2.2.0",sep=""))
  Sys.setenv("GTK_BASEPATH"=paste(aa,"R/GTK",sep=""))
  Sys.setenv("R_HOME"=paste(aa,"R/R-2.12.1",sep=""))
  Sys.setenv("R_LIBS"=paste(aa,"R/Cran",sep=""))
  Sys.setenv("DYLD_LIBRARY_PATH"=paste(aa,"R/R-2.12.1/bin/i386",sep=""))
  PathNew<- 
paste(aa,"R/GTK/bin;",aa,"R/ImageMagick;",Sys.getenv("PATH"),";",aa,"PortableUSB/PortableApps/MikeTex/miktex/bin",sep="")
  Sys.setenv("DirWinBugs"=paste(aa,"PortableUSB/PortableApps/WinBUGS14",sep=""))
  Sys.setenv("PATH"=PathNew)
#
  })


Any ideas what can be wrong?
Thanks

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


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


[R] [R-pkgs] Package animation update 2.0-1

2011-01-11 Thread Yihui Xie
Hi,

The package animation 2.0-1 is on CRAN now
(http://cran.r-project.org/package=animation).


                 CHANGES IN animation VERSION 2.0-1


NEW FEATURES

   o demo('Xmas_card') contributed by Yuan Huang

   o demo('flowers') to show how to download images from the Internet
 and create an animation

   o a new function pdftk() as a wrapper to call the Pdftk toolkit
 (mainly for compressing PDF graphics output)

   o saveLatex(), saveSWF() and saveMovie() can compress the PDF
 graphics using Pdftk (if available) now; if ani.options('ani.type')
 is 'pdf' and the 'pdftk' option is set, these functions will try the
 compression first

   o new functions ani.record() and ani.replay() to record and replay
 the animations; they can be used to capture the changes in the
 graphics made by low-level plotting commands (see ?ani.record for
 examples)


SIGNIFICANT CHANGES

   o the function tidy.source() was completely removed from this
 package; users should go to the formatR package (tidy.source() is
 there now)


MINOR CHANGES

   o the argument 'expr' in saveHTML(), saveLatex(), saveMovie() and
 saveSWF() will be evaluated by eval(), so we may pass a real R
 expression (see ?expression) to 'expr', e.g. saveHTML(expression(for
 (i in 1:10) plot(runif(30; but note this argument does not *have
 to* be a real R expression -- you can still use an arbitrary R code
 chunk, e.g. saveHTML(for (i in 1:10) plot(runif(30))) (thanks,
 Aquery)


Regards,
Yihui
--
Yihui Xie 
Phone: 515-294-2465 Web: http://yihui.name
Department of Statistics, Iowa State University
2215 Snedecor Hall, Ames, IA

___
R-packages mailing list
r-packa...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-packages

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


Re: [R] ThinkCell type waterfall charts in R?

2011-01-11 Thread Jim Lemon

On 01/11/2011 09:20 AM, ang wrote:


I was actually looking to create a similar type of graph in R.
But I guess I am going to try to approach it by using a stacked column chart
and hiding the 'net' series, while only showing the increases/decreases in
value.

I'll post an update later on what I come up with.



Hi ang,
You might want to look at the waterfall plot function in the plotrix 
package.


Jim

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


[R] sorting question

2011-01-11 Thread Albert-Jan Roskam
Hi,

I have a data frame with variables a, b, c (character vars) and t (time var, 
could be represented as POSIXct or character, depending on which is most 
useful. 
The format is "-mm-dd hh:mm:ss CET"). Now, I want to sort the data frame in 
ascending order by a, b, c and then in descending order by t.

Here's what I've got, but I'm not sure how to put the 'descending' bit in the 
code. The trick with the minus sign doesn't work as the var isn't numeric. 


idvars <-- c("a", "b", "c")
myData["t"] <- as.POSIXct(myData[,"t"])
result <- myData[ do.call(order, myData[c(idvars, "t")]), ]

This should be simple, I guess, but I've been staring at this a bit too long 
now. Suggestions anyone? TIA!
Cheers!!
Albert-Jan 


~~
All right, but apart from the sanitation, the medicine, education, wine, public 
order, irrigation, roads, a fresh water system, and public health, what have 
the 
Romans ever done for us?
~~ 


  
[[alternative HTML version deleted]]

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


Re: [R] Rectangle height in lattice xyplot key

2011-01-11 Thread Duncan Mackay

Hi Deepayan

Thank you very much for the update to lattice and the solution.

For the current problem the texts are all single lines and the solution 
works well.
The example was to accentuate the problem and I thought that it would be 
useful in the future.
I was able to make all the heights look equal in the example by creating a 
vector and reducing the height of the 2-lined rows.


xyplot(1~1,
   key = list(corner = c(0.8,0.8),
  cex   = 0.6,
  title = "Functional Groups",
  cex.title = 0.7,
  columns = 1,
  padding.text = 4,
  text  = list(label =
c("Tree","Shrub\n(low)","Herb","Grass","Shrub\n(small)")),
  points = list(pch = c(1,3,4,20,16),
col = c(2,3,4,5,6)),
  lines = list(col = c(2,3,4,5,6),
   size = 2),
  rectangles = list(col = c(2,3,4,5,6),
size = 1,
height = c(0.7,0.45,0.7,0.7,0.45),
border = FALSE)) )


Regards

Duncan

At 17:14 11/01/2011, you wrote:

On Sun, Jan 9, 2011 at 9:31 AM, Duncan Mackay  wrote:
> Dear All
>
> I have a problem with the height of the boxes in the key in the following.
> (The text is over 2 lines to  accentuate the problem of no space 
between the

> rectangles.)
> Is there an easy way to put a space between the rectangles; size controls
> the width but there appears to be nothing for the height?

There is now (in the last update of lattice, released last week). So
the following should work:

xyplot(1~1,
   key = list(corner = c(0.8,0.8),
  cex   = 0.6,
  title = "Functional Groups",
  cex.title = 0.7,
  columns = 1,
  padding.text = 4,
  text  = list(label =
c("Tree","Shrub\n(low)","Herb","Grass","Shrub\n(small)")),
  points = list(pch = c(1,3,4,20,16),
col = c(2,3,4,5,6)),
  lines = list(col = c(2,3,4,5,6),
   size = 2),
  rectangles = list(col = c(2,3,4,5,6),
size = 1,
height = 0.7, ## newly added
border = FALSE)) )

Of course 'height' is just a multiplier, so the heights are still not
all the same because of differing row heights, but I assume that's not
going to be a problem in your real example.

-Deepayan

>
> Â  Â  Â  Â xyplot(1~1,
> Â  Â  Â  Â  Â  key = list(corner = c(0.8,0.8),
>                      cex   = 0.6,
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â title = "Functional Groups",
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â cex.title = 0.7,
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â columns = 1,
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â padding.text = 4,
>                      text  = list(label =
> c("Tree","Shrub\n(low)","Herb","Grass","Shrub\n(small)")),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â points = list(pch = c(1,3,4,20,16),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â col = c(2,3,4,5,6)),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â lines = list(col = c(2,3,4,5,6),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â size = 2),),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â rectangles = list(col = c(2,3,4,5,6),
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â size = 1,
> Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â border = 
FALSE)) )

>
> I could get things to work when the legend was on top using
> Â http://finzi.psych.upenn.edu/R/Rhelp02/archive/46654.html on something
> similar
> but could not get it to work in the last panel of a multipanel plot.
>
> Regards
>
> Duncan Mackay
> Department of Agronomy and Soil Science
> University of New England
> Armidale NSW 2351
> Email: home mac...@northnet.com.au
>
> R version 2.12.1 (2010-12-16)
> Platform: i386-pc-mingw32/i386 (32-bit)
>
> locale:
> [1] LC_COLLATE=English_Australia.1252 Â LC_CTYPE=English_Australia.1252
> LC_MONETARY=English_Australia.1252
> [4] LC_NUMERIC=C 
                      LC_TIME=English_Australia.1252

>
> attached base packages:
> [1] splines   datasets  utils     stats     graphics  grDevices grid
>  methods   base
>
> other attached packages:
> Â [1] locfit_1.5-6 Â  Â  Â  Â akima_0.5-4 Â  Â  Â  Â  lme4_0.999375-37
> Matrix_0.999375-46 Â mgcv_1.7-2
> Â [6] geepack_1.0-17 Â  Â  Â doBy_4.1.2 Â  Â  Â  Â  Â contrast_0.13 
Design_2.3-0

> Â  Â Hmisc_3.8-3
> [11] survival_2.36-2 Â  Â  gee_4.13-16 Â  Â  Â  Â  R2HTML_2.2 
         zoo_1.6-4

> Â  Â  Â  Â  gsubfn_0.5-5
> [16] proto_0.3-8 Â  Â  Â  Â  RODBC_1.3-2 Â  Â  Â  Â  gtools_2.6.2 
xtable_1.5-6

> Â R.oo_1.7.4
> [21] R.methodsS3_1.2.1 Â  foreign_0.8-41 Â  Â  Â chron_2.3-39 MASS_7.3-9
> Â latticeExtra_0.6-14
> [26] RColorBrewer_1.0-2 Â lattice_0.19-13
>
> loaded via a namespace (

Re: [R] sorting question

2011-01-11 Thread Albert-Jan Roskam
Got it, after rtmíng I realized I had to use xtfrm: 

ord <- order(c(idvars, -xtfrm("t")))
myData[ord, ]
 
It's downright ugly, as it confronts the user with some implementation detail 
(cf SAS or SPSS), but well, it works.
Cheers!!
Albert-Jan 


~~
All right, but apart from the sanitation, the medicine, education, wine, public 
order, irrigation, roads, a fresh water system, and public health, what have 
the 
Romans ever done for us?
~~ 






To: R Mailing List 
Sent: Tue, January 11, 2011 12:11:01 PM
Subject: [R] sorting question

Hi,

I have a data frame with variables a, b, c (character vars) and t (time var, 
could be represented as POSIXct or character, depending on which is most 
useful. 

The format is "-mm-dd hh:mm:ss CET"). Now, I want to sort the data frame in 
ascending order by a, b, c and then in descending order by t.

Here's what I've got, but I'm not sure how to put the 'descending' bit in the 
code. The trick with the minus sign doesn't work as the var isn't numeric.


idvars <-- c("a", "b", "c")
myData["t"] <- as.POSIXct(myData[,"t"])
result <- myData[ do.call(order, myData[c(idvars, "t")]), ]

This should be simple, I guess, but I've been staring at this a bit too long 
[[elided Yahoo spam]]
Cheers!!
Albert-Jan 


~~
All right, but apart from the sanitation, the medicine, education, wine, public 
order, irrigation, roads, a fresh water system, and public health, what have 
the 

Romans ever done for us?
~~ 


      
    [[alternative HTML version deleted]]


  
[[alternative HTML version deleted]]

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


[R] Add line to plot from stl decomposed time series?

2011-01-11 Thread Rainer M Krug
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Hi

I would like to add a line to the plot of the data panel of a stl
decomposed time series.

Example:

plot(stl(nottem, "per"))

plots a 4(or is t 8?) panel graph. I thought that I might be able to use
mfg to plot in the top (left?) panel ("data"), but it is not working.

Any help appreciated,

Rainer






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

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Tel:+33 - (0)9 53 10 27 44
Cell:   +27 - (0)8 39 47 90 42
Fax (SA):   +27 - (0)8 65 16 27 82
Fax (D) :   +49 - (0)3 21 21 25 22 44
Fax (FR):   +33 - (0)9 58 10 27 44
email:  rai...@krugs.de

Skype:  RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk0sRFQACgkQoYgNqgF2egpFHgCggCR2bxRmc+2rt1DILLoW3a04
tkcAni8QAfaCK6ulKcvSKggRZKcoe00Q
=18GE
-END PGP SIGNATURE-

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

2011-01-11 Thread Ed Keith
This is not exactly an R specific question, but I think the people on this list 
can probably help.

I'm working on a simulation. In the model I have the first three moments of the 
distributions of the variables. I know how to generate a random number from a 
distribution given the first two moments assuming the third moment is 0. But I 
do not know how to generate a number drawn from a distribution with a nonzero 
third monument.

If someone could point me to a good reference I would appreciate it.

Thank you in advance,

   -EdK

Ed Keith
e_...@yahoo.com

Blog: edkeith.blogspot.com

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


Re: [R] how to calculate this natural logarithm

2011-01-11 Thread Mike Marchywka



> Date: Sun, 9 Jan 2011 12:57:49 +
> From: c...@witthoft.com
> To: r-help@r-project.org
> Subject: Re: [R] how to calculate this natural logarithm
>
>
>
> You could re-do part of your code to run with mpfr-class variables, and use
> this function:
> # mpfr choose(n,k)
> rmpfac<-function(n,k,prec=50)
> factorial(mpfr(n,prec))/factorial(mpfr(k,prec))/factorial(mpfr(n-k,prec))
> Converting into and out of mpfr may not be worth it, but you will get all
> the precision you want without any nasty "Inf" results

Is the OP still here? I was curious what end result you were
after and if you tried any of these approaches. Apparently the quantity
here is something like 5000 which could probably itself
be an intermediate. If you took the approach of using algebra, as suggested
by earlier posters, and your
end result was supposed to be between 0 and 1, there would probably be
more simplifications to make. IIRC, you had terms like 
n!(2n!)
---
k!(n-k)!(600-k)!(2n-600+k)!
and you could cancel a lot ahead of time. You'd think
somewhere, like in statistical mechanics, this would be a well known
issue. 





> :-)
> Carl

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


Re: [R] From vector to a step function

2011-01-11 Thread Alaios
Dear Jim Holtman,
I would like to thank you for your help.
Just to update also the community that worked fine :)

Best Regards
Alex

--- On Mon, 1/10/11, jim holtman  wrote:

> From: jim holtman 
> Subject: Re: [R] From vector to a step function
> To: "Alaios" 
> Cc: r-help@r-project.org
> Date: Monday, January 10, 2011, 5:44 PM
> try this:
> 
> > f.main <- function(vec){
> +     breaks <- seq(-3, by = 1,
> length = length(vec) + 1L)
> +     function(x){
> +         indx <-
> findInterval(x, breaks)
> +         vec[indx]
> +     }
> + }
> > f1 <- f.main(c(3,4,5,1,2,3,4,5))
> > f2 <- f.main(c(5,6,2,4,7,3,2,5))
> > f3 <- f.main(c(1,2,4,7,3,1,3,5))
> >
> > f1(c(-2.5,-2,-.5,0,.5,1.5,2.5))
> [1] 3 4 5 1 1 2 3
> > f2(c(-2.5,-2,-.5,0,.5,1.5,2.5))
> [1] 5 6 2 4 4 7 3
> > f3(c(-2.5,-2,-.5,0,.5,1.5,2.5))
> [1] 1 2 4 7 7 3 1
> >
> 
> 
> On Mon, Jan 10, 2011 at 11:42 AM, Alaios 
> wrote:
> > Greetings R members.
> >
> > I have a few vectors that denote the 'steps' of
> different step functions
> > v1=c(3,4,5,1,2,3,4,5)
> > v2=c(5,6,2,4,7,3,2,5)
> > v3=c(1,2,4,7,3,1,3,5)
> >
> > Here v1,v2,v3 are considered as the steps for the
> f1,f2,f3 step functions.
> >
> > For example f1 looks like that (step size is always
> same and fixed)
> >
> > f1= 3 (x>=-3,x<-2)
> > f1= 4 (x>=-2,x<-1)
> > f1= 5 (x>=-1,x<0)
> > f1= 1 (x>=0, x<1)
> > and so on.
> >
> > What I only have are these vectors that are
> interpreted as functions (as shown above, x (step is 1
> here). I would like to ask your help of how I can create
> some function that reads one of the vector v1,v2,v3 and
> returns results that are acceptable by integrate()
> function.
> >
> > Usually integrate wants a pre-defined function like
> > myfunc(x)<-function{ x^2 }
> > but this is not the case here.
> >
> > Could you please give me some hints how I can
> proceed?
> >
> > I would like to thank u in advance for your help
> >
> > Regards
> > Alex
> >
> > __
> > R-help@r-project.org
> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained,
> reproducible code.
> >
> 
> 
> 
> -- 
> Jim Holtman
> Data Munger Guru
> 
> What is the problem that you are trying to solve?
> 




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


[R] Postscript function Bug at R x64 2.12.1?

2011-01-11 Thread Rodrigo Aluizio
Hi list. I was saving my modified pairs graphic using a custom panel from
the R Graphics site, and I got an interesting difference in the final image
when I save it as eps or png.

This custom panel make possible to show at the left side of the pairs plot
the p-value as symbol and the correlation r value with its cex proportional
to the r value itself.

 

Well, when saving as a png file everything is pretty fine. But, if I save it
as an eps, one of the r values (0.0056) gets extremely bigger than it should
, and it's the only one presenting this behavior.

 

P.S.: I use GIMP 2.6.11 (x86) and GhostScript 9.00 (x86) to render the eps
file, I don't know if it could be an issue of one of these other softwares.

 

Below are the script that generates everything, a link to an example matrix
(csv) and a Sys.info result.

 

The link to an example data

 

http://dl.dropbox.com/u/2613625/Groups.csv

 

The script without the part used to import the data

 

# Analysis

cor.prob <- function(X, dfr = nrow(X) - 2) {

R <- cor(X)

above <- row(R) < col(R)

r2 <- R[above]^2

Fstat <- r2 * dfr / (1 - r2)

R[above] <- 1 - pf(Fstat, 1, dfr)

R

}

cor.prob(Grupos[,-7])

 

# Custom upper pairs panel

panel.cor <- function(x, y, digits=2, prefix="", cex.cor) 

{

usr <- par("usr"); on.exit(par(usr)) 

par(usr = c(0, 1, 0, 1)) 

r <- abs(cor(x, y)) 

txt <- format(c(r, 0.123456789), digits=digits)[1] 

txt <- paste(prefix, txt, sep="") 

if(missing(cex.cor)) cex <- 0.8/strwidth(txt) 

 

test <- cor.test(x,y) 

# borrowed from printCoefmat

Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, 

  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),

  symbols = c("***", "**", "*", ".", " ")) 

 

text(0.5, 0.5, txt, cex = cex * r) 

text(.8, .8, Signif, cex=cex, col=2) 

}

 

# Saving the Graphics

postscript('CorMatrix.eps',paper='special',onefile=F,horizontal=F,bg='white'
,width=7,height=7,pointsize=9) # Here the issue pops up

par(pty='m')

pairs(Grupos[,-7],lower.panel=panel.smooth,upper.panel=panel.cor)

dev.off()

 

png('CorMatrix.png',bg='white',width=2800,height=2800,res=300) # Here
everything goes fine

par(pty='m')

pairs(Grupos[,-7],lower.panel=panel.smooth,upper.panel=panel.cor)

dev.off()

 

# System Information

Sys.info()

 sysname  release  version nodename  machine
login 

   "Windows"  "7 x64" "build 7600" "RODRIGO-PC" "x86-64"
"Rodrigo" 

user 

   "Rodrigo"

 

Hope you can help me find out what is going on!

Regards

---

MSc.   Rodrigo Aluizio


[[alternative HTML version deleted]]

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


[R] Postscript function Bug at R x64 2.12.1?

2011-01-11 Thread Rodrigo Aluizio
Hi, a correction in the previous script. Anyone who tries to use it must
ignore (remove) the [,-7] after the Grupos object.

 

Sorry.

 

De: Rodrigo Aluizio [mailto:r.alui...@gmail.com] 
Enviada em: terça-feira, 11 de janeiro de 2011 11:01
Para: R Help
Assunto: Postscript function Bug at R x64 2.12.1?

 

Hi list. I was saving my modified pairs graphic using a custom panel from
the R Graphics site, and I got an interesting difference in the final image
when I save it as eps or png.

This custom panel make possible to show at the left side of the pairs plot
the p-value as symbol and the correlation r value with its cex proportional
to the r value itself.

 

Well, when saving as a png file everything is pretty fine. But, if I save it
as an eps, one of the r values (0.0056) gets extremely bigger than it should
, and it’s the only one presenting this behavior.

 

P.S.: I use GIMP 2.6.11 (x86) and GhostScript 9.00 (x86) to render the eps
file, I don’t know if it could be an issue of one of these other softwares.

 

Below are the script that generates everything, a link to an example matrix
(csv) and a Sys.info result.

 

The link to an example data

 

http://dl.dropbox.com/u/2613625/Groups.csv

 

The script without the part used to import the data

 

# Analysis

cor.prob <- function(X, dfr = nrow(X) - 2) {

R <- cor(X)

above <- row(R) < col(R)

r2 <- R[above]^2

Fstat <- r2 * dfr / (1 - r2)

R[above] <- 1 - pf(Fstat, 1, dfr)

R

}

cor.prob(Grupos[,-7])

 

# Custom upper pairs panel

panel.cor <- function(x, y, digits=2, prefix="", cex.cor) 

{

usr <- par("usr"); on.exit(par(usr)) 

par(usr = c(0, 1, 0, 1)) 

r <- abs(cor(x, y)) 

txt <- format(c(r, 0.123456789), digits=digits)[1] 

txt <- paste(prefix, txt, sep="") 

if(missing(cex.cor)) cex <- 0.8/strwidth(txt) 

 

test <- cor.test(x,y) 

# borrowed from printCoefmat

Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, 

  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),

  symbols = c("***", "**", "*", ".", " ")) 

 

text(0.5, 0.5, txt, cex = cex * r) 

text(.8, .8, Signif, cex=cex, col=2) 

}

 

# Saving the Graphics

postscript('CorMatrix.eps',paper='special',onefile=F,horizontal=F,bg='white'
,width=7,height=7,pointsize=9) # Here the issue pops up

par(pty='m')

pairs(Grupos[,-7],lower.panel=panel.smooth,upper.panel=panel.cor)

dev.off()

 

png('CorMatrix.png',bg='white',width=2800,height=2800,res=300) # Here
everything goes fine

par(pty='m')

pairs(Grupos[,-7],lower.panel=panel.smooth,upper.panel=panel.cor)

dev.off()

 

# System Information

Sys.info()

 sysname  release  version nodename  machine
login 

   "Windows"  "7 x64" "build 7600" "RODRIGO-PC" "x86-64"
"Rodrigo" 

user 

   "Rodrigo"

 

Hope you can help me find out what is going on!

Regards

---

MSc. Rodrigo Aluizio  


[[alternative HTML version deleted]]

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


[R] using R on an php webapp

2011-01-11 Thread mihai

i write an php R web app on a linux server that takes an user input matrix
(20X20) and makes an Multidimensional Scaling calculation. Now my answer is,
can R support (let's say) 1000 requests per second? (my app is intended to
be use on facebook) What's the best way, use R as a server, deployng more
than one server with R? i must say that i'ts the first time that i use R on
an app. 
thanks 
mihai
-- 
View this message in context: 
http://r.789695.n4.nabble.com/using-R-on-an-php-webapp-tp3208868p3208868.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Confidence interval on quantile regression predictions

2011-01-11 Thread Davey, Andrew
I am using the quantreg package to build a quantile regression model and
wish to generate confidence intervals for the fitted values.

After fitting the model, I have tried running predict() and
predict.rq(), but in each case I obtain a vector of the fitted values
only.

For example:

library(quantreg)
y<-rnorm(50,10,2)
x<-seq(1,50,1)
model<-rq(y~x,tau=0.5,method="br",model=TRUE)
model$fitted
predict(model,type=c("percentile"),interval=c("confidence"),level=0.95)
#returns same output as model$fitted

Am I missing something obvious? I have tried altering the settings for
the type and interval arguments but without success.

(I am running R v 2.10.1).

Thank you for your help.

Regards,


Andrew Davey
Senior Consultant

Tel:+ 44 (0) 1793 865023
Mob:+44 (0) 7747 863190
Fax:+ 44 (0) 1793 865001
E-mail: andrew.da...@wrcplc.co.uk
Website:www.wrcplc.co.uk
P Before printing, please think about the environment.

---
To read our latest newsletter visit 
http://www.wrcplc.co.uk/default.aspx?item=835 - Keeping our clients up-to-date 
with WRc's business developments
---
Visit our websites www.wrcplc.co.uk and www.wrcnsf.com, as well as 
www.waterportfolio.com for collaborative research projects.
---
The Information in this e-mail is confidential and is intended solely for the 
addressee. Access to this e-mail by any other party is unauthorised. If you are 
not the intended recipient, any disclosure, copying, distribution or any action 
taken in reliance on the information contained in this e-mail is prohibited and 
maybe unlawful. When addressed to WRc Clients, any opinions or advice contained 
in this e-mail are subject to the terms and conditions expressed in the 
governing WRc Client Agreement.
---
WRc plc is a company registered in England and Wales. Registered office 
address: Frankland Road, Blagrove, Swindon, Wiltshire SN5 8YF. Company 
registration number 2262098. VAT number 527 1804 53.
---
[[alternative HTML version deleted]]

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


Re: [R] Real time dataset

2011-01-11 Thread Daróczi Gergely
Hello,

this little example my work for demo purposes:


x1 = runif(45)
x2 = runif(45)
x3 = runif(45)
df = data.frame(x1=x1,x2=x2,x3=x3)

test <- function(df, values, n) {
if (nrow(df) < n) {
 df <- rbind(df, values)
} else {
df[1:(nrow(df)-1),] <- df[2:nrow(df),]
 df[n,] <- values
}
return(df)
}

for (i in 1:10) {
 df <- test(df, rep(1, ncol(df)), 50)
}

The demo is not exactly what you specified, as a dataframe cannot have
different number of observations in columns/variables, so I modified it to
update it with a "full" observation (=row) each run.

Gergely

On Tue, Jan 11, 2011 at 6:14 AM, antujsrv  wrote:

>
> my dataset looks like :
>
> >df
>
>   VIX   GLD   FAS
> 12 4  5
> 28 9  10
> 356   9   98 .. continued
>
> the dataset has n observations which is fixed and i need to create a
> function :
>
> test <- function(variable name, value) -- this function has to insert the
> value under the respective variable and if the column length exceeds "n"
> then automatically the first value from the top of the variable should get
> deleted, the way we have in stacking-first in first out.
>
> Any help writing this function would be highly appreciated
> Thanks
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Real-time-dataset-tp3208473p3208473.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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


Re: [R] Meaning of pterms in survreg object?

2011-01-11 Thread Terry Therneau
I'm having a hard time fathoming your question, since "pterms" appears
nowhere in my source code to residuals.survreg.  Can we start at the
beginning, per the posting guide?

 1. What version of R are you running?
 2. What exactly are the statements you typed which give the message?

Note that residuals are not saved in a fit, but are computed later.
  fit <- survreg(Surv(time, status) ~ age + height, data=mydata)
  resid(fit, type='response')

Terry Therneau

--- begin included message 
Residual and predict methods are not carried out on the
survreg
object, although the component linear.predictors exists for
the
survreg object. Looking at the code I see that the residual
method 
refuses to run if any of the components of the pterms vector
is
equal to
2.  
  
--- end inclusion ---

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


Re: [R] Global variables

2011-01-11 Thread Sebastien Bihorel

Thanks,

I will have a look at it.

Sebastien

Michael Bedward wrote:

Hi Sebastian,

You might also find the proto package useful as a way of restricting
the scope of variables. It provides a more intuitive (at least to me)
way of packaging variables and functions up into environments that can
be related in a hierarchy.

Michael

On 10 January 2011 23:48, Sebastien Bihorel
 wrote:
  

Thank Gabor and Duncan,

That will be helpful.

Gabor Grothendieck wrote:


On Thu, Jan 6, 2011 at 4:59 PM, Duncan Murdoch  wrote:

  

On 06/01/2011 4:45 PM, Sebastien Bihorel wrote:



Dear R-users,

Is there a way I can prevent global variables to be visible within my
functions?

  

Yes, but you probably shouldn't.  You would do it by setting the environment
of the function to something that doesn't have the global environment as a
parent, or grandparent, etc.  The only common examples of that are baseenv()
and emptyenv().  For example,

x <- 1
f <- function() print(x)

Then f() will work, and print the 1.  But if I do

environment(f) <- baseenv()

then it won't work:




f()

  

Error in print(x) : object 'x' not found

The problem with doing this is that it is not the way users expect functions
to work, and it will probably have weird side effects.  It is not the way
things work in packages (even packages with namespaces will eventually
search the global environment, the namespace just comes first).  There's no
simple way to do it and yet get access to functions in other packages
besides base without explicitly specifying them (e.g. you'd need to use
stats::lm(), not just lm(), etc.)




A variation of this would be:

environment(f) <- as.environment(2)

which would skip over the global environment, .GlobEnv, but would
still search the loaded packages.  In the example above x would not be
found but it still could find lm, etc.



  

   [[alternative HTML version deleted]]

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




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


Re: [R] Confidence interval on quantile regression predictions

2011-01-11 Thread Roger Koenker
You need to add explicitly newdata = list(x=x) and if you want percentile 
method you also
need se = "boot".  


Roger Koenker
rkoen...@illinois.edu

On Jan 11, 2011, at 3:54 AM, Davey, Andrew wrote:

> I am using the quantreg package to build a quantile regression model and
> wish to generate confidence intervals for the fitted values.
> 
> After fitting the model, I have tried running predict() and
> predict.rq(), but in each case I obtain a vector of the fitted values
> only.
> 
> For example:
> 
> library(quantreg)
> y<-rnorm(50,10,2)
> x<-seq(1,50,1)
> model<-rq(y~x,tau=0.5,method="br",model=TRUE)
> model$fitted
> predict(model,type=c("percentile"),interval=c("confidence"),level=0.95)
> #returns same output as model$fitted
> 
> Am I missing something obvious? I have tried altering the settings for
> the type and interval arguments but without success.
> 
> (I am running R v 2.10.1).
> 
> Thank you for your help.
> 
> Regards,
> 
> 
> Andrew Davey
> Senior Consultant
> 
> Tel:  + 44 (0) 1793 865023
> Mob:+44 (0) 7747 863190
> Fax:  + 44 (0) 1793 865001
> E-mail:   andrew.da...@wrcplc.co.uk
> Website:  www.wrcplc.co.uk
> P Before printing, please think about the environment.
> 
> ---
> To read our latest newsletter visit 
> http://www.wrcplc.co.uk/default.aspx?item=835 - Keeping our clients 
> up-to-date with WRc's business developments
> ---
> Visit our websites www.wrcplc.co.uk and www.wrcnsf.com, as well as 
> www.waterportfolio.com for collaborative research projects.
> ---
> The Information in this e-mail is confidential and is intended solely for the 
> addressee. Access to this e-mail by any other party is unauthorised. If you 
> are not the intended recipient, any disclosure, copying, distribution or any 
> action taken in reliance on the information contained in this e-mail is 
> prohibited and maybe unlawful. When addressed to WRc Clients, any opinions or 
> advice contained in this e-mail are subject to the terms and conditions 
> expressed in the governing WRc Client Agreement.
> ---
> WRc plc is a company registered in England and Wales. Registered office 
> address: Frankland Road, Blagrove, Swindon, Wiltshire SN5 8YF. Company 
> registration number 2262098. VAT number 527 1804 53.
> ---
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] Interpolate xts

2011-01-11 Thread Joshua Ulrich
Use na.approx:

set.seed(21)
x <- xts(rnorm(10), Sys.time()-10:1)
is.na(x) <- 2:4
is.na(x) <- 8:9
na.approx(x)
na.spline(x)

Best,
--
Joshua Ulrich  |  FOSS Trading: www.fosstrading.com



On Tue, Jan 11, 2011 at 12:08 AM, Rustamali Manesiya
 wrote:
> Hello,
>
>       I have a xts object, I would like to fill the NA with linear
> interpolated data. Can anyone please help.
>
>> str(zz)
> An ‘xts’ object from 2010-11-24 15:59:29 to 2010-11-24 16:00:00 containing:
>  Data: num [1:23401, 1] 312 312 312 312 312 ...
>  Indexed by objects of class: [POSIXct,POSIXt] TZ:
>  xts Attributes:
> List of 2
>  $ src    : chr "datafeed"
>  $ updated: POSIXct[1:1], format: "2011-01-08 00:33:05"
>
>>zz
> 2010-11-24 15:59:29 315.0300
> 2010-11-24 15:59:30       NA
> 2010-11-24 15:59:31       NA
> 2010-11-24 15:59:32       NA
> 2010-11-24 15:59:33       NA
> 2010-11-24 16:00:00 314.7000
>
>
> Rusty
>
>        [[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

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


[R] Catch separability problems in a (multinom) regression

2011-01-11 Thread sovo0815
So far I have "learned" that a binary logit model can either not 
converge or show linear separability (two different warnings). And 
so far I assume that this is also the case for multinomial models. 
For these models, if such errors occur, the coefficients are not 
very useful. In a loop (bootstrapping XYZ repetitions), I can 
catch the convergence problem via the $convergence variable 
(set to 1). But how can I catch cases of separability?


library(nnet)
y <- factor(rep(LETTERS[1:3], each=20))
a <- c(rep(1:4, each=10), rep(5:6, 10))
b <- c(rep(1:2, 20), rep(3:4, each=10))
mod <- multinom(y ~ a + b, model=T)


--
Sören Vogel, sovo0...@gmail.com, http://sovo0815.wordpress.com/

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


Re: [R] rprofile.site

2011-01-11 Thread PtitBleu

Hello Julian,

I'm not sure my answer will be very helpful to you but it seems that Tinn-R
looks for Rprofile.site at the old place (R/bin/etc). With the new R
configuration, this file is in now R/bin. 
Maybe you can try to copy the Rprofil.site file to the place expected by
Tinn-R.

I also found this post  http://r.789695.n4.nabble.com/Tinn-R-td2952577.html
http://r.789695.n4.nabble.com/Tinn-R-td2952577.html 

Personaly, after a lot of tests, I have a running Tinn-R 1.19.4.7 + R 2.12.1
(Rgui) configuration.

I don't use the last version of Tinn-R because the re-call of the command
doesn't work for me (I have the same message than you have, that is
source(.trPaths[5], echo=TRUE, max.deparse.length=150)).
 
At the beginning, I had problems because again Tinn-R was looking Rgui at
the previous place and I had to give the path each time I started Tinn-R.
But changing the 'InstallPath' value in the
HKEY_LOCAL_MACHINE/SOFTWARE/R-core/R32/2.12.1 entry of the register solved
my problem (I put R\bin\i386 instead of R\bin)

Tinn-R is now able to keep the path in the .ini file and then is able to
directly connect to Rgui.
My Tinn-R 1.19.4.7 + Rgui 2.12.1 configuration is now ok (it matches my
expectations at least).  

I hope you will find a way (mine or another one) for your problem.
Ptit Bleu.
  
-- 
View this message in context: 
http://r.789695.n4.nabble.com/rprofile-site-tp3193599p3209134.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Help with Data Transformation

2011-01-11 Thread Hadley Wickham
> The data is initially extracted from an SQL database into Excel, then saved 
> as a tab-delimited text file for use in R.

You might also want to look at the SQL packages for R so you can skip
this manual step. I'd recommend starting with
http://cran.r-project.org/doc/manuals/R-data.html#Relational-databases

Hadley

-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

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


[R] list concatenation

2011-01-11 Thread Georg Otto
Dear R gurus,


first let me apologize for a question that might hve been answered
before. I was not able to find the solution yet. I want to concatenate
two lists of lists at their lowest level.

Suppose I have two lists of lists:

list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),

   "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))


list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),
   
   "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))



> list.1
$I
$I$A
[1] "a" "b" "c"

$I$B
[1] "d" "e" "f"


$II
$II$A
[1] "g" "h" "i"

$II$B
[1] "j" "k" "l"


> list.2
$I
$I$A
[1] "A" "B" "C"

$I$B
[1] "D" "E" "F"


$II
$II$A
[1] "G" "H" "I"

$II$B
[1] "J" "K" "L"


Now I want to concatenate list elements of the lowest levels, so the
result looks like this:


$I
$I$A
[1] "a" "b" "c" "A" "B" "C"

$I$B
[1] "d" "e" "f" "D" "E" "F"


$II
$II$A
[1] "g" "h" "i" "G" "H" "I"

$II$B
[1] "j" "k" "l" "J" "K" "L"


Has anybody a good solution for that?

Best,

Georg

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

2011-01-11 Thread Jason Edgecombe

Hello,

I have logging information for multiple machines, which I am trying to 
summarize and graph. So far, I process each host individually, but I 
would like to summarize the user count across multiple hosts. I want to 
answer the question "how many unique users logged in on a certain day 
across a group of machines"?


I'm not quite sure how to scale the data frame and analysis to summarize 
multiple hosts, though. I'm still getting a feel for using R.


Here is a snippet of data for one host. the user_count column is 
generated from the users column using my custom function "usercount()". 
the samples are taken roughly once per minute and only unique samples 
are recorded. (i.e. use na.locf() to uncompress the data). Samples may 
occur twice in the same minute and are rarely aligned on the same time.


Here is the original data before I turn t into a zoo series and run 
na.locf() over it so I can aggregate a single host by day. I'm open to a 
better way.

> foo
  usersdatetime user_count
1 user1 & user2 2007-03-29 19:16:30  2
2 user1 & user2 2007-03-31 00:04:46  2
3 user1 & user2 2007-04-02 11:49:20  2
4 user1 & user2 2007-04-02 12:02:04  2
5 user1 & user2 2007-04-02 12:44:02  2
6 user1 & user2 & user3 2007-04-02 16:34:05  3

> dput(foo)
structure(list(users = c("user1 & user2", "user1 & user2", "user1 & user2",
"user1 & user2", "user1 & user2", "user1 & user2 & user3"), datetime = 
structure(c(1175210190,
1175313886, 1175528960, 1175529724, 1175532242, 1175546045), class = 
c("POSIXt",

"POSIXct"), tzone = "US/Eastern"), user_count = c(2, 2, 2, 2,
2, 3)), .Names = c("users", "datetime", "user_count"), row.names = c(NA,
6L), class = "data.frame")


Thanks,
Jason

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


Re: [R] list concatenation

2011-01-11 Thread Joshua Wiley
Hi Georg,

This was an interesting challenge to me.  Here's what I came up with.
The first option meets your desired result, but could get messy with
deeper nesting.  The second is less code, but is not quite what you
want and requires as.data.frame() to give a reasonable result for each
list.  Calling either option a "good solution" would be rather
generous.  I'm iinterested to see what other people do.

Josh

## Data
list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),
  "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))
list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),
  "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))

## Try 1
list.t1 <- list.1

for(i in length(list.1)) {
  for(j in length(list.1[[i]])) {
list.t1[[c(i, j)]] <- c(list.1[[c(i, j)]], list.2[[c(i, j)]])
  }
}

## Try 2
list.t2 <- as.list(do.call("rbind", lapply(list(list.1, list.2),
  as.data.frame, stringsAsFactors = FALSE)))

## Results

list.t1
list.t2


On Tue, Jan 11, 2011 at 7:44 AM, Georg Otto  wrote:
> Dear R gurus,
>
>
> first let me apologize for a question that might hve been answered
> before. I was not able to find the solution yet. I want to concatenate
> two lists of lists at their lowest level.
>
> Suppose I have two lists of lists:
>
> list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),
>
>               "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))
>
>
> list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),
>
>               "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))
>
>
>
>> list.1
> $I
> $I$A
> [1] "a" "b" "c"
>
> $I$B
> [1] "d" "e" "f"
>
>
> $II
> $II$A
> [1] "g" "h" "i"
>
> $II$B
> [1] "j" "k" "l"
>
>
>> list.2
> $I
> $I$A
> [1] "A" "B" "C"
>
> $I$B
> [1] "D" "E" "F"
>
>
> $II
> $II$A
> [1] "G" "H" "I"
>
> $II$B
> [1] "J" "K" "L"
>
>
> Now I want to concatenate list elements of the lowest levels, so the
> result looks like this:
>
>
> $I
> $I$A
> [1] "a" "b" "c" "A" "B" "C"
>
> $I$B
> [1] "d" "e" "f" "D" "E" "F"
>
>
> $II
> $II$A
> [1] "g" "h" "i" "G" "H" "I"
>
> $II$B
> [1] "j" "k" "l" "J" "K" "L"
>
>
> Has anybody a good solution for that?
>
> Best,
>
> Georg
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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


Re: [R] scaling to multiple data files

2011-01-11 Thread jim holtman
I am not sure exactly what your data represents.  For example, from
looking at the data it appears that user1 and user2 have been logged
on for about 4 days; is that what the data is saying?  If you are
keeping track of users, why not write out a file that has the
start/end time for each user's session.  The first time you see them,
put an entry in a table and as soon as they don't show up in your
sample, write out a record for them.  With that information is it easy
to create a report of the number of unique people over time.

On Tue, Jan 11, 2011 at 10:47 AM, Jason Edgecombe
 wrote:
> Hello,
>
> I have logging information for multiple machines, which I am trying to
> summarize and graph. So far, I process each host individually, but I would
> like to summarize the user count across multiple hosts. I want to answer the
> question "how many unique users logged in on a certain day across a group of
> machines"?
>
> I'm not quite sure how to scale the data frame and analysis to summarize
> multiple hosts, though. I'm still getting a feel for using R.
>
> Here is a snippet of data for one host. the user_count column is generated
> from the users column using my custom function "usercount()". the samples
> are taken roughly once per minute and only unique samples are recorded.
> (i.e. use na.locf() to uncompress the data). Samples may occur twice in the
> same minute and are rarely aligned on the same time.
>
> Here is the original data before I turn t into a zoo series and run
> na.locf() over it so I can aggregate a single host by day. I'm open to a
> better way.
>> foo
>                  users            datetime user_count
> 1         user1 & user2 2007-03-29 19:16:30          2
> 2         user1 & user2 2007-03-31 00:04:46          2
> 3         user1 & user2 2007-04-02 11:49:20          2
> 4         user1 & user2 2007-04-02 12:02:04          2
> 5         user1 & user2 2007-04-02 12:44:02          2
> 6 user1 & user2 & user3 2007-04-02 16:34:05          3
>
>> dput(foo)
> structure(list(users = c("user1 & user2", "user1 & user2", "user1 & user2",
> "user1 & user2", "user1 & user2", "user1 & user2 & user3"), datetime =
> structure(c(1175210190,
> 1175313886, 1175528960, 1175529724, 1175532242, 1175546045), class =
> c("POSIXt",
> "POSIXct"), tzone = "US/Eastern"), user_count = c(2, 2, 2, 2,
> 2, 3)), .Names = c("users", "datetime", "user_count"), row.names = c(NA,
> 6L), class = "data.frame")
>
>
> Thanks,
> Jason
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Jim Holtman
Data Munger Guru

What is the problem that you are trying to solve?

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


Re: [R] Match numeric vector against rows in a matrix?

2011-01-11 Thread Kevin Ummel
Hi Petr,

Sorry for the delay. I ended up implementing a simple and fast but 
not-universal or elegant solution, since it turned out that I only needed to 
test a rather small subset of combinations for my purposes.

Nevertheless, I did go back and try the solutions posted. Bill's 'binary 
expansion' approach was the fastest by about a factor of 8.

Cheers,
Kevin

On Jan 11, 2011, at 9:20 AM, Petr Savicky wrote:

> Dear Kevin Ummel:
> 
> There were several suggestions on R-help concering your question below.
> None of them suggests a base function. I would also expect that
> there is a base function for matching the rows. A related function is
> dist(), but it takes only one matrix as input and computes the distances
> between all pairs of its rows. Matching only some rows would require
> a modification of dist(), which would take two matrices and compare
> rows of one of them to the rows of the other. However, i do not know
> such a function in R.
> 
> Is some of the suggestions on R-help suitable for your purposes?
> 
> Thank you in advance for your kind reply.
> 
> Best regards, Petr Savicky.
> 
> On Wed, Jan 05, 2011 at 07:16:47PM +, Kevin Ummel wrote:
>> Two posts in one day is not a good day...and this question seems like it 
>> should have an obvious answer:
>> 
>> I have a matrix where rows are unique combinations of 1's and 0's:
>> 
>>> combs=as.matrix(expand.grid(c(0,1),c(0,1)))
>>> combs
>> Var1 Var2
>> [1,]00
>> [2,]10
>> [3,]01
>> [4,]11
>> 
>> I want a single function that will give the row index containing an exact 
>> match with vector x:
>> 
>>> x=c(0,1)
>> 
>> The solution needs to be applied many times, so I need something quick -- I 
>> was hoping a base function would do it, but I'm drawing a blank.
>> 
>> Thanks!
>> Kevin
>> 
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>> 

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


Re: [R] list concatenation

2011-01-11 Thread Charles C. Berry

On Tue, 11 Jan 2011, Georg Otto wrote:


Dear R gurus,


first let me apologize for a question that might hve been answered
before. I was not able to find the solution yet. I want to concatenate
two lists of lists at their lowest level.

Suppose I have two lists of lists:

list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),

  "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))


list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),

  "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))




list.1

$I
$I$A
[1] "a" "b" "c"

$I$B
[1] "d" "e" "f"


$II
$II$A
[1] "g" "h" "i"

$II$B
[1] "j" "k" "l"



list.2

$I
$I$A
[1] "A" "B" "C"

$I$B
[1] "D" "E" "F"


$II
$II$A
[1] "G" "H" "I"

$II$B
[1] "J" "K" "L"


Now I want to concatenate list elements of the lowest levels, so the
result looks like this:


$I
$I$A
[1] "a" "b" "c" "A" "B" "C"

$I$B
[1] "d" "e" "f" "D" "E" "F"


$II
$II$A
[1] "g" "h" "i" "G" "H" "I"

$II$B
[1] "j" "k" "l" "J" "K" "L"


Has anybody a good solution for that?




mapply( function(x,y) mapply(c, x, y, SIMPLIFY=FALSE),

+ list.1, list.2, SIMPLIFY=FALSE )
$I
$I$A
[1] "a" "b" "c" "A" "B" "C"

$I$B
[1] "d" "e" "f" "D" "E" "F"


$II
$II$A
[1] "g" "h" "i" "G" "H" "I"

$II$B
[1] "j" "k" "l" "J" "K" "L"

HTH,

Chuck








Best,

Georg

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



Charles C. BerryDept of Family/Preventive Medicine
cbe...@tajo.ucsd.eduUC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

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


Re: [R] list concatenation

2011-01-11 Thread Bert Gunter
Lists are (isomorphic to) trees with (possibly) labelled nodes. A
completely general solution in which two trees have possibly different
topologies and different labels would therefore involve identifying
the paths to leaves on each tree, e.g. via depth first search using
recursion, and unioning leaves with the same paths (which could be
quickly found in R via match() on the paths). This is a standard
exercise in a data structures course.

Considerable simplification could be effected if tree topologies
and/or labels are identical or have other restrictions on them.
However, you have not made it clear in your post whether this is the
case (it is in your example).

-- Bert

On Tue, Jan 11, 2011 at 8:04 AM, Joshua Wiley  wrote:
> Hi Georg,
>
> This was an interesting challenge to me.  Here's what I came up with.
> The first option meets your desired result, but could get messy with
> deeper nesting.  The second is less code, but is not quite what you
> want and requires as.data.frame() to give a reasonable result for each
> list.  Calling either option a "good solution" would be rather
> generous.  I'm iinterested to see what other people do.
>
> Josh
>
> ## Data
> list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),
>              "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))
> list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),
>              "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))
>
> ## Try 1
> list.t1 <- list.1
>
> for(i in length(list.1)) {
>  for(j in length(list.1[[i]])) {
>    list.t1[[c(i, j)]] <- c(list.1[[c(i, j)]], list.2[[c(i, j)]])
>  }
> }
>
> ## Try 2
> list.t2 <- as.list(do.call("rbind", lapply(list(list.1, list.2),
>  as.data.frame, stringsAsFactors = FALSE)))
>
> ## Results
>
> list.t1
> list.t2
>
>
> On Tue, Jan 11, 2011 at 7:44 AM, Georg Otto  wrote:
>> Dear R gurus,
>>
>>
>> first let me apologize for a question that might hve been answered
>> before. I was not able to find the solution yet. I want to concatenate
>> two lists of lists at their lowest level.
>>
>> Suppose I have two lists of lists:
>>
>> list.1 <- list("I"=list("A"=c("a", "b", "c"), "B"=c("d", "e", "f")),
>>
>>               "II"=list("A"=c("g", "h", "i"), "B"=c("j", "k", "l")))
>>
>>
>> list.2 <- list("I"=list("A"=c("A", "B", "C"), "B"=c("D", "E", "F")),
>>
>>               "II"=list("A"=c("G", "H", "I"), "B"=c("J", "K", "L")))
>>
>>
>>
>>> list.1
>> $I
>> $I$A
>> [1] "a" "b" "c"
>>
>> $I$B
>> [1] "d" "e" "f"
>>
>>
>> $II
>> $II$A
>> [1] "g" "h" "i"
>>
>> $II$B
>> [1] "j" "k" "l"
>>
>>
>>> list.2
>> $I
>> $I$A
>> [1] "A" "B" "C"
>>
>> $I$B
>> [1] "D" "E" "F"
>>
>>
>> $II
>> $II$A
>> [1] "G" "H" "I"
>>
>> $II$B
>> [1] "J" "K" "L"
>>
>>
>> Now I want to concatenate list elements of the lowest levels, so the
>> result looks like this:
>>
>>
>> $I
>> $I$A
>> [1] "a" "b" "c" "A" "B" "C"
>>
>> $I$B
>> [1] "d" "e" "f" "D" "E" "F"
>>
>>
>> $II
>> $II$A
>> [1] "g" "h" "i" "G" "H" "I"
>>
>> $II$B
>> [1] "j" "k" "l" "J" "K" "L"
>>
>>
>> Has anybody a good solution for that?
>>
>> Best,
>>
>> Georg
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
>
> --
> Joshua Wiley
> Ph.D. Student, Health Psychology
> University of California, Los Angeles
> http://www.joshuawiley.com/
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Bert Gunter
Genentech Nonclinical Biostatistics

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 specification where response is a 2col matrix

2011-01-11 Thread Uwe Ligges

Hi,

when I apply a glm() model in two ways,
first with the response in a two column matrix specification with 
successes and failures


y <- matrix(c(
5, 1,
3, 3,
2, 2,
0, 4), ncol=2, byrow=TRUE)

X <- data.frame(x1 = factor(c(1,1,0,0)),
x2 = factor(c(0,1,0,1)))

glm(y ~ x1 + x2, data = X, family="binomial")


second with a model matrix that full rows (i.e. has as many rows as real 
observations) and represents identical data:



Xf <- data.frame(x1 = factor(rep(c(1,1,0,0), rowSums(y))),
 x2 = factor(rep(c(0,1,0,1), rowSums(y
yf <- factor(rep(rep(0:1, 4), t(y)))

glm(yf ~ x1 + x2, data = Xf, family="binomial")


we will find that the number of degrees of freedom and the AIC etc. 
differ --  I'd expect them to be identical (as the coefficient estimates 
and such things are).


maybe I am confused tonight, hence I do not file it as a bug report 
right away and wait to be enlightened ...



Thanks and best wishes,
Uwe

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 specification where response is a 2col matrix

2011-01-11 Thread Prof Brian Ripley
Your first model is a binomial glm witb 4 observations of 6,6,4,4 
trials.


Your second model is a Bernoulli glm with 20 observations of one trial 
each.


The saturated models are different, as are the likelihoods 
(unsurprising given the data is different): the binomial model has 
comnbinarial factors (e.g. choose(10,5)*choose(6,3)*choose(4,2)) that 
the Bernoulli model does not have, so the AICs differ.


I am not sure where these issues of aggregating Bernoulli trials is 
explained (nor am I near my books), but this is a common question.


On Tue, 11 Jan 2011, Uwe Ligges wrote:


Hi,

when I apply a glm() model in two ways,
first with the response in a two column matrix specification with successes 
and failures


y <- matrix(c(
   5, 1,
   3, 3,
   2, 2,
   0, 4), ncol=2, byrow=TRUE)

X <- data.frame(x1 = factor(c(1,1,0,0)),
   x2 = factor(c(0,1,0,1)))

glm(y ~ x1 + x2, data = X, family="binomial")


second with a model matrix that full rows (i.e. has as many rows as real 
observations) and represents identical data:



Xf <- data.frame(x1 = factor(rep(c(1,1,0,0), rowSums(y))),
x2 = factor(rep(c(0,1,0,1), rowSums(y
yf <- factor(rep(rep(0:1, 4), t(y)))

glm(yf ~ x1 + x2, data = Xf, family="binomial")


we will find that the number of degrees of freedom and the AIC etc. differ -- 
I'd expect them to be identical (as the coefficient estimates and such things 
are).


maybe I am confused tonight, hence I do not file it as a bug report right 
away and wait to be enlightened ...



Thanks and best wishes,
Uwe


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

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


[R] adding plots to existent multiuple plots

2011-01-11 Thread Jannis
Dear list members,


does anybody know a way how to add plots to already existing subplots created 
with layout? I have created two subplots via

layout(matrix(1:2))

and would like to add different plots to each of them. As I need to load data 
in between I need to add a plot to the first subplot, then the second, load 
data, and add a plot to the first again, followed by the second. Therefore I 
can not plot all plots in subplot 1 first, change the active plot and fill 
subplot 2.

Is there a way to set the subplots to be the active plotting region repeatedly?


Thanks in advance
Jannis



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


Re: [R] list concatenation

2011-01-11 Thread David Katz

Without any error checking, I'd try this:

recurse <-
function(l1,l2){
if(is(l1[[1]],"list"))
mapply(recurse,l1,l2,SIMPLIFY=F) else
{mapply(c,l1,l2,SIMPLIFY=F)}
}

recurse(list.1,list.2)

which recursively traverses each tree (list) in tandem until one is not
composed of lists and then concatenates. If the structure of the lists does
not agree, this will fail.

Regards,

David Katz
da...@davidkatzconsulting.com
-- 
View this message in context: 
http://r.789695.n4.nabble.com/list-concatenation-tp3209182p3209324.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Writing diagonal matrix in opposite direction

2011-01-11 Thread Ron Michael
Hi, is there any direct R function to write an diagonal matrix in an opposite 
way? for example I want to get like:
 
> diag(rnorm(5))[,5:1]
   [,1]   [,2]   [,3]  [,4]   [,5]
[1,]  0.000  0.000  0.000  0.00 -0.1504687
[2,]  0.000  0.000  0.000 -2.139669  0.000
[3,]  0.000  0.000 -0.2102133  0.00  0.000
[4,]  0.000 -0.2609686  0.000  0.00  0.000
[5,] -0.6818889  0.000  0.000  0.00  0.000

Is there any typical name for this type of matrix?
 
Thanks,


[[alternative HTML version deleted]]

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


[R] how to sort new data frame based on the original data frame

2011-01-11 Thread wangwallace

I have a really simple question

I have a data frame of 8 variables (the first column is the subjects' id): 

SubID G1G2 G3 G4W1W2  W3W4 
  1  6  5   6   2  6  22   4 
  2  6  4   7   2  6  62   3 
  3  5  5   5   5  5  54   5 
  4  5  4   3   4  4  45   2 
  5  5  6   7   5  6  44   1 
  6  5  4   3   6  4  37   3 
  7  3  6   6   3  6  52   1 
  8  3  6   6   3  6  54   7 


this data frame have two sets of variables. each set simply represent one
scale. as shown above, the first scale, say G, consists of four items: G1,
G2, G3, and G4, whereas the second scale, say W, also has four items: W1,
W2, W3, W4. 
the leftmost column lists the subjects' ID. 

I drew 100 new random samples based on the data frame. here is the structure
of each new random sample:
   varvar   var var 
g  g  g   g   
g  g  g   g   
w  w w   w   
w  w w   w   
w  w w   w 
w  w w   w 
w  w w   w 
w  w w   w  

each random sample satisfies the following rules:

###the top two rows have to be filled with 2 random rows of the 8 rows of G
numbers. the rest should be filled with 6 random rows of the 8 rows of W
numbers. At the same time, the SubIDs of all eight rows should be different
among each other.

here below is the syntax I've used:

> fff<-function(dat,g=2,w=6){
+ sel1<-sample(1:8,g)
+ sel2<-sample((1:8)[-sel1],w)
+ M=dat[sel1,2:5]
+ N=dat[sel2,6:9]
+ colnames(N)<-colnames(M)
+ rbind(M,N)
+}

> result<-vector("list",100)
> for(i in 1:100)result[[i]]<-fff(data,2,6)
> result

here is the first random sample:

> result[[1]]

   G1G2 G3 G4
3  5  5   5   5
6  5  4   3   6
4  4  4   5   2  
1  6  2   2   4 
7  6  5   2   1
2  6  6   2   3
5  6  4   4   1
8  6  5   4   7

I am wondering how can I sort the rows of each new random samples in the
same order that is corresponding to the SubID in the original data.
Specifically, what kind of syntax should I've added into the one I've used
above to make the random sample, say, the first random sample, look like:

> result[[1]]

   G1G2 G3 G4
1  6  2   2   4 
2  6  6   2   3
3  5  5   5   5
4  4  4   5   2  
5  6  4   4   1
6  5  4   3   6
7  6  5   2   1
8  6  5   4   7

Many thanks!!
-- 
View this message in context: 
http://r.789695.n4.nabble.com/how-to-sort-new-data-frame-based-on-the-original-data-frame-tp3209353p3209353.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Writing diagonal matrix in opposite direction

2011-01-11 Thread RICHARD M. HEIBERGER
> x <- 1:5
> diag(x)[rev(seq(length(x))),]


On Tue, Jan 11, 2011 at 12:56 PM, Ron Michael wrote:

> Hi, is there any direct R function to write an diagonal matrix in an
> opposite way? for example I want to get like:
>
> > diag(rnorm(5))[,5:1]
>[,1]   [,2]   [,3]  [,4]   [,5]
> [1,]  0.000  0.000  0.000  0.00 -0.1504687
> [2,]  0.000  0.000  0.000 -2.139669  0.000
> [3,]  0.000  0.000 -0.2102133  0.00  0.000
> [4,]  0.000 -0.2609686  0.000  0.00  0.000
> [5,] -0.6818889  0.000  0.000  0.00  0.000
>
> Is there any typical name for this type of matrix?
>
> Thanks,
>
>
>[[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

[[alternative HTML version deleted]]

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


[R] Alphabetic labels on multi-plot graphics

2011-01-11 Thread Eduardo de Oliveira Horta
Is there a way to achieve

lbl=c("a", "b", "c", "d")

par(mfrow=c(2,2), ann=FALSE)
for (t in 1:4){
  plot(seq(from=1,to=2*pi,length=100),
sin(t*seq(from=1,to=2*pi,length=100)), type="l")
  title(main=paste("(", lbl[t], ")", sep=""))
}

without having to use an object like 'lbl'?

More generally: is it possible to iteratively (as in a loop) add
alphabetic titles to multi-plot graphics when the range over which 't'
above varies is of an arbitrary length? It is important that
ann=FALSE, because I don't want the axes labels.

Thanks in advance, and best regards.

Eduardo Horta

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


[R] Displaying a PLS biplot under a 3d surface plot

2011-01-11 Thread Lorenzen, Ted {Quaker}
All:

Anyone have an idea about how to plot a PLS biplot on the x-y plane under a 3-D 
surface?  Preferably a 3-D surface such as one created with rgl.surface() 
(meaning, I would love an interactive plot when all is said and done).

I have seen examples which combine a 3-d surface with contour plot (e.g. 
http://stackoverflow.com/questions/1896419/plotting-a-3d-surface-plot-with-contour-map-overlay-using-r)
  . . .  and I imagine I could construct the biplot 'by hand' via plot3d, but 
it'd be awfully handy to be able to take the output of the biplot and tuck it 
under the surface plot.

Thanks for your suggestions,
Ted

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


Re: [R] ThinkCell type waterfall charts in R?

2011-01-11 Thread ang

Hi Jim,

I looked through the plotrix documentation, and the waterfall plot comes
from the stackpoly function right?  I'm not sure if I can modify the
stackpoly to create the plot I want, since stackpoly is a line plot and
fills the area under with color.

I haven't played with all the options yet, but what I was looking for was
more similar to staircase.plot, but instead of horizontal bars, they would
be vertical columns.  Would you happen to know any packages or existing
plots that could be easily modified to do this?

Thanks in advance for your help.

Adrian
-- 
View this message in context: 
http://r.789695.n4.nabble.com/ThinkCell-type-waterfall-charts-in-R-tp881466p3209490.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] how to sort new data frame based on the original data frame

2011-01-11 Thread Joshua Wiley
Hi,

On Tue, Jan 11, 2011 at 9:19 AM, wangwallace  wrote:
>
> I have a really simple question
>
> I have a data frame of 8 variables (the first column is the subjects' id):
>
>    SubID     G1    G2     G3     G4    W1    W2      W3    W4
>      1          6      5       6       2      6      2        2       4
>      2          6      4       7       2      6      6        2       3
>      3          5      5       5       5      5      5        4       5
>      4          5      4       3       4      4      4        5       2
>      5          5      6       7       5      6      4        4       1
>      6          5      4       3       6      4      3        7       3
>      7          3      6       6       3      6      5        2       1
>      8          3      6       6       3      6      5        4       7
>
>
> this data frame have two sets of variables. each set simply represent one
> scale. as shown above, the first scale, say G, consists of four items: G1,
> G2, G3, and G4, whereas the second scale, say W, also has four items: W1,
> W2, W3, W4.
> the leftmost column lists the subjects' ID.
>
> I drew 100 new random samples based on the data frame. here is the structure
> of each new random sample:
>               var    var   var     var
>                g      g      g       g
>                g      g      g       g
>                w      w     w       w
>                w      w     w       w
>                w      w     w       w
>                w      w     w       w
>                w      w     w       w
>                w      w     w       w
>
> each random sample satisfies the following rules:
>
> ###the top two rows have to be filled with 2 random rows of the 8 rows of G
> numbers. the rest should be filled with 6 random rows of the 8 rows of W
> numbers. At the same time, the SubIDs of all eight rows should be different
> among each other.
>
> here below is the syntax I've used:
>
>> fff<-function(dat,g=2,w=6){
> + sel1<-sample(1:8,g)
> + sel2<-sample((1:8)[-sel1],w)
> + M=dat[sel1,2:5]
> + N=dat[sel2,6:9]
> + colnames(N)<-colnames(M)
## just add the order to rbind()
> + rbind(M, N)[order(c(sel1, sel2)), ]
> +}

>> result<-vector("list",100)
>> for(i in 1:100)result[[i]]<-fff(data,2,6)
>> result

For convenience and speed, consider building this (the for loop) into
your function.  The only part that you actually need looped is the
sample(), so you could get some performance gains, if you are
interested/that is an issue.

HTH,

Josh

-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] adding plots to existent multiuple plots

2011-01-11 Thread Jannis
Okseems as I found a way. For anybody looking for a similar solution:

?split.screen


Cheers
Jannis

--- Jannis  schrieb am Di, 11.1.2011:

> Von: Jannis 
> Betreff: [R] adding plots to existent multiuple plots
> An: r-help@r-project.org
> Datum: Dienstag, 11. Januar, 2011 18:13 Uhr
> Dear list members,
> 
> 
> does anybody know a way how to add plots to already
> existing subplots created with layout? I have created two
> subplots via
> 
> layout(matrix(1:2))
> 
> and would like to add different plots to each of them. As I
> need to load data in between I need to add a plot to the
> first subplot, then the second, load data, and add a plot to
> the first again, followed by the second. Therefore I can not
> plot all plots in subplot 1 first, change the active plot
> and fill subplot 2.
> 
> Is there a way to set the subplots to be the active
> plotting region repeatedly?
> 
> 
> Thanks in advance
> Jannis
> 
> 
> 
> __
> R-help@r-project.org
> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
> 



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


Re: [R] Alphabetic labels on multi-plot graphics

2011-01-11 Thread David Winsemius


On Jan 11, 2011, at 2:06 PM, Eduardo de Oliveira Horta wrote:


Is there a way to achieve

lbl=c("a", "b", "c", "d")

opar <- par(mfrow=c(2,2), ann=FALSE)
for (t in 1:4){
 plot(seq(from=1,to=2*pi,length=100),
sin(t*seq(from=1,to=2*pi,length=100)), type="l")
 title(main=paste("(", lbl[t], ")", sep=""))
}
par(opar)


opar <- par(mfrow=c(2,2), ann=FALSE)
 for (t in 1:4){
  plot(seq(from=1,to=2*pi,length=100),
 sin(t*seq(from=1,to=2*pi,length=100)), type="l")
  title(main=bquote("("*.(letters[t])*")") )
 }
par(opar)


without having to use an object like 'lbl'?


(You obviously need something from which to construct the titles.)



More generally: is it possible to iteratively (as in a loop) add
alphabetic titles to multi-plot graphics when the range over which 't'
above varies is of an arbitrary length? It is important that
ann=FALSE, because I don't want the axes labels.

Thanks in advance, and best regards.

Eduardo Horta

__


David Winsemius, MD
West Hartford, CT

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

2011-01-11 Thread Christofer Bogaso
Dear all, I would like to ask one question related to statistics, for
specifically on defining dummy variables. As of now, I have come across 3
different kind of dummy variables (assuming I am working with Seasonal
dummy, and number of season is 4):

> dummy1 <- diag(4)
> for(i in 1:3) dummy1 <- rbind(dummy1, diag(4))
> dummy1 <- dummy1[,-4]
>
> dummy2 <- dummy1
> dummy2[dummy2 == 0] = -1/(4-1)
>
> dummy3 <- dummy1 - 1/4
>
> head(dummy1)
 [,1] [,2] [,3]
[1,]100
[2,]010
[3,]001
[4,]000
[5,]100
[6,]010
> head(dummy2)
   [,1]   [,2]   [,3]
[1,]  1.000 -0.333 -0.333
[2,] -0.333  1.000 -0.333
[3,] -0.333 -0.333  1.000
[4,] -0.333 -0.333 -0.333
[5,]  1.000 -0.333 -0.333
[6,] -0.333  1.000 -0.333
> head(dummy3)
  [,1]  [,2]  [,3]
[1,]  0.75 -0.25 -0.25
[2,] -0.25  0.75 -0.25
[3,] -0.25 -0.25  0.75
[4,] -0.25 -0.25 -0.25
[5,]  0.75 -0.25 -0.25
[6,] -0.25  0.75 -0.25
Now I want to know which type of dummy definition is called Centered dummy
and why it is called so? Is it equivalent to use any of the above
definitions (atleast 2nd and 3rd?) It would really be very helpful if
somebody point any suggestion and clarification.

Thanks and regards,

[[alternative HTML version deleted]]

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


Re: [R] scaling to multiple data files

2011-01-11 Thread Jason Edgecombe
That's correct, those users have been logged in or had processes running 
on this machine for four days. The machines in question are time-sharing 
Linux servers for college students and professors to use. multi-day jobs 
are common.


The "last" command does what you suggest, but it doesn't capture 
processes left running in the background when a user logs out.  This 
data is simpler and collected across Windows and Linux hosts. Sessions 
are somewhat ambiguous. We just care who is running processes on a 
machine at a certain time. We also collect the process name of processes 
that each user is running so that we can gauge how often applications 
are being used, and by whom. For this analysis, I'm not worried about 
which processes are running, only unique users per day. i have almost 
four years of historical data for some machines in this format.


We have multiple tools written in different languages that parse this 
data. I'm writing one that does better graphing.


On 01/11/2011 11:39 AM, jim holtman wrote:

I am not sure exactly what your data represents.  For example, from
looking at the data it appears that user1 and user2 have been logged
on for about 4 days; is that what the data is saying?  If you are
keeping track of users, why not write out a file that has the
start/end time for each user's session.  The first time you see them,
put an entry in a table and as soon as they don't show up in your
sample, write out a record for them.  With that information is it easy
to create a report of the number of unique people over time.

On Tue, Jan 11, 2011 at 10:47 AM, Jason Edgecombe
  wrote:

Hello,

I have logging information for multiple machines, which I am trying to
summarize and graph. So far, I process each host individually, but I would
like to summarize the user count across multiple hosts. I want to answer the
question "how many unique users logged in on a certain day across a group of
machines"?

I'm not quite sure how to scale the data frame and analysis to summarize
multiple hosts, though. I'm still getting a feel for using R.

Here is a snippet of data for one host. the user_count column is generated
from the users column using my custom function "usercount()". the samples
are taken roughly once per minute and only unique samples are recorded.
(i.e. use na.locf() to uncompress the data). Samples may occur twice in the
same minute and are rarely aligned on the same time.

Here is the original data before I turn t into a zoo series and run
na.locf() over it so I can aggregate a single host by day. I'm open to a
better way.

foo

  usersdatetime user_count
1 user1&  user2 2007-03-29 19:16:30  2
2 user1&  user2 2007-03-31 00:04:46  2
3 user1&  user2 2007-04-02 11:49:20  2
4 user1&  user2 2007-04-02 12:02:04  2
5 user1&  user2 2007-04-02 12:44:02  2
6 user1&  user2&  user3 2007-04-02 16:34:05  3


dput(foo)

structure(list(users = c("user1&  user2", "user1&  user2", "user1&  user2",
"user1&  user2", "user1&  user2", "user1&  user2&  user3"), datetime =
structure(c(1175210190,
1175313886, 1175528960, 1175529724, 1175532242, 1175546045), class =
c("POSIXt",
"POSIXct"), tzone = "US/Eastern"), user_count = c(2, 2, 2, 2,
2, 3)), .Names = c("users", "datetime", "user_count"), row.names = c(NA,
6L), class = "data.frame")


Thanks,
Jason

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






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

2011-01-11 Thread Bert Gunter
R does not use dummy variables. Models and contrasts are specified in
more natural, model formula based ways. See

?arima

 and/or CRAN's Time Series task view for numerous packages that fit time series.

-- Bert

On Tue, Jan 11, 2011 at 12:18 PM, Christofer Bogaso
 wrote:
> Dear all, I would like to ask one question related to statistics, for
> specifically on defining dummy variables. As of now, I have come across 3
> different kind of dummy variables (assuming I am working with Seasonal
> dummy, and number of season is 4):
>
>> dummy1 <- diag(4)
>> for(i in 1:3) dummy1 <- rbind(dummy1, diag(4))
>> dummy1 <- dummy1[,-4]
>>
>> dummy2 <- dummy1
>> dummy2[dummy2 == 0] = -1/(4-1)
>>
>> dummy3 <- dummy1 - 1/4
>>
>> head(dummy1)
>     [,1] [,2] [,3]
> [1,]    1    0    0
> [2,]    0    1    0
> [3,]    0    0    1
> [4,]    0    0    0
> [5,]    1    0    0
> [6,]    0    1    0
>> head(dummy2)
>           [,1]       [,2]       [,3]
> [1,]  1.000 -0.333 -0.333
> [2,] -0.333  1.000 -0.333
> [3,] -0.333 -0.333  1.000
> [4,] -0.333 -0.333 -0.333
> [5,]  1.000 -0.333 -0.333
> [6,] -0.333  1.000 -0.333
>> head(dummy3)
>      [,1]  [,2]  [,3]
> [1,]  0.75 -0.25 -0.25
> [2,] -0.25  0.75 -0.25
> [3,] -0.25 -0.25  0.75
> [4,] -0.25 -0.25 -0.25
> [5,]  0.75 -0.25 -0.25
> [6,] -0.25  0.75 -0.25
> Now I want to know which type of dummy definition is called Centered dummy
> and why it is called so? Is it equivalent to use any of the above
> definitions (atleast 2nd and 3rd?) It would really be very helpful if
> somebody point any suggestion and clarification.
>
> Thanks and regards,
>
>        [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Bert Gunter
Genentech Nonclinical Biostatistics
467-7374
http://devo.gene.com/groups/devo/depts/ncb/home.shtml

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


Re: [R] Assumptions for ANOVA: the right way to check the normality

2011-01-11 Thread Greg Snow
> From: Frodo Jedi [mailto:frodo.j...@yahoo.com] 
> Sent: Monday, January 10, 2011 5:44 PM
> To: Greg Snow
> Cc: r-help@r-project.org
> Subject: Re: [R] Assumptions for ANOVA: the right way to check the normality
> 
> Dear Greg,
> first of all thanks for your reply. And I add also many thanks to all of you 
> guys who are helping me, sorry for the amount of questions I recently posted 
> ;-) 
> 
> I don´t have a solid statistics background (I am not a statician) and I am 
> basically learning everything by myself. 
> So my first goal is TO UNDERSTAND. I need to have general guidelines because 
> for my PhD I am doing and I will do several psycophysic experiments.
> I am totally alone in this challenge, so I am asking some help to you guys as 
> I think that here is the best place to exchange the thing that I miss
> and that will never found in any book: the experience.

Isn't there a single statistician anywhere in the University?  Does your 
committee have any experience with any of this?

> >What is the question you are really trying to find the answer for?  Knowing 
> >that may help us give more meaningful answers.
> 
> Concerning your question I thought to have been clear. I want to understand 
> which analysis I have to use in order to understand if 
> the differences I am having are statistically significant or not. Now, as in 
> all the books I read there is written that to apply ANOVA 
> I must respect the assumption of normality then I am try to find a way to 
> understand this.

A general run of anova procedures will produce multiple p-values addressing 
multiple null hypotheses addressing many different questions (often many of 
which are uninteresting).  Which terms are you really trying to test and which 
are included because you already know that they have an effect.

Are you including interactions because you find them actually interesting? Or 
just because that is what everyone else does?

[snip]
 
> >Also remember that the normality of the data/residuals/etc. is not as 
> >important as the CLT for your sample size.  The main things that make the 
> >CLT not work (for samples that are >not large enough) are outliers and 
> >strong skewness, since your outcome is limited to the numbers 1-7, I don’t 
> >see outliers or skewness being a real problem.  So you are probably >fine 
> >for fixed effects style models (though checking with experts in your area or 
> >doing simulations can support/counter this).  
> 
> As far as I have seen everyone in my field does ANOVA.

[imagine best Mom voice] and if everyone in your field jumped off a cliff . . .

Do you want to do what everyone else is doing, or something new and different?

What does your committee chair say about this?

> >But when you add in random effects then there is a lot of uncertainty about 
> >if the normal theory still holds, the latest lme code uses mcmc sampling 
> >rather than depending on >normal theory and is still being developed.
> 
> 
> For "random effects" do you mean the repeated measures right? So why 
> staticians developed the ANOVA with repeated measure if there is so much 
> uncertainty?

Repeated measures are one type of random effect analysis, but random and mixed 
effects is more general than just repeated measures.

Statisticians developed those methods because they worked for simple cases, 
made some sense for more complicated cases, and they did not have anything that 
was both better and practical.  Now with modern computers we can see when those 
do work (unfortunately not as often as had been hoped) and what was once 
impractical is now much simpler (but inertia is to do it the old way, even 
though the people who developed the old way would have preferred to do it our 
way).  The article: 

Why Permutation Tests Are Superior to t and F Tests in Biomedical Research
John Ludbrook and Hugh Dudley
The American Statistician
Vol. 52, No. 2 (May, 1998), pp. 127-132

May be enlightening here (and give possible alternatives).

Also see: 
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q1/001819.html

for some simulation involving mixed models.  One shows that the normal theory 
works fine for that particular case, the next one shows a case where the normal 
theory does not work, then shows how to use simulation (parametric bootstrap) 
to get a more appropriate p-value.  You can adapt those examples for your own 
situation.

>  
> >This now comes back to my first question: what are you trying to find out?
> 
> My ultimate goal is to find the p-values in order to understand if my results 
> are significative or not. So I can write them on the paper ;-)

There is a function in the TeachingDemos package that will produce p-values if 
that is all your want, these are independent of any normality assumptions, 
independent of any data in fact.  However they don't really help with 
understanding.

Graphing the data (I think you have done this already) is the best route to 
understanding.  If you need more than that, the

Re: [R] A question on dummy variable

2011-01-11 Thread David Winsemius
You are not offering example of real codings but are rather showing  
something1 that you think looks like something2 (in R) that looks like  
something3 (in a textbook?). My guess is that the something2 might be  
contrast matrices or model matrices.


If you want a contrast matrix whose columns sum to zero (which is one  
possible situation that some people might call "centered" then look at  
the documentation for sum and poly contrasts. If you want to see a  
situation where model matrices are constructed which compare to the  
overall mean (another possible interpretation of "centered"), then  
look at the documentation for model.matrix and run the examples on  
that page with a -1 in the formula.


--
David.
On Jan 11, 2011, at 3:18 PM, Christofer Bogaso wrote:


Dear all, I would like to ask one question related to statistics, for
specifically on defining dummy variables. As of now, I have come  
across 3

different kind of dummy variables (assuming I am working with Seasonal
dummy, and number of season is 4):


dummy1 <- diag(4)
for(i in 1:3) dummy1 <- rbind(dummy1, diag(4))
dummy1 <- dummy1[,-4]

dummy2 <- dummy1
dummy2[dummy2 == 0] = -1/(4-1)

dummy3 <- dummy1 - 1/4

head(dummy1)

[,1] [,2] [,3]
[1,]100
[2,]010
[3,]001
[4,]000
[5,]100
[6,]010

head(dummy2)

  [,1]   [,2]   [,3]
[1,]  1.000 -0.333 -0.333
[2,] -0.333  1.000 -0.333
[3,] -0.333 -0.333  1.000
[4,] -0.333 -0.333 -0.333
[5,]  1.000 -0.333 -0.333
[6,] -0.333  1.000 -0.333

head(dummy3)

 [,1]  [,2]  [,3]
[1,]  0.75 -0.25 -0.25
[2,] -0.25  0.75 -0.25
[3,] -0.25 -0.25  0.75
[4,] -0.25 -0.25 -0.25
[5,]  0.75 -0.25 -0.25
[6,] -0.25  0.75 -0.25
Now I want to know which type of dummy definition is called Centered  
dummy

and why it is called so? Is it equivalent to use any of the above
definitions (atleast 2nd and 3rd?) It would really be very helpful if
somebody point any suggestion and clarification.

Thanks and regards,

[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] Assumptions for ANOVA: the right way to check the normality

2011-01-11 Thread Frodo Jedi
Many many thanks for your feedback Greg.
You have been very enlightening for me.

Now is time for me to study the material you kindly provided me. Thanks.








From: Greg Snow 

Cc: "r-help@r-project.org" 
Sent: Tue, January 11, 2011 10:13:34 PM
Subject: RE: [R] Assumptions for ANOVA: the right way to check the normality


> Sent: Monday, January 10, 2011 5:44 PM
> To: Greg Snow
> Cc: r-help@r-project.org
> Subject: Re: [R] Assumptions for ANOVA: the right way to check the normality
> 
> Dear Greg,
> first of all thanks for your reply. And I add also many thanks to all of you 
>guys who are helping me, sorry for the amount of questions I recently posted 
>;-) 
>
> 
> I don´t have a solid statistics background (I am not a statician) and I am 
>basically learning everything by myself. 
>
> So my first goal is TO UNDERSTAND. I need to have general guidelines because 
>for my PhD I am doing and I will do several psycophysic experiments.
> I am totally alone in this challenge, so I am asking some help to you guys as 
> I 
>think that here is the best place to exchange the thing that I miss
> and that will never found in any book: the experience.

Isn't there a single statistician anywhere in the University?  Does your 
committee have any experience with any of this?

> >What is the question you are really trying to find the answer for?  Knowing 
>that may help us give more meaningful answers.
> 
> Concerning your question I thought to have been clear. I want to understand 
>which analysis I have to use in order to understand if 
>
> the differences I am having are statistically significant or not. Now, as in 
>all the books I read there is written that to apply ANOVA 
>
> I must respect the assumption of normality then I am try to find a way to
>understand this.

A general run of anova procedures will produce multiple p-values addressing
multiple null hypotheses addressing many different questions (often many of
which are uninteresting).  Which terms are you really trying to test and which 
are included because you already know that they have an effect.

Are you including interactions because you find them actually interesting? Or 
just because that is what everyone else does?

[snip]

> >Also remember that the normality of the data/residuals/etc. is not as 
>important as the CLT for your sample size.  The main things that make the CLT 
>not work (for samples that are >not large enough) are outliers and strong
>skewness, since your outcome is limited to the numbers 1-7, I don’t see 
>outliers 
>or skewness being a real problem.  So you are probably >fine for fixed effects 
>style models (though checking with experts in your area or doing simulations 
>can 
>support/counter this).  
>
> 
> As far as I have seen everyone in my field does ANOVA.

[imagine best Mom voice] and if everyone in your field jumped off a cliff . . .

Do you want to do what everyone else is doing, or something new and different?

What does your committee chair say about this?

> >But when you add in random effects then there is a lot of uncertainty about 
> >if 
>the normal theory still holds, the latest lme code uses mcmc sampling rather 
>than depending on >normal theory and is still being developed.
> 
> 
> For "random effects" do you mean the repeated measures right? So why 
> staticians 
>developed the ANOVA with repeated measure if there is so much uncertainty?

Repeated measures are one type of random effect analysis, but random and mixed 
effects is more general than just repeated measures.

Statisticians developed those methods because they worked for simple cases, 
made 
some sense for more complicated cases, and they did not have anything that was 
both better and practical.  Now with modern computers we can see when those do 
work (unfortunately not as often as had been hoped) and what was once 
impractical is now much simpler (but inertia is to do it the old way, even
though the people who developed the old way would have preferred to do it our 
way).  The article: 


Why Permutation Tests Are Superior to t and F Tests in Biomedical Research
John Ludbrook and Hugh Dudley
The American Statistician
Vol. 52, No. 2 (May, 1998), pp. 127-132

May be enlightening here (and give possible alternatives).

Also see: 
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q1/001819.html 

for some simulation involving mixed models.  One shows that the normal theory 
works fine for that particular case, the next one shows a case where the normal 
theory does not work, then shows how to use simulation (parametric bootstrap) 
to 
get a more appropriate p-value.  You can adapt those examples for your own
situation.

>  
> >This now comes back to my first question: what are you trying to find out?
> 
> My ultimate goal is to find the p-values in order to understand if my results 
>are significative or not. So I can write them on the paper ;-)

There is a function in the TeachingDemos package that will produ

Re: [R] Alphabetic labels on multi-plot graphics

2011-01-11 Thread Eduardo de Oliveira Horta
Thanks,

I wasn't even aware that 'letters' was there.

Best regards

On Tue, Jan 11, 2011 at 6:08 PM, David Winsemius  wrote:
>
> On Jan 11, 2011, at 2:06 PM, Eduardo de Oliveira Horta wrote:
>
>> Is there a way to achieve
>>
>> lbl=c("a", "b", "c", "d")
>>
>> opar <- par(mfrow=c(2,2), ann=FALSE)
>> for (t in 1:4){
>>  plot(seq(from=1,to=2*pi,length=100),
>> sin(t*seq(from=1,to=2*pi,length=100)), type="l")
>>  title(main=paste("(", lbl[t], ")", sep=""))
>> }
>> par(opar)
>
> opar <- par(mfrow=c(2,2), ann=FALSE)
>  for (t in 1:4){
>  plot(seq(from=1,to=2*pi,length=100),
>  sin(t*seq(from=1,to=2*pi,length=100)), type="l")
>  title(main=bquote("("*.(letters[t])*")") )
>  }
> par(opar)
>
>> without having to use an object like 'lbl'?
>
> (You obviously need something from which to construct the titles.)
>
>>
>> More generally: is it possible to iteratively (as in a loop) add
>> alphabetic titles to multi-plot graphics when the range over which 't'
>> above varies is of an arbitrary length? It is important that
>> ann=FALSE, because I don't want the axes labels.
>>
>> Thanks in advance, and best regards.
>>
>> Eduardo Horta
>>
>> __
>
> David Winsemius, MD
> West Hartford, CT
>
>

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

2011-01-11 Thread Dickison, Daniel
The documentation for `aggregate` makes it sound like aggregate.formula should 
behave identically to aggregate.data.frame (apart from the way the parameters 
are passed).  But it looks like aggregate.formula is quietly removing rows 
where any of the "output" variables (those on the LHS of the formula) are NA.  
This differs from how aggregate.data.frame works.  Is this expected behavior?

Here are a couple of examples:

> d <- data.frame(a=rep(1:2, each=2),
+ b=c(1,2,NA,3))
> aggregate(d["b"], d["a"], mean)
  a   b
1 1 1.5
2 2  NA
> aggregate(b ~ a, d, mean)
  a   b
1 1 1.5
2 2 3.0

It's removing whole rows even if just one of the columns is NA, i.e.:

> d <- data.frame(a=rep(1:2, each=2),
+ b=c(1,2,NA,3),
+ c=c(NA,2,3,NA))
> aggregate(cbind(b,c) ~ a, d, mean)
  a b c
1 1 2 2

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


Re: [R] Strange behavior when trying to piggyback off of "fitdistr"

2011-01-11 Thread emorway

Avi, 

Its been nearly a year since you made this post, but I'm curious if you were
able to find a solution to your problem?  Your inquiry is closely related to
a problem I'm having.

Thanks, 
Eric
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Strange-behavior-when-trying-to-piggyback-off-of-fitdistr-tp1012536p3209887.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] gennerating skewed random numbers

2011-01-11 Thread Michael Bedward
Hi Ed,

The sn package can generate random values from skewed Normal and t
distributions.

Also see here: http://azzalini.stat.unipd.it/SN/faq-r.html

Michael

On 11 January 2011 23:37, Ed Keith  wrote:
> This is not exactly an R specific question, but I think the people on this 
> list can probably help.
>
> I'm working on a simulation. In the model I have the first three moments of 
> the distributions of the variables. I know how to generate a random number 
> from a distribution given the first two moments assuming the third moment is 
> 0. But I do not know how to generate a number drawn from a distribution with 
> a nonzero third monument.
>
> If someone could point me to a good reference I would appreciate it.
>
> Thank you in advance,
>
>   -EdK
>
> Ed Keith
> e_...@yahoo.com
>
> Blog: edkeith.blogspot.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

2011-01-11 Thread David Winsemius


On Jan 11, 2011, at 5:41 PM, Dickison, Daniel wrote:

The documentation for `aggregate` makes it sound like  
aggregate.formula should behave identically to aggregate.data.frame  
(apart from the way the parameters are passed).  But it looks like  
aggregate.formula is quietly removing rows where any of the "output"  
variables (those on the LHS of the formula) are NA.  This differs  
from how aggregate.data.frame works.  Is this expected behavior?


Here are a couple of examples:


d <- data.frame(a=rep(1:2, each=2),

+ b=c(1,2,NA,3))

aggregate(d["b"], d["a"], mean)

 a   b
1 1 1.5
2 2  NA

aggregate(b ~ a, d, mean)

 a   b
1 1 1.5
2 2 3.0

It's removing whole rows even if just one of the columns is NA, i.e.:


d <- data.frame(a=rep(1:2, each=2),

+ b=c(1,2,NA,3),
+ c=c(NA,2,3,NA))

aggregate(cbind(b,c) ~ a, d, mean)

 a b c
1 1 2 2



The help page for aggregate gives the calling defaults for  
aggregate.formula as:
## S3 method for class 'formula' aggregate(formula, data, FUN, ...,  
subset, na.action = na.omit)
So the description you give seems to be adhering to what I would have  
expected (had I initially read the help page.)

--
David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Problems creating a PNG file for a dendrogram: "Error in plot.window(...) : need finite 'xlim' values"

2011-01-11 Thread Richard Vlasimsky

Has anyone successfully created a PNG file for a dendrogram?

I am able to successfully launch and view a dendrogram in Quartz.  However, the 
dendrogram is quite large (too large to read on a computer screen), so I am 
trying to save it to a file (1000x4000 pixels) for viewing in other apps.  
However, whenever I try to initiate a PNG device, I get a "need finitite 'xlim' 
values" error. 



Here is some example code to illustrate my point:

cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#This works!  Plot is generated in quartz no problem.


#Now, try this:
png(filename="delme.png",width=4000,height=1000);
cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#Error in plot.window(...) : need finite 'xlim' values
#In addition: Warning messages:
#1: In min(x) : no non-missing arguments to min; returning Inf
#2: In max(x) : no non-missing arguments to max; returning -Inf
#3: In min(x) : no non-missing arguments to min; returning Inf
#4: In max(x) : no non-missing arguments to max; returning -Inf 

This is the exact same code, only a prior call to png() causes the seemingly 
unrelated xlim to fail.  Why is this?

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


Re: [R] Problems creating a PNG file for a dendrogram: "Error in plot.window(...) : need finite 'xlim' values"

2011-01-11 Thread Peter Langfelder
I'm guessing that your code actually generates two plots, one with the
command p <- plot(hc) and one with plot(p), which doesn't work for a
png. Try getting rid of the plot(p).

Peter

On Tue, Jan 11, 2011 at 4:01 PM, Richard Vlasimsky
 wrote:
>
> Has anyone successfully created a PNG file for a dendrogram?
>
> I am able to successfully launch and view a dendrogram in Quartz.  However, 
> the dendrogram is quite large (too large to read on a computer screen), so I 
> am trying to save it to a file (1000x4000 pixels) for viewing in other apps.  
> However, whenever I try to initiate a PNG device, I get a "need finitite 
> 'xlim' values" error.
>
>
>
> Here is some example code to illustrate my point:
>
> cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
> distance <- as.dist(1.0-cor.matrix);
> hc <- hclust(distance);
> p <- plot(hc);
> plot(p);
> #This works!  Plot is generated in quartz no problem.
>
>
> #Now, try this:
> png(filename="delme.png",width=4000,height=1000);
> cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
> distance <- as.dist(1.0-cor.matrix);
> hc <- hclust(distance);
> p <- plot(hc);
> plot(p);
> #Error in plot.window(...) : need finite 'xlim' values
> #In addition: Warning messages:
> #1: In min(x) : no non-missing arguments to min; returning Inf
> #2: In max(x) : no non-missing arguments to max; returning -Inf
> #3: In min(x) : no non-missing arguments to min; returning Inf
> #4: In max(x) : no non-missing arguments to max; returning -Inf
>
> This is the exact same code, only a prior call to png() causes the seemingly 
> unrelated xlim to fail.  Why is this?
>
> Thanks,
> Richard Vlasimsky
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] A question on dummy variable

2011-01-11 Thread Gabor Grothendieck
On Tue, Jan 11, 2011 at 3:18 PM, Christofer Bogaso
 wrote:
> Dear all, I would like to ask one question related to statistics, for
> specifically on defining dummy variables. As of now, I have come across 3
> different kind of dummy variables (assuming I am working with Seasonal
> dummy, and number of season is 4):
>
>> dummy1 <- diag(4)
>> for(i in 1:3) dummy1 <- rbind(dummy1, diag(4))
>> dummy1 <- dummy1[,-4]
>>
>> dummy2 <- dummy1
>> dummy2[dummy2 == 0] = -1/(4-1)
>>
>> dummy3 <- dummy1 - 1/4
>>
>> head(dummy1)
>     [,1] [,2] [,3]
> [1,]    1    0    0
> [2,]    0    1    0
> [3,]    0    0    1
> [4,]    0    0    0
> [5,]    1    0    0
> [6,]    0    1    0
>> head(dummy2)
>           [,1]       [,2]       [,3]
> [1,]  1.000 -0.333 -0.333
> [2,] -0.333  1.000 -0.333
> [3,] -0.333 -0.333  1.000
> [4,] -0.333 -0.333 -0.333
> [5,]  1.000 -0.333 -0.333
> [6,] -0.333  1.000 -0.333
>> head(dummy3)
>      [,1]  [,2]  [,3]
> [1,]  0.75 -0.25 -0.25
> [2,] -0.25  0.75 -0.25
> [3,] -0.25 -0.25  0.75
> [4,] -0.25 -0.25 -0.25
> [5,]  0.75 -0.25 -0.25
> [6,] -0.25  0.75 -0.25
> Now I want to know which type of dummy definition is called Centered dummy
> and why it is called so? Is it equivalent to use any of the above
> definitions (atleast 2nd and 3rd?) It would really be very helpful if
> somebody point any suggestion and clarification.
>

The contrasts of your dummy1 matrix are contr.SAS contrasts in R.
(The default contrasts in R are contr.treatment which are the same as
contr.SAS except contr.SAS uses the last level as the base whereas
treatment contrasts use the first level as the base.)

   options(contrasts = c("contr.SAS", "contr.poly"))
   f <- gl(4, 1, 16)
   M <- model.matrix( ~ f )
   all( M[, -1] == dummy1) # TRUE

Centered contrasts are ones which have been centered -- i.e. the mean
of each column has been subtracted from that column.  This is
equivalent to saying that the column sums are zero.

The means of the three columns of dummy1 are c(1/4, 1/4, 1/4) so if we
subtract 1/4 from dummy1 we get a centered contrasts matrix. That is
precisely what you did to get dummy3.  We can check that dummy3 is
centered:

   colSums(dummy3) # 0 0 0

dummy2 is just a scaled version of dummy3.  In fact dummy2 equals
dummy3 / .75 so its not fundamentally different.  Its columns still
sum to zero so its still centered.

   all( dummy2 == dummy3 / .75) # TRUE
   colSums(dummy2) # 0 0 0 except for floating point error

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

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


Re: [R] Problems creating a PNG file for a dendrogram: "Error in plot.window(...) : need finite 'xlim' values"

2011-01-11 Thread Bill.Venables
I very much doubt your first example does work.  the value of plot() is NULL 
which if you plot again will give the error message you see in your second 
example.

What where you trying to achieve doing

p <- plot(hc)
plot(p) ### this one is trying to plot NULL

?

Here is an example (such as you might have given, according to the posting 
guide):

> x <- matrix(rnorm(500*5), 500, 5)
> dc <- as.dist(1-cor(x))
> hc <- hclust(dc)
> p <- plot(hc)
> plot(p)
Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf
> 

Look familiar? This is why:

> p
NULL
> 


-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Richard Vlasimsky
Sent: Wednesday, 12 January 2011 10:01 AM
To: r-help@r-project.org
Subject: [R] Problems creating a PNG file for a dendrogram: "Error in 
plot.window(...) : need finite 'xlim' values"


Has anyone successfully created a PNG file for a dendrogram?

I am able to successfully launch and view a dendrogram in Quartz.  However, the 
dendrogram is quite large (too large to read on a computer screen), so I am 
trying to save it to a file (1000x4000 pixels) for viewing in other apps.  
However, whenever I try to initiate a PNG device, I get a "need finitite 'xlim' 
values" error. 



Here is some example code to illustrate my point:

cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#This works!  Plot is generated in quartz no problem.


#Now, try this:
png(filename="delme.png",width=4000,height=1000);
cor.matrix <- cor(mydata,method="pearson",use="pairwise.complete.obs");
distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#Error in plot.window(...) : need finite 'xlim' values
#In addition: Warning messages:
#1: In min(x) : no non-missing arguments to min; returning Inf
#2: In max(x) : no non-missing arguments to max; returning -Inf
#3: In min(x) : no non-missing arguments to min; returning Inf
#4: In max(x) : no non-missing arguments to max; returning -Inf 

This is the exact same code, only a prior call to png() causes the seemingly 
unrelated xlim to fail.  Why is this?

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

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


Re: [R] aggregate.formula implicitly removes rows containing NA

2011-01-11 Thread Peter Ehlers

On 2011-01-11 14:41, Dickison, Daniel wrote:

The documentation for `aggregate` makes it sound like aggregate.formula should behave 
identically to aggregate.data.frame (apart from the way the parameters are passed).  But 
it looks like aggregate.formula is quietly removing rows where any of the 
"output" variables (those on the LHS of the formula) are NA.  This differs from 
how aggregate.data.frame works.  Is this expected behavior?

Here are a couple of examples:


d<- data.frame(a=rep(1:2, each=2),

+ b=c(1,2,NA,3))

aggregate(d["b"], d["a"], mean)

   a   b
1 1 1.5
2 2  NA

aggregate(b ~ a, d, mean)

   a   b
1 1 1.5
2 2 3.0

It's removing whole rows even if just one of the columns is NA, i.e.:


d<- data.frame(a=rep(1:2, each=2),

+ b=c(1,2,NA,3),
+ c=c(NA,2,3,NA))

aggregate(cbind(b,c) ~ a, d, mean)

   a b c
1 1 2 2

Daniel


Try setting na.acton = na.pass.

Peter Ehlers

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

2011-01-11 Thread Dickison, Daniel
Oh wow, that would be it. Not sure how I missed that. Thanks for the tip.

Sent from my iPhone

On Jan 11, 2011, at 18:56, "David Winsemius"  wrote:

> 
> On Jan 11, 2011, at 5:41 PM, Dickison, Daniel wrote:
> 
>> The documentation for `aggregate` makes it sound like  
>> aggregate.formula should behave identically to aggregate.data.frame  
>> (apart from the way the parameters are passed).  But it looks like  
>> aggregate.formula is quietly removing rows where any of the "output"  
>> variables (those on the LHS of the formula) are NA.  This differs  
>> from how aggregate.data.frame works.  Is this expected behavior?
>> 
>> Here are a couple of examples:
>> 
>>> d <- data.frame(a=rep(1:2, each=2),
>> + b=c(1,2,NA,3))
>>> aggregate(d["b"], d["a"], mean)
>> a   b
>> 1 1 1.5
>> 2 2  NA
>>> aggregate(b ~ a, d, mean)
>> a   b
>> 1 1 1.5
>> 2 2 3.0
>> 
>> It's removing whole rows even if just one of the columns is NA, i.e.:
>> 
>>> d <- data.frame(a=rep(1:2, each=2),
>> + b=c(1,2,NA,3),
>> + c=c(NA,2,3,NA))
>>> aggregate(cbind(b,c) ~ a, d, mean)
>> a b c
>> 1 1 2 2
>> 
> 
> The help page for aggregate gives the calling defaults for  
> aggregate.formula as:
> ## S3 method for class 'formula' aggregate(formula, data, FUN, ...,  
> subset, na.action = na.omit)
> So the description you give seems to be adhering to what I would have  
> expected (had I initially read the help page.)
> -- 
> David Winsemius, MD
> West Hartford, CT
> 

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


Re: [R] Problems creating a PNG file for a dendrogram: "Error in plot.window(...) : need finite 'xlim' values"

2011-01-11 Thread David Winsemius


On Jan 11, 2011, at 7:01 PM, Richard Vlasimsky wrote:



Has anyone successfully created a PNG file for a dendrogram?

I am able to successfully launch and view a dendrogram in Quartz.   
However, the dendrogram is quite large (too large to read on a  
computer screen), so I am trying to save it to a file (1000x4000  
pixels) for viewing in other apps.  However, whenever I try to  
initiate a PNG device, I get a "need finitite 'xlim' values" error.




Here is some example code to illustrate my point:

cor.matrix <-  
cor(mydata,method="pearson",use="pairwise.complete.obs");

distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#This works!  Plot is generated in quartz no problem.


#Now, try this:
png(filename="delme.png",width=4000,height=1000);
cor.matrix <-  
cor(mydata,method="pearson",use="pairwise.complete.obs");

distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#Error in plot.window(...) : need finite 'xlim' values
#In addition: Warning messages:
#1: In min(x) : no non-missing arguments to min; returning Inf
#2: In max(x) : no non-missing arguments to max; returning -Inf
#3: In min(x) : no non-missing arguments to min; returning Inf
#4: In max(x) : no non-missing arguments to max; returning -Inf


I'm not sure the other two answers address the problems I found. When  
I try to set up a png file with the parameters width=4000,height=1000,  
on a Mac I intially got no plot with what is an otherwise valid  
command. But after successfully getting plotting to a png device the  
logjam appear broken. Try:


 graphics.off()
 dev.list()
#NULL
 png(filename="delme.png",width=4000,height=1000);
 plot(hc)
 dev.off()


(Of course I used dev.off() which you did not, but even adding  
dev.off() was not enough to get success, at least initially. I don't  
understand the suggestion to get rid of plot(hc) or the suggestion  
that hclust() returns NULL. That's certainly not how I read the help  
page and examples for hclust.)


This is the exact same code, only a prior call to png() causes the  
seemingly unrelated xlim to fail.  Why is this?


Thanks,
Richard Vlasimsky



David Winsemius, MD
West Hartford, CT

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

2011-01-11 Thread Wincent
Dear r heper,

How can I change the strip text, for example (16,23] in the following
example, to other more informative text such as "high level" on the
fly?

library(effects)
Cowles$ex2 <- cut(Cowles$extraversion,3)
mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, family=binomial)
eff.cowles <- allEffects(mod.cowles)
plot(eff.cowles, 'neuroticism:ex2',factor.names=F)

Thank you.

Ronggui


-- 
Wincent Ronggui HUANG (Ph.D.)
City University of Hong Kong
http://asrr.r-forge.r-project.org/rghuang.html

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


Re: [R] Problems creating a PNG file for a dendrogram: "Error in plot.window(...) : need finite 'xlim' values"

2011-01-11 Thread David Winsemius


On Jan 11, 2011, at 9:27 PM, David Winsemius wrote:



On Jan 11, 2011, at 7:01 PM, Richard Vlasimsky wrote:



Has anyone successfully created a PNG file for a dendrogram?

I am able to successfully launch and view a dendrogram in Quartz.   
However, the dendrogram is quite large (too large to read on a  
computer screen), so I am trying to save it to a file (1000x4000  
pixels) for viewing in other apps.  However, whenever I try to  
initiate a PNG device, I get a "need finitite 'xlim' values" error.




Here is some example code to illustrate my point:

cor.matrix <-  
cor(mydata,method="pearson",use="pairwise.complete.obs");

distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#This works!  Plot is generated in quartz no problem.


#Now, try this:
png(filename="delme.png",width=4000,height=1000);
cor.matrix <-  
cor(mydata,method="pearson",use="pairwise.complete.obs");

distance <- as.dist(1.0-cor.matrix);
hc <- hclust(distance);
p <- plot(hc);
plot(p);
#Error in plot.window(...) : need finite 'xlim' values
#In addition: Warning messages:
#1: In min(x) : no non-missing arguments to min; returning Inf
#2: In max(x) : no non-missing arguments to max; returning -Inf
#3: In min(x) : no non-missing arguments to min; returning Inf
#4: In max(x) : no non-missing arguments to max; returning -Inf


I'm not sure the other two answers address the problems I found.  
When I try to set up a png file with the parameters  
width=4000,height=1000, on a Mac I intially got no plot with what is  
an otherwise valid command. But after successfully getting plotting  
to a png device the logjam appear broken. Try:


graphics.off()
dev.list()
#NULL
png(filename="delme.png",width=4000,height=1000);
plot(hc)
dev.off()


(Of course I used dev.off() which you did not, but even adding  
dev.off() was not enough to get success, at least initially. I don't  
understand the suggestion to get rid of plot(hc) or the suggestion  
that hclust() returns NULL. That's certainly not how I read the help  
page and examples for hclust.)


I guess it's true to say that I misunderstood when Venables and  
Langfelder _didn't_ say either of those things.  They said that  
plot(plot(hc)) was not needed and I have now seen the light. The  
missing ingredient in the OP's frustrations is still dev.off().




This is the exact same code, only a prior call to png() causes the  
seemingly unrelated xlim to fail.  Why is this?


Thanks,
Richard Vlasimsky



David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT

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

2011-01-11 Thread John Sorkin
Christofer,
I am not sure I understand how you are using your dummy variables. Generally if 
you have n categories you need n-1 dummy variables. Thus if you have three 
categories, low, medium, high and want to compare two of the levels to a 
reference level (a coding scheme sometimes called reference cell coding) you 
could use the following coding which medium and high to the reference level, 
low:

level   dummy1 dummy2
low0 0
medium 0 1
high   1 0

You will notice that for three categories, my dummy variables from  an 3 by 2 
matrix. In general the dummy variable matrix for n categories will be an n by 
n-1 matrix. You say your have four seasons. I would expect your dummy variable 
matrix to be of size 4 by 3. Your matrices are 6 by 3. Am I not understanding 
what you are trying to do?
John

John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
Baltimore VA Medical Center GRECC,
University of Maryland School of Medicine Claude D. Pepper OAIC,
University of Maryland Clinical Nutrition Research Unit, and
Baltimore VA Center Stroke of Excellence

University of Maryland School of Medicine
Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524

(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
jsor...@grecc.umaryland.edu

>>> Christofer Bogaso  1/11/2011 3:18 PM >>>
Dear all, I would like to ask one question related to statistics, for
specifically on defining dummy variables. As of now, I have come across 3
different kind of dummy variables (assuming I am working with Seasonal
dummy, and number of season is 4):

> dummy1 <- diag(4)
> for(i in 1:3) dummy1 <- rbind(dummy1, diag(4))
> dummy1 <- dummy1[,-4]
>
> dummy2 <- dummy1
> dummy2[dummy2 == 0] = -1/(4-1)
>
> dummy3 <- dummy1 - 1/4
>
> head(dummy1)
 [,1] [,2] [,3]
[1,]100
[2,]010
[3,]001
[4,]000
[5,]100
[6,]010
> head(dummy2)
   [,1]   [,2]   [,3]
[1,]  1.000 -0.333 -0.333
[2,] -0.333  1.000 -0.333
[3,] -0.333 -0.333  1.000
[4,] -0.333 -0.333 -0.333
[5,]  1.000 -0.333 -0.333
[6,] -0.333  1.000 -0.333
> head(dummy3)
  [,1]  [,2]  [,3]
[1,]  0.75 -0.25 -0.25
[2,] -0.25  0.75 -0.25
[3,] -0.25 -0.25  0.75
[4,] -0.25 -0.25 -0.25
[5,]  0.75 -0.25 -0.25
[6,] -0.25  0.75 -0.25
Now I want to know which type of dummy definition is called Centered dummy
and why it is called so? Is it equivalent to use any of the above
definitions (atleast 2nd and 3rd?) It would really be very helpful if
somebody point any suggestion and clarification.

Thanks and regards,

[[alternative HTML version deleted]]

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

Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}

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


[R] GLMM with lme4 and octopus behaviour

2011-01-11 Thread Samaritan

Hi all,

First time poster and a relatively new R user, I'm beginning analysis for my
masters degree. I'm doing a bit of work on octopus behaviour, and while it's
been fascinating, the stats behind it are a bit beyond my grasp at the
moment. I was hoping that somebody with more experience my be able to look
at my example and offer their wisdom, much to my appreciation :-)

At the most basic level, I'm testing the effect of sleep deprivation on
various behaviours (e.g. amount of time spent awake, amount of time spend
expressing difference textures, patterns, colours etc). I take a video
sample of each octopus every hour for 72 hours and score their behaviour in
Jwatcher.

I'm using GLMM so that I can nest Individual as a random factor within Time,
which I'm told will reduce the problem of making repeated measures (is this
effectively blocking by time and by octopus?). Currently my model looks like
this:

octopus.lmer<-lmer(awake~as.factor(Treatment)+Sex+Weight+(1|Time/Octopus))

where "Treatment" is "0" or "1" representing sleep-deprived or
sleep-allowed. When I try to fit the model I get the following error
messages:

Error: length(f1) == length(f2) is not TRUE
In addition: Warning messages:
1: In Octopus:Time :
  numerical expression has 190 elements: only the first used
2: In Octopus:Time :
  numerical expression has 190 elements: only the first used

This is the point at which I become lost. What does this mean? Clearly I'm
not doing something right, perhaps in my data preparation? So far as I can
see the length of each of the variables is the same (although I'm no certain
as to what f1 and f2 refer to).

If anyone could offer me some kind of advice about this I would appreciate
it very, VERY much. Both of my supervisors have no experience with R and so
have kind of washed their hands, so I'm alone in this and your expertise
would be a big help.

Dean


-- 
View this message in context: 
http://r.789695.n4.nabble.com/GLMM-with-lme4-and-octopus-behaviour-tp3212716p3212716.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] how to change strip text of effect plot

2011-01-11 Thread Joshua Wiley
Hi,

I am guessing this is not what you meant by "on the fly", but I think
it will be by far the easiest way.  Plotting an effects object is a
high level plot with a lot of defaults and automation built in to make
your life simple.  The cost is that it is less flexible---you work its
way, not vice versa.  If you want the factor named high, just label it
that way to begin with.  If you think it makes the graphs more
interpretable/meaningful, then it will make model summaries, etc.
better also.  Worst case, you fit the model twice (one with fancy
names, one with numbers), which unless you have a massive dataset will
not take long or be an onerous burden anyways.  Here's how you can
include labels directly in cut():

Cowles$ex2 <- cut(Cowles$extraversion, 3, c("low", "medium", "high"))
mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, family=binomial)
eff.cowles <- allEffects(mod.cowles)
plot(eff.cowles, 'neuroticism:ex2',factor.names=F)

Cheers,

Josh

On Tue, Jan 11, 2011 at 7:21 PM, Wincent  wrote:
> Dear r heper,
>
> How can I change the strip text, for example (16,23] in the following
> example, to other more informative text such as "high level" on the
> fly?
>
> library(effects)
> Cowles$ex2 <- cut(Cowles$extraversion,3)
> mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, 
> family=binomial)
> eff.cowles <- allEffects(mod.cowles)
> plot(eff.cowles, 'neuroticism:ex2',factor.names=F)
>
> Thank you.
>
> Ronggui
>
>
> --
> Wincent Ronggui HUANG (Ph.D.)
> City University of Hong Kong
> http://asrr.r-forge.r-project.org/rghuang.html
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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


Re: [R] how to change strip text of effect plot

2011-01-11 Thread Wincent
Sure that is one way to go.
Since it is possible to pass lattice arguments to that high level
function, and there should be ways to relabel/custom the strip text in
the lattice plotting system, I would think such an indirect method a
last resort.

Thank you.

Ronggui

On 12 January 2011 11:51, Joshua Wiley  wrote:
> Hi,
>
> I am guessing this is not what you meant by "on the fly", but I think
> it will be by far the easiest way.  Plotting an effects object is a
> high level plot with a lot of defaults and automation built in to make
> your life simple.  The cost is that it is less flexible---you work its
> way, not vice versa.  If you want the factor named high, just label it
> that way to begin with.  If you think it makes the graphs more
> interpretable/meaningful, then it will make model summaries, etc.
> better also.  Worst case, you fit the model twice (one with fancy
> names, one with numbers), which unless you have a massive dataset will
> not take long or be an onerous burden anyways.  Here's how you can
> include labels directly in cut():
>
> Cowles$ex2 <- cut(Cowles$extraversion, 3, c("low", "medium", "high"))
> mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, 
> family=binomial)
> eff.cowles <- allEffects(mod.cowles)
> plot(eff.cowles, 'neuroticism:ex2',factor.names=F)
>
> Cheers,
>
> Josh
>
> On Tue, Jan 11, 2011 at 7:21 PM, Wincent  wrote:
>> Dear r heper,
>>
>> How can I change the strip text, for example (16,23] in the following
>> example, to other more informative text such as "high level" on the
>> fly?
>>
>> library(effects)
>> Cowles$ex2 <- cut(Cowles$extraversion,3)
>> mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, 
>> family=binomial)
>> eff.cowles <- allEffects(mod.cowles)
>> plot(eff.cowles, 'neuroticism:ex2',factor.names=F)
>>
>> Thank you.
>>
>> Ronggui
>>
>>
>> --
>> Wincent Ronggui HUANG (Ph.D.)
>> City University of Hong Kong
>> http://asrr.r-forge.r-project.org/rghuang.html
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
>
> --
> Joshua Wiley
> Ph.D. Student, Health Psychology
> University of California, Los Angeles
> http://www.joshuawiley.com/
>



-- 
Wincent Ronggui HUANG (Ph.D.)
City University of Hong Kong
http://asrr.r-forge.r-project.org/rghuang.html

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


Re: [R] Aggragating subsets of data in larger vector with sapply

2011-01-11 Thread Joshua Ulrich
Hi Chris,

This seems to work on the sample data you provided.

FUN <- function(x) {
  x <- xts(as.numeric(x),index(x))
  period.apply(x, endpoints(x,"secs"), sum)
}
lapply(split.default(xSym$Size,xSym$Direction), FUN)

Best,
--
Joshua Ulrich  |  FOSS Trading: www.fosstrading.com



On Sun, Jan 9, 2011 at 6:10 PM, rivercode  wrote:
>
>
> Have 40,000 rows of buy/sell trade data and am trying to add up the buys for
> each second, the code works but it is very slow.  Any suggestions how to
> improve the sapply function ?
>
> secEP = endpoints(xSym$Direction, "secs")  # vector of last second on an XTS
> timeseries object with multiple entries for each second.
> d = xSym$Direction
> s = xSym$Size
> buySize = sapply(1:(length(secEP)-1), function(y) {
>        i =  (secEP[y]+ 1):secEP[y+1]; # index of vectors between each secEP
>        return(sum(as.numeric(s[i][d[i] == "buy"])));
> } )
>
> Object details:
>
> secEP = numeric Vector of one second Endpoints in xSym$Direction.
>
>> head(xSym$Direction)
>                    Direction
> 2011-01-05 09:30:00 "unkn"
> 2011-01-05 09:30:02 "sell"
> 2011-01-05 09:30:02 "buy"
> 2011-01-05 09:30:04 "buy"
> 2011-01-05 09:30:04 "buy"
> 2011-01-05 09:30:04 "buy"
>
>> head(xSym$Size)
>                    Size
> 2011-01-05 09:30:00 " 865"
> 2011-01-05 09:30:02 " 100"
> 2011-01-05 09:30:02 " 100"
> 2011-01-05 09:30:04 " 100"
> 2011-01-05 09:30:04 " 100"
> 2011-01-05 09:30:04 "  41"
>
> Thanks,
> Chris
>
>
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Aggragating-subsets-of-data-in-larger-vector-with-sapply-tp3206445p3206445.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] how to change strip text of effect plot

2011-01-11 Thread Joshua Wiley
Hmm, I felt like it was the clearest, most direct way.  Then again,
things of this nature (overide defaults using arguments rather than
changing what you feed the functions) do seem to be a common request
both for lattice and ggplot2.  In any case, my memory of the lattice
book and a quick search of the archives suggest that the typical way
to do this when dealing with lattice functions such as xyplot() would
be to define a custom strip function.  For example:

strip = function(..., factor.levels) strip.default(..., factor.levels
= c("low", "medium", "high"))

however, if I am not mistaken, you are plotting an object of class
"efflist", which dispatches plot.efflist which eventually dispatches
plot.eff, in the bowels of plot.eff is the call to xyplot(), with its
own strip function already defined.  I did not see a way to override
this, nor did I see an option to pass an argument to it.

It is trivial to edit the function's code to work with such a change.
You could probably create a local copy of it, make the necessary
changes, and then just ensure that your object gets dispatched to your
revised function rather than the packages original.  Even easier, I
suppose:

debug(plot)
plot(eff.cowles, 'neuroticism:ex2',factor.names=F)
## follow along until it has just created the object "plot"
## then overwrite that with your desired changes (factor.levels)
## before it is printed to the screen

plot <- xyplot(eval(parse(text = paste("fit ~", predictors[x.var],
"|", paste(predictors[-x.var], collapse = "*", strip =
function(..., factor.levels) strip.default(..., factor.levels =
c("low", "medium", "high"),
strip.names = c(factor.names, TRUE)), panel = function(x,
y, subscripts, x.vals, rug, lower, upper, has.se, ...) {
llines(x, y, lwd = 2, col = colors[1], ...)
if (rug)
lrug(x.vals)
if (has.se) {
llines(x, lower[subscripts], lty = 2, col = colors[2])
llines(x, upper[subscripts], lty = 2, col = colors[2])
}
if (has.thresholds) {
panel.abline(h = thresholds, lty = 3)
panel.text(rep(current.panel.limits()$xlim[1], length(thresholds)),
thresholds, threshold.labels, adj = c(0, 0),
cex = 0.75)
panel.text(rep(current.panel.limits()$xlim[2], length(thresholds)),
thresholds, threshold.labels, adj = c(1, 0),
cex = 0.75)
}
}, ylim = ylim, ylab = ylab, xlab = if (missing(xlab))
predictors[x.var]
else xlab, x.vals = x.vals, rug = rug, main = main, lower = x$lower,
upper = x$upper, has.se = has.se, data = x, scales = list(y =
list(at = tickmarks$at,
labels = tickmarks$labels), alternating = alternating),
...)

alternately still, you could rename the factor levels in the object
"eff.cowles", which is sort of inbetween changing your data and
changing the strip label defaults.

I can understand why it seems like there should be a simpler solution,
but I honestly do not see one.  The factor.levels argument does not
even work inside xyplot(), it needs to be in the strip function, and
that is nested far away from plot(effobject).  I do not see any
documentation that suggests a way, nor do I see any formal arguments
to facilitate it.  Then again, I'm not an expert in lattice or the
effects package.

Best regards,

Josh

On Tue, Jan 11, 2011 at 8:06 PM, Wincent  wrote:
> Sure that is one way to go.
> Since it is possible to pass lattice arguments to that high level
> function, and there should be ways to relabel/custom the strip text in
> the lattice plotting system, I would think such an indirect method a
> last resort.
>
> Thank you.
>
> Ronggui
>
> On 12 January 2011 11:51, Joshua Wiley  wrote:
>> Hi,
>>
>> I am guessing this is not what you meant by "on the fly", but I think
>> it will be by far the easiest way.  Plotting an effects object is a
>> high level plot with a lot of defaults and automation built in to make
>> your life simple.  The cost is that it is less flexible---you work its
>> way, not vice versa.  If you want the factor named high, just label it
>> that way to begin with.  If you think it makes the graphs more
>> interpretable/meaningful, then it will make model summaries, etc.
>> better also.  Worst case, you fit the model twice (one with fancy
>> names, one with numbers), which unless you have a massive dataset will
>> not take long or be an onerous burden anyways.  Here's how you can
>> include labels directly in cut():
>>
>> Cowles$ex2 <- cut(Cowles$extraversion, 3, c("low", "medium", "high"))
>> mod.cowles <- glm(volunteer ~ sex+neuroticism*ex2,data=Cowles, 
>> family=binomial)
>> eff.cowles <- allEffects(mod.cowles)
>> plot(eff.cowles, 'neuroticism:ex2',factor.names=F)
>>
>> Cheers,
>>
>> Josh
>>
>> On Tue, Jan 11, 2011 at 7:21 PM, Wincent  wrote:
>>> Dear r heper,
>>>
>>> How can I change the strip text, for example (16,23] in the following
>>> example,

[R] plot: skip a range of axis

2011-01-11 Thread Yuan Jian
Hi,
I am using plot to show scatter points in 2_D.
in my data, there is no data between -1 and +1 in x-axis.
I want to skip this region, i.e. x axis becomes [-Inf:-1, 1:Inf].
can any one tell me how to do?

YU




  
[[alternative HTML version deleted]]

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


Re: [R] GLMM with lme4 and octopus behaviour

2011-01-11 Thread Ben Bolker
Samaritan  gmail.com> writes:
 
> [snip]

  You might want to ask follow-up questions on the R-sig-mixed-models list

> At the most basic level, I'm testing the effect of sleep deprivation on
> various behaviours (e.g. amount of time spent awake, amount of time spend
> expressing difference textures, patterns, colours etc). I take a video
> sample of each octopus every hour for 72 hours and score their behaviour in
> Jwatcher.
> 
> I'm using GLMM so that I can nest Individual as a random factor within Time,
> which I'm told will reduce the problem of making repeated measures (is this
> effectively blocking by time and by octopus?). Currently my model looks like
> this:
> 
> octopus.lmer<-lmer(awake~as.factor(Treatment)+Sex+Weight+(1|Time/Octopus))

  You say you want to use GLMM -- presumably awake is a binary variable
that you want to treat as such? If so, you need the argument
'family=binomial' in your model. (You might want to use the 'glmer'
function instead, for clarity, although in practice R takes care of
this for you.)

> 
> where "Treatment" is "0" or "1" representing sleep-deprived or
> sleep-allowed. 

  As a general practice you should probably code your data as
"deprived"/"allowed" rather than 0/1: you won't have to use as.factor
and you will automatically be able to keep track of the coding.

When I try to fit the model I get the following error
> messages:
> 
> Error: length(f1) == length(f2) is not TRUE
> In addition: Warning messages:
> 1: In Octopus:Time :
>   numerical expression has 190 elements: only the first used
> 2: In Octopus:Time :
>   numerical expression has 190 elements: only the first used
> 
> This is the point at which I become lost. What does this mean? Clearly I'm
> not doing something right, perhaps in my data preparation? So far as I can
> see the length of each of the variables is the same (although I'm no certain
> as to what f1 and f2 refer to).

  This is coming from inside lmer.  My guess is that you want
to make sure Time is a factor.

dat <- expand.grid(Time=1:5,Octopus=1:5)
dat$awake <- sample(0:1,prob=c(0.5,0.5),replace=TRUE,size=nrow(dat))

dat$Time <- factor(dat$Time)
dat$Octopus <- factor(dat$Octopus)

library(lme4)

lmer(awake~(1|Time/Octopus),data=dat)

  I got a similar error when both Time and Octopus were numeric.
When I turned one or the other but not both into a factor I got
different errors.  When they are both factors I get an error which
is related to the fact that I don't have enough data (only one
observation per block, so the block effects are confounded with
the residual variation)

lmer(awake~(1|Time)+(1|Octopus),data=dat)

which represents a crossed effect, does work.

> If anyone could offer me some kind of advice about this I would appreciate
> it very, VERY much. Both of my supervisors have no experience with R and so
> have kind of washed their hands, so I'm alone in this and your expertise
> would be a big help.
> 
> Dean

  Even though R is wonderful, if your supervisors use a different
statistical package, I would strongly recommend you stick to what 
they use so that you can get properly trained,  unless you are *very*
self-sufficient.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] extracting more information from optim in R?

2011-01-11 Thread Steve Su
Dear All,

I am in the process of translating some of my functions into C and one such 
function sometimes get stuck when I applied optim. The pure R codes work fine 
but the translated C codes occasionally get stuck and I like to find out why. 

Is it possible to extract the parameters of function calculated in each step of 
the Nelder-Simplex algorithm, the trace option only produce the value of the 
objective function in each iteration? Is there a simple way to extract this 
information?

Any help greatly appreciated.

###

Assistant Professor Steve Su
School of Mathematics and Statistics
Faculty of Engineering, Computing and Mathematics

M019, 35 Stirling Highway
Crawley, 6009, WA, Australia

Phone:+6164883369

http://www.uwa.edu.au/people/steve.su
CRICOS Provider Code: 00126G
[[alternative HTML version deleted]]

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


[R] Degrees of freedom

2011-01-11 Thread Umit Tokac
Hello,

I have a little problem about degree of freedom in R.
if you can help me, I will be happy.
I used nlme function to analyze my data and run the linear mixed 
effects model in R.
I did the linear mixed effect analysis in SAS and SPSS as well.
However, R gave the different degrees of freedom than SAS and SPSS did.
Can you help me to learn what the reason is to obtain different 
degrees of freedom from R?

Thanks 

Umit Tokac
Graduate Student 
Measurement and Statistics 
Florida State University
Tel:(850)345-7487  

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