Re: [R] playing with R: make a animated GIF file...
klebyn a écrit : my objective: 1) to save PNG files; - i don't know the best way to make this; ?png (also bmp(), jpg() available) hih __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem in installing a package
Claire Lee wrote: I'm using R in Windows XP. I created a package myself. I've used R CMD check to check it. Everything seems OK except the latex. I get the error message: * checking bbHist-manual.tex ... ERROR LaTeX errors when creating DVI version. This typically indicates Rd problems. I ignored it because I didn't want to submit it to CRAN. Then I tried to use R CMD INSTALL to install it. First I get: mv: cannot move `c:/PROGRA~1/R/rw2011/library/bbHist' to `c:/PROGRA~1/R/rw2011/library/00LOCK/bbHist ': Permission denied and a bunch of making DLL errors. Then when I tried a second time, I get: open(c:/progra~1/r/rw2011/library/bbHist/DESCRIPTION): No such file or directory I can see a 00LOCK directory is created in the c:/PROGRA~1/R/rw2011/library directory. Any idea why this is happening? No, information is still too sparse, unfortunately. Do you have full write access? The 00LOCK directory is used to save the older package in order to be able to restore it if a new installtion fails. Something went wrong and you have to remove it manually now, I guess. Uwe Ligges Thanks. Claire __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] how to use tune.knn() for dataset with missing values
Hi Everybody, i again have the problem in using tune.knn(), its giving an error saying missing values are not allowed again here is the script for BreastCancer Data, library(e1071) library(mda) trdata-data.frame(train,row.names=NULL) attach(trdata) xtr - subset(trdata, select = -Class) ytr - Class bestpara -tune.knn(xtr,ytr, k = 1:25, tunecontrol = tune.control(sampling = cross)) and here i got the mentioned error. can anybody help me in this regard... Thanks Regards, Uttam Phulwale Tata Consultancy Services Limited Mailto: [EMAIL PROTECTED] Website: http://www.tcs.com [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Changing the value of variables passed to functions as arguments
Paul Baer [EMAIL PROTECTED] writes: Is it possible to write functions in such a way that, rather than having to write a=function(a), one can just write function(a) and have the variable passed as the argument be modified? My real interest here is being able to invoke the editor by writing ed(filename) rather than filename=edit(filename). You might want to check out fix()... -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Testing strata by covariate interactions in coxph
Dear list members, I am working with a Cox ph model for the duration of unemployment. The event of interest in my analysis is getting employed. I have various background variables explaining this event: age, sex, education etc. I have multiple unemployment spells per person. I use a model with person-specific frailty terms in order to take into account the correlation of spells by the same person. The persons can be divided into 3 groups, say A, B and C. I am interested to see whether there are differences between estimated covariate effects between the groups. Therefore I specify a model with strata by covariate interactions. I would like to conduct a Wald test for the null hypothesis no differences between any covariate effects in the 3 groups. This is similar to the example by Therneau Grambsch in their book Modeling survival data. Extending the Cox model, p. 47, except that I have interactions with many covariates that I would like to test jointly (and a frailty term). In S-plus there seems to be a function waldtest that does the job. Is there anything similar in R that I could use? Here is the code for my model. The covariates are a subset of all the covariates that I use. fit=coxph(Surv(tykesto, tyoll)~nainen*strata(group)+I(ika-mean(ika))*strata(group) +I(ika2-mean(ika2))*strata(group)+keski*strata(group)+korkea*strata(group) +frailty(hnro),data=coxdata) tyoll2 =1, if event getting employed has occurred, 0 otherwise nainen = 1, if female 0 otherwise ika = age in years at the start of the spell ika2= age squared keski =1, if secondary education, 0 otherwise korkea=1, if higher education, 0 otherwise hnro= person identifier group = group identifier Thank you in advance for any help. Marjo Pyy-Martikainen __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] multiple line plots
Marc Schwartz wrote:
On Wed, 2005-10-05 at 22:19 +1000, sosman wrote:
I have some data in a CSV file:
time,pos,t,tl
15:23:44:350,M1_01,4511,1127
15:23:44:350,M1_02,4514,1128
15:23:44:350,M1_03,4503,1125
...
15:23:44:491,M2_01,4500,1125
15:23:44:491,M2_02,4496,1124
15:23:44:491,M2_03,4516,1129
...
15:23:44:710,M3_01,4504,1126
15:23:44:710,M3_02,4516,1129
15:23:44:710,M3_03,4498,1124
...
Each pos (eg M1_01) is an independent time series. I would like to plot
each time series as lines on a single plot and I wondered if there was
something more straight forward than I was attempting.
I got as far as:
fname = 't100.csv'
t = read.csv(fname)
tpos = split(t, t$pos)
plot(tpos[[M1_01]]$t, type='l')
for (p in names(tpos)) {
lines(tpos[[p]]$t)
}
which seems to work but then I got stuck on how to make each line a
different colour and figured that there might a be a one liner R command
to do what I want.
Any tips would be appreciated.
See the examples in ?plot.ts for some approaches.
You will need to review ?ts to create time series objects from your data
to be used in plot.ts().
Another approach, which is not specific to time series, is ?matplot.
The matplot example looks like the go.
The example data didn't really show the grouping and even though I
mentioned time series, simply plotting the t values as an ordered
sequence is fine for this application (sorry about the red herring).
The dataset below is what I should have shown:
post
1 M1_01 4511
2 M1_02 4514
3 M1_03 4503
4 M1_01 4500
5 M1_02 4496
6 M1_03 4516
7 M1_01 4504
8 M1_02 4516
9 M1_03 4498
So what I ended up with was:
# Make a wide data set
tw = unstack(t, t ~ pos)
# Results in a list since not all series the same length
# Find the shortest dataset
len = min(sapply(tw, length))
setlen = function(l, newlen) { length(l) = newlen }
# Not sure why this did not work
#sapply(tw, setlen, len)
for (n in names(tw)) {
length(tw[[n]]) = len
}
matplot(data.frame(tw), type='l')
Apart from flying a bit blind, I obtained the plot I was after.
thanks
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] --gui none problem
Hi, As part of a batch process I run R with option --gui none. This works nicely on one maching running Version 1.9.0 Under development (unstable) (2004-03-06) but on another machine using the newer Version 2.1.0 (2005-04-18) I get the following fatal error message: ERROR: unknown GUI none Does anyone know if this is a bug or a problem with the installation? -- Morten Lindow, Bioinformatics Centre, University of Copenhagen phone: +45 3532 1348 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] --gui none problem
Version 2.x no longer has a -gui=none option. Running BATCH jobs supposedly disables a gui though. Luc Vereecken At 10:43 AM 10/6/2005, Morten Lindow wrote: Hi, As part of a batch process I run R with option --gui none. This works nicely on one maching running Version 1.9.0 Under development (unstable) (2004-03-06) but on another machine using the newer Version 2.1.0 (2005-04-18) I get the following fatal error message: ERROR: unknown GUI none Does anyone know if this is a bug or a problem with the installation? -- Morten Lindow, Bioinformatics Centre, University of Copenhagen phone: +45 3532 1348 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] --gui none problem
On Thu, 6 Oct 2005, Luc Vereecken wrote: Version 2.x no longer has a -gui=none option. Running BATCH jobs supposedly disables a gui though. There is no --gui=none option in 2.1.x (but there is in 2.0.x). The gui is not disabled in 2.1.x either according to the help page nor in practice (although using an actual GUI without an operator in attendance is not in general sensible, and any error will terminate R). Here is the NEWS item for R 2.1.0 o BATCH on Unix no longer sets --gui=none as the X11 module is only loaded if needed. Try R CMD BATCH --gui=tk (Should be Tk, but that was a buglet in 2.1.x). Luc Vereecken At 10:43 AM 10/6/2005, Morten Lindow wrote: Hi, As part of a batch process I run R with option --gui none. This works nicely on one maching running Version 1.9.0 Under development (unstable) (2004-03-06) but on another machine using the newer Version 2.1.0 (2005-04-18) I get the following fatal error message: ERROR: unknown GUI none Does anyone know if this is a bug or a problem with the installation? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] multiple line plots
# Not sure why this did not work
#sapply(tw, setlen, len)
It probably did, but you discarded the result. Try
tw - sapply(tw, setlen, len)
On Thu, 6 Oct 2005, sosman wrote:
Marc Schwartz wrote:
On Wed, 2005-10-05 at 22:19 +1000, sosman wrote:
I have some data in a CSV file:
time,pos,t,tl
15:23:44:350,M1_01,4511,1127
15:23:44:350,M1_02,4514,1128
15:23:44:350,M1_03,4503,1125
...
15:23:44:491,M2_01,4500,1125
15:23:44:491,M2_02,4496,1124
15:23:44:491,M2_03,4516,1129
...
15:23:44:710,M3_01,4504,1126
15:23:44:710,M3_02,4516,1129
15:23:44:710,M3_03,4498,1124
...
Each pos (eg M1_01) is an independent time series. I would like to plot
each time series as lines on a single plot and I wondered if there was
something more straight forward than I was attempting.
I got as far as:
fname = 't100.csv'
t = read.csv(fname)
tpos = split(t, t$pos)
plot(tpos[[M1_01]]$t, type='l')
for (p in names(tpos)) {
lines(tpos[[p]]$t)
}
which seems to work but then I got stuck on how to make each line a
different colour and figured that there might a be a one liner R command
to do what I want.
Any tips would be appreciated.
See the examples in ?plot.ts for some approaches.
You will need to review ?ts to create time series objects from your data
to be used in plot.ts().
Another approach, which is not specific to time series, is ?matplot.
The matplot example looks like the go.
The example data didn't really show the grouping and even though I
mentioned time series, simply plotting the t values as an ordered
sequence is fine for this application (sorry about the red herring).
The dataset below is what I should have shown:
post
1 M1_01 4511
2 M1_02 4514
3 M1_03 4503
4 M1_01 4500
5 M1_02 4496
6 M1_03 4516
7 M1_01 4504
8 M1_02 4516
9 M1_03 4498
So what I ended up with was:
# Make a wide data set
tw = unstack(t, t ~ pos)
# Results in a list since not all series the same length
# Find the shortest dataset
len = min(sapply(tw, length))
setlen = function(l, newlen) { length(l) = newlen }
# Not sure why this did not work
#sapply(tw, setlen, len)
for (n in names(tw)) {
length(tw[[n]]) = len
}
matplot(data.frame(tw), type='l')
Apart from flying a bit blind, I obtained the plot I was after.
thanks
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Compare two distance matrices
Hi all, I am trying to compare two distance matrices with R. I would like to create a XY plot of these matrices and do some linear regression on it. But, I am a bit new to R, so i have a few questions (I searched in the documentation with no success). The first problem is loading a distance matrix into R. This matrix is the output of a the Phylip program Protdist and lookes like this: 5 n_crassa0.00 0.690737 0.895257 0.882576 2.365386 c_neufor0.690737 0.00 0.956910 0.979988 2.103041 a_thaliana 0.895257 0.956910 0.00 1.003668 2.724847 pompep 0.882576 0.979988 1.003668 0.00 2.065202 s_cerevis 2.365386 2.103041 2.724847 2.065202 0.00 5 n_crassa0.00 0.739560 0.933986 0.861644 2.207467 c_neufor0.739560 0.00 0.988779 0.925168 1.941141 a_thaliana 0.933986 0.988779 0.00 1.007803 2.415320 pompep 0.861644 0.925168 1.007803 0.00 2.394490 s_cerevis 2.207467 1.941141 2.415320 2.394490 0.00 . I tried with the scan() function to load the files, but with no success. How should i load in these files? Second i saw there is package for distance matrices (http://stat.ethz.ch/R-manual/R-devel/library/stats/html/dist.html). I thought of using as.dist() to convert the files to a R dist matrix. I have not been able to try this, because the first step didn't succeed. Is this the right approach, or is there any other way of comparing distance matrices with R? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] playing with R: make a animated GIF file...
you could first make the frames to put together in R, then outside of R
glue the frames into a gif.
eg:
frames-10
for(i in 1:frames) {
jpeg(paste(ani_, i, .jpg, sep = ))
plot(1:10,1:10, col = i)
dev.off()
}
then use an image editing program to glue the jpgs together -- eg free and
crossplatform imagemagick (www.imagemagick.org).
under linux and once imagemagick is installed, the command
$ convert -delay 10 ani_*.jpg animation.gif
makes a looping animated gif `animation.gif' with 10/100 sec. between
frame flips.
Katharine Mullen
Department of Physics and Astronomy
Faculty of Sciences
Vrije Universiteit
de Boelelaan 1081
1081 HV Amsterdam
The Netherlands
room: T.1.06
tel: +31 205987870
fax: +31 205987992
e-mail: [EMAIL PROTECTED]
http://www.nat.vu.nl/~kate/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Interpolation in time
Can anybody help me write a code on the following data example, which fills out all NA values by using a linear interpolation with the two closest values? Doy is day of year (%j). Code example: yr-c(rep(2000,14)) doy-c(16:29) dat-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) ta-cbind(yr,doy,dat) ta yr doy dat [1,] 2000 16 3.2 [2,] 2000 17 NA [3,] 2000 18 NA [4,] 2000 19 NA [5,] 2000 20 NA [6,] 2000 21 NA [7,] 2000 22 NA [8,] 2000 23 5.1 [9,] 2000 24 NA [10,] 2000 25 NA [11,] 2000 26 NA [12,] 2000 27 NA [13,] 2000 28 NA [14,] 2000 29 4.6 Anette Norgaard [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] eliminate t() and %*% using crossprod() and solve(A,b)
you might want to see Douglas Bates. Least squares calculations in R. R News, 4(1):17-20, June 2004. http://CRAN.R-project.org/doc/Rnews/ he gives some rules of thumb, eg use solve(A,b) not solve(A) %*% b use crossprod(X) not t(X) %*% X use crossprod(X,y) not t(X) y Katharine Mullen Department of Physics and Astronomy Faculty of Sciences Vrije Universiteit de Boelelaan 1081 1081 HV Amsterdam The Netherlands room: T.1.06 tel: +31 205987870 fax: +31 205987992 e-mail: [EMAIL PROTECTED] http://www.nat.vu.nl/~kate/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Animation of Mandelbrot Set
Very great job
On 10/4/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote:
Hi,
I was playing with Mandelbrot sets and come up with the following code, I
thought I would share:
library(fields) # for tim.colors
library(caTools) # for write.gif
m = 400 # grid size
C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ),
imag=rep(seq(-1.2,1.2, length.out=m), m ) )
C = matrix(C,m,m)
Z = 0
X = array(0, c(m,m,20))
for (k in 1:20) {
Z = Z^2+C
X[,,k] = exp(-abs(Z))
}
image(X[,,k], col=tim.colors(256)) # show final image in R
write.gif(X, Mandelbrot.gif, col=tim.colors(256), delay=100)
# drop Mandelbrot.gif file from current directory on any web brouser to
see the animation
Jarek
\
Jarek Tuszynski, PhD. o / \
Science Applications International Corporation \__,|
(703) 676-4192 \
[EMAIL PROTECTED] ` \
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
--
Xiaohua Dai, Dr.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] R-2.2.0 is released
I've rolled up R-2.2.0.tar.gz a short while ago. This version contains
several changes and additions, mostly incremental. See the full list
of changes below.
You can get it from
http://cran.r-project.org/src/base/R-2/R-2.2.0.tar.gz
or wait for it to be mirrored at a CRAN site nearer to you. Binaries
for various platforms will appear in due course.
There is also a version split for floppies.
For the R Core Team
Peter Dalgaard
These are the md5sums for the freshly created files, in case you wish
to check that they are uncorrupted:
94d55d512a9ba36caa9b7df079bae19f COPYING
d8045f3b8f929c1cb29a1e3fd737b499 COPYING.LIB
043a28ec5378bfaba88e4fb34f805980 FAQ
70447ae7f2c35233d3065b004aa4f331 INSTALL
5209c94d85a195fb92cdf796912a732b NEWS
88bbd6781faedc788a1cbd434194480c ONEWS
4f004de59e24a52d0f500063b4603bcb OONEWS
6bddf439ae417a48bd31892996ea111c R-2.2.0.tar.gz
f8763b77147796b3adf52045183ee0c3 R-2.2.0.tar.gz-split.aa
ba00cb5ff9c3e82038c3b3abcce60855 R-2.2.0.tar.gz-split.ab
9668413beca51736390b63afa489b2f1 R-2.2.0.tar.gz-split.ac
b19e3a225a66b50e14671f1bb36e1d07 R-2.2.0.tar.gz-split.ad
2465e208aab735e20d1efee7c72f6c23 R-2.2.0.tar.gz-split.ae
903f37e74de637e71ef619c5801f719e R-2.2.0.tar.gz-split.af
46502602ec014ba2f261ed6c81811ea6 R-2.2.0.tar.gz-split.ag
58d5e7d99ec15388687f2a7dca78b647 R-2.2.0.tar.gz-split.ah
d8c2356d0e3e650b5bfc92e5ee22a91d R-2.2.0.tar.gz-split.ai
cabdf55568d9f90115faaaf18cddfa07 R-2.2.0.tar.gz-split.aj
56a780cdec835c5c598f8dfc0738f7f3 README
020479f381d5f9038dcb18708997f5da RESOURCES
Here is the relevant bit of the NEWS file:
CHANGES IN R VERSION 2.2.0
USER-VISIBLE CHANGES
o plot(lm object) uses a new default 'which = 5'
for the fourth panel when 'which' is not specified.
o The SVN revision number will appear after the date in the
welcome message. The date shown is now the date of the last
change to the sources rather than the date the sources were
prepared.
o is.null(expression()) now returns FALSE. Only NULL gives TRUE
in is.null().
o graphics::xy.coords, xyz.coords and n2mfrow have been moved to
the grDevices name space (to be available for grid as well).
graphics::boxplot.stats, contourLines, nclass.*, and chull
have been moved to the grDevices name space. The C code
underlying chull() has been moved to package grDevices.
o split(x, f), split-() and unsplit() now by default split by all
levels of a factor f, even when some are empty.
Use split(x, f, drop = TRUE) if you want the old behavior of
dropping empty levels. split() and split-() are S3 generic
functions with new arguments 'drop' and '...' and all methods now
should have 'drop' and '...' arguments as well.
o The default for 'allowEscapes' in both read.table() and scan()
has been changed to FALSE.
o The default for 'gcFirst' in system.time() is now TRUE.
NEW FEATURES
o .Platform has a new component 'path.sep', the separator used
between paths in environment variables such as PATH and TEXINPUTS.
o anova.mlm() now handles the single-model case.
o Hexadecimal values are now allowed for as.numeric() and
as.integer() on all platforms, and as integer constants in R code.
o attach() now prints an information message when objects are
masked on the search path by or from a newly attached database.
o axis() now returns 'at' positions.
o axis() has a new argument 'hadj' to control horizontal
adjustment of labels.
o axis.Date() and axis.POSIXct() now accept a 'labels' argument
(contributed by Gavin Simpson).
o barplot() now has arguments 'log = ' and 'add = FALSE'
(as in barplot2() from package 'gplots').
o baseenv() has been added, to return the base environment. This
is currently still NULL, but will change in a future release.
o boxplot() now responds to supplying 'yaxs' (via bxp()).
(Wish of PR#8072.)
o capabilities() has a new component 'NLS'.
o cbind() and rbind() now react to 'deparse.level' = {0,1,2}
(as in another system not unlike R).
o Experimental versions of cbind() and rbind() in methods package,
based on new generic function cbind2(x,y) and rbind2(). This will
allow the equivalent of S4 methods for cbind() and rbind() ---
currently only after an explicit activation call, see ?cbind2.
o New functions cdplot() and spineplot() for conditional density
plots and spine plots or spinograms. Spine plots are now used
instead of bar plots for x-y scatterplots where y is a factor.
o checkDocFiles() in package 'tools' now checks for bad \usage
lines (syntactically invalid R code).
o The nonparametric variants of cor.test() now behave better in
the presence of ties. The
Re: [R] Animation of Mandelbrot Set
Anyone know why I would get an Error: couldn't find function write.gif
despite loading library(caTools) with no errors in R 2.1.1 under XP?
Thanks,
Roger
library(fields) # for tim.colors
fields is loaded use help(fields) for an overview of this library
library(caTools) # for write.gif
Loading required package: bitops
m = 400 # grid size
C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ),
+ imag=rep(seq(-1.2,1.2, length.out=m), m ) )
C = matrix(C,m,m)
Z = 0
X = array(0, c(m,m,20))
for (k in 1:20) {
+ Z = Z^2+C
+ X[,,k] = exp(-abs(Z))
+ }
image(X[,,k], col=tim.colors(256)) # show final image in R
write.gif(X, Mandelbrot.gif, col=tim.colors(256), delay=100)
Error: couldn't find function write.gif
On 10/4/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote:
Hi,
I was playing with Mandelbrot sets and come up with the following code, I
thought I would share:
library(fields) # for tim.colors
library(caTools) # for write.gif
m = 400 # grid size
C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ),
imag=rep(seq(-1.2,1.2, length.out=m), m ) )
C = matrix(C,m,m)
Z = 0
X = array(0, c(m,m,20))
for (k in 1:20) {
Z = Z^2+C
X[,,k] = exp(-abs(Z))
}
image(X[,,k], col=tim.colors(256)) # show final image in R
write.gif(X, Mandelbrot.gif, col=tim.colors(256), delay=100)
# drop Mandelbrot.gif file from current directory on any web brouser to
see the animation
Jarek
\
Jarek Tuszynski, PhD. o / \
Science Applications International Corporation \__,|
(703) 676-4192 \
[EMAIL PROTECTED] ` \
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Animation of Mandelbrot Set
The binary version of caTools on CRAN (1.0) does not have the write.gif function but the source version (1.4) does ... hth, ingmar From: roger bos [EMAIL PROTECTED] Reply-To: roger bos [EMAIL PROTECTED] Date: Thu, 6 Oct 2005 08:14:42 -0400 To: Tuszynski, Jaroslaw W. [EMAIL PROTECTED] Cc: \([EMAIL PROTECTED]) r-help@stat.math.ethz.ch Subject: Re: [R] Animation of Mandelbrot Set Anyone know why I would get an Error: couldn't find function write.gif despite loading library(caTools) with no errors in R 2.1.1 under XP? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] playing with R: make a animated GIF file...
See write.gif function in caTools.
Jarek
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of klebyn
Sent: Wednesday, October 05, 2005 8:05 PM
To: r-help@stat.math.ethz.ch
Subject: [R] playing with R: make a animated GIF file...
Hello all
I am playing with R for to make a animated GIF.
any suggestions, improvements are welcome :-)
case somebody could help me, i thanks!
Cleber N. Borges ( klebyn )
my objective:
(steps TODO)
---
1) to save PNG files;
- i don't know the best way to make this;
2) transform the PNG files into GIF files (easy! no problem! ... i think
...)
3) reload the GiF files in R and use the caTools package to make a
animated GIF.
--
the code
reverse the STRING
strReverse - function(x) sapply(lapply(strsplit(x, NULL), rev), paste,
collapse=)
logotype to animate
yourLogo =Is Nice to play with R-package
logoWidth = 1.5
logoHeight = 2.5
L = nchar(yourLogo)
TrigSplit = 360 / L
yourLogo = strReverse(yourLogo)
posx = numeric(L)
posy = numeric(L)
for( i in 0:L){
posx[i] = logoHeight * sin(i * TrigSplit * pi / 180)
posy[i] = logoWidth * cos(i * TrigSplit * pi / 180)
}
max_x = max(posx)*1.1
max_y = max(posy)*3
min_x = min(posx)*1.1
min_y = min(posy)*3
cex = 2/(posy + 2)
idx = 1:L
for(j in 1:L-1) {
###file = paste(CQM_,j,.png,sep=)
###png(filename=file, bg=transparent)
plot(0,t='n', xlim=c(min_x,max_x), ylim=c(min_y,max_y), axes=FALSE,
ann=FALSE, font=3 )
for( i in 1:L){text(x=posx[i], y=posy[i],
labels=substr(yourLogo,idx[i],idx[i]), col='blue', cex=cex[i] ) }
idx = (append(idx[L],idx))[1:L]
Sys.sleep(0.2)
###dev.off()
}
## final code
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Interpolation in time
Is this what you want? yr-c(rep(2000,14)) doy-c(16:29) dat-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) ta-cbind(yr,doy,dat) ta yr doy dat [1,] 2000 16 3.2 [2,] 2000 17 NA [3,] 2000 18 NA [4,] 2000 19 NA [5,] 2000 20 NA [6,] 2000 21 NA [7,] 2000 22 NA [8,] 2000 23 5.1 [9,] 2000 24 NA [10,] 2000 25 NA [11,] 2000 26 NA [12,] 2000 27 NA [13,] 2000 28 NA [14,] 2000 29 4.6 good - !is.na(ta[,'dat']) x.f - approxfun(ta[good,'doy'], ta[good,'dat'], rule=2) ta[!good, 'dat'] - x.f(ta[!good, 'doy']) ta yr doy dat [1,] 2000 16 3.20 [2,] 2000 17 3.471429 [3,] 2000 18 3.742857 [4,] 2000 19 4.014286 [5,] 2000 20 4.285714 [6,] 2000 21 4.557143 [7,] 2000 22 4.828571 [8,] 2000 23 5.10 [9,] 2000 24 5.016667 [10,] 2000 25 4.93 [11,] 2000 26 4.85 [12,] 2000 27 4.77 [13,] 2000 28 4.68 [14,] 2000 29 4.60 On 10/6/05, Anette Nørgaard [EMAIL PROTECTED] wrote: Can anybody help me write a code on the following data example, which fills out all NA values by using a linear interpolation with the two closest values? Doy is day of year (%j). Code example: yr-c(rep(2000,14)) doy-c(16:29) dat-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) ta-cbind(yr,doy,dat) ta yr doy dat [1,] 2000 16 3.2 [2,] 2000 17 NA [3,] 2000 18 NA [4,] 2000 19 NA [5,] 2000 20 NA [6,] 2000 21 NA [7,] 2000 22 NA [8,] 2000 23 5.1 [9,] 2000 24 NA [10,] 2000 25 NA [11,] 2000 26 NA [12,] 2000 27 NA [13,] 2000 28 NA [14,] 2000 29 4.6 Anette Norgaard [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Jim Holtman Cincinnati, OH +1 513 247 0281 What the problem you are trying to solve? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Animation of Mandelbrot Set
http://cran.r-project.org/src/contrib/Descriptions/caTools.html has current versions (1.4) of both source and binary. Jarek \ Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] ` \ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Sent: Thursday, October 06, 2005 8:19 AM To: roger bos; Tuszynski, Jaroslaw W. Cc: ([EMAIL PROTECTED]) Subject: Re: [R] Animation of Mandelbrot Set The binary version of caTools on CRAN (1.0) does not have the write.gif function but the source version (1.4) does ... hth, ingmar From: roger bos [EMAIL PROTECTED] Reply-To: roger bos [EMAIL PROTECTED] Date: Thu, 6 Oct 2005 08:14:42 -0400 To: Tuszynski, Jaroslaw W. [EMAIL PROTECTED] Cc: \([EMAIL PROTECTED]) r-help@stat.math.ethz.ch Subject: Re: [R] Animation of Mandelbrot Set Anyone know why I would get an Error: couldn't find function write.gif despite loading library(caTools) with no errors in R 2.1.1 under XP? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] modeling language for optimization problems
Prof, My apologies if you believe I was contradicting you. I though I was concurring with you that 1) R does not have anything that is as powerful as NuOpt and 2) there is no 'one' optimization package in R. I did, however, want to point out that R does have separate packages that offer a fairy broad range of the simpler types of optimzation. I mentioned that R can do QPs , LPs, and integer programming, among others. Thanks, Roger On 10/3/05, Prof Brian Ripley [EMAIL PROTECTED] wrote: On Mon, 3 Oct 2005, roger bos wrote: As others have alluded, R does not have any one package that is as versatile and powerful as NuOpt, but R does have many different optimization packages, so you can do LP, QP, Integer programming, and many more types of optimization, all without having to learn a new language. But being free, they are admittedly not as powerful as NuOpt. The main thing I know you can do in NuOpt that you cannot do in R (to my knowledge) is mixed integer quadratic programming, which would be nice, but there are many work arounds. Please tell us which R packages provide - primal-dual interior point method based on line search for general Convex Programming (CP) models including convex Quadratic Programming (CQP) models. - primal-dual interior point method based on trust region method for general Non-Linear Programming (NLP) models. - primal-dual interior point method based on quasi-Newton method for general Non-Linear Programming (NLP) models. - active set method for convex Quadratic Programming (CQP) models and mixed integer Quadratic Programming(MIQP) models.' and how you found them? I get help.search(interior point) ... rq.fit.fn(quantreg) Quantile Regression Fitting via Interior Point Methods rq.fit.fnb(quantreg) Quantile Regression Fitting via Interior Point Methods rq.fit.fnc(quantreg) Quantile Regression Fitting via Interior Point Methods help.search(convex programming) No help files found with alias or concept or title matching 'convex programming' using fuzzy matching. ... Thanks, Roger On 10/3/05, Prof Brian Ripley [EMAIL PROTECTED] wrote: On Mon, 3 Oct 2005, Huntsinger, Reid wrote: Have you looked at the R interface to GLPK (the GNU Linear Programming Kit)? http://cran.r-project.org/src/contrib/Descriptions/glpk.html NUOPT is not just about LP: the subject was `language for optimization'. Its manual says `NUOPT is a collection of powerful optimization methods, including: - primal-dual interior point method with higher order correction for Linear Programming (LP) models. - simplex method for Linear Programming (LP) and mixed integer programming (MILP) models. - primal-dual interior point method based on line search for general Convex Programming (CP) models including convex Quadratic Programming (CQP) models. - primal-dual interior point method based on trust region method for general Non-Linear Programming (NLP) models. - primal-dual interior point method based on quasi-Newton method for general Non-Linear Programming (NLP) models. - active set method for convex Quadratic Programming (CQP) models and mixed integer Quadratic Programming(MIQP) models.' In any case, I don't see GLPK as a `language' and other LP solvers are available in R. From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Prof Brian Ripley Sent: Sunday, October 02, 2005 10:38 AM To: Paolo Cavatore Cc: r-help@stat.math.ethz.ch Subject: Re: [R] modeling language for optimization problems On Sun, 2 Oct 2005, Paolo Cavatore wrote: Does anyone know whether R has its own modeling language for optimization problems (like SIMPLE in NuOPT for S-plus)? No. Note that SIMPLE is the language of NUOPT, not of S-PLUS. There is an (extra-cost) interface module S+NUOPT, but it is an interface to NUOPT's engine. As far as I am aware R itself covers almost none of the ground of S+NUOPT, and available packages cover only a small part of it. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865
[R] how to use tune.knn() for dataset with missing values
Well, since knn() can't handle incomplete data as it says, you can choose to either omit incomplete observations (e.g., using na.omit()), or to impute the data if the conditions are met (missingness at random, ...); see, e.g., packages cat, mix, norm, and e1071 for that. HTH, David Hi Everybody, i again have the problem in using tune.knn(), its giving an error saying missing values are not allowed again here is the script for BreastCancer Data, library(e1071) library(mda) trdata-data.frame(train,row.names=NULL) attach(trdata) xtr - subset(trdata, select = -Class) ytr - Class bestpara -tune.knn(xtr,ytr, k = 1:25, tunecontrol = tune.control(sampling = cross)) and here i got the mentioned error. can anybody help me in this regard... Thanks Regards, Uttam Phulwale Tata Consultancy Services Limited Mailto: [EMAIL PROTECTED] Website: http://www.tcs.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Compare two distance matrices
Hi, hi all, I am trying to compare two distance matrices with R. I would like to create a XY plot of these matrices and do some linear regression on it. But, I am a bit new to R, so i have a few questions (I searched in the documentation with no success). The first problem is loading a distance matrix into R. This matrix is the output of a the Phylip program Protdist and lookes like this: I tried with the scan() function to load the files, but with no success. How should i load in these files? you can separately load each matrix with two text files. require(ade4) mat1 - read.table(mat1.txt) nam1 - mat1[,1] mat1 - mat1[,-1] row.names(mat1) - names(mat1) - nam1 mat2 - read.table(mat2.txt) nam2 - mat2[,1] mat2 - mat2[,-1] row.names(mat2) - names(mat2) - nam2 dist1 - mat2dist(mat1) dist2 - mat2dist(mat2) # with mat1: # n_crassa0.00 0.690737 0.895257 0.882576 2.365386 # c_neufor0.690737 0.00 0.956910 0.979988 2.103041 # a_thaliana 0.895257 0.956910 0.00 1.003668 2.724847 # pompep 0.882576 0.979988 1.003668 0.00 2.065202 # s_cerevis 2.365386 2.103041 2.724847 2.065202 0.00 # and mat2: # n_crassa0.00 0.739560 0.933986 0.861644 2.207467 # c_neufor0.739560 0.00 0.988779 0.925168 1.941141 # a_thaliana 0.933986 0.988779 0.00 1.007803 2.415320 # pompep 0.861644 0.925168 1.007803 0.00 2.394490 # s_cerevis 2.207467 1.941141 2.415320 2.394490 0.00 I think that its the more simple solution, NOT the optimal solution. Maybe, there is an interface to read Phylip file in Bioconductor project(?). To transform matrix to dist, you can use the function mat2dist of the library ade4 (see example) To compare 2 distances matrices, there are many possibility ! For example, if the distances matrices are Euclidian, you can used directly principal coordinate analyses (dudi.pco), and co-inertia analysis or procuste analysis,etc ... # examples: # Data preparation require(ade4) mat1 - read.table(mat1.txt) nam1 - mat1[,1] mat1 - mat1[,-1] row.names(mat1) - names(mat1) - nam1 mat2 - read.table(mat2.txt) nam2 - mat2[,1] mat2 - mat2[,-1] row.names(mat2) - names(mat2) - nam2 ? mat2dist dist1 - mat2dist(mat1) dist2 - mat2dist(mat2) dist1 dist2 is.euclid(dist1) is.euclid(dist2) # cool the distances matrices are euclidian :) # to compare the 2 matrices # example 1 : mantel test ? mantel.randtest mt1 - mantel.randtest(dist1,dist2,nrepet=10) plot(mt1) # Example 2: coefficient of vectorial correlation (Escoufier, 1973) ? RVdist.randtest RV1 - RVdist.randtest(dist1,dist2,nrepet=10) plot(RV1) # Example 3: coinertia analysis ?dudi.pco # ?cmdscale ?coinertia ?randtest.coinertia pco1 - dudi.pco(dist1,nf=3,scannf=F) pco2 - dudi.pco(dist2,nf=3,scannf=F) co1 - coinertia(pco1,pco2,nf=3,scannf=F) plot(co1) testco1 -randtest.coinertia(co1,nrepet=10) plot(testco1) # Example 4: procuste analysis ? procuste ? procuste.randtest pco1 - dudi.pco(dist1,nf=3,scannf=F) pco2 - dudi.pco(dist2,nf=3,scannf=F) proc1 - procuste(pco1$tab,pco2$tab) plot(proc1) testproc1 -procuste.randtest(pco1$tab,pco2$tab,nrepet=10) plot(testproc1) the choice depend to your outcome. hopes this help P.BADY __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Animation of Mandelbrot Set
I use the following code to keep my packages up to date. Is there some reason why I am not getting full function availablility with this updating code: x - packageStatus(repositories=http://cran.r-project.org/src/contrib;) st - x$avai[Status] install.packages(rownames(st)[which(st$Status==not installed)]) Thanks, Roger On 10/6/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote: http://cran.r-project.org/src/contrib/Descriptions/caTools.html has current versions (1.4) of both source and binary. Jarek \ Jarek Tuszynski, PhD. o / \ Science Applications International Corporation \__,| (703) 676-4192 \ [EMAIL PROTECTED] ` \ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Sent: Thursday, October 06, 2005 8:19 AM To: roger bos; Tuszynski, Jaroslaw W. Cc: ([EMAIL PROTECTED]) Subject: Re: [R] Animation of Mandelbrot Set The binary version of caTools on CRAN (1.0) does not have the write.gif function but the source version (1.4) does ... hth, ingmar From: roger bos [EMAIL PROTECTED] Reply-To: roger bos [EMAIL PROTECTED] Date: Thu, 6 Oct 2005 08:14:42 -0400 To: Tuszynski, Jaroslaw W. [EMAIL PROTECTED] Cc: \([EMAIL PROTECTED]) r-help@stat.math.ethz.ch Subject: Re: [R] Animation of Mandelbrot Set Anyone know why I would get an Error: couldn't find function write.gif despite loading library(caTools) with no errors in R 2.1.1 under XP? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R/S-Plus equivalent to Genstat predict: predictions over averages of covariates
Dear Peter, See the effects package, described in http://www.jstatsoft.org/counter.php?id=75url=v08/i15/effect-displays-revi sed.pdf. I hope this helps, John John Fox Department of Sociology McMaster University Hamilton, Ontario Canada L8S 4M4 905-525-9140x23604 http://socserv.mcmaster.ca/jfox -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Peter Dunn Sent: Wednesday, October 05, 2005 9:06 PM To: R-help mailing list Subject: [R] R/S-Plus equivalent to Genstat predict: predictions over averages of covariates Hi all I'm doing some things with a colleague comparing different sorts of models. My colleague has fitted a number of glms in Genstat (which I have never used), while the glm I have been using is only available for R. He has a spreadsheet of fitted means from each of his models obtained from using the Genstat predict function. For example, suppose we fit the model of the type glm.out - glm( y ~ factor(F1) + factor(F2) + X1 + poly(X2,2) + poly(X3,2), family=...) Then he produces a table like this (made up, but similar): F1(level1)12.2 F1(level2)14.2 F1(level3)15.3 F2(level1)10.3 F2(level2)9.1 X1=0 10.2 X1=0.510.4 X1=1 10.4 X1=1.510.5 X1=2 10.9 X1=2.511.9 X1=3 11.8 X2=0 12.0 X2=0.512.2 X2=1 12.5 X2=1.512.9 X2=2 13.0 X2=2.513.1 X2=3 13.5 Each of the numbers are a predicted mean. So when X1=0, on average we predict an outcome of 10.2. To obtain these figures in Genstat, he uses the Genstat predict function. When I asked for an explanation of how it was done (ie to make the predictions, what values of the other covariates were used) I was told: So, for a one-dimensional table of fitted means for any factor (or variate), all other variates are set to their average values; and the factor constants (including the first, at zero) are given a weighted average depending on their respective numbers of observations. So for quantitative variables (such as pH), one uses the mean pH in the data set when making the predictions. Reasonable anmd easy. But for categorical variables (like Month), he implies we use a weighted average of the fitted coefficients for all the months, depending on the proportion of times those factor levels appear in the data. (I hope I explained that OK...) Is there an equivalent way in R or S-Plus of doing this? I have to do it for a number of sites and species, so an automated way would be useful. I have tried searching to no avail (but may not be searching on the correct terms), and tried hard-coding something myself as yet unsuccessfully: The poly terms and the use of the weighted averaging over the factor levels are proving a bit too much for my limited skills. Any assistance appreciated. (Any clarification of what I mean can be provided if I have not been clear.) Thanks, as always. P. version _ platform i386-pc-linux-gnu arch i386 os linux-gnu system i386, linux-gnu status major2 minor1.0 year 2005 month04 day 18 language R -- Dr Peter Dunn | Senior Lecturer in Statistics Faculty of Sciences, University of Southern Queensland Web:http://www.sci.usq.edu.au/staff/dunn Email: dunn at usq.edu.au CRICOS: QLD 00244B | NSW 02225M | VIC 02387D | WA 02521C __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R/S-Plus equivalent to Genstat predict: predictions over averages of covariates
check out the effects package on CRAN. Kjetil Peter Dunn wrote: Hi all I'm doing some things with a colleague comparing different sorts of models. My colleague has fitted a number of glms in Genstat (which I have never used), while the glm I have been using is only available for R. He has a spreadsheet of fitted means from each of his models obtained from using the Genstat predict function. For example, suppose we fit the model of the type glm.out - glm( y ~ factor(F1) + factor(F2) + X1 + poly(X2,2) + poly(X3,2), family=...) Then he produces a table like this (made up, but similar): F1(level1)12.2 F1(level2)14.2 F1(level3)15.3 F2(level1)10.3 F2(level2)9.1 X1=0 10.2 X1=0.510.4 X1=1 10.4 X1=1.510.5 X1=2 10.9 X1=2.511.9 X1=3 11.8 X2=0 12.0 X2=0.512.2 X2=1 12.5 X2=1.512.9 X2=2 13.0 X2=2.513.1 X2=3 13.5 Each of the numbers are a predicted mean. So when X1=0, on average we predict an outcome of 10.2. To obtain these figures in Genstat, he uses the Genstat predict function. When I asked for an explanation of how it was done (ie to make the predictions, what values of the other covariates were used) I was told: So, for a one-dimensional table of fitted means for any factor (or variate), all other variates are set to their average values; and the factor constants (including the first, at zero) are given a weighted average depending on their respective numbers of observations. So for quantitative variables (such as pH), one uses the mean pH in the data set when making the predictions. Reasonable anmd easy. But for categorical variables (like Month), he implies we use a weighted average of the fitted coefficients for all the months, depending on the proportion of times those factor levels appear in the data. (I hope I explained that OK...) Is there an equivalent way in R or S-Plus of doing this? I have to do it for a number of sites and species, so an automated way would be useful. I have tried searching to no avail (but may not be searching on the correct terms), and tried hard-coding something myself as yet unsuccessfully: The poly terms and the use of the weighted averaging over the factor levels are proving a bit too much for my limited skills. Any assistance appreciated. (Any clarification of what I mean can be provided if I have not been clear.) Thanks, as always. P. version _ platform i386-pc-linux-gnu arch i386 os linux-gnu system i386, linux-gnu status major2 minor1.0 year 2005 month04 day 18 language R -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Singular matrix
Dear All,
I have written the following programs to find a non-singular (10*10)
covariance matrix.
Here is the program:
nitems - 10
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
while(inherits(inverse, try-error))
{
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
}
The loop doesn't stop ... This means that no non-singular matrix was
found!!!
some thing wrong !!
Thanks a lot for any reply
Bernard
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] factor : how does it work ?
Dear all, I try for long to understand exactly what is the factor type and especially how it works, but it seems too difficult for me I read paragraphs about it, and I understand quite well what it is (I think) but I still can't figure how to deal with. Especially these 2 mysteries (for me) : 1st when I make a dataframe (with the as.data.frame() or the data.frame() commands) from vectors, it seems that some columns of the dataframe (which where vectors) are factors and some not, but I didn't find an explanation for which become factor and which don't. (I know I can use I() to avoid the factor transformaton but I think it is not an optimal solution to avoid the factor type just because I don't kno how to deal with) 2d I can't manage to deal with factors, so when I have some, I transform them in vectors (with levels()), but I think I miss the power and utility of the factor type ? Any help GREATLY appreciated, best regards, Florence. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
On 10/6/2005 9:14 AM, Florence Combes wrote:
Dear all,
I try for long to understand exactly what is the factor type and especially
how it works, but it seems too difficult for me
I read paragraphs about it, and I understand quite well what it is (I think)
but I still can't figure how to deal with.
Especially these 2 mysteries (for me) :
1st when I make a dataframe (with the as.data.frame() or the data.frame()
commands) from vectors, it seems that some columns of the dataframe (which
where vectors) are factors and some not, but I didn't find an explanation
for which become factor and which don't.
(I know I can use I() to avoid the factor transformaton but I think it is
not an optimal solution to avoid the factor type just because I don't kno
how to deal with)
This is described in the ?data.frame man page: Character variables
passed to 'data.frame' are converted to factor columns unless protected
by 'I'.
2d I can't manage to deal with factors, so when I have some, I transform
them in vectors (with levels()), but I think I miss the power and utility of
the factor type ?
levels() is not the conversion you want. That lists all the levels, but
it doesn't tell you how they correspond to individual observations. For
example,
df - data.frame(x=1:3, y=c('a','b','a'))
df
x y
1 1 a
2 2 b
3 3 a
levels(df$y)
[1] a b
If you need to convert back to character values, use as.character():
as.character(df$y)
[1] a b a
For many purposes, you can ignore the fact that your data is stored as a
factor instead of a character vector. There are a few differences:
1. You can't compare the levels of a factor unless you declared it to
be ordered:
df$y[1] df$y[2]
[1] NA
Warning message:
not meaningful for factors in: Ops.factor(df$y[1], df$y[2])
but
df$y - ordered(df$y)
df$y[1] df$y[2]
[1] FALSE
However, you need to watch out here: the comparison is done by the order
of the factors, not an alphabetic comparison of their names:
levels(df$y) - c(before, after)
df
x y
1 1 before
2 2 after
3 3 before
df$y[1] df$y[2]
[1] FALSE
2. as.integer() works differently on factors: it gets the position in
the levels vector. For example,
as.integer(df$y)
[1] 1 2 1
as.integer(as.character(df$y))
[1] NA NA NA
Warning message:
NAs introduced by coercion
There are other differences, but these are the two main ones that are
likely to cause you trouble.
Duncan Murdoch
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Singular matrix
Don't you have the dimension of x backward? Try:
set.seed(1)
x - matrix(rnorm(50, 3, 3), 10, 5)
vinv - solve(crossprod(x))
vinv
[,1] [,2] [,3] [,4] [,5]
[1,] 0.019918251 -0.006247646 0.006600209 0.003687249 -0.018670806
[2,] -0.006247646 0.018121025 -0.014815905 -0.005647350 0.003434065
[3,] 0.006600209 -0.014815905 0.023411617 -0.002250342 -0.003258960
[4,] 0.003687249 -0.005647350 -0.002250342 0.025168959 -0.020070844
[5,] -0.018670806 0.003434065 -0.003258960 -0.020070844 0.039593016
If you really have 5 cases and 10 variables, the covariance matrix will
have to be singular.
Andy
From: Marc Bernard
Dear All,
I have written the following programs to find a
non-singular (10*10) covariance matrix.
Here is the program:
nitems - 10
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
while(inherits(inverse, try-error))
{
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
}
The loop doesn't stop ... This means that no non-singular
matrix was found!!!
some thing wrong !!
Thanks a lot for any reply
Bernard
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Singular matrix
'x' is a rank 5 matrix from which you want to create a rank 10
crossproduct! Try the following instead:
x - array(rnorm(100 * nitems, 3, 3), c(100, nitems))
sigma - crossprod(x)
I hope it helps.
Best,
Dimitris
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://www.med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
- Original Message -
From: Marc Bernard [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Thursday, October 06, 2005 3:11 PM
Subject: [R] Singular matrix
Dear All,
I have written the following programs to find a non-singular
(10*10) covariance matrix.
Here is the program:
nitems - 10
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
while(inherits(inverse, try-error))
{
x - array(rnorm(5*nitems,3,3), c(5,nitems))
sigma - t(x)%*%x
inverse - try(solve(sigma), TRUE)
}
The loop doesn't stop ... This means that no non-singular matrix
was found!!!
some thing wrong !!
Thanks a lot for any reply
Bernard
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Kouros Owzar is out
I will be out of the office starting 10/05/2005 and will not return until 10/07/2005. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Kouros Owzar is out
Well, cool, why do you need to tell this to more than 3000 readers of R-help ?? Kouros I will be out of the office starting 10/05/2005 and will not return until Kouros 10/07/2005. Kouros __ Kouros R-help@stat.math.ethz.ch mailing list Kouros https://stat.ethz.ch/mailman/listinfo/r-help Kouros PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html PLEASE __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Interpolation in time
On Thu, 6 Oct 2005 16:10:15 +0200 Anette Nørgaard wrote: This is exactly what I requested, thank you!! However I do actually have several columns in my data sheet where I need to do the same thing, then how do I come about that? Look at na.approx() in package zoo. Best, Z e.g. yr-c(rep(2000,14)) doy-c(16:29) dat1-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) dat2-c(2.2,NA,NA,NA,NA,NA,NA,6.1,NA,NA,NA,NA,NA,4.2) dat3-c(3.4,NA,NA,NA,NA,NA,NA,4.1,NA,NA,NA,NA,NA,4.7) ta-cbind(yr,doy,dat1,dat2,dat3) [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Interpolation in time
Is doy intended to represent the number of days since the beginning of the year? In that case convert the first two columns to class Date and interpolate using approx. See ?approx for variations: tt - as.Date(paste(yr, 1, 1, sep = -)) + doy - 1 ta[,dat] - approx(tt, dat, tt)$y Even better would be to create an irregular time series object. library(zoo) tt - as.Date(paste(yr, 1, 1, sep = -)) + doy - 1 ta.z - na.approx(zoo(dat, tt)) Now ta.z is a zoo object representing your time series. coredata(ta.z) is the data and time(ta.z) are the dates. See: library(zoo) vignette(zoo) for more info. On 10/6/05, Anette Nørgaard [EMAIL PROTECTED] wrote: Can anybody help me write a code on the following data example, which fills out all NA values by using a linear interpolation with the two closest values? Doy is day of year (%j). Code example: yr-c(rep(2000,14)) doy-c(16:29) dat-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) ta-cbind(yr,doy,dat) ta yr doy dat [1,] 2000 16 3.2 [2,] 2000 17 NA [3,] 2000 18 NA [4,] 2000 19 NA [5,] 2000 20 NA [6,] 2000 21 NA [7,] 2000 22 NA [8,] 2000 23 5.1 [9,] 2000 24 NA [10,] 2000 25 NA [11,] 2000 26 NA [12,] 2000 27 NA [13,] 2000 28 NA [14,] 2000 29 4.6 Anette Norgaard [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
2d I can't manage to deal with factors, so when I have some, I
transform
them in vectors (with levels()), but I think I miss the power and
utility
of
the factor type ?
levels() is not the conversion you want.
in fact I use
'as.numeric(levels(f))[f]'
(from the ?factor description)
That lists all the levels, but
it doesn't tell you how they correspond to individual observations. For
example,
df - data.frame(x=1:3, y=c('a','b','a'))
df
x y
1 1 a
2 2 b
3 3 a
levels(df$y)
[1] a b
If you need to convert back to character values, use as.character():
as.character(df$y)
[1] a b a
got it.
1. You can't compare the levels of a factor unless you declared it to
be ordered:
df$y[1] df$y[2]
[1] NA
Warning message:
not meaningful for factors in: Ops.factor(df$y[1], df$y[2])
but
df$y - ordered(df$y)
df$y[1] df$y[2]
[1] FALSE
However, you need to watch out here: the comparison is done by the order
of the factors
I am sorry I don't understand this.
here you compare the position of a in the factor and the position of b in
the factor ?
, not an alphabetic comparison of their names:
levels(df$y) - c(before, after)
df
x y
1 1 before
2 2 after
3 3 before
df$y[1] df$y[2]
[1] FALSE
best regards,
florence.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
On 10/6/2005 10:20 AM, Florence Combes wrote:
2d I can't manage to deal with factors, so when I have some, I
transform
them in vectors (with levels()), but I think I miss the power and
utility
of
the factor type ?
levels() is not the conversion you want.
in fact I use
'as.numeric(levels(f))[f]'
(from the ?factor description)
That will only work if the levels have names that can be converted to
numbers. In the example below, the levels are a and b, so you'll
get NA values if you try this.
That lists all the levels, but
it doesn't tell you how they correspond to individual observations. For
example,
df - data.frame(x=1:3, y=c('a','b','a'))
df
x y
1 1 a
2 2 b
3 3 a
levels(df$y)
[1] a b
If you need to convert back to character values, use as.character():
as.character(df$y)
[1] a b a
got it.
1. You can't compare the levels of a factor unless you declared it to
be ordered:
df$y[1] df$y[2]
[1] NA
Warning message:
not meaningful for factors in: Ops.factor(df$y[1], df$y[2])
but
df$y - ordered(df$y)
df$y[1] df$y[2]
[1] FALSE
However, you need to watch out here: the comparison is done by the order
of the factors
I am sorry I don't understand this.
here you compare the position of a in the factor and the position of b in
the factor ?
It's the position of a in the levels() vector that is being compared.
I declared that the factor had ordered levels, and R interprets that
to mean that the first level is less than the second level, etc. This
is useful if you want to use meaningful names for ordered categories.
Comparison will be by the order of the categories, not by the name you
chose.
Duncan Murdoch
, not an alphabetic comparison of their names:
levels(df$y) - c(before, after)
df
x y
1 1 before
2 2 after
3 3 before
df$y[1] df$y[2]
[1] FALSE
best regards,
florence.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
a last question, and thanks a million for your patience and your
explanations ...
I tried with a df called merged and a column named Pcc_0h_A (which is
numeric values):
length(as.vector(merged$Pcc_0h_A))
[1] 12202
as.numeric(as.vector(merged$Pcc_0h_A)[1:10])
[1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA
ord-ordered(merged$Pcc_0h_A)
length(ord)
[1] 12202
ord[1:10]
[1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA
5386 Levels: 10.001 10.002 10.003 10.005 10.006 10.010 ...
9.999
here I have NA instead of NA because ord is a factor and the notation is
different ?
length(as.numeric(merged$Pcc_0h_A))
[1] 12202
as.numeric(merged$Pcc_0h_A[1:10])
[1] 1812 1547 3308 3114 1960 1370 NA NA NA NA
are these the levels names converted into numbers ? I don't think because
levels are like 10.001, 10.002 etc and 1812, 1547 etc are not in this form.
thanks a million
florence;
On 10/6/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
On 10/6/2005 10:20 AM, Florence Combes wrote:
2d I can't manage to deal with factors, so when I have some, I
transform
them in vectors (with levels()), but I think I miss the power and
utility
of
the factor type ?
levels() is not the conversion you want.
in fact I use
'as.numeric(levels(f))[f]'
(from the ?factor description)
That will only work if the levels have names that can be converted to
numbers. In the example below, the levels are a and b, so you'll
get NA values if you try this.
That lists all the levels, but
it doesn't tell you how they correspond to individual observations.
For
example,
df - data.frame(x=1:3, y=c('a','b','a'))
df
x y
1 1 a
2 2 b
3 3 a
levels(df$y)
[1] a b
If you need to convert back to character values, use as.character():
as.character(df$y)
[1] a b a
got it.
1. You can't compare the levels of a factor unless you declared it to
be ordered:
df$y[1] df$y[2]
[1] NA
Warning message:
not meaningful for factors in: Ops.factor(df$y[1], df$y[2])
but
df$y - ordered(df$y)
df$y[1] df$y[2]
[1] FALSE
However, you need to watch out here: the comparison is done by the
order
of the factors
I am sorry I don't understand this.
here you compare the position of a in the factor and the position of b
in
the factor ?
It's the position of a in the levels() vector that is being compared.
I declared that the factor had ordered levels, and R interprets that
to mean that the first level is less than the second level, etc. This
is useful if you want to use meaningful names for ordered categories.
Comparison will be by the order of the categories, not by the name you
chose.
Duncan Murdoch
, not an alphabetic comparison of their names:
levels(df$y) - c(before, after)
df
x y
1 1 before
2 2 after
3 3 before
df$y[1] df$y[2]
[1] FALSE
best regards,
florence.
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Interpolation in time
na.approx(zoo(ta[,-seq(2)], tt)) where tt is as before. On 10/6/05, Anette Nørgaard [EMAIL PROTECTED] wrote: This is exactly what I requested, thank you!! However I do actually have several columns in my data sheet where I need to do the same thing, then how do I come about that? e.g. yr-c(rep(2000,14)) doy-c(16:29) dat1-c(3.2,NA,NA,NA,NA,NA,NA,5.1,NA,NA,NA,NA,NA,4.6) dat2-c(2.2,NA,NA,NA,NA,NA,NA,6.1,NA,NA,NA,NA,NA,4.2) dat3-c(3.4,NA,NA,NA,NA,NA,NA,4.1,NA,NA,NA,NA,NA,4.7) ta-cbind(yr,doy,dat1,dat2,dat3) [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
On 10/6/2005 10:50 AM, Florence Combes wrote:
a last question, and thanks a million for your patience and your
explanations ...
I tried with a df called merged and a column named Pcc_0h_A (which is
numeric values):
length(as.vector(merged$Pcc_0h_A))
[1] 12202
as.numeric(as.vector(merged$Pcc_0h_A)[1:10])
[1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA
ord-ordered(merged$Pcc_0h_A)
length(ord)
[1] 12202
ord[1:10]
[1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA
5386 Levels: 10.001 10.002 10.003 10.005 10.006 10.010 ...
9.999
here I have NA instead of NA because ord is a factor and the notation is
different ?
I can't tell what's going on here. Since you are only showing me
converted values of each column (as.vector(), as.numeric(), ordered(),
etc.) I can't tell what the original looked like.
A useful way to get an overview of a dataframe is to look at the results
of three function calls:
head(merged)# list the first few rows
str(merged) # describe the structure of the dataframe
summary(merged) # summarize the data in each of the columns.
Duncan Murdoch
length(as.numeric(merged$Pcc_0h_A))
[1] 12202
as.numeric(merged$Pcc_0h_A[1:10])
[1] 1812 1547 3308 3114 1960 1370 NA NA NA NA
are these the levels names converted into numbers ? I don't think because
levels are like 10.001, 10.002 etc and 1812, 1547 etc are not in this form.
thanks a million
florence;
On 10/6/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
On 10/6/2005 10:20 AM, Florence Combes wrote:
2d I can't manage to deal with factors, so when I have some, I
transform
them in vectors (with levels()), but I think I miss the power and
utility
of
the factor type ?
levels() is not the conversion you want.
in fact I use
'as.numeric(levels(f))[f]'
(from the ?factor description)
That will only work if the levels have names that can be converted to
numbers. In the example below, the levels are a and b, so you'll
get NA values if you try this.
That lists all the levels, but
it doesn't tell you how they correspond to individual observations.
For
example,
df - data.frame(x=1:3, y=c('a','b','a'))
df
x y
1 1 a
2 2 b
3 3 a
levels(df$y)
[1] a b
If you need to convert back to character values, use as.character():
as.character(df$y)
[1] a b a
got it.
1. You can't compare the levels of a factor unless you declared it to
be ordered:
df$y[1] df$y[2]
[1] NA
Warning message:
not meaningful for factors in: Ops.factor(df$y[1], df$y[2])
but
df$y - ordered(df$y)
df$y[1] df$y[2]
[1] FALSE
However, you need to watch out here: the comparison is done by the
order
of the factors
I am sorry I don't understand this.
here you compare the position of a in the factor and the position of b
in
the factor ?
It's the position of a in the levels() vector that is being compared.
I declared that the factor had ordered levels, and R interprets that
to mean that the first level is less than the second level, etc. This
is useful if you want to use meaningful names for ordered categories.
Comparison will be by the order of the categories, not by the name you
chose.
Duncan Murdoch
, not an alphabetic comparison of their names:
levels(df$y) - c(before, after)
df
x y
1 1 before
2 2 after
3 3 before
df$y[1] df$y[2]
[1] FALSE
best regards,
florence.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] factor : how does it work ?
head(merged) ID Name Pcc_0h_A Pcc_0h_swapped_A 3302 301495 Q0010_01 |Q0010||Hypothetical ORF 12.276 11.716 6943 309175 Q0010_01 |Q0010||Hypothetical ORF 11.958 11.271 14065 298935 Q0017_01 |Q0017||Hypothetical ORF 14.098 13.122 6420 306615 Q0017_01 |Q0017||Hypothetical ORF 13.843 13.061 5066 296375 Q0032_01 |Q0032||Hypothetical ORF 12.451 11.467 12707 304055 Q0032_01 |Q0032||Hypothetical ORF 11.745 11.482 Pcc_0h_M Pcc_0h_swapped_M 3302 -0.249 0.316 6943 -0.115 0.780 14065 -0.053 0.263 6420 0.009 0.323 5066 0.015 0.687 12707 0.074 0.768 str(merged) `data.frame': 12202 obs. of 6 variables: $ ID : Factor w/ 12202 levels 295080,295081,..: 5076 11177 3046 9147 1009 7110 5136 11237 3106 9207 ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... $ Name : Factor w/ 6101 levels Q0010_01 ..,..: 1 1 2 2 3 3 4 4 5 5 ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... $ Pcc_0h_A : Factor w/ 5386 levels 10.001,10.002,..: 1812 1547 3308 3114 1960 1370 NA NA NA NA ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... $ Pcc_0h_swapped_A: Factor w/ 5082 levels 10.001,10.002,..: 1256 885 2533 2477 1051 1064 NA NA NA NA ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... $ Pcc_0h_M : Factor w/ 1940 levels 0.000, 0.001,..: 499 231 107 18 30 148 NA NA NA NA ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... $ Pcc_0h_swapped_M: Factor w/ 2343 levels 0.000, 0.001,..: 632 1453 526 646 1319 1434 NA NA NA NA ... ..- attr(*, names)= chr 3302 6943 14065 6420 ... a last question, and thanks a million for your patience and your explanations ... I tried with a df called merged and a column named Pcc_0h_A (which is numeric values): length(as.vector(merged$Pcc_0h_A)) [1] 12202 as.numeric(as.vector(merged$Pcc_0h_A)[1:10]) [1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA ord-ordered(merged$Pcc_0h_A) length(ord) [1] 12202 ord[1:10] [1] 12.276 11.958 14.098 13.843 12.451 11.745 NA NA NA NA 5386 Levels: 10.001 10.002 10.003 10.005 10.006 10.010 ... 9.999 here I have NA instead of NA because ord is a factor and the notation is different ? length(as.numeric(merged$Pcc_0h_A)) [1] 12202 as.numeric(merged$Pcc_0h_A[1:10]) [1] 1812 1547 3308 3114 1960 1370 NA NA NA NA are these the levels names converted into numbers ? I don't think because levels are like 10.001, 10.002 etc and 1812, 1547 etc are not in this form. with the str(merged) value I guess that 1812, 1547 etc are a sort of rank , am I right ? thanks a million florence; [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Use crossprod() instead of t(x) %*% x !
Marc == Marc Bernard [EMAIL PROTECTED] on Thu, 6 Oct 2005 15:11:09 +0200 (CEST) writes: ... Marc Here is the program: Marc nitems - 10 Marc x - array(rnorm(5*nitems,3,3), c(5,nitems)) Marc sigma - t(x)%*%x Marc inverse - try(solve(sigma), TRUE) . Just a side remark on your code above: You should learn about and use crossprod(x) rather than t(x) %*% x because of a 1) more efficient implementation 2) more accurate implementation The exact details depends on the (accelerated/optimized or not) version BLAS/Lapack your version of R is using and also on the kind of matrices. Further note, that for crossprod(t(X)) == X %*% t(X) LAPACK also provides a direc version to which the 'Matrix' package interfaces via function tcrossprod() which in particular also works *fast* for some of the sparse matrix classes. Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Simple question.....
Hi all user R,
My simple question is...I have a vector of names of predictors,
text-c(datem,cola,eslom)...I try to plot the model with this predictor
in sequence loop,
for(i in 1:3){
png(paste(fig_,i,sep=))
plot(preplot.gam(mod9)[[i]],se=T,rug=F,main=,xaxt=n,ylab=,xlab=)
axis(1,as.numeric(text[i]),as.character(text[i]),cex.axis=.9)
dev.off()
}
But the line with function axis get error
Error in axis(side, at, labels, tick, line, pos, outer, font, vfont, lty, :
no locations are finite
I put in shell text[i] give datem, I try to erase the , can not search what
is the function to erase this character (). Samebody can help me to erase ,
when put the predictor without , datem, there not problem
Thank for all
Fernando Espindola R.
Division Investigacion Pesquera
Instituto de Fomento Pesquero
Blanco 839
Valparaiso - CHILE
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] boxplot statistics
From: Graham Williams Received Tue 04 Oct 2005 7:26pm +1000 from Karin Lagesen: I have read and reread the boxplot and the boxplot stats page, and I still cannot understand how and what boxplot shows. I realize that this might be due to me not knowing enough statistics, but anyway... First, how does boxplot determine the size of the box? And is the line inside the box the mean or the median (or something completely different?) And how does it determine how long out the whiskers should go? Also, the boxplot.stats page talks about hinges, what are those? The two hinges are versions of the first and third quartile, i.e., close to 'quantile(x, c(1,3)/4)'. Wikipedia has a reasonable description http://en.wikipedia.org/wiki/Boxplot ... but not quite accurate. If I'm not mistaken, boxplots are based on Tukey's letter values. Here's one description of what they are: http://www.math.yorku.ca/SCS/Courses/eda/eda1.html#H2_32:1.2 Andy __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Simple question.....
See ?get.
Andy
From: Fernando Espíndola
Hi all user R,
My simple question is...I have a vector of names of predictors,
text-c(datem,cola,eslom)...I try to plot the model
with this predictor in sequence loop,
for(i in 1:3){
png(paste(fig_,i,sep=))
plot(preplot.gam(mod9)[[i]],se=T,rug=F,main=,xaxt=n,ylab=
,xlab=)
axis(1,as.numeric(text[i]),as.character(text[i]),cex.axis=.9)
dev.off()
}
But the line with function axis get error
Error in axis(side, at, labels, tick, line, pos, outer, font,
vfont, lty, :
no locations are finite
I put in shell text[i] give datem, I try to erase the ,
can not search what is the function to erase this character
(). Samebody can help me to erase , when put the
predictor without , datem, there not problem
Thank for all
Fernando Espindola R.
Division Investigacion Pesquera
Instituto de Fomento Pesquero
Blanco 839
Valparaiso - CHILE
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem with installing a package
I do have full access to that directory. I have the bbHist package in c:/PROGRA~1/R/rw2011/library/bbHist directory. Then under the library directory I did check and build. Here's what I got: $ R CMD check bbHist * checking for working latex ... OK * using log directory 'c:/progra~1/r/rw2011/library/bbHist.Rcheck' * using R version 2.1.1, 2005-06-20 * checking for file 'bbHist/DESCRIPTION' ... OK * this is package 'bbHist' version '0.1-1' * checking if this is a source package ... WARNING Subdirectory 'bbHist/src' contains object files. installing R.css in c:/progra~1/r/rw2011/library/bbHist.Rcheck -- Making package bbHist adding build stamp to DESCRIPTION making DLL ... ... DLL made installing DLL installing R files installing man source files installing indices installing help Building/Updating help pages for package 'bbHist' Formats: text html latex example chm bbHisttexthtml latex example chm Microsoft HTML Help Compiler 4.74.8702 Compiling c:\progra~1\r\rw2011\library\bbHist\chm\bbHist.chm Compile time: 0 minutes, 0 seconds 2 Topics 2 Local links 0 Internet links 1 Graphic Created c:\progra~1\r\rw2011\library\bbHist\chm\bbHist.chm, 21,741 bytes Compression increased file by 9,099 bytes. adding MD5 sums * DONE (bbHist) * checking package directory ... OK * checking for portable file names ... OK * checking DESCRIPTION meta-information ... OK * checking package dependencies ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for syntax errors ... OK * checking R files for library.dynam ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... WARNING Foreign function calls without 'PACKAGE' argument: .C(getBBData, ...) See section 'System and foreign language interfaces' of the 'Writing R Extensions' manual. * checking Rd files ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking for CRLF line endings in C/C++/Fortran sources/headers ... WARNING Found the following sources/headers with CRLF line endings: src/bbapi.h src/bbunix.h src/rbb.c Some Unix compilers require LF line endings. * creating bbHist-Ex.R ... OK * checking examples ... OK * creating bbHist-manual.tex ... OK * checking bbHist-manual.tex ... ERROR LaTeX errors when creating DVI version. This typically indicates Rd problems. $ R CMD build bbHist * checking for file 'bbHist/DESCRIPTION' ... OK * preparing 'bbHist': * checking DESCRIPTION meta-information ... OK * cleaning src * removing junk files * checking for LF line-endings in source files file 'bbHist/src/bbapi.h' had CRLF line endings file 'bbHist/src/bbunix.h' had CRLF line endings file 'bbHist/src/rbb.c' had CRLF line endings * checking for empty directories * building 'bbHist_0.1-1.tar.gz' And then when I do R CMD INSTALL bbHist, I get open(c:/progra~1/r/rw2011/library/bbHist/DESCRIPTION): No such file or directory The bbHist directory is removed and a 00LOCK directory is created. What I don't understand is if INSTALL removed the bbHist directory, why is it stilling looking for the DESCRIPTION file in that diretory? I'm really puzzled. If anyone has any idea, please let me know. Thanks. Date: Thu, 06 Oct 2005 08:47:23 +0200 From: Uwe Ligges [EMAIL PROTECTED] Subject: Re: [R] problem in installing a package To: Claire Lee [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=us-ascii; format=flowed Claire Lee wrote: I'm using R in Windows XP. I created a package myself. I've used R CMD check to check it. Everything seems OK except the latex. I get the error message: * checking bbHist-manual.tex ... ERROR LaTeX errors when creating DVI version. This typically indicates Rd problems. I ignored it because I didn't want to submit it to CRAN. Then I tried to use R CMD INSTALL to install it. First I get: mv: cannot move `c:/PROGRA~1/R/rw2011/library/bbHist' to `c:/PROGRA~1/R/rw2011/library/00LOCK/bbHist ': Permission denied and a bunch of making DLL errors. Then when I tried a second time, I get: open(c:/progra~1/r/rw2011/library/bbHist/DESCRIPTION): No such file or directory I can see a 00LOCK directory is created in the c:/PROGRA~1/R/rw2011/library directory. Any idea why this is happening? No, information is still too sparse, unfortunately. Do you have full write access? The 00LOCK directory is used to save the older package in order to be able to restore it if a new installtion fails. Something went wrong and you have to remove it manually now, I guess. Uwe Ligges Thanks. Claire __ R-help@stat.math.ethz.ch
Re: [R] problem with installing a package
Claire Lee wrote: I do have full access to that directory. I have the bbHist package in c:/PROGRA~1/R/rw2011/library/bbHist directory. Then under the library directory I did check and build. Here's what I got: Now I understand: You *must not* have the sources in that library! You want to install from some *other arbitrary* source directory *into* the library mentioned above. Uwe Ligges $ R CMD check bbHist * checking for working latex ... OK * using log directory 'c:/progra~1/r/rw2011/library/bbHist.Rcheck' * using R version 2.1.1, 2005-06-20 * checking for file 'bbHist/DESCRIPTION' ... OK * this is package 'bbHist' version '0.1-1' * checking if this is a source package ... WARNING Subdirectory 'bbHist/src' contains object files. installing R.css in c:/progra~1/r/rw2011/library/bbHist.Rcheck -- Making package bbHist adding build stamp to DESCRIPTION making DLL ... ... DLL made installing DLL installing R files installing man source files installing indices installing help Building/Updating help pages for package 'bbHist' Formats: text html latex example chm bbHisttexthtml latex example chm Microsoft HTML Help Compiler 4.74.8702 Compiling c:\progra~1\r\rw2011\library\bbHist\chm\bbHist.chm Compile time: 0 minutes, 0 seconds 2 Topics 2 Local links 0 Internet links 1 Graphic Created c:\progra~1\r\rw2011\library\bbHist\chm\bbHist.chm, 21,741 bytes Compression increased file by 9,099 bytes. adding MD5 sums * DONE (bbHist) * checking package directory ... OK * checking for portable file names ... OK * checking DESCRIPTION meta-information ... OK * checking package dependencies ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for syntax errors ... OK * checking R files for library.dynam ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... WARNING Foreign function calls without 'PACKAGE' argument: .C(getBBData, ...) See section 'System and foreign language interfaces' of the 'Writing R Extensions' manual. * checking Rd files ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking for CRLF line endings in C/C++/Fortran sources/headers ... WARNING Found the following sources/headers with CRLF line endings: src/bbapi.h src/bbunix.h src/rbb.c Some Unix compilers require LF line endings. * creating bbHist-Ex.R ... OK * checking examples ... OK * creating bbHist-manual.tex ... OK * checking bbHist-manual.tex ... ERROR LaTeX errors when creating DVI version. This typically indicates Rd problems. $ R CMD build bbHist * checking for file 'bbHist/DESCRIPTION' ... OK * preparing 'bbHist': * checking DESCRIPTION meta-information ... OK * cleaning src * removing junk files * checking for LF line-endings in source files file 'bbHist/src/bbapi.h' had CRLF line endings file 'bbHist/src/bbunix.h' had CRLF line endings file 'bbHist/src/rbb.c' had CRLF line endings * checking for empty directories * building 'bbHist_0.1-1.tar.gz' And then when I do R CMD INSTALL bbHist, I get open(c:/progra~1/r/rw2011/library/bbHist/DESCRIPTION): No such file or directory The bbHist directory is removed and a 00LOCK directory is created. What I don't understand is if INSTALL removed the bbHist directory, why is it stilling looking for the DESCRIPTION file in that diretory? I'm really puzzled. If anyone has any idea, please let me know. Thanks. Date: Thu, 06 Oct 2005 08:47:23 +0200 From: Uwe Ligges [EMAIL PROTECTED] Subject: Re: [R] problem in installing a package To: Claire Lee [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=us-ascii; format=flowed Claire Lee wrote: I'm using R in Windows XP. I created a package myself. I've used R CMD check to check it. Everything seems OK except the latex. I get the error message: * checking bbHist-manual.tex ... ERROR LaTeX errors when creating DVI version. This typically indicates Rd problems. I ignored it because I didn't want to submit it to CRAN. Then I tried to use R CMD INSTALL to install it. First I get: mv: cannot move `c:/PROGRA~1/R/rw2011/library/bbHist' to `c:/PROGRA~1/R/rw2011/library/00LOCK/bbHist ': Permission denied and a bunch of making DLL errors. Then when I tried a second time, I get: open(c:/progra~1/r/rw2011/library/bbHist/DESCRIPTION): No such file or directory I can see a 00LOCK directory is created in the c:/PROGRA~1/R/rw2011/library directory. Any idea why this is happening? No, information is still too sparse, unfortunately. Do you have full write access? The 00LOCK directory is used to save the
Re: [R] Compare two distance matrices
bady == bady [EMAIL PROTECTED]
on Thu, 06 Oct 2005 14:39:27 +0200 writes:
bady Hi, hi all,
I am trying to compare two distance matrices with R. I would like to
create a XY plot of these matrices and do some linear regression on
it. But, I am a bit new to R, so i have a few questions (I searched in
the documentation with no success).
The first problem is loading a distance matrix into R. This matrix is
the output of a the Phylip program Protdist and lookes like this:
I tried with the scan() function to load the files, but with no
success. How should i load in these files?
bady you can separately load each matrix with two text files.
bady require(ade4)
bady mat1 - read.table(mat1.txt)
bady nam1 - mat1[,1]
bady mat1 - mat1[,-1]
bady row.names(mat1) - names(mat1) - nam1
bady mat2 - read.table(mat2.txt)
bady nam2 - mat2[,1]
bady mat2 - mat2[,-1]
bady row.names(mat2) - names(mat2) - nam2
bady dist1 - mat2dist(mat1)
bady dist2 - mat2dist(mat2)
but I don't see why you would need an extra package ade4 and
its extra - function mat2dist().
when the 'stats' package already provides the function
as.dist(.) {the help page of which was mentioned by the
original poster}.
Here is a reproducible example showing how I think as.dist()
works sufficiently:
(m - toeplitz(round(rnorm(6),2)))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.42 -0.78 -0.42 -2.24 0.74 1.31
[2,] -0.78 -0.42 -0.78 -0.42 -2.24 0.74
[3,] -0.42 -0.78 -0.42 -0.78 -0.42 -2.24
[4,] -2.24 -0.42 -0.78 -0.42 -0.78 -0.42
[5,] 0.74 -2.24 -0.42 -0.78 -0.42 -0.78
[6,] 1.31 0.74 -2.24 -0.42 -0.78 -0.42
as.dist(m)
1 2 3 4 5
2 -0.78
3 -0.42 -0.78
4 -2.24 -0.42 -0.78
5 0.74 -2.24 -0.42 -0.78
6 1.31 0.74 -2.24 -0.42 -0.78
## it also works for data frames {if really needed}:
dm - as.data.frame(m)
as.dist(dm)
1 2 3 4 5
2 -0.78
3 -0.42 -0.78
4 -2.24 -0.42 -0.78
5 0.74 -2.24 -0.42 -0.78
6 1.31 0.74 -2.24 -0.42 -0.78
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] isdst
Can someone, please, explain the difference is results below (notice
the isdst value)
unlist(as.POSIXlt('2005-7-1'))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 1
unlist(as.POSIXlt(as.Date('2005-7-1')))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 0
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Testing strata by covariate interactions in coxph
On Thu, 6 Oct 2005, Pyy-Martikainen Marjo wrote: I specify a model with strata by covariate interactions. I would like to conduct a Wald test for the null hypothesis no differences between any covariate effects in the 3 groups. The survey package has a function regTermTest (which isn't specific to survey models -- it works on anything with vcov() and coef() methods). -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] isdst
See the Avoiding Errors section of the Help Desk article
in R News 4/1.
On 10/6/05, Omar Lakkis [EMAIL PROTECTED] wrote:
Can someone, please, explain the difference is results below (notice
the isdst value)
unlist(as.POSIXlt('2005-7-1'))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 1
unlist(as.POSIXlt(as.Date('2005-7-1')))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 0
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] blank graphics window while R is working
Dear R-help, When R starts to execute some code, all the graphics windows associated with that R session go blank, and remain blank until the R prompt returns. I am using the ion window manager (http://modeemi.cs.tut.fi/~tuomov/ion/) under Debian linux. I notice that under the gnome desktop environment with Fedora Core 4 linux, the graphics windows retain their images while R evaluates. Does anyone else experience this problem or know of any solution, or is this just an incompatibility between R and my (presumably somewhat uncommon) window manager? Thanks a lot, Dan p.s. the problem remains whether R is run via ESS or not. version _ platform i386-pc-linux-gnu arch i386 os linux-gnu system i386, linux-gnu status major2 minor1.0 year 2005 month04 day 18 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] isdst
On Thu, 6 Oct 2005, Omar Lakkis wrote:
Can someone, please, explain the difference is results below (notice
the isdst value)
Timezones. (That is carefully documented on the help page for
as.POSIXlt.)
unlist(as.POSIXlt('2005-7-1'))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 1
This is midnight on 2005-07-01 in your timezone, which presumably does
have DST and it was in force then.
unlist(as.POSIXlt(as.Date('2005-7-1')))
sec min hour mday mon year wday yday isdst
0 0 0 1 6 105 5 181 0
This is midnight on 2005-07-01 in UTC (which does not have DST).
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] blank graphics window while R is working
Your R is an obselete version. There are some changes in 2.2.0 which affect the hinting given to X11 WMs, so please try the current version of R. (I know it had only recently been released, but the posting guide refers you to the development versions.) I don't know if this will solve it (I have never met the problem) but do know that we will not patch obselete versions. On Thu, 6 Oct 2005, Dan Davison wrote: Dear R-help, When R starts to execute some code, all the graphics windows associated with that R session go blank, and remain blank until the R prompt returns. I am using the ion window manager (http://modeemi.cs.tut.fi/~tuomov/ion/) under Debian linux. I notice that under the gnome desktop environment with Fedora Core 4 linux, the graphics windows retain their images while R evaluates. Does anyone else experience this problem or know of any solution, or is this just an incompatibility between R and my (presumably somewhat uncommon) window manager? Thanks a lot, Dan p.s. the problem remains whether R is run via ESS or not. version _ platform i386-pc-linux-gnu arch i386 os linux-gnu system i386, linux-gnu status major2 minor1.0 year 2005 month04 day 18 language R -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] R Cocoa GUI syntax color bug?
I am running R 2.1.1 and R Cocoa GUI 1.12(1622) on a Mac Dual G5 with OS 10.4.2 Anytime I try to use the syntax color preference for the built-in editor, R crashes within a few minutes of working on a file. What information can I provide that might facilitate investigation of this? Thank you, Hank Stevens Dr. Martin Henry H. Stevens, Assistant Professor 338 Pearson Hall Botany Department Miami University Oxford, OH 45056 Office: (513) 529-4206 Lab: (513) 529-4262 FAX: (513) 529-4243 http://www.cas.muohio.edu/~stevenmh/ http://www.muohio.edu/ecology/ http://www.muohio.edu/botany/ E Pluribus Unum __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] boxplot statistics
A related comment - don't rely (too much) on boxplots. They show only a few things, which may be limiting in many cases and completely misleading in others. Here are a couple of suggestions for plots which you may find more useful than the standard box plots: - figure 3.27 from http://www.stat.auckland.ac.nz/~paul/RGraphics/chapter3.html - violin plots (see package vioplot) - density plots - histograms - box-percentile plots (bpplot from Hmisc) - quantile plots - if comparing 2 distributions, qq plots, quantile-difference plots, mean-difference plots etc. -Original Message- From: Karin Lagesen [mailto:[EMAIL PROTECTED] Sent: Tuesday, October 04, 2005 5:24 AM To: [EMAIL PROTECTED] Subject: [R] boxplot statistics I have read and reread the boxplot and the boxplot stats page, and I still cannot understand how and what boxplot shows. I realize that this might be due to me not knowing enough statistics, but anyway... First, how does boxplot determine the size of the box? And is the line inside the box the mean or the median (or something completely different?) And how does it determine how long out the whiskers should go? Also, the boxplot.stats page talks about hinges, what are those? The two hinges are versions of the first and third quartile, i.e., close to 'quantile(x, c(1,3)/4)'. Thankyou very much. Karin -- Karin Lagesen, PhD student [EMAIL PROTECTED] http://www.cmbn.no/rognes/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] circular statistics plotting
Hi all, I'm new to the list here, and I have what I think is a simple question. Using the circular package, is there a way to plot the mean and variance on top of a rose diagram or other plot of the data? Thanks in advance... - Jason Jason Horn Boston University Department of Biology 5 Cumington Street Boston, MA 02215 [EMAIL PROTECTED] office: 617 353 6987 cell: 401 588 2766 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Compare two distance matrices
You are right,
in the help page of mat2dist, dist2mat we found:
mat2dist and dist2mat are local copies of as.dist and as.matrix.dist of mva.
now these functions are in stats. We have to check it.
Thanks.
Quoting Martin Maechler [EMAIL PROTECTED]:
bady == bady [EMAIL PROTECTED]
on Thu, 06 Oct 2005 14:39:27 +0200 writes:
bady Hi, hi all,
I am trying to compare two distance matrices with R. I would like to
create a XY plot of these matrices and do some linear regression on
it. But, I am a bit new to R, so i have a few questions (I searched in
the documentation with no success).
The first problem is loading a distance matrix into R. This matrix is
the output of a the Phylip program Protdist and lookes like this:
I tried with the scan() function to load the files, but with no
success. How should i load in these files?
bady you can separately load each matrix with two text files.
bady require(ade4)
bady mat1 - read.table(mat1.txt)
bady nam1 - mat1[,1]
bady mat1 - mat1[,-1]
bady row.names(mat1) - names(mat1) - nam1
bady mat2 - read.table(mat2.txt)
bady nam2 - mat2[,1]
bady mat2 - mat2[,-1]
bady row.names(mat2) - names(mat2) - nam2
bady dist1 - mat2dist(mat1)
bady dist2 - mat2dist(mat2)
but I don't see why you would need an extra package ade4 and
its extra - function mat2dist().
when the 'stats' package already provides the function
as.dist(.) {the help page of which was mentioned by the
original poster}.
Here is a reproducible example showing how I think as.dist()
works sufficiently:
(m - toeplitz(round(rnorm(6),2)))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.42 -0.78 -0.42 -2.24 0.74 1.31
[2,] -0.78 -0.42 -0.78 -0.42 -2.24 0.74
[3,] -0.42 -0.78 -0.42 -0.78 -0.42 -2.24
[4,] -2.24 -0.42 -0.78 -0.42 -0.78 -0.42
[5,] 0.74 -2.24 -0.42 -0.78 -0.42 -0.78
[6,] 1.31 0.74 -2.24 -0.42 -0.78 -0.42
as.dist(m)
1 2 3 4 5
2 -0.78
3 -0.42 -0.78
4 -2.24 -0.42 -0.78
5 0.74 -2.24 -0.42 -0.78
6 1.31 0.74 -2.24 -0.42 -0.78
## it also works for data frames {if really needed}:
dm - as.data.frame(m)
as.dist(dm)
1 2 3 4 5
2 -0.78
3 -0.42 -0.78
4 -2.24 -0.42 -0.78
5 0.74 -2.24 -0.42 -0.78
6 1.31 0.74 -2.24 -0.42 -0.78
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Stephane DRAY
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] [BioC] (no subject)
I got the 'hgu133a' for R 2.0.0, would it be useful ? On 10/6/05, Francesca Buffa [EMAIL PROTECTED] wrote: Dear all, does anybody have, or know where could I find, a version of the hgu133b library compatible with R 2.0.1? (I need it sooner than I can update to the new R version) thank you so much in advance francesca ___ Bioconductor mailing list [EMAIL PROTECTED] https://stat.ethz.ch/mailman/listinfo/bioconductor [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] a question about LMS and what constitutes outliers
Hi, I have been using the lqs function with method='lms'. However the results I get are a little different from the results noted by Rousseeuw Leroy (Robust Regression and Outlier Detection) and I was wondering how to use these results for outlier detection. I'm using the stackloss dataset, for which the original Rousseeuw et al. program points out that observations 1,2,3,4 and 21 are outliers. This conslusion is arrived at by testing whether the residual is greater than 2.5 * standard error Netx I ran lqs as: m - lqs(stackloss[,-4], stackloss[,4], method='lms', control=list (psamp=4, nsamp='exact', adjust=TRUE)) (I ran it exhaustively since that was how I ran the original program from Rousseeuw) The coefficients obtained from lqs() are more or less identical to that obtained by the original program. However the scale estimates do not match. I assume that this would be becuase of the per sample adjustments. Now if I want to decide whether an observation is an outlier I use the condition which( abs(m$resid) 2.5 * m$scale[1] ) and this gives me 1 2 3 4 8 13 14 20 21 1 2 3 4 8 13 14 20 21 Now, it includes the original outliers as noted by Rousseuw, but also 4 extra ones. From a plot of the residuals I can see obs 13,14,20 possibly being regarded as outliers but 8 seems a stretch. I tried evaluating the above condition with m$scale[2] but I get the same result. I also tried running lqs() with adjust=FALSE in which case using the above condition obs 1,2,3,4,13,20,21 are regarded as outliers. So my questions are 1) Am I correct in using the above condition to determine whether an observation is an outlier? 2) If so, is it correct that lqs() will detect more outliers than noted by the original book/program? Thanks, --- Rajarshi Guha [EMAIL PROTECTED] http://jijo.cjb.net GPG Fingerprint: 0CCA 8EE2 2EEB 25E2 AB04 06F7 1BB9 E634 9B87 56EE --- After an instrument has been assembled, extra components will be found on the bench. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] data.frame error using sem package
I keep getting this error when I try to use the sem package. I and another person who has successfully used the sem package for similar analysis (fMRI effective connectivity) cannot figure out what is wrong with my code. I would appreciate any suggestions. The error message: Error in data.frame(object$coeff, se, z, 2 * (1 - pnorm(abs(z))), par.code) : arguments imply differing number of rows: 6, 0 In addition: Warning message: Could not compute QR decomposition of Hessian. Optimization probably did not converge. in: sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars, Thank you, Suzanne __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] data.frame error using sem package
Suzanne Witt wrote: I keep getting this error when I try to use the sem package. I and another person who has successfully used the sem package for similar analysis (fMRI effective connectivity) cannot figure out what is wrong with my code. I would appreciate any suggestions. It is almost impossible to help if you do not specify a toy example that shows how you produced that error message. Please read the psoting guide which tells you how to specify such examples most appropriately in order to get a good answer. Uwe Ligges The error message: Error in data.frame(object$coeff, se, z, 2 * (1 - pnorm(abs(z))), par.code) : arguments imply differing number of rows: 6, 0 In addition: Warning message: Could not compute QR decomposition of Hessian. Optimization probably did not converge. in: sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars, Thank you, Suzanne __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] data.frame error using sem package
I am trying to use sem to measure effective connectivity among four
brain regions. I have pasted the code that I am trying to run since
that seems easier than trying to come up with another example.
The input data is time series data taken from SPM; they are each
1x121 columns of numbers. I get the error either when I source the
whole code, or if I enter it line by line when I go to get the summary.
Thanks,
Suzanne
library(sem)
# Load the region timecourses.
lsma1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_LSMA.dat)
rsma1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_RSMA.dat)
lmc1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_LM1.dat)
rmc1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_RM1.dat)
# Combine all the timecourses from each session into a single data
frame and name the columns appropriately.
lcomf - cbind(lsma1, rsma1, lmc1, rmc1)
names(lcomf) - c(LSMA, RSMA, LM1, RM1)
# Type this at the command line to see a summary of your data
str(lcomf)
# Set up the structural equation model.
p.ram - matrix(c( 'LM1 - LSMA', 'LM1 - LSMA', NA,
'LSMA - RSMA', 'LSMA - RSMA', NA,
'RSMA - RM1', 'RSMA - RM1', NA,
'LSMA - LSMA', 'LSMA - LSMA', NA,
'RSMA - RSMA', 'RSMA - RSMA', NA,
'RM1 - RM1', 'RM1 - RM1', NA),
ncol = 3, byrow = TRUE)
# Tell which variables are exogenous ('fixed').
p.fixed - c('LM1')
# Do the fitting for session 1.
C - cor(lcomf)
nt - dim(lcomf)[1]
#attach(lcomf)
lcomf.results - sem(p.ram, C, nt, obs.variables = rownames(C),
fixed.x = p.fixed)
# Check out the results using the summary function
summary(lcomf.results)
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] data.frame error using sem package
Dear Suzanne,
Take a look at your model specification:
p.ram
[,1][,2][,3]
[1,] LM1 - LSMA LM1 - LSMA NA
[2,] LSMA - RSMA LSMA - RSMA NA
[3,] RSMA - RM1 RSMA - RM1 NA
[4,] LSMA - LSMA LSMA - LSMA NA
[5,] RSMA - RSMA RSMA - RSMA NA
[6,] RM1 - RM1 RM1 - RM1 NA
This matrix should have three columns, the first giving the path, the second
the name of the corresponding parameter, and the third the start value for
the parameter (or NA if you want sem() to compute a start value). You've
apparently left out the parameter names. Please see the sem examples for
details and the paper at
http://socserv.socsci.mcmaster.ca/jfox/sem-package.pdf.
I hope this helps,
John
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Suzanne Witt
Sent: Thursday, October 06, 2005 2:55 PM
To: r-help@stat.math.ethz.ch
Subject: [R] data.frame error using sem package
I am trying to use sem to measure effective connectivity
among four brain regions. I have pasted the code that I am
trying to run since that seems easier than trying to come up
with another example.
The input data is time series data taken from SPM; they are each
1x121 columns of numbers. I get the error either when I
source the whole code, or if I enter it line by line when I
go to get the summary.
Thanks,
Suzanne
library(sem)
# Load the region timecourses.
lsma1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_LSMA.dat)
rsma1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_RSMA.dat)
lmc1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_LM1.dat)
rmc1 - read.table(/Users/witt/parkinsons/rmrkm010905/R_files/
010905_lcomf_RM1.dat)
# Combine all the timecourses from each session into a single
data frame and name the columns appropriately.
lcomf - cbind(lsma1, rsma1, lmc1, rmc1)
names(lcomf) - c(LSMA, RSMA, LM1, RM1)
# Type this at the command line to see a summary of your data
str(lcomf)
# Set up the structural equation model.
p.ram - matrix(c( 'LM1 - LSMA', 'LM1 - LSMA', NA,
'LSMA - RSMA', 'LSMA - RSMA', NA,
'RSMA - RM1', 'RSMA - RM1', NA,
'LSMA - LSMA', 'LSMA - LSMA', NA,
'RSMA - RSMA', 'RSMA - RSMA', NA,
'RM1 - RM1', 'RM1 - RM1', NA),
ncol = 3, byrow = TRUE)
# Tell which variables are exogenous ('fixed').
p.fixed - c('LM1')
# Do the fitting for session 1.
C - cor(lcomf)
nt - dim(lcomf)[1]
#attach(lcomf)
lcomf.results - sem(p.ram, C, nt, obs.variables =
rownames(C), fixed.x = p.fixed)
# Check out the results using the summary function
summary(lcomf.results)
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] circular statistics plotting
Jason Horn wrote: Hi all, I'm new to the list here, and I have what I think is a simple question. Using the circular package, is there a way to plot the mean and variance on top of a rose diagram or other plot of the data? The package seems to be very well designed: I did not know anything about this package at the time reading your question. Reading the help page and trying out the first example, it took me less than a minute to figure it out (hence I wonder why you haven't?). From the example: library(circular) x - circular(runif(50, 0, 2*pi)) rose.diag(x, bins = 18, main = 'Uniform Data') ## plot the mean: points(mean(x), col = red) Uwe Ligges Thanks in advance... - Jason Jason Horn Boston University Department of Biology 5 Cumington Street Boston, MA 02215 [EMAIL PROTECTED] office: 617 353 6987 cell: 401 588 2766 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Replicated goodness-of-fit test
Dear all, I am facing a problem that seems to me more tricky now that it did at first sight. I have a collection of data which consist in frequency distributions: 6 patches had been proposed to female insects (for oviposition), 3 of them corresponding to one treatment (A), the other 3 to another treatment (B). The results are the number of patches that have been chosen for oviposition by each female (0 to 3 patches can be used for oviposition, for each treatment). Thus, the experiment is replicated (each female's choice is a replication). I am interested in testing if females do prefer one treatment to the other (if they oviposit preferently on one sort of patch), i.e. testing the goodness of fit of my data to an equal ratio of oviposition on patches A and B. I intended to use the replicated goodness-of-fit test (G-statistic) described in Sokal Rohlf (1981): this test seemed accurate for such data... But a first problem is that some frequencies in my data set are = 0... I wonder if an x-( x + 1) transformation could be reasonnably performed to permit the taking of ln(X) in calculations. Moreover, low frequencies (less than 5) are usually reported as a problem for G-tests, and because my females had only the choice to oviposit on 6 patches (3 for treatment A, 3 for treatment B), I have indeed low frequencies! Do you know another way to perform such an analysis, more accurate in my case (or do you thing a G-test could nevertheless be used)? I hope my explanations were comprehensible, and if it isn't the case, don't hesitate to ask for more. Thanks for your reply Kind regards, Karine __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] boxplot statistics
Selon bogdan romocea [EMAIL PROTECTED]: A related comment - don't rely (too much) on boxplots. They show only a few things, which may be limiting in many cases and completely misleading in others. Here are a couple of suggestions for plots which you may find more useful than the standard box plots: - figure 3.27 from http://www.stat.auckland.ac.nz/~paul/RGraphics/chapter3.html - violin plots (see package vioplot) - density plots - histograms - box-percentile plots (bpplot from Hmisc) - quantile plots - if comparing 2 distributions, qq plots, quantile-difference plots, mean-difference plots etc. Hi, HDR (highest density regions) boxplots are interresting. See http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=102 Romain -Original Message- From: Karin Lagesen [mailto:[EMAIL PROTECTED] Sent: Tuesday, October 04, 2005 5:24 AM To: [EMAIL PROTECTED] Subject: [R] boxplot statistics I have read and reread the boxplot and the boxplot stats page, and I still cannot understand how and what boxplot shows. I realize that this might be due to me not knowing enough statistics, but anyway... First, how does boxplot determine the size of the box? And is the line inside the box the mean or the median (or something completely different?) And how does it determine how long out the whiskers should go? Also, the boxplot.stats page talks about hinges, what are those? The two hinges are versions of the first and third quartile, i.e., close to 'quantile(x, c(1,3)/4)'. Thankyou very much. Karin -- Karin Lagesen, PhD student [EMAIL PROTECTED] http://www.cmbn.no/rognes/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] error bars on xy plots and barplots and histograms
Dear all I have yet to find code for making and showing error bars. Is there some idea that the simple presentation of confidence intervals is not good statistics? Any leads appreciated. Chris Buddenhagen, Botany Department, Charles Darwin Research Station, Santa Cruz,Galapagos. Mail: Charles Darwin Foundation, Casilla 17-01-3891 Avenida 6 de Diciembre N36-109 y Pasaje California Quito, ECUADOR __ EL CONTENIDO DE ESTE MENSAJE ES DE ABSOLUTA RESPONSABILIDAD DEL AUTOR. FUNDACION CHARLES DARWIN WWW.DARWINFOUNDATION.ORG __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] error bars on xy plots and barplots and histograms
Chris, computation of confidence intervals is left to the user. Showing them can be achieved using the arrows() function. I hope that this helps, Andrew On Thu, Oct 06, 2005 at 08:19:27AM -0600, Chris Buddenhagen wrote: Dear all I have yet to find code for making and showing error bars. Is there some idea that the simple presentation of confidence intervals is not good statistics? Any leads appreciated. Chris Buddenhagen, Botany Department, Charles Darwin Research Station, Santa Cruz,Galapagos. Mail: Charles Darwin Foundation, Casilla 17-01-3891 Avenida 6 de Diciembre N36-109 y Pasaje California Quito, ECUADOR __ EL CONTENIDO DE ESTE MENSAJE ES DE ABSOLUTA RESPONSABILIDAD DEL AUTOR. FUNDACION CHARLES DARWIN WWW.DARWINFOUNDATION.ORG __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Andrew Robinson Senior Lecturer in Statistics Tel: +61-3-8344-9763 Department of Mathematics and StatisticsFax: +61-3-8344-4599 University of Melbourne, VIC 3010 Australia Email: [EMAIL PROTECTED]Website: http://www.ms.unimelb.edu.au __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Different sings for correlations in OLS and TSA
I haven't seen a reply to this, so I will offer a comment:
It looks to me like the correlation you question is the correlation
between the estimated intercept and slope of the gls regression line.
This is different from the correlation between ts.mar and ts.anr. To
understand this, consider the following addition to your script:
plot(ts.mar, ts.anr)
# The intercept and slope are negative
# because the mean of ts.mar is positive.
fit3a - gls(ts.anr ~ I(ts.mar-1000),correlation =
corARMA(value=c(mod3$coef[1],mod3$coef[2]),p=2))
summary(fit3a)
snip
Correlation:
(Intr)
I(ts.mar - 1000) 0.988
# The estimated slope is the same
# in fit3 and fit3a
# but the sign of the correlation between
# intercept and slope has changed.
I hope this response still interests you, even though it comes well
over a month after your post. You are to commended for providing a
good, reasonably complete example. If it had been shorter, it might
have received more comments sooner.
If this does not answer your question or you have another, please let
us know.
spencer graves
Michael Tiemann wrote:
Dear list,
I am trying to re-analyse something. I do have two time series, one
of which (ts.mar) might help explaining the other (ts.anr). In the
original analysis, no-one seems to have cared about the data being
time-series and they just did OLS. This yielded a strong positive
correlation.
I want to know if this correlation is still as strong when the
autocorrelations are taken into account. There are autocorrelations, so
I model the data with arima() to get the parameters and fit it with
gls(). So far, the code seems to work fine, but what puzzles me is that
I get different sings: the gls-fit yields a strong negative correlation.
This shouldn't be so, so I suspect I am doing something wrong.
Here is my code:
# this is my data
ts.mar-ts(c(431.3,438,389.7,353.3,354.6,371.8,397.7,438.5,467.9,505.7,574.7,644.7,667.8,616.4,509.6,447,413.1,384.1),start=1980,freq=1)
ts.anr-ts(c(104.1,102.4,97.9,96.2,95.1,95.1,97.9,101.6,105.9,111.1,117.9,121.3,121.8,114.2,107.6,105.1,101.9,98.6),start=1980,freq=1)
# to find autocorrelations via (p)acf's and mle I do:
fun.tsa.mle-function(x){
par(mfrow=c(3,1))
acf(x)
pacf(x)
# AR model is estimated
m1- ar.mle(x)
# An estimation of the unexplained portion of variance
m1.1-m1$var.pred
# plot the function
plot(x)
# Give a printout
print(m1)
print(unexplained portion of variance:)
print(m1.1)
print(Mean:)
print(m1$x.mean)
par(mfrow=c(1,1))
}
#now, the autocorrelations should be consistent with following processes:
fun.tsa.mle(ts.mar) #following DAAG a p=2 AR
fun.tsa.mle(ts.anr) #following DAAG a p=2 AR
#I need to know, wether ts.anr can be explained with ts.mar, so
#according to ar.mle:
mod3-arima(ts.anr,order=c(2,0,0),xreg=ts.mar,transform.pars=TRUE)
fit3 - gls(ts.anr ~ ts.mar,correlation =
corARMA(value=c(mod3$coef[1],mod3$coef[2]),p=2))
summary(fit3)
ts.plot(ts.anr,fit3$fitted,col=1:2)
#the puzzling bit is the negative correlation. It ought to be positive,
I think.
#a simple OLS (this is what the people before me have done) yields
test3-ols(ts.anr~ts.mar)
test3 #with a positive correlation. Why?
Where is the mistake? Up to now, I just thought time-series analyses
would correct parameters and estimations, but simply changing signs?
Appreciating your help and suggestions,
Michael.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA
[EMAIL PROTECTED]
www.pdf.com http://www.pdf.com
Tel: 408-938-4420
Fax: 408-280-7915
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Inverse autocorrelation function?
Dear Dr. Hartley: RSiteSearch(Inverse autocorrelation function) produced 11 hits, none of which seemed to relate to your question. If you have a reasonable algorithm for computing the IACF, it might not be difficult to program in R. If you have compiled code, e.g., C++ or Fortran, it might not be difficult to link to it. Googling for inverse autocorrelation function led me to an article On Autoregressive Model Identification Ette Harrison Etuk (www.jos.nu/Articles/abstract.asp?article=42113) which compares IACF and PACF, concluding (a) neither is consistently more powerful than the other, but On the whole the partial autocorrelation function exhibits better performance. Conclusions: (1) R does not seem to have an IACF function, unless the IACF is identical to something else that R has under a different name. (2) The R Project for Stastical Computing is NOT a completed product but a project perpetually under renewal and extension. As such, it has many contributed packages and is happy to accept more. Spencer Graves David Hartley wrote: In time series analysis it is helpful to plot the autocorrelation function (ACF), partial autocorrelation function (PACF), and the inverse autocorrelation function (IACF). The stats library provides the ability to compute and plot the ACF and PACF, but I cannot find an [R] procedure to compute and plot the IACF. Is there one? Best regards, David Hartley, PhD __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Error in integrate
I'm using R 2.0.1 under Windows XP. I get the following message when I
run the code listed below. I don't seem to have any problems using the
function slice outside integrate.
Error in integrate(slice, 0, Inf, , , , , , , a, b) * delta :
non-numeric argument to binary operator
[ By the way, I'm trying to evaluate a two-dimensional integral by
slicing it up into one-dimensional bits which I will loop over to evaluate.]
Here's the code:
mu - 0.3
sig2 - 0.01
alpha - 20
beta - 15.75
rho - 0.1
k - 1/(rho^2.5*gamma(rho)^2*sqrt(2*pi*sig2))
slice - function(w,a,b)
{
g - w^(1/rho)
g1 - w1^(1/rho)
h - g1^a*g^b
E - -(y-rho*mu -g1/alpha + g/beta)^2/(2*sig2*rho)
k*h*exp(E- g1 - g)
}
hi - 10
delta - 0.05
grid - seq(delta/2,hi,delta)
y - -0.3
a - 0
b - 0
m - length(grid)
A - rep(0,m)
j - 0
for (w1 in grid)
{
j - j+1
A[j] - integrate(slice,0,Inf,,,a,b)*delta
cat(A[j],\n)
}
--
Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: [EMAIL PROTECTED]Fax 7 838 4155
Phone +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Error in integrate
Murray Jorgensen [EMAIL PROTECTED] wrote: I'm using R 2.0.1 under Windows XP. I get the following message when I run the code listed below. I don't seem to have any problems using the function slice outside integrate. Error in integrate(slice, 0, Inf, , , , , , , a, b) * delta : non-numeric argument to binary operator RTFM! Integrate returns: Value: A list of class 'integrate' with components : : You need: A[j] - integrate(slice,0,Inf,,,a,b)$value*delta Ray __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] error bars on xy plots and barplots and histograms
Andrew Robinson wrote: Chris, computation of confidence intervals is left to the user. Showing them can be achieved using the arrows() function. I hope that this helps, Andrew Also do ?xYplot after installing the Hmisc package. -Frank On Thu, Oct 06, 2005 at 08:19:27AM -0600, Chris Buddenhagen wrote: Dear all I have yet to find code for making and showing error bars. Is there some idea that the simple presentation of confidence intervals is not good statistics? Any leads appreciated. Chris Buddenhagen, Botany Department, Charles Darwin Research Station, Santa Cruz,Galapagos. Mail: Charles Darwin Foundation, Casilla 17-01-3891 Avenida 6 de Diciembre N36-109 y Pasaje California Quito, ECUADOR __ EL CONTENIDO DE ESTE MENSAJE ES DE ABSOLUTA RESPONSABILIDAD DEL AUTOR. FUNDACION CHARLES DARWIN WWW.DARWINFOUNDATION.ORG __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
