[R] More questions about gam and plot
When I use a gam command c - gam(depvar~var1+var2+s(var3)+s(var4)) Followed by a plot command plot(c,) The standard errors cross in the middle of the graph, (This variable is almost linearly related with the outcome variable) why? When I use the persp command after the gam fit, it does not graph the smooth terms but usually two of the other terms in the model, why? Dr. Tor A Strand Centre for International Health Haukeland Hospital University of Bergen 5021 Bergen Norway Phone: (country prefix 47) Residence:56 51 10 88, office: 55 97 49 80, fax: 55 97 49 79, cellular: 90 97 10 86 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Keeping track of occurrence of warning message
On Tue, Jul 15, 2003 at 02:58:10PM +0800, Tan Chuen Seng wrote:
Hi there,
I am interested if there is anyway to keep track of the occurrence of
warning message.
I know that warnings will only be printed out at the end of the program
if warn=0. However I am also interested at which particular set of data
does the warnings occur too. This is because I am running 1000 data, so
if there are 2 or 3 data that give warnings, I would like to know which
are the ones out of the 1000 data.
I tried using the following code in the program to indicate where the
warning occur but was unable to get anything recorded although the
warnings() gave me 12 messages.
track.warning-NULL
if(options(warn)$warn=0){
track.warning-c(track.warning,data.no)
}
Your help is greatly appreciated. Thanks.
From chuen seng
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
One way to do things is to have a list to store warnings has you
hit them in you loop: you can access what is the list last.warning.
Hopin' it helps,
L.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] R, geochemistry, ternary diagrams
Are there enough geochemists using R already that he'd find like-minded people to discuss technical issues with if he _did_ switch to R? Is there a package somewhere already that does ternary and other geochemistry diagrams? Another possibility for a ternary plot was mentioned by Prof Ripley in http://maths.newcastle.edu.au/~rking/R/help/02b/3637.html library(MASS) example(Skye) gives code and an example HTH, Tobias __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] bug?
Before everyone re-invents the wheel: All current versions of the S language (including R for several years) have the generic function all.equal() with a method for simple numeric vectors (but with many other methods too). The point is that you need a *combination* of relative and absolute difference. Look at all.equal.numeric() if you want to learn .. Martin Maechler [EMAIL PROTECTED] http://stat.ethz.ch/~maechler/ Seminar fuer Statistik, ETH-Zentrum LEO C16Leonhardstr. 27 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-1-632-3408 fax: ...-1228 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] dbApply (R newbee)
I am trying to use R interfaced with MySQL. Present goal is that R should calculate the 85% quantile of AvgSpeed for each LinieID. Looking through documentation of the RMySQL Package, I guessed that dbApply would do the trick due to this example ## compute quanitiles for each network agent con - dbConnect(MySQL(), group=vitalAnalysis) res - dbSendQuery(con, select Agent, ip_addr, DATA from pseudo_data order by Agent) out - dbApply(res, INDEX = Agent, FUN = function(x, grp) quantile(x$DATA, names=FALSE)) But when I try I get this: con = dbConnect(MySQL(),group=Speed) res - dbSendQuery(con, Select LinieID, AvgSpeed from speedlinietur order by LinieID, AvgSpeed) out - dbApply(res, INDEX = LinieID,FUN = function(x, grp) quantile(x$AvgSpeed, names=FALSE)) Error in mysqlDBApply(res, ...) : couldn't find function dbGetConnection I saw in an earlier posting that: There has been a change in function names in the new RMySQL 0.5-0 version (also in the new ROracle and RSQLite). The reason is that the R-SIG-DB has agreed on a common interface to all databases, and as part of this new common interface, most functions have been renamed. The following simple name mapping may help you upgrade existing code to the new DBI: pre-DBI DBI 0.1-4 .. dbGetConnection = getConnection .. when I do this I get this error message dbGetConnection = getConnection out - dbApply(res, INDEX = LinieID,FUN = function(x, grp) quantile(x$AvgSpeed, names=FALSE)) Error in if (what %in% set) structure(what, class = connection) else NULL : argument is of length zero Can anyone please help me Jesper Runge Madsen Ph.D. Trafikforskningsgruppen Institut for Samfundsudvikling og Planlægning Aalborg Universitet mailto:[EMAIL PROTECTED] [EMAIL PROTECTED] [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] More questions about gam and plot
When I use a gam command c - gam(depvar~var1+var2+s(var3)+s(var4)) Followed by a plot command plot(c,) The standard errors cross in the middle of the graph, (This variable is almost linearly related with the outcome variable) why? - Smooth terms are centred - which means that they are constrained to sum to zero over the covariate values. For a straight line this means that there must be one point that is known exactly. When I use the persp command after the gam fit, it does not graph the smooth terms but usually two of the other terms in the model, why? - By default, persp picks up the first two covariates in your model. If you want to plot against other variables you need to tell persp this... see ?persp.gam. best, Simon _ Simon Wood [EMAIL PROTECTED]www.stats.gla.ac.uk/~simon/ Department of Statistics, University of Glasgow, Glasgow, G12 8QQ Direct telephone: (0)141 330 4530 Fax: (0)141 330 4814 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] printf and friends in R?
Hi folks
Does R have anything straightforwardly resembling the commands
fprintf, sprintf, printf (derived from C, and present in octave
and matlab)?
As in printf(format_string, ... )
where format_string defines the print format (including any
fixed text) and ... is a list of variables whose values are
to be inserted into the line.
Example:
printf(Case %d#%.2f#%.2f%.4f\n, n,x1,x2,x3 )
which would output a line like
Case 10#3.21#7.65#0.4321
this particular case being a line in the right format for processing
by groff's 'tbl' table-formatter. In fact you could write the whole
table definition using printf:
printf(.TS\ntab(#);\nr n n n.\n
for(i in (1:nrow(X))){
printf(Case %d#%.2f#%.2f%.4f\n, n,X[i,1],X[i,2],X[i,3] )
}
printf(.TE\n)
Other variants could be used for similar purposes.
(If something like this is not already around in R, I suppose it
wouldn't be too difficult for me to write it; but it's worth asking
first!)
Thanks,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 11:32:39
-- XFMail --
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Hypothesis testing after optim
Hi Peter, Here are two links that may be useful: ML hypothesis tests: http://statgen.iop.kcl.ac.uk/bgim/mle/sslike_5.html ML confidence bounds: http://www.weibull.com/LifeDataWeb/likelihood_ratio_confidence_bounds.htm Yours, Gábor -- Gabor BORGULYA MD MSc Semmelweis University of Budapest, 2nd Dept of Paediatrics Hungarian Paediatric Cancer Registry phone: +36 - 1 - 4591500 / 2834 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Subsetting a matrix [Again!]
On Tuesday 15 July 2003 12:14, Ted Harding wrote: Hi Folks, People's suggestion of drop=FALSE seemed to do the trick (preserving matrix character when subestting to a row, i.e. creating 1xk matrix). However, I seem to have encountered a case where even this does not work: mu-c(1,2,3) mu-matrix(mu,nrow=1) mu [,1] [,2] [,3] [1,]123 iX1-c(T,F,F); iX2- !iX1 mu1-mu[iX1,drop=FALSE]; mu2-mu[iX2,drop=FALSE]; mu1 [1] 1 mu2 [1] 2 3 So now I still don't get my 1xk matrices, even though mu is a matrix and I've used drop=FALSE. Why? Because you are subsetting mu like a vector not like a matrix. The following produces the desired output: R mu1-mu[,iX1,drop=FALSE]; mu2-mu[,iX2,drop=FALSE]; R mu1 [,1] [1,]1 R mu2 [,1] [,2] [1,]23 Thus, iX1 and iX2 may only be used as the column index. Z (I'm getting a bit bewildered by all this, and must get it pinned down: the code this will go into is too complicated to allow easy debugging if the subsetting does unpredicted things.) [BTW: Just in case anyone gets the thought that it might work if you matrify iX1, iX2 e.g. iX1-matrix(iX1,nrow=1) -- well, it doesn't!] Best wishes to all, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 167 1972 Date: 15-Jul-03 Time: 11:14:05 -- XFMail -- __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] R, geochemistry, ternary diagrams
And in package vcd function ternaryplot(). g., -d On 2003.07.15 09:33, Tobias Verbeke wrote: Are there enough geochemists using R already that he'd find like-minded people to discuss technical issues with if he _did_ switch to R? Is there a package somewhere already that does ternary and other geochemistry diagrams? Another possibility for a ternary plot was mentioned by Prof Ripley in http://maths.newcastle.edu.au/~rking/R/help/02b/3637.html library(MASS) example(Skye) gives code and an example HTH, Tobias __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Subsetting a matrix [Again!]
On 15-Jul-03 Achim Zeileis wrote: [...] mu1-mu[iX1,drop=FALSE]; mu2-mu[iX2,drop=FALSE]; mu1 [1] 1 mu2 [1] 2 3 So now I still don't get my 1xk matrices, even though mu is a matrix and I've used drop=FALSE. Why? Because you are subsetting mu like a vector not like a matrix. The following produces the desired output: R mu1-mu[,iX1,drop=FALSE]; mu2-mu[,iX2,drop=FALSE]; Thanks!! (Must have some coffee). Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 167 1972 Date: 15-Jul-03 Time: 12:01:02 -- XFMail -- __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Multinomial Logit with multiple groups
Hi, I am inexperienced with ML and R so please be tolerant :) I am trying to replicate the method I have seen in a paper without success. If my understanding is correct (a big 'IF') it seems to use Multinomial Logit on multiple groups of various sizes, with 'nature' selecting the choice (the winner) - then uses Maximum Likelihood to optimise the parameters to produce a model for prediction. I have not found any examples which use this technique. What is worse the paper only really provides a summary of the method. So I am stuck! Here is an short summary extract from the paper describing the method: ** Suppose horse h* is observed to win a race. The multinomial logit model gives: exp(Vh*) Ph*= --for h* = 1,2,...,H. H 'Sigma'exp(Vh) h=1 A linear-in-parameters specification leads to: N Vh = 'Sigma' An*Zhn n=1 where Zhn=Zhn(Xh,Yn) is the measured value of attribute n for horse h in a race. The 'A' values in the equation are the parameters of the stochastic utility model that must be estimated from a sample of races. The likelihiood function can be written: j exp(L) = 'Product' Pjh* j=1 where j denetes a race, h* is the horse observed to win race j, and L is the log-likelihood function. * In an ideal world I would hope for the R code to solve a toy problem using the above method. I can provide a jpg of the paper and a dataset if required. But really, *any* help you could give to help me get to grips with it woul be great. Thanks, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] package announcement: Generalized Boosted Models (gbm)
Generalized Boosted Models (gbm) This package implements extensions to Y. Freund and R. Schapire's AdaBoost algorithm and J. Friedman's gradient boosting machine (aka multivariate adaptive regression trees, MART). It includes regression methods for least squares, absolute loss, logistic, Poisson, Cox proportional hazards/partial likelihood, and the AdaBoost exponential loss. It handles continuous, nominal, ordinal covariates as well as those containing missing values. This package also includes a preliminary out-of-bag estimator for the optimal number of iterations, graphical tools for lower dimensional projections of the fitted surface, and a few demos of example gbm sessions. gbm 1.0 will soon appear on CRAN. Earlier versions have been up for a few months and the latest includes many of the suggestions and fixes sent to me by the early adopters. Enjoy! Greg ___ Greg Ridgeway, Ph.D. Statistician RAND http://www.rand.org/methodology/stat/ ___ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-announce __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] eigen vector sign reversal
I've just installed R 1.7.1 under linux red hat I noticed sign reversal of eigen vectors ,some of them not all, upon using diag function relative to those obtained using R 1.4.1 this is gonna miss up lots of my previous scripts I wonder if there is a way to avoid this. best regards karim __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Multivariate regression method
Hi Folks,
Thanks to several people's suggestions and clarifications,
I think I have implemented a function which computes the
conditional mean and covariance matrix of a subset of the
dimensions of an MV-normal variable, given the values on the
other dimensions (these conditioning value can be presented
as a matrix, to deal with several cases at once).
The code is below, for anyone who would like to use it.
Comments will be welcome.
Two auxiliary functions .+. and ixX are defined as well as the main
function MV.regn
Example:
U1-rnorm(10);U2-rnorm(10);U3-rnorm(10);
X-cbind(U2+U3.U1+U3,U1+U2); mu-matrix(c(0,0,0),nrow=1);
S-matrix(c(2,1,1, 1,2,1, 1,1,2),nrow=3);
#Ex 1
MV.regn(S,mu,X[,1,drop=FALSE],1)
#Ex 2
MV.regn(S,mu,X[c(1,3,5,7),1:2],1,2)
==
%.+%-function(X,x){sweep(X,2,x,+)} ## Adds x to rows of X
ixX-function(X,...){(1:ncol(X))%in%c( ... )}
MV.regn-function(S,mu,x1,...){
### NB NB The k-variate MV variables etc are ROW vectors throughout
### (as in nxk matrix of data on n cases with k variables observed).
### NB NB S,mu,x1 MUST be arrays (matrices): create with drop=FALSE
### or specify when entering arguments, e.g.
### MV.regn(S,mu,X[,1,drop=FALSE],1)
### Arguments: S is the covariance matrix of MV X
###mu is the ROW (1xk matrix) expectation(X)
###x1 is matrix: rows are conditioning values for selected
### columns of X (NB if a single column make sure it's
### a matrix).
###... is an indexing vector or comma-list of numbers
### selecting the conditioning columns of X for the
### conditioning variable X1 (implies complementary set of
### columns of X for the variable X2 whose conditional
### distribution (X2 | X1=x1) is to be found).
iX1-ixX(S, ... ); iX2-!iX1;
s11-solve(S[iX1,iX1,drop=FALSE]); s12-S[iX1,iX2,drop=FALSE];
s21-S[iX2,iX1,drop=FALSE]; s22-S[iX2,iX2,drop=FALSE];
mu1-mu[,iX1,drop=FALSE]; mu2-mu[,iX2,drop=FALSE];
Cmu - (x1%.+%(-mu1))%*%s11%*%s12 %.+% mu2;
Cvar - s22 - s21%*%s11%*%s12;
list(Cmu=Cmu,Cvar=Cvar,iX1=iX1,iX2=iX2)
}
=
Best wishes to all,
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 13:23:25
-- XFMail --
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] printf and friends in R?
(Ted Harding) wrote: On 15-Jul-03 Barry Rowlingson wrote: help.search(printf) reveals the sprintf function (amongst others). Thanks! (For some reason this drew a blank when I tried it before ... ). perhaps because help.search does approximate matching and found everything with 'print' in as well? I get about 12 screens of matches to that, but because I knew sprintf was in there _somewhere_ I found it. You can turn off approximate matching with: help.search(printf,agrep=F) which gives the one true match. Baz __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] clearing some part of a graph
Hi R lovers
2 questions:
1) I'd like to know how to clean the title, the sub title or the labels of
a graph. I know how to redefine it with the function title() but it
overwrites the previous title and do not replace it
2) How could I clear a whole plot (for example in a multiple figure
environment)
thanks for your help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] function acf in package ts
Hi R lovers!
I'd like to know if the 2 blue lines in the graph produced by the function
acf in the package ts represents the level for the test of significance of
the autocorrelation
thanks for help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] printf and friends in R?
On Tue, 2003-07-15 at 08:27, Barry Rowlingson wrote: (Ted Harding) wrote: On 15-Jul-03 Barry Rowlingson wrote: help.search(printf) reveals the sprintf function (amongst others). Thanks! (For some reason this drew a blank when I tried it before ... ). perhaps because help.search does approximate matching and found everything with 'print' in as well? I get about 12 screens of matches to that, but because I knew sprintf was in there _somewhere_ I found it. You can turn off approximate matching with: help.search(printf,agrep=F) which gives the one true match. Baz Yet another option, which I forget about myself sometimes, is apropos(). For example: apropos(printf) [1] sprintf See ?apropos for more info. HTH, Marc Schwartz __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] clearing some part of a graph
Dear Vincent,
You can't clear a title only. The idea is to plot without one, like
plot(...,main=,sub=,xlab=,ylab=)
and add the title by hand (as you mentioned).
Moreover, AFAIK you can't clear one plot in a multiple figure environment,
since you can only clear the graphics device (not the plot),
but you can put a new plot onto a specific part of your m.f.e. using
par(mfg=c(i,j))
see ?par and An Introduction to R for more details.
(I have to admit that I am unsure whether that will remove a plot there,
or create the new one on top of it - hiding the old one. That matters
with respect to dev.copy() I assume.)
HTH
Thomas
---
Thomas Hotz
Research Associate in Medical Statistics
University of Leicester
United Kingdom
Department of Epidemiology and Public Health
22-28 Princess Road West
Leicester
LE1 6TP
Tel +44 116 252-5410
Fax +44 116 252-5423
Division of Medicine for the Elderly
Department of Medicine
The Glenfield Hospital
Leicester
LE3 9QP
Tel +44 116 256-3643
Fax +44 116 232-2976
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
Sent: 15 July 2003 14:31
To: [EMAIL PROTECTED]
Subject: [R] clearing some part of a graph
Hi R lovers
2 questions:
1) I'd like to know how to clean the title, the sub title or
the labels of
a graph. I know how to redefine it with the function title() but it
overwrites the previous title and do not replace it
2) How could I clear a whole plot (for example in a multiple figure
environment)
thanks for your help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are
confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] (no subject)
Hi I got a problem with creating a textfile: how can I create a textfile, which has a headline like: #data1 1 2 3 4 5 6 7 8 the problem is, how to bind text and matrix, so that the read.table-function will ignore the text and read the numbers. Thank you for help Michael -- 2 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] clearing some part of a graph
Many thanks
I don't see any other solution than what you've proposed.
I have tried the par(mfg=c(i,j,k,l))
But it overwrites rather than replaces the previous graph
take care
vincnet
|-+
| | [EMAIL PROTECTED]|
| | .uk |
| ||
| | 07/15/03 03:53 PM|
| ||
|-+
--|
|
|
| To: Vincent STOLIAROFF/fr/[EMAIL PROTECTED], [EMAIL PROTECTED]
|
| cc:
|
| Subject: RE: [R] clearing some part of a graph
|
--|
Dear Vincent,
You can't clear a title only. The idea is to plot without one, like
plot(...,main=,sub=,xlab=,ylab=)
and add the title by hand (as you mentioned).
Moreover, AFAIK you can't clear one plot in a multiple figure environment,
since you can only clear the graphics device (not the plot),
but you can put a new plot onto a specific part of your m.f.e. using
par(mfg=c(i,j))
see ?par and An Introduction to R for more details.
(I have to admit that I am unsure whether that will remove a plot there,
or create the new one on top of it - hiding the old one. That matters
with respect to dev.copy() I assume.)
HTH
Thomas
---
Thomas Hotz
Research Associate in Medical Statistics
University of Leicester
United Kingdom
Department of Epidemiology and Public Health
22-28 Princess Road West
Leicester
LE1 6TP
Tel +44 116 252-5410
Fax +44 116 252-5423
Division of Medicine for the Elderly
Department of Medicine
The Glenfield Hospital
Leicester
LE3 9QP
Tel +44 116 256-3643
Fax +44 116 232-2976
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
Sent: 15 July 2003 14:31
To: [EMAIL PROTECTED]
Subject: [R] clearing some part of a graph
Hi R lovers
2 questions:
1) I'd like to know how to clean the title, the sub title or
the labels of
a graph. I know how to redefine it with the function title() but it
overwrites the previous title and do not replace it
2) How could I clear a whole plot (for example in a multiple figure
environment)
thanks for your help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are
confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] skip first line with read.table (was: no subject)
Dear Michael, There was an email thread midst June 2003 on a related issue. The suggestions made there will certainly help you Have a look for thread Programcode and data in the same textfile (just type it into the R Site Search, and the thread will come up). I recall that read.table() has an argument skip which allows the first lines to be skipped, see ?read.table for details. The manual R Data Import/Export is very useful as well. Please provide a meaningful subject that makes searching the archives easier. HTH Thomas --- Thomas Hotz Research Associate in Medical Statistics University of Leicester United Kingdom Department of Epidemiology and Public Health 22-28 Princess Road West Leicester LE1 6TP Tel +44 116 252-5410 Fax +44 116 252-5423 Division of Medicine for the Elderly Department of Medicine The Glenfield Hospital Leicester LE3 9QP Tel +44 116 256-3643 Fax +44 116 232-2976 -Original Message- From: michael kirschbaum [mailto:[EMAIL PROTECTED] Sent: 15 July 2003 14:56 To: [EMAIL PROTECTED] Subject: [R] (no subject) Hi I got a problem with creating a textfile: how can I create a textfile, which has a headline like: #data1 1 2 3 4 5 6 7 8 the problem is, how to bind text and matrix, so that the read.table-function will ignore the text and read the numbers. Thank you for help Michael -- 2 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] How to use fortrane compiler to install unix R???
Hi: I am a R beginner,and I had a problem about installing R on unix(Dec-Tru64 alpha) my system had been installed fortrane complier,and I also follow manual document to do following procedure. ./configure make make install but ,when I do 'make' there are something wrong,and 'make' procedure will be stopped,and also 'make install'. please help me to solve this problem.thanks. [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] (no subject)
sink may be the simplest but won't give you the control you may want. For full control, have you considered cat and write.table? Both have append arguments plus others. hope this helps. spencer graves michael kirschbaum wrote: Hi I got a problem with creating a textfile: how can I create a textfile, which has a headline like: #data1 1 2 3 4 5 6 7 8 the problem is, how to bind text and matrix, so that the read.table-function will ignore the text and read the numbers. Thank you for help Michael __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] (no subject)
OK [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] (no subject)
[read text file :] #data1 1 2 3 4 5 6 7 8 Is this what you mean ? b - read.table(filewithyourdata, header=F, sep= ) b V1 V2 V3 V4 1 1 2 3 4 2 5 6 7 8 HTH, Tobias __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] (no subject)
Sorry Michael, I should read more carefully. You asked to create the file, not to read it. Tobias __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] How to read in data
Hello, I'm new to R and in the process of testing it My first question: I fail to read in my data (ANSI toto.txt file, tab separated) test -read.table(toto.txt) Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file `toto.txt' test - scan(C:\\toto.txt) Error in scan(C:\\toto.txt) : scan expected a real, got No_D test -scan(test.dat) Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file `toto.txt (and no, it is not read only or locked or whatever) I use Windows 2000/XP second question...what are the size limits of statistical files I can handle? I plan to analize plant datas (up to 500'000 records, from which I will analize a restrictive set of variates ) Even when broken down by some chracteristics, the data to analize can have 50'000-100'000 records Well thank for the help Anne [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] How to read in data
You probably need to change the directory to the one that contains the text file you are reading in. This is done under the File menu. Jim James W. MacDonald Affymetrix and cDNA Microarray Core University of Michigan Cancer Center 1500 E. Medical Center Drive 7410 CCGC Ann Arbor MI 48109 734-647-5623 Anne Piotet [EMAIL PROTECTED] 08/15/03 10:42AM Hello, I'm new to R and in the process of testing it My first question: I fail to read in my data (ANSI toto.txt file, tab separated) test -read.table(toto.txt) Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file `toto.txt' test - scan(C:\\toto.txt) Error in scan(C:\\toto.txt) : scan expected a real, got No_D test -scan(test.dat) Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file `toto.txt (and no, it is not read only or locked or whatever) I use Windows 2000/XP second question...what are the size limits of statistical files I can handle? I plan to analize plant datas (up to 500'000 records, from which I will analize a restrictive set of variates ) Even when broken down by some chracteristics, the data to analize can have 50'000-100'000 records Well thank for the help Anne [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] How to read in data
AP == Anne Piotet [EMAIL PROTECTED] disait: AP Hello, I'm new to R and in the process of testing it My first AP question: I fail to read in my data (ANSI toto.txt file, tab AP separated) test -read.table(toto.txt) AP Error in file(file, r) : unable to open AP connection In addition: Warning message: cannot open file AP `toto.txt' test - scan(C:\\toto.txt) AP Error in scan(C:\\toto.txt) : scan expected a AP real, got No_D test -scan(test.dat) AP Error in file(file, r) : unable to open AP connection In addition: Warning message: cannot open file AP `toto.txt (and no, it is not read only or locked or whatever) AP I use Windows 2000/XP I think read.table(C:\\toto.txt,header=TRUE) will do the job : the message you got on your first and third attempts means that you gave a wrong path to your file. otherwise, read the help for read.table carefully (header and skip parameters). AP second question...what are the size limits of statistical AP files I can handle? I plan to analize plant datas (up to AP 500'000 records, from which I will analize a restrictive set AP of variates ) Even when broken down by some chracteristics, AP the data to analize can have 50'000-100'000 records AP Well thank for the help Anne de rien ;) regards, --Mathieu __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] How to read in data
Anne Piotet wrote:
Hello,
I'm new to R and in the process of testing it
My first question: I fail to read in my data (ANSI toto.txt file, tab separated)
test -read.table(toto.txt)
Error in file(file, r) : unable to open connection
In addition: Warning message:
cannot open file `toto.txt'
- that's because it didnt find the file in that location.
test - scan(C:\\toto.txt)
Error in scan(C:\\toto.txt) : scan expected a real, got No_D
- that's because it did find the file, but there was the text No_D
in it. scan() will only read numbers unless you tell it otherwise.
test -scan(test.dat)
Error in file(file, r) : unable to open connection
In addition: Warning message:
cannot open file `toto.txt
again, its not looking in c:\, so it doesn't find it. Funny how
scan(test.dat) brings up an error about toto.txt :)
R has a working directory which is where scan() and read.file() will
start looking for files without a full path - type getwd() to see where
that is at any time.
You didnt try the other option:
test - read.table(c:\\toto.txt, sep='\t')
- I give a full path to toto.txt and tell it the columns are separated
with tabs ('\t'). You may need other options - popular ones are as.is=T
which keeps character variables as text rather than converting to
categorical data (factors), and head=T if the first line of the file is
a header with column names.
If this works, then do names(test) and summary(test) to see what
you've got.
second question...what are the size limits of statistical files I can handle? I plan to analize plant datas (up to 500'000 records, from which I will analize a restrictive set of variates ) Even when broken down by some chracteristics, the data to analize can have 50'000-100'000 records
Depends - whats the size of the machine you are using (and dont say
its a small box that fits under my monitor). How much RAM and disk space
does it have?
Baz
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Specifying an lme model
Ross Darnell [EMAIL PROTECTED] writes: I would like some advice on how if possible, to test the following I have subjects each measured several times. The subjects are sampled from 3 subpopulations (groups). The question is Is the between-subject variance the same for the three groups? The null model is lme0 - lme(y~group,random=~1|subject) I did think that the model that defined a specific between-subject variance for each group was update(lme0,.~., weights=varIdent(form=~1|group)) but I am not sure. I think you have it right. You should then compare the two fitted models using the anova generic, which will provide a likelihood ratio test statistic and a p-value based on a chi-squared reference distribution. Regard the p-value as an approximation. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] How to use fortrane compiler to install unix R???
On 15-Jul-03 kny wrote: Hi: I am a R beginner,and I had a problem about installing R on unix(Dec-Tru64 alpha) my system had been installed fortrane complier,and I also follow manual document to do following procedure. ./configure make make install but ,when I do 'make' there are something wrong,and 'make' procedure will be stopped,and also 'make install'. please help me to solve this problem.thanks. Did you get any apparent error messages from './configure'? If you did, this will indicate whether your system lacks something essential, or is not properly set up for the job. If not: An indication of where things went wrong, and possibly of what went wrong, should be available in the output from 'make'. However, this is often voluminous, so I suggest you try make 2somewhere/make.errors where somewhere is a convenient directory to store the stuff. Then page through this, and try to identify indications of problems. Good luck, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 167 1972 Date: 15-Jul-03 Time: 15:41:28 -- XFMail -- __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] function acf in package ts
From ?plot.acf:
Note:
The confidence interval plotted in `plot.acf' is based on an
uncorrelated series and should be treated with appropriate
caution. Using `ci.type = ma' may be less potentially
misleading.
Andy
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
Sent: Tuesday, July 15, 2003 9:36 AM
To: [EMAIL PROTECTED]
Subject: [R] function acf in package ts
Hi R lovers!
I'd like to know if the 2 blue lines in the graph produced by
the function acf in the package ts represents the level for
the test of significance of the autocorrelation
thanks for help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are
confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo /r-help
--
Notice: This e-mail message, together with any attachments, ...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] clearing some part of a graph
From: [EMAIL PROTECTED]
Many thanks
I don't see any other solution than what you've proposed.
I have tried the par(mfg=c(i,j,k,l))
But it overwrites rather than replaces the previous graph
Set the background color of the new plot to something like white so the
previous graph will be completely covered.
HTH,
Andy
take care
vincnet
Dear Vincent,
You can't clear a title only. The idea is to plot without one, like
plot(...,main=,sub=,xlab=,ylab=)
and add the title by hand (as you mentioned).
Moreover, AFAIK you can't clear one plot in a multiple figure
environment, since you can only clear the graphics device
(not the plot), but you can put a new plot onto a specific
part of your m.f.e. using
par(mfg=c(i,j))
see ?par and An Introduction to R for more details.
(I have to admit that I am unsure whether that will remove a
plot there, or create the new one on top of it - hiding the
old one. That matters with respect to dev.copy() I assume.)
HTH
Thomas
---
Thomas Hotz
Research Associate in Medical Statistics
University of Leicester
United Kingdom
Department of Epidemiology and Public Health
22-28 Princess Road West
Leicester
LE1 6TP
Tel +44 116 252-5410
Fax +44 116 252-5423
Division of Medicine for the Elderly
Department of Medicine
The Glenfield Hospital
Leicester
LE3 9QP
Tel +44 116 256-3643
Fax +44 116 232-2976
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
Sent: 15 July 2003 14:31
To: [EMAIL PROTECTED]
Subject: [R] clearing some part of a graph
Hi R lovers
2 questions:
1) I'd like to know how to clean the title, the sub title or the
labels of a graph. I know how to redefine it with the
function title()
but it overwrites the previous title and do not replace it
2) How could I clear a whole plot (for example in a multiple figure
environment)
thanks for your help
vincent
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your personal
address book. Please see http://www.sgcib.com for more information.
**
This message and any attachments (the message) are
confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
**
The sender's email address has changed to
firstname.lastname@ sgcib.com. You may want to update your
personal address book. Please see http://www.sgcib.com for more
information.
**
This message and any attachments (the message) are
confide...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo /r-help
--
Notice: This e-mail message, together with any attachments, ...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] predict
Good afternoon, Can you please tell me how R computes : predict(ar.x)$pred in : #let x be a vector ar.x-ar(x) predict(ar.x)$pred Thank you Bye __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] matrix manipulations
Hi cor(x,apply(x,1,sum)) gives me the correlations of each column with the sums of each row (correct me if I'm wrong, please). What I need are the correlations of each column with the sums of each row except the entry in the given column. It seems that for any one column i I get it by doing: cor(x[,i],apply(x[,-i],1,sum)) But I struggle to get it for all the columns. I was trying things like: for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum)) which doesn't generate any output at all, and rbind(for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum))) [,1] [1,] 0.1880237 outputs just the result of the very last column. I know that it shouldn't be necessary to use for(), but I couldn't figure out a way how to do the task using e.g. apply(). How do you get the results of all columns? Thank you, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] eigen vector sign reversal
I think at version 1.7.0 R started using LAPACK for its eigen/svd routines. I think using `eigen(x, EISPACK = TRUE)' uses the previous version. -roger Karim Elsawy wrote: I've just installed R 1.7.1 under linux red hat I noticed sign reversal of eigen vectors ,some of them not all, upon using diag function relative to those obtained using R 1.4.1 this is gonna miss up lots of my previous scripts I wonder if there is a way to avoid this. best regards karim __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] matrix manipulations
I don't think for loop is so bad here, but if you insist on not using it,
try:
x-matrix(rnorm(25), 5, 5)
sapply(1:5, function(i) cor(x[,i], rowSums(x[,-i])))
[1] -0.04179336 -0.08613796 0.48194936 0.38317629 -0.22081706
HTH,
Andy
-Original Message-
From: David Andel [mailto:[EMAIL PROTECTED]
Sent: Tuesday, July 15, 2003 12:52 PM
To: [EMAIL PROTECTED]
Subject: [R] matrix manipulations
Hi
cor(x,apply(x,1,sum)) gives me the correlations of each
column with the sums of each row (correct me if I'm wrong, please).
What I need are the correlations of each column with the sums
of each row except the entry in the given column. It seems
that for any one column i I get it by doing:
cor(x[,i],apply(x[,-i],1,sum))
But I struggle to get it for all the columns. I was trying
things like:
for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum))
which doesn't generate any output at all, and
rbind(for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum)))
[,1]
[1,] 0.1880237
outputs just the result of the very last column.
I know that it shouldn't be necessary to use for(), but I
couldn't figure out a way how to do the task using e.g. apply().
How do you get the results of all columns?
Thank you,
David
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo /r-help
--
Notice: This e-mail message, together with any attachments, ...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] (no subject)
One way might be to use a connection and `writeLines'. For example: a - matrix(1:8, byrow = TRUE, ncol = 4) a [,1] [,2] [,3] [,4] [1,]1234 [2,]5678 con - file(testfile.txt, w) writeLines(#data, con) write(a, con, ncol = 4) close(con) -roger michael kirschbaum wrote: Hi I got a problem with creating a textfile: how can I create a textfile, which has a headline like: #data1 1 2 3 4 5 6 7 8 the problem is, how to bind text and matrix, so that the read.table-function will ignore the text and read the numbers. Thank you for help Michael __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] matrix manipulations
David, I am not sure if it can be done in a vectorized form, but this should work y - NULL for(i in 1:ncol(x)) y - c(y, cor(x[,i],apply(x[,-i],1,sum))) Cheers, Andy __ Andy Jaworski Engineering Systems Technology Center 3M Center, 518-1-01 St. Paul, MN 55144-1000 - E-mail: [EMAIL PROTECTED] Tel: (651) 733-6092 Fax: (651) 736-3122 |-+ | | David Andel| | | [EMAIL PROTECTED] | | | Sent by: | | | [EMAIL PROTECTED]| | | ath.ethz.ch | | || | || | | 07/15/2003 11:51 | | || |-+ -| | | | To: [EMAIL PROTECTED] | | cc: | | Subject: [R] matrix manipulations | -| Hi cor(x,apply(x,1,sum)) gives me the correlations of each column with the sums of each row (correct me if I'm wrong, please). What I need are the correlations of each column with the sums of each row except the entry in the given column. It seems that for any one column i I get it by doing: cor(x[,i],apply(x[,-i],1,sum)) But I struggle to get it for all the columns. I was trying things like: for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum)) which doesn't generate any output at all, and rbind(for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum))) [,1] [1,] 0.1880237 outputs just the result of the very last column. I know that it shouldn't be necessary to use for(), but I couldn't figure out a way how to do the task using e.g. apply(). How do you get the results of all columns? Thank you, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Excel can do what R can't?????
Hi there
I thought this would be of particular interest to people using 'optim'
functions and perhaps people involved with R development.
I've been beaten down by R trying to get it to perform an optimization on a
mass-balance model. I've written the same program in excel, and using the
'solver' function, it comes up with an answer for my variables (p, ACT,
which I've assigned to q in R) that gives a solution to the function f in
about 3 seconds, with a value of the function around 0.0004. R, on the
other hand, appears to get stuck in local minima, and spits back an
approximation that is close the the p, ACT values excel does, but not
nearly precise enough for my needs, and not nearly as precise as excel, and
it takes about 3 minutes. Also, the solution for the value it returns for
the function is about 8 orders of magnitude greater than the excel version,
so I can't really say the function is approximating zero. I was able to
determine this using a trace command on function f, which is listed below.
This is very likely due to the fact that I've made some coding error along
the way, or have done something else wrong, but I don't know. Either way,
I am shocked and surprised that a program like excel is outperforming
R. I've attached my command file and the dataset temp.dat at the bottom
of this e-mail for anyone who would like to fiddle around with it, and if
you come up with something, PLEASE let me know- In the meantime, I've got
to start fiddling with excel and figuring out how to automate the solver
calculation.
Briefly, the point of the program is to approximate the model output from
an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt
and Hgt, by seeking the optimal values of p, ACT involved in the iterative
process.
Also, if your interested in recent correspondence that explains the point
of the program a bit, and how the function ties in to the iterative
process, search the R help forum for e-mails entitled [R] problem with
coding for 'optim' in R. Thanks also to Roger Peng and numerous others
for helping me get this far.
The whole point of me doing this in R was because it's supposed to be
spectacularly fast at automating complex loops, but seems to be falling
short for this application. Hopefully it's something wrong with my coding
and not with R itself.
Mike
R COMMAND FILE:
#perch.R #
# Hewett and Johnson bioenergetics #
# model combined with #
# Trudel MMBM to estimate #
# Consumption in perch in R code #
# Execute with #
# R --vanilla perch.R perch.out#
#USER INPUT BELOW
#Weight at time 0
Wo- 9.2
#Hg concentration at time 0 (ugHg/g wet weight)
Hgo- 0.08
#Weight at time t
Wt- 32.2
#Hg concentration at time t (ugHg/g wet weight)
Hgt- 0.110
#Prey methylmercury concentration (as constant)
Hgp- 0.033
#Prey caloric value (as constant)
Pc- 800
#Energy density of fish (as constant, calories)
Ef - 1000
#Maturity status, 0=immature, 1=mature
Mat- 0
#Sex, 1=male, 2=female
Sex- 1
#USER INPUT ABOVE
#Bioenergetics parameters for perch
CA - 0.25
CB - 0.73 #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d
CQ - 2.3
CTO - 23
CTM - 28
Zc- (log(CQ))*(CTM-CTO)
Yc- (log(CQ))*(CTM-CTO+2)
Xc- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400
RA - 34.992 #0.0108*3240 cal/g 02, converting weight of 02 to cal
RB - 0.8 #same as 1+(-0.2) see above...
RQ - 2.1
RTO - 28
RTM - 33
Za - (log(RQ))*(RTM-RTO)
Ya- (log(RQ))*(RTM-RTO+2)
Xa- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400
S - 0.172
FA - 0.158
FB - -0.222
FG - 0.631
UA- 0.0253
UB- 0.58
UG- -0.299
#Mass balance model parameters
EA - 0.002938
EB - -0.2
EQ - 0.066
a - 0.8
#Specifying sex-specific parameters
GSI- NULL
if (Sex==1) GSI-0.05 else
if (Sex==2) GSI-0.17
# Define margin of error functions
#merror - function(phat,M,alpha) # (1-alpha)*100% merror for a proportion
#{
#z - qnorm(1-alpha/2)
#merror - z * sqrt(phat*(1-phat)/M) # M is (Monte Carlo) sample size
#merror
#}
#Bring in temp file
temper - scan(temp.dat, na.strings = ., list(Day=0, jday=0, Temp=0))
Day-temper$Day ; jday-temper$jday ; Temp-temper$Temp ;
temp- cbind (Day, jday, Temp)
#Day = number of days modelled, jday=julian day, Temp = daily avg. temp.
#temp [,2]
Vc-(CTM-(temp[,3]))/(CTM-CTO)
Vr-(RTM-(temp[,3]))/(RTM-RTO)
comp- cbind (Day, jday, Temp, Vc, Vr)
#comp
bio-matrix(NA, ncol=13, nrow=length(Day))
W-NULL
C-NULL
ASMR-NULL
SMR-NULL
A-NULL
F-NULL
U-NULL
SDA-NULL
Gr-NULL
Hg-NULL
Ed-NULL
GHg-NULL
K-NULL
Expegk-NULL
EGK-NULL
p-NULL
ACT-NULL
#starting values for p, ACT
p - 1 # 0.558626306252032 #solution set for p, ACT from excel 'solver' f'n
ACT - 2 # 1.66764519286918
q-c(p,ACT)
#specify sttarting values
#q0-c(p = 1, ACT = 1)
#introduce function to solve
f - function (q)
{
M- length(Day) #number of days iterated
for (i in 1:M)
{
#Bioenergetics model
[R] Problem with my simulation
Dear R-users, thanks for your help last time! But now I've got a new problem with my simulation program. In order to save time I decided to divide my program in different parts. So I can use several computers (at the moment 3). On all these computers I installed the R-version 1.7.1. One works without problems. The other 2 don't work. During program`s running Rgui.exe causes an error, so the computer interrupts and closes R. Sometimes I don't get any error message and in the other case it seems useless. Could anybody give me hints? Thanks a lot! -- Peggy Ramlau __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Problem with my simulation
Seems like you need to get help here:
http://www.catb.org/~esr/faqs/smart-questions.html first before anyone here
can help you.
Andy
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Sent: Tuesday, July 15, 2003 1:51 PM
To: R-Project
Subject: [R] Problem with my simulation
Dear R-users,
thanks for your help last time!
But now I've got a new problem with my simulation program.
In order to save time I decided to divide my program in
different parts. So I can use several computers (at the
moment 3). On all these computers I installed the R-version
1.7.1. One works without problems. The other 2 don't work.
During program`s running Rgui.exe causes an error, so the
computer interrupts and closes R. Sometimes I don't get any
error message and in the other case it seems useless. Could
anybody give me hints?
Thanks a lot!
--
Peggy Ramlau
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo /r-help
--
Notice: This e-mail message, together with any attachments, ...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] matrix manipulations
What about the following: x - array(1:15, dim=c(5,3)) k - dim(x)[2] (P.k - outer(rep(1,k), rep(1,k))-diag(k)) [,1] [,2] [,3] [1,]011 [2,]101 [3,]110 (x.1 - (x %*%P.k)) [,1] [,2] [,3] [1,] 17 127 [2,] 19 149 [3,] 21 16 11 [4,] 23 18 13 [5,] 25 20 15 cor(x, x.1) [,1] [,2] [,3] [1,]111 [2,]111 [3,]111 This has the right pieces and the correct answer in this case and presumably in others. hope this helps. spencer graves [EMAIL PROTECTED] wrote: David, I am not sure if it can be done in a vectorized form, but this should work y - NULL for(i in 1:ncol(x)) y - c(y, cor(x[,i],apply(x[,-i],1,sum))) Cheers, Andy __ Andy Jaworski Engineering Systems Technology Center 3M Center, 518-1-01 St. Paul, MN 55144-1000 - E-mail: [EMAIL PROTECTED] Tel: (651) 733-6092 Fax: (651) 736-3122 |-+ | | David Andel| | | [EMAIL PROTECTED] | | | Sent by: | | | [EMAIL PROTECTED]| | | ath.ethz.ch | | || | || | | 07/15/2003 11:51 | | || |-+ -| | | | To: [EMAIL PROTECTED] | | cc: | | Subject: [R] matrix manipulations | -| Hi cor(x,apply(x,1,sum)) gives me the correlations of each column with the sums of each row (correct me if I'm wrong, please). What I need are the correlations of each column with the sums of each row except the entry in the given column. It seems that for any one column i I get it by doing: cor(x[,i],apply(x[,-i],1,sum)) But I struggle to get it for all the columns. I was trying things like: for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum)) which doesn't generate any output at all, and rbind(for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum))) [,1] [1,] 0.1880237 outputs just the result of the very last column. I know that it shouldn't be necessary to use for(), but I couldn't figure out a way how to do the task using e.g. apply(). How do you get the results of all columns? Thank you, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
I've programmed many things like this in both Excel and R. When I did not get the same answer from both, it was because I had an error in one (or both). I do this routinely as part of debugging: I catch many mistakes this way, and I often feel I can not trust my answers without this level of checking. I use Excel with Solver if I only need one solution or if I'm working with someone who doesn't have R or S-Plus. Otherwise, I prefer S-Plus of R. First forget about optim: Do you get the same numbers from your function f and from Excel? Have you plotted the function to be sure you even have local minima? Naively, I would expect Excel to be more likely to get stuck in local minima than optim. I'm sorry you've had such a frustrating experience with R. The S language is very powerful but does have a steep learning curve. I had S-Plus for 3-5 years before I was finally able to do anything useful with it. Now, it is an integral part of how I do much of what I do. hope this helps. spencer graves Michael Rennie wrote: Hi there I thought this would be of particular interest to people using 'optim' functions and perhaps people involved with R development. I've been beaten down by R trying to get it to perform an optimization on a mass-balance model. I've written the same program in excel, and using the 'solver' function, it comes up with an answer for my variables (p, ACT, which I've assigned to q in R) that gives a solution to the function f in about 3 seconds, with a value of the function around 0.0004. R, on the other hand, appears to get stuck in local minima, and spits back an approximation that is close the the p, ACT values excel does, but not nearly precise enough for my needs, and not nearly as precise as excel, and it takes about 3 minutes. Also, the solution for the value it returns for the function is about 8 orders of magnitude greater than the excel version, so I can't really say the function is approximating zero. I was able to determine this using a trace command on function f, which is listed below. This is very likely due to the fact that I've made some coding error along the way, or have done something else wrong, but I don't know. Either way, I am shocked and surprised that a program like excel is outperforming R. I've attached my command file and the dataset temp.dat at the bottom of this e-mail for anyone who would like to fiddle around with it, and if you come up with something, PLEASE let me know- In the meantime, I've got to start fiddling with excel and figuring out how to automate the solver calculation. Briefly, the point of the program is to approximate the model output from an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt and Hgt, by seeking the optimal values of p, ACT involved in the iterative process. Also, if your interested in recent correspondence that explains the point of the program a bit, and how the function ties in to the iterative process, search the R help forum for e-mails entitled [R] problem with coding for 'optim' in R. Thanks also to Roger Peng and numerous others for helping me get this far. The whole point of me doing this in R was because it's supposed to be spectacularly fast at automating complex loops, but seems to be falling short for this application. Hopefully it's something wrong with my coding and not with R itself. Mike R COMMAND FILE: #perch.R # # Hewett and Johnson bioenergetics # # model combined with # # Trudel MMBM to estimate # # Consumption in perch in R code # # Execute with # # R --vanilla perch.R perch.out# #USER INPUT BELOW #Weight at time 0 Wo- 9.2 #Hg concentration at time 0 (ugHg/g wet weight) Hgo- 0.08 #Weight at time t Wt- 32.2 #Hg concentration at time t (ugHg/g wet weight) Hgt- 0.110 #Prey methylmercury concentration (as constant) Hgp- 0.033 #Prey caloric value (as constant) Pc- 800 #Energy density of fish (as constant, calories) Ef - 1000 #Maturity status, 0=immature, 1=mature Mat- 0 #Sex, 1=male, 2=female Sex- 1 #USER INPUT ABOVE #Bioenergetics parameters for perch CA - 0.25 CB - 0.73 #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d CQ - 2.3 CTO - 23 CTM - 28 Zc- (log(CQ))*(CTM-CTO) Yc- (log(CQ))*(CTM-CTO+2) Xc- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400 RA - 34.992 #0.0108*3240 cal/g 02, converting weight of 02 to cal RB - 0.8 #same as 1+(-0.2) see above... RQ - 2.1 RTO - 28 RTM - 33 Za - (log(RQ))*(RTM-RTO) Ya- (log(RQ))*(RTM-RTO+2) Xa- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400 S - 0.172 FA - 0.158 FB - -0.222 FG - 0.631 UA- 0.0253 UB- 0.58 UG- -0.299 #Mass balance model parameters EA - 0.002938 EB - -0.2 EQ - 0.066 a - 0.8 #Specifying sex-specific parameters GSI- NULL if (Sex==1) GSI-0.05 else if (Sex==2) GSI-0.17 # Define margin of error functions #merror -
[R] Tree question
I was under the impression that the tree method (e.g. as implemented in rpart) was insensitive to monotonic transformations of the dependent variable. e.g. Breiman Olshen et al. Classification and Regression Trees state In a standard data structure [a tree] is invariant under all monotone transformations of individual ordered varaibles (p. 57) However, I get very different results from tr.hh.pri - rpart((log(YPRISX+1)~AGE+DRUGUSEY+SEX+OBSXNUM)) and tr.hh.pri - rpart(YPRISX~AGE+DRUGUSEY+SEX+OBSXNUM) the former gives more splits and different splits. Some notes: The DV is a count variable, and highly skew, with some 0s, many 1s, and a long right tail out to 99. AGE ranges from 18-25 DRUGUSEY is ordered (hardest drug used) and OBSXNUM is also ordered (proportion of your friends who object to your having 'casual sex') printing the first tree gives 1) root 307 23.472040 0.7114605 2) AGE=19.5 196 13.811070 0.6857971 4) OBSXNUM 2.5 69 5.712526 0.6338252 8) DRUGUSEY=1.5 15 2.261203 0.5161601 * 9) DRUGUSEY 1.5 54 3.185960 0.6665100 * 5) OBSXNUM=2.5 127 7.810911 0.7140339 * 3) AGE 19.5 111 9.303947 0.7567761 6) DRUGUSEY 0.5 48 1.105266 0.6727132 * 7) DRUGUSEY=0.5 63 7.601052 0.8208239 14) SEX=1.5 21 1.258395 0.7317629 * 15) SEX 1.5 42 6.092803 0.8653544 * printing the second tree gives 1) root 307 144.540700 1.1205210 2) AGE=19.5 196 68.382650 1.0561220 * 3) AGE 19.5 111 73.909910 1.2342340 6) DRUGUSEY 0.5 48 2.979167 0.9791667 * 7) DRUGUSEY=0.5 63 65.428570 1.4285710 14) SEX=1.5 21 6.571429 1.1428570 * 15) SEX 1.5 42 56.285710 1.5714290 * So, is this the 'exception that proves the rule'? Have I done something wrong? Or what? Any ideas or thoughts? Thanks in advance Peter Peter L. Flom, PhD Assistant Director, Statistics and Data Analysis Core Center for Drug Use and HIV Research National Development and Research Institutes 71 W. 23rd St www.peterflom.com New York, NY 10010 (212) 845-4485 (voice) (917) 438-0894 (fax) __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
Mike, The definition of your function f() seems quite inefficient. You could vectorize it, which would shorten and speed up your code, especially if M is large. See the R introduction file available online to learn how to do it if you don't already know how. Also, you have to return only one argument. Unless I'm wrong, your function wants to return Wtmod, Hgtmod, q and f. I'm don't think this would change anything in this case, but you should definitely clean this up! Another advice... If you can simplify your example into a few lines of ready-to-execute code with a toy dataset, then it's easy for everyone to try it out and you can get more feedback. The code you've included is quite large and cumbersome. For one thing, you could easily have removed the lines of code that were commented out. Meanwhile, I would suggest that you go back to the basics of R to clean up your code. Sorry I can't be more helpful. Jerome On July 15, 2003 10:46 am, Michael Rennie wrote: Hi there I thought this would be of particular interest to people using 'optim' functions and perhaps people involved with R development. I've been beaten down by R trying to get it to perform an optimization on a mass-balance model. I've written the same program in excel, and using the 'solver' function, it comes up with an answer for my variables (p, ACT, which I've assigned to q in R) that gives a solution to the function f in about 3 seconds, with a value of the function around 0.0004. R, on the other hand, appears to get stuck in local minima, and spits back an approximation that is close the the p, ACT values excel does, but not nearly precise enough for my needs, and not nearly as precise as excel, and it takes about 3 minutes. Also, the solution for the value it returns for the function is about 8 orders of magnitude greater than the excel version, so I can't really say the function is approximating zero. I was able to determine this using a trace command on function f, which is listed below. This is very likely due to the fact that I've made some coding error along the way, or have done something else wrong, but I don't know. Either way, I am shocked and surprised that a program like excel is outperforming R. I've attached my command file and the dataset temp.dat at the bottom of this e-mail for anyone who would like to fiddle around with it, and if you come up with something, PLEASE let me know- In the meantime, I've got to start fiddling with excel and figuring out how to automate the solver calculation. Briefly, the point of the program is to approximate the model output from an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt and Hgt, by seeking the optimal values of p, ACT involved in the iterative process. Also, if your interested in recent correspondence that explains the point of the program a bit, and how the function ties in to the iterative process, search the R help forum for e-mails entitled [R] problem with coding for 'optim' in R. Thanks also to Roger Peng and numerous others for helping me get this far. The whole point of me doing this in R was because it's supposed to be spectacularly fast at automating complex loops, but seems to be falling short for this application. Hopefully it's something wrong with my coding and not with R itself. Mike R COMMAND FILE: #perch.R # # Hewett and Johnson bioenergetics # # model combined with # # Trudel MMBM to estimate # # Consumption in perch in R code # # Execute with # # R --vanilla perch.R perch.out# #USER INPUT BELOW #Weight at time 0 Wo- 9.2 #Hg concentration at time 0 (ugHg/g wet weight) Hgo- 0.08 #Weight at time t Wt- 32.2 #Hg concentration at time t (ugHg/g wet weight) Hgt- 0.110 #Prey methylmercury concentration (as constant) Hgp- 0.033 #Prey caloric value (as constant) Pc- 800 #Energy density of fish (as constant, calories) Ef - 1000 #Maturity status, 0=immature, 1=mature Mat- 0 #Sex, 1=male, 2=female Sex- 1 #USER INPUT ABOVE #Bioenergetics parameters for perch CA - 0.25 CB - 0.73 #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d CQ - 2.3 CTO - 23 CTM - 28 Zc- (log(CQ))*(CTM-CTO) Yc- (log(CQ))*(CTM-CTO+2) Xc- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400 RA - 34.992 #0.0108*3240 cal/g 02, converting weight of 02 to cal RB - 0.8 #same as 1+(-0.2) see above... RQ - 2.1 RTO - 28 RTM - 33 Za - (log(RQ))*(RTM-RTO) Ya- (log(RQ))*(RTM-RTO+2) Xa- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400 S - 0.172 FA - 0.158 FB - -0.222 FG - 0.631 UA- 0.0253 UB- 0.58 UG- -0.299 #Mass balance model parameters EA - 0.002938 EB - -0.2 EQ - 0.066 a - 0.8 #Specifying sex-specific parameters GSI- NULL if (Sex==1) GSI-0.05 else if (Sex==2) GSI-0.17 # Define margin of error functions #merror -
[R] tree problem solved
You guys are fast! Several people pointed out that the tree method is insenstive to monotone transforms of the independent variables, not the dependent variable Thanks Peter __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Why two chisq.test p values differ when the contingency tableis transposed?
I'm using R1.7.0 runing with Win XP. Thanks, ...Tao x [,1] [,2] [1,] 149 151 [2,]18 t(x) [,1] [,2] [1,] 1491 [2,] 1518 chisq.test(x, simulate.p.value=T, B=10) Pearson's Chi-squared test with simulated p-value (based on 1e+05 replicates) data: x X-squared = 5.2001, df = NA, p-value = 0.03774 chisq.test(t(x), simulate.p.value=T, B=10) Pearson's Chi-squared test with simulated p-value (based on 1e+05 replicates) data: t(x) X-squared = 5.2001, df = NA, p-value = 0.01642 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Plotting a graph of many lines between groups of points...
I have a data file read into a data frame. For example, V1 V2 V3 V4 1 1 1 3 4 2 2 3 5 10 . . . . . . . . . . n V1[n] V2[n] V3[n] V4[n] to n=many thousand I want to plot a graph with many line segments, where v1[i]=x1, v2[i]=y1, v3[i]=x2, v4[i]=y2 for i=1,n. This seems relatively simple in theory but I've spent quite a bit of time trying to make it happen with plot(type=), points(x,y), or lines(x,y) to no avail. Do I need to turn these into vectors before plotting them? Any help would be greatly appreciated. Thanks, Andrew __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Why two chisq.test p values differ when the contingency
On 15-Jul-03 Tao Shi wrote:
x
[,1] [,2]
[1,] 149 151
[2,]18
t(x)
[,1] [,2]
[1,] 1491
[2,] 1518
chisq.test(x, simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: x
X-squared = 5.2001, df = NA, p-value = 0.03774
chisq.test(t(x), simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: t(x)
X-squared = 5.2001, df = NA, p-value = 0.01642
Possibly you may just have been unlucky, though the 0.03774 seems large:
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
B=10)$p.value)}
c2tx-chisq.test(tx, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2tx-c(c2tx,chisq.test(tx, simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx)
c2xc2tx
[1,] 0.01627 0.01720
[2,] 0.01672 0.01690
[3,] 0.01662 0.01669
[4,] 0.01733 0.01656
[5,] 0.01679 0.01777
[6,] 0.01715 0.01769
[7,] 0.01765 0.01769
[8,] 0.01703 0.01740
[9,] 0.01704 0.01708
[10,] 0.01669 0.01655
sd(c2x)
[1] 0.0003946715
sd(c2tx)
[1] 0.0004737099
Ted.
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 21:00:04
-- XFMail --
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Plotting a graph of many lines between groups ofpoints...
[EMAIL PROTECTED] wrote: I have a data file read into a data frame. For example, V1 V2 V3 V4 1 1 1 3 4 2 2 3 5 10 . . . . . . . . . . n V1[n] V2[n] V3[n] V4[n] to n=many thousand is this stored as a data.frame? a matrix? In any case there seems to be a problem with your indexing ... I want to plot a graph with many line segments, where v1[i]=x1, v2[i]=y1, v3[i]=x2, v4[i]=y2 for i=1,n. indexing is ok if you have df$v3[i] This seems relatively simple in theory but I've spent quite a bit of time trying to make it happen with plot(type=), points(x,y), or lines(x,y) to no avail. could you send the actual code for your plot? Perhaps some way to construct the data. at this point things are a bit murky. Do I need to turn these into vectors before plotting them? Thus my question about how the data are stored... Any help would be greatly appreciated. Thanks, Andrew __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help Cheers, Joel Joel F. Kincaid, Ph. D. Assistant Professor Department of Economics and Finance Franklin P. Perdue School of Business Salisbury University Salisbury Maryland, 21801 Phone: (410) 548-4416 Email: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
At 11:47 AM 7/15/03 -0700, Jerome Asselin wrote: Mike, The definition of your function f() seems quite inefficient. You could vectorize it, which would shorten and speed up your code, especially if M is large. Hi, Jerome I don;t think I can vectorize it, since in the iteration loop, the value for each [i] is dependent on the value of [i-1], so I require the loop to go through each [i] before I can get my values for any particular vector (variable). I actually had my program operating this way in the first place, but I get all sorts of warnings and the 'optim' function especially doesn't seem to appreciate it. See the R introduction file available online to learn how to do it if you don't already know how. Also, you have to return only one argument. Unless I'm wrong, your function wants to return Wtmod, Hgtmod, q and f. I'm don't think this would change anything in this case, but you should definitely clean this up! The calls to Wtmod, q, and Hgtmod are all just residual from the development of the loop inside function f. I would like to get the last line of 'bioday' reported from within the loop, had I figured out the optimization, but that point is rather moot unless I can get the optimization functioning. Another advice... If you can simplify your example into a few lines of ready-to-execute code with a toy dataset, then it's easy for everyone to try it out and you can get more feedback. The code you've included is quite large and cumbersome. For one thing, you could easily have removed the lines of code that were commented out. Meanwhile, I would suggest that you go back to the basics of R to clean up your code. Thanks for the advice- every bit helps if I eventually get this thing to work. Mike Sorry I can't be more helpful. Jerome On July 15, 2003 10:46 am, Michael Rennie wrote: Hi there I thought this would be of particular interest to people using 'optim' functions and perhaps people involved with R development. I've been beaten down by R trying to get it to perform an optimization on a mass-balance model. I've written the same program in excel, and using the 'solver' function, it comes up with an answer for my variables (p, ACT, which I've assigned to q in R) that gives a solution to the function f in about 3 seconds, with a value of the function around 0.0004. R, on the other hand, appears to get stuck in local minima, and spits back an approximation that is close the the p, ACT values excel does, but not nearly precise enough for my needs, and not nearly as precise as excel, and it takes about 3 minutes. Also, the solution for the value it returns for the function is about 8 orders of magnitude greater than the excel version, so I can't really say the function is approximating zero. I was able to determine this using a trace command on function f, which is listed below. This is very likely due to the fact that I've made some coding error along the way, or have done something else wrong, but I don't know. Either way, I am shocked and surprised that a program like excel is outperforming R. I've attached my command file and the dataset temp.dat at the bottom of this e-mail for anyone who would like to fiddle around with it, and if you come up with something, PLEASE let me know- In the meantime, I've got to start fiddling with excel and figuring out how to automate the solver calculation. Briefly, the point of the program is to approximate the model output from an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt and Hgt, by seeking the optimal values of p, ACT involved in the iterative process. Also, if your interested in recent correspondence that explains the point of the program a bit, and how the function ties in to the iterative process, search the R help forum for e-mails entitled [R] problem with coding for 'optim' in R. Thanks also to Roger Peng and numerous others for helping me get this far. The whole point of me doing this in R was because it's supposed to be spectacularly fast at automating complex loops, but seems to be falling short for this application. Hopefully it's something wrong with my coding and not with R itself. Mike R COMMAND FILE: #perch.R # # Hewett and Johnson bioenergetics # # model combined with # # Trudel MMBM to estimate # # Consumption in perch in R code # # Execute with # # R --vanilla perch.R perch.out# #USER INPUT BELOW #Weight at time 0 Wo- 9.2 #Hg concentration at time 0 (ugHg/g wet weight) Hgo- 0.08 #Weight at time t Wt- 32.2 #Hg concentration at time t (ugHg/g wet weight) Hgt- 0.110 #Prey methylmercury concentration (as constant) Hgp- 0.033 #Prey caloric value (as constant) Pc- 800 #Energy density of fish (as constant, calories) Ef - 1000 #Maturity status, 0=immature, 1=mature
Re: [R] Plotting a graph of many lines between groups of points...
[EMAIL PROTECTED] wrote: I have a data file read into a data frame. For example, V1 V2 V3 V4 1 1 1 3 4 2 2 3 5 10 . . . . . . . . . . n V1[n] V2[n] V3[n] V4[n] to n=many thousand is this stored as a data.frame? a matrix? In any case there seems to be a problem with your indexing ... I want to plot a graph with many line segments, where v1[i]=x1, v2[i]=y1, v3[i]=x2, v4[i]=y2 for i=1,n. indexing is ok if you have df$v3[i] This seems relatively simple in theory but I've spent quite a bit of time trying to make it happen with plot(type=), points(x,y), or lines(x,y) to no avail. could you send the actual code for your plot? Perhaps some way to construct the data. at this point things are a bit murky. Do I need to turn these into vectors before plotting them? Thus my question about how the data are stored... Any help would be greatly appreciated. Thanks, Andrew __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help Cheers, Joel Joel F. Kincaid, Ph. D. Assistant Professor Department of Economics and Finance Franklin P. Perdue School of Business Salisbury University Salisbury Maryland, 21801 Phone: (410) 548-4416 Email: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Why two chisq.test p values differ when the contingency
Hi, Ted and Dennis:
Thanks for your speedy replies! I don't think this happens just randomly, rather, I'm
thinking it may be due to the way chisq.test function handles simulation. Here shows
why: (Ted, I think there is an error in your code, tx should be t(x) )
x
[,1] [,2]
[1,] 149 151
[2,]18
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
+B=10)$p.value)}
c2tx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2tx-c(c2tx,chisq.test(t(x), simulate.p.value=T,
+ B=10)$p.value)}
cbind(c2x,c2tx)
c2xc2tx
[1,] 0.03727 0.01629
[2,] 0.03682 0.01662
[3,] 0.03671 0.01665
[4,] 0.03788 0.01745
[5,] 0.03706 0.01646
[6,] 0.03715 0.01728
[7,] 0.03664 0.01683
[8,] 0.03681 0.01720
[9,] 0.03742 0.01758
[10,] 0.03712 0.01685
[11,] 0.03739 0.01615
[12,] 0.03811 0.01653
[13,] 0.03711 0.01673
[14,] 0.03639 0.01678
[15,] 0.03714 0.01719
[16,] 0.03774 0.01780
[17,] 0.03574 0.01707
[18,] 0.03661 0.01705
[19,] 0.03751 0.01711
[20,] 0.03683 0.01718
[21,] 0.03678 0.01653
...Tao
[EMAIL PROTECTED] wrote:
On 15-Jul-03 Tao Shi wrote:
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
t(x)
[,1] [,2]
[1,] 149 1
[2,] 151 8
chisq.test(x, simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: x
X-squared = 5.2001, df = NA, p-value = 0.03774
chisq.test(t(x), simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: t(x)
X-squared = 5.2001, df = NA, p-value = 0.01642
Possibly you may just have been unlucky, though the 0.03774 seems large:
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
B=10)$p.value)}
c2tx-chisq.test(tx, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2tx-c(c2tx,chisq.test(tx, simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx)
c2x c2tx
[1,] 0.01627 0.01720
[2,] 0.01672 0.01690
[3,] 0.01662 0.01669
[4,] 0.01733 0.01656
[5,] 0.01679 0.01777
[6,] 0.01715 0.01769
[7,] 0.01765 0.01769
[8,] 0.01703 0.01740
[9,] 0.01704 0.01708
[10,] 0.01669 0.01655
sd(c2x)
[1] 0.0003946715
sd(c2tx)
[1] 0.0004737099
Ted.
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 21:00:04
-- XFMail --
-
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Excel can do what R can't?????
I thought that you could simplify your code by using something like c(0,W[-length(W)]) as opposed to W[i-1] in a loop, but now I understand it's not that easy. Unless you can analytically simplify the calculation of W in order to vectorize it, it's going to be slow. However, many of the lines don't depend on [i] and not on [i-1]. Therefore you could simplify those as they don't need to be calculated within the loop. HTH, Jerome On July 15, 2003 01:24 pm, Michael Rennie wrote: At 11:47 AM 7/15/03 -0700, Jerome Asselin wrote: Mike, The definition of your function f() seems quite inefficient. You could vectorize it, which would shorten and speed up your code, especially if M is large. Hi, Jerome I don;t think I can vectorize it, since in the iteration loop, the value for each [i] is dependent on the value of [i-1], so I require the loop to go through each [i] before I can get my values for any particular vector (variable). I actually had my program operating this way in the first place, but I get all sorts of warnings and the 'optim' function especially doesn't seem to appreciate it. See the R introduction file available online to learn how to do it if you don't already know how. Also, you have to return only one argument. Unless I'm wrong, your function wants to return Wtmod, Hgtmod, q and f. I'm don't think this would change anything in this case, but you should definitely clean this up! The calls to Wtmod, q, and Hgtmod are all just residual from the development of the loop inside function f. I would like to get the last line of 'bioday' reported from within the loop, had I figured out the optimization, but that point is rather moot unless I can get the optimization functioning. Another advice... If you can simplify your example into a few lines of ready-to-execute code with a toy dataset, then it's easy for everyone to try it out and you can get more feedback. The code you've included is quite large and cumbersome. For one thing, you could easily have removed the lines of code that were commented out. Meanwhile, I would suggest that you go back to the basics of R to clean up your code. Thanks for the advice- every bit helps if I eventually get this thing to work. Mike Sorry I can't be more helpful. Jerome On July 15, 2003 10:46 am, Michael Rennie wrote: Hi there I thought this would be of particular interest to people using 'optim' functions and perhaps people involved with R development. I've been beaten down by R trying to get it to perform an optimization on a mass-balance model. I've written the same program in excel, and using the 'solver' function, it comes up with an answer for my variables (p, ACT, which I've assigned to q in R) that gives a solution to the function f in about 3 seconds, with a value of the function around 0.0004. R, on the other hand, appears to get stuck in local minima, and spits back an approximation that is close the the p, ACT values excel does, but not nearly precise enough for my needs, and not nearly as precise as excel, and it takes about 3 minutes. Also, the solution for the value it returns for the function is about 8 orders of magnitude greater than the excel version, so I can't really say the function is approximating zero. I was able to determine this using a trace command on function f, which is listed below. This is very likely due to the fact that I've made some coding error along the way, or have done something else wrong, but I don't know. Either way, I am shocked and surprised that a program like excel is outperforming R. I've attached my command file and the dataset temp.dat at the bottom of this e-mail for anyone who would like to fiddle around with it, and if you come up with something, PLEASE let me know- In the meantime, I've got to start fiddling with excel and figuring out how to automate the solver calculation. Briefly, the point of the program is to approximate the model output from an iterative calculation, Wtmod and Hgtmod, to user-specified endpoints Wt and Hgt, by seeking the optimal values of p, ACT involved in the iterative process. Also, if your interested in recent correspondence that explains the point of the program a bit, and how the function ties in to the iterative process, search the R help forum for e-mails entitled [R] problem with coding for 'optim' in R. Thanks also to Roger Peng and numerous others for helping me get this far. The whole point of me doing this in R was because it's supposed to be spectacularly fast at automating complex loops, but seems to be falling short for this application. Hopefully it's something wrong with my coding and not with R itself. Mike R COMMAND FILE: #perch.R # # Hewett and Johnson bioenergetics # #
RE: [R] Plotting a graph of many lines between groups ofpoints.. .
From: Andrew Johnson [mailto:[EMAIL PROTECTED]
I have a data file read into a data frame.
For example,
V1 V2 V3 V4
1 1 1 3 4
2 2 3 5 10
. . . . .
. . . . .
n V1[n] V2[n] V3[n] V4[n]
to n=many thousand
I want to plot a graph with many line segments, where
v1[i]=x1, v2[i]=y1,
v3[i]=x2, v4[i]=y2 for i=1,n.
This seems relatively simple in theory but I've spent quite a
bit of time
trying to make it happen with plot(type=), points(x,y), or
lines(x,y) to no
avail.
Do I need to turn these into vectors before plotting them?
See if matplot(mydata[,c(1,3)], mydata[,c(2,4)]) does what you want. If so,
check the help page for matplot for other options.
HTH,
Andy
Any help would be greatly appreciated.
Thanks,
Andrew
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo /r-help
--
Notice: This e-mail message, together with any attachments, ...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: RE: [R] Why two chisq.test p values differ when the contingency
Hi Tao:
The P-values for 2x2 table are generated based on a random (discrete
uniform distribution) sampling of all possible 2x2 tables, conditioning
on the observed margin totals. If one of the cells is extremely small,
as in your case, you get a big difference in P-values. Suppose, you
changed the cell with value 1 to, say, 5 or 6, then the two P-values
are nearly the same. However, I don't understand why they should be so
different, since the set of all possible 2x2 tables will be the same in
both cases. I would be interested in knowing how this happens.
Ravi.
- Original Message -
From: Shi, Tao [EMAIL PROTECTED]
Date: Tuesday, July 15, 2003 4:37 pm
Subject: RE: [R] Why two chisq.test p values differ when the contingency
Hi, Ted and Dennis:
Thanks for your speedy replies! I don't think this happens just
randomly, rather, I'm thinking it may be due to the way chisq.test
function handles simulation. Here shows why: (Ted, I think there
is an error in your code, tx should be t(x) )
x
[,1] [,2]
[1,] 149 151
[2,]18
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
+B=10)$p.value)}
c2tx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2tx-c(c2tx,chisq.test(t(x), simulate.p.value=T,
+ B=10)$p.value)}
cbind(c2x,c2tx)
c2xc2tx
[1,] 0.03727 0.01629
[2,] 0.03682 0.01662
[3,] 0.03671 0.01665
[4,] 0.03788 0.01745
[5,] 0.03706 0.01646
[6,] 0.03715 0.01728
[7,] 0.03664 0.01683
[8,] 0.03681 0.01720
[9,] 0.03742 0.01758
[10,] 0.03712 0.01685
[11,] 0.03739 0.01615
[12,] 0.03811 0.01653
[13,] 0.03711 0.01673
[14,] 0.03639 0.01678
[15,] 0.03714 0.01719
[16,] 0.03774 0.01780
[17,] 0.03574 0.01707
[18,] 0.03661 0.01705
[19,] 0.03751 0.01711
[20,] 0.03683 0.01718
[21,] 0.03678 0.01653
...Tao
[EMAIL PROTECTED] wrote:
On 15-Jul-03 Tao Shi wrote:
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
t(x)
[,1] [,2]
[1,] 149 1
[2,] 151 8
chisq.test(x, simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: x
X-squared = 5.2001, df = NA, p-value = 0.03774
chisq.test(t(x), simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: t(x)
X-squared = 5.2001, df = NA, p-value = 0.01642
Possibly you may just have been unlucky, though the 0.03774 seems
large:
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
B=10)$p.value)}
c2tx-chisq.test(tx, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2tx-c(c2tx,chisq.test(tx, simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx)
c2x c2tx
[1,] 0.01627 0.01720
[2,] 0.01672 0.01690
[3,] 0.01662 0.01669
[4,] 0.01733 0.01656
[5,] 0.01679 0.01777
[6,] 0.01715 0.01769
[7,] 0.01765 0.01769
[8,] 0.01703 0.01740
[9,] 0.01704 0.01708
[10,] 0.01669 0.01655
sd(c2x)
[1] 0.0003946715
sd(c2tx)
[1] 0.0004737099
Ted.
---
-
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 21:00:04
-- XFMail -
-
-
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] greek problem
Hello everybody.I have a problem with R.At the print or similar commands, I am writing at greek but R Gui doesn't understant much of the it! Perhaps if I could change the font to a compatible with the greek's one, the problem would be solved.Although I don't know how I could change font.I would be pleased, if someone could help me.I will be looking forward for an answer. My e-mail is [EMAIL PROTECTED] - ÁðïêôÞóôå ôçí äùñåÜí [EMAIL PROTECTED] [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] matrix manipulations
Thanks a lot, this does exactly what I was looking for. On Tue, 15 Jul 2003, Liaw, Andy wrote: I don't think for loop is so bad here, but if you insist on not using it, try: x-matrix(rnorm(25), 5, 5) sapply(1:5, function(i) cor(x[,i], rowSums(x[,-i]))) [1] -0.04179336 -0.08613796 0.48194936 0.38317629 -0.22081706 HTH, Andy -Original Message- From: David Andel [mailto:[EMAIL PROTECTED] Sent: Tuesday, July 15, 2003 12:52 PM To: [EMAIL PROTECTED] Subject: [R] matrix manipulations Hi cor(x,apply(x,1,sum)) gives me the correlations of each column with the sums of each row (correct me if I'm wrong, please). What I need are the correlations of each column with the sums of each row except the entry in the given column. It seems that for any one column i I get it by doing: cor(x[,i],apply(x[,-i],1,sum)) But I struggle to get it for all the columns. I was trying things like: for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum)) which doesn't generate any output at all, and rbind(for(i in 1:ncol(x)) cor(x[,i],apply(x[,-i],1,sum))) [,1] [1,] 0.1880237 outputs just the result of the very last column. I know that it shouldn't be necessary to use for(), but I couldn't figure out a way how to do the task using e.g. apply(). How do you get the results of all columns? Thank you, David __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo /r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Specifying an lme model
Thanks for the help. Sorry but now comes the second question. What we have now assumes equal within-subject variances. How can I fit separate (pooled) within-subject variances for each group (similar to an unequal variance t-test) and test if this is better than the existing constant within-subject variance. Many thanks Ross Douglas Bates [EMAIL PROTECTED] writes: Ross Darnell [EMAIL PROTECTED] writes: I would like some advice on how if possible, to test the following I have subjects each measured several times. The subjects are sampled from 3 subpopulations (groups). The question is Is the between-subject variance the same for the three groups? The null model is lme0 - lme(y~group,random=~1|subject) I did think that the model that defined a specific between-subject variance for each group was update(lme0,.~., weights=varIdent(form=~1|group)) but I am not sure. I think you have it right. You should then compare the two fitted models using the anova generic, which will provide a likelihood ratio test statistic and a p-value based on a chi-squared reference distribution. Regard the p-value as an approximation. -- Ross Darnell School of Health and Rehabilitation Sciences University of Queensland, Brisbane QLD 4067 AUSTRALIA Email: [EMAIL PROTECTED] Phone +61 7 3365 6087 http://www.uq.edu.au/~uqrdarne/ __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why two chisq.test p values differ when the contingency
Shi, Tao [EMAIL PROTECTED] writes:
Hi, Ted and Dennis:
Thanks for your speedy replies! I don't think this happens just randomly, rather,
I'm thinking it may be due to the way chisq.test function handles simulation. Here
shows why: (Ted, I think there is an error in your code, tx should be t(x) )
x
[,1] [,2]
[1,] 149 151
[2,]18
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
+B=10)$p.value)}
c2tx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2tx-c(c2tx,chisq.test(t(x), simulate.p.value=T,
+ B=10)$p.value)}
cbind(c2x,c2tx)
c2xc2tx
[1,] 0.03727 0.01629
[2,] 0.03682 0.01662
I agree that this looks dodgy. The simulation is by taking samples of
tables consistent with the given marginals, so should be invariant
under transpose operations. I venture a guess that the algorithm is
somehow forgetting to count tables that are identical to the current
table. (Notice that there are really only ten tables consistent with
those marginals, with probabilities
dhyper(0:9,150,159,9)
[1] 0.002258834 0.020194876 0.079185170 0.178727311 0.255905014
0.241046013
[7] 0.149366119 0.058713525 0.013284864 0.001318274
and the differences between c2x and c2tx look suspiciously close to
0.020194876...)
You might want to file this as a bug report.
--
O__ Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Why two chisq.test p values differ when the contingency
On 15-Jul-03 Shi, Tao wrote:
Hi, Ted and Dennis:
Thanks for your speedy replies! I don't think this happens just
randomly, rather, I'm thinking it may be due to the way chisq.test
function handles simulation. Here shows why: (Ted, I think there is an
error in your code, tx should be t(x) )
Not so -- I had already transposed them:
x
[,1] [,2]
[1,] 149 151
[2,]18
tx
[,1] [,2]
[1,] 1491
[2,] 1518
Anyway, just to check, a third run (which as it happens follows straight
on from the ones previously reported since I had not used that instance
of R since) amd putting the new results alongside the previous ones:
c2trx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2trx-c(c2trx,chisq.test(t(x), simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx,c2trx)
c2xc2tx c2trx
[1,] 0.01627 0.01720 0.01628
[2,] 0.01672 0.01690 0.01686
[3,] 0.01662 0.01669 0.01706
[4,] 0.01733 0.01656 0.01705
[5,] 0.01679 0.01777 0.01633
[6,] 0.01715 0.01769 0.01782
[7,] 0.01765 0.01769 0.01688
[8,] 0.01703 0.01740 0.01683
[9,] 0.01704 0.01708 0.01689
[10,] 0.01669 0.01655 0.01721
Maybe Peter Dalgaard's suspicion is true, that one case is not being
counted. But this must be R-implementation/version-dependent. In my
case:
version
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major1
minor6.1
year 2002
month11
day 01
language R
Ted.
x
[,1] [,2]
[1,] 149 151
[2,]18
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
+B=10)$p.value)}
c2tx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2tx-c(c2tx,chisq.test(t(x), simulate.p.value=T,
+ B=10)$p.value)}
cbind(c2x,c2tx)
c2xc2tx
[1,] 0.03727 0.01629
[2,] 0.03682 0.01662
[3,] 0.03671 0.01665
[4,] 0.03788 0.01745
[5,] 0.03706 0.01646
[6,] 0.03715 0.01728
[7,] 0.03664 0.01683
[8,] 0.03681 0.01720
[9,] 0.03742 0.01758
[10,] 0.03712 0.01685
[11,] 0.03739 0.01615
[12,] 0.03811 0.01653
[13,] 0.03711 0.01673
[14,] 0.03639 0.01678
[15,] 0.03714 0.01719
[16,] 0.03774 0.01780
[17,] 0.03574 0.01707
[18,] 0.03661 0.01705
[19,] 0.03751 0.01711
[20,] 0.03683 0.01718
[21,] 0.03678 0.01653
...Tao
[EMAIL PROTECTED] wrote:
On 15-Jul-03 Tao Shi wrote:
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
t(x)
[,1] [,2]
[1,] 149 1
[2,] 151 8
chisq.test(x, simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: x
X-squared = 5.2001, df = NA, p-value = 0.03774
chisq.test(t(x), simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: t(x)
X-squared = 5.2001, df = NA, p-value = 0.01642
Possibly you may just have been unlucky, though the 0.03774 seems
large:
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
B=10)$p.value)}
c2tx-chisq.test(tx, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2tx-c(c2tx,chisq.test(tx, simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx)
c2x c2tx
[1,] 0.01627 0.01720
[2,] 0.01672 0.01690
[3,] 0.01662 0.01669
[4,] 0.01733 0.01656
[5,] 0.01679 0.01777
[6,] 0.01715 0.01769
[7,] 0.01765 0.01769
[8,] 0.01703 0.01740
[9,] 0.01704 0.01708
[10,] 0.01669 0.01655
sd(c2x)
[1] 0.0003946715
sd(c2tx)
[1] 0.0004737099
Ted.
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 21:00:04
-- XFMail --
-
Do you Yahoo!?
E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 167 1972
Date: 16-Jul-03 Time: 01:30:59
-- XFMail --
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Why two chisq.test p values differ when the contingency
Hi, Ted:
I guess this problem is platform-dependent. I just tied it on a R 1.6.1 runing on
Win2K, it gave me two different p values. But when I tried it on R1.7.0 on a Linux
Server, I got the similar result as you did. I have filed a bug-report as Peter
suggested.
...Tao
[EMAIL PROTECTED] wrote:
On 15-Jul-03 Shi, Tao wrote:
Hi, Ted and Dennis:
Thanks for your speedy replies! I don't think this happens just
randomly, rather, I'm thinking it may be due to the way chisq.test
function handles simulation. Here shows why: (Ted, I think there is an
error in your code, tx should be t(x) )
Not so -- I had already transposed them:
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
tx
[,1] [,2]
[1,] 149 1
[2,] 151 8
Anyway, just to check, a third run (which as it happens follows straight
on from the ones previously reported since I had not used that instance
of R since) amd putting the new results alongside the previous ones:
c2trx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2trx-c(c2trx,chisq.test(t(x), simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx,c2trx)
c2x c2tx c2trx
[1,] 0.01627 0.01720 0.01628
[2,] 0.01672 0.01690 0.01686
[3,] 0.01662 0.01669 0.01706
[4,] 0.01733 0.01656 0.01705
[5,] 0.01679 0.01777 0.01633
[6,] 0.01715 0.01769 0.01782
[7,] 0.01765 0.01769 0.01688
[8,] 0.01703 0.01740 0.01683
[9,] 0.01704 0.01708 0.01689
[10,] 0.01669 0.01655 0.01721
Maybe Peter Dalgaard's suspicion is true, that one case is not being
counted. But this must be R-implementation/version-dependent. In my
case:
version
_
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major 1
minor 6.1
year 2002
month 11
day 01
language R
Ted.
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
+ B=10)$p.value)}
c2tx-chisq.test(t(x), simulate.p.value=T, B=10)$p.value
for(i in (1:20)){c2tx-c(c2tx,chisq.test(t(x), simulate.p.value=T,
+ B=10)$p.value)}
cbind(c2x,c2tx)
c2x c2tx
[1,] 0.03727 0.01629
[2,] 0.03682 0.01662
[3,] 0.03671 0.01665
[4,] 0.03788 0.01745
[5,] 0.03706 0.01646
[6,] 0.03715 0.01728
[7,] 0.03664 0.01683
[8,] 0.03681 0.01720
[9,] 0.03742 0.01758
[10,] 0.03712 0.01685
[11,] 0.03739 0.01615
[12,] 0.03811 0.01653
[13,] 0.03711 0.01673
[14,] 0.03639 0.01678
[15,] 0.03714 0.01719
[16,] 0.03774 0.01780
[17,] 0.03574 0.01707
[18,] 0.03661 0.01705
[19,] 0.03751 0.01711
[20,] 0.03683 0.01718
[21,] 0.03678 0.01653
...Tao
[EMAIL PROTECTED] wrote:
On 15-Jul-03 Tao Shi wrote:
x
[,1] [,2]
[1,] 149 151
[2,] 1 8
t(x)
[,1] [,2]
[1,] 149 1
[2,] 151 8
chisq.test(x, simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: x
X-squared = 5.2001, df = NA, p-value = 0.03774
chisq.test(t(x), simulate.p.value=T, B=10)
Pearson's Chi-squared test with simulated p-value (based on
1e+05 replicates)
data: t(x)
X-squared = 5.2001, df = NA, p-value = 0.01642
Possibly you may just have been unlucky, though the 0.03774 seems
large:
c2x-chisq.test(x, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2x-c(c2x,chisq.test(x, simulate.p.value=T,
B=10)$p.value)}
c2tx-chisq.test(tx, simulate.p.value=T, B=10)$p.value
for(i in (1:9)){c2tx-c(c2tx,chisq.test(tx, simulate.p.value=T,
B=10)$p.value)}
cbind(c2x,c2tx)
c2x c2tx
[1,] 0.01627 0.01720
[2,] 0.01672 0.01690
[3,] 0.01662 0.01669
[4,] 0.01733 0.01656
[5,] 0.01679 0.01777
[6,] 0.01715 0.01769
[7,] 0.01765 0.01769
[8,] 0.01703 0.01740
[9,] 0.01704 0.01708
[10,] 0.01669 0.01655
sd(c2x)
[1] 0.0003946715
sd(c2tx)
[1] 0.0004737099
Ted.
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03 Time: 21:00:04
-- XFMail --
-
Do you Yahoo!?
E-Mail: (Ted Harding)
Fax-to-email: +44 (0)870 167 1972
Date: 16-Jul-03 Time: 01:30:59
-- XFMail --
-
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
