Re: [R] R-1.8.0 memory.limit()
James MacDonald wrote: Thanks Uwe, What is the theoretical limit for 32 bit Windows anyway? BTW, --max-mem-size=2000M works (although maybe not all usable?) I don't know the usable amount exactly, one has to look into Microsoft's docs. It must be a bit smaller than 2GB (differently from the problem of ~3.7GB usable out of 4GB, because of other address spaces). I've just looked into the sources: R calculates with integers and of course (too stupid I had not realized it at once) has the same 32bit limit! So try 2047M (last working one) and 2049M which becomes negative. memory.limit() [1] 2097152000 In R-1.7.1, --max-mem-size=2G *did* work (I get the same results as above). In R-1.7.1 (2003-06-16) on WinNT4.0 it did *not*. Uwe 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 Uwe Ligges [EMAIL PROTECTED] 10/07/03 12:28PM James MacDonald wrote: Using R-1.8.0 (d/l and compiled on 2003-10-01) on WinXP, I seem to be unable to determine the maximum memory allocated to R. The help still says to use memory.limit(size=NA), but this returns the value NA. In addition, I have set --max-mem-size=2G but I run out of memory somewhere around 500Mb (which is why I am trying to find out how much memory is allocated). I don't have any other programs open, so I should certainly be able to address more memory than that. Any ideas? Yes. 2GB is a bit more than the possible theoretical maximal value for a single process on 32-bit Windows. Instead, try starting with, e.g., --max-mem-size=1700M and you'll get memory.limit() [1] 1782579200 Uwe Ligges TIA 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 __ [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: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Richard A. O'Keefe wrote: I am puzzled by the advice to use is.na(x) - TRUE instead of x - NA. ?NA says Function `is.na-' may provide a safer way to set missingness. It behaves differently for factors, for example. However, MAY provide is a bit scary, and it doesn't say WHAT the difference in behaviour is. I must say that is.na(x) - ... is rather repugnant, because it doesn't work. What do I mean? Well, as the designers of SETL who many years ago coined the term sinister function call to talk about f(...)-..., pointed out, if you do f(x) - y then afterwards you expect f(x) == y to be true. So let's try it: x - c(1,NA,3) is.na(x) - c(FALSE,FALSE,TRUE) x [1] 1 NA NA is.na(x) [1] FALSE TRUE TRUE v So I _assigned_ c(FALSE,FALSE,TRUE) to is.na(x), but I _got_ c(FALSE,TRUE, TRUE) instead. ^ That is not how a well behaved sinister function call should work, and it's enough to scare someone off is.na()- forever. The obvious way to set elements of a variable to missing is ... - NA. Wouldn't it be better if that just plain worked? Can someone give an example of is.na()- and -NA working differently with a factor? I just tried it: x - factor(c(3,1,4,1,5,9)) y - x is.na(x) - x==1 y[y==1] - NA x [1] 3NA 4NA 59 Levels: 1 3 4 5 9 y [1] 3NA 4NA 59 Levels: 1 3 4 5 9 Both approaches seem to have given the same answer. What did I miss? As mentioned in another mail to R-help. I'm pretty sure there was (is?) a problem with character (and/or factor) and assignment of NAs, but I cannot (re)produce an example. I think something for the x - NA case has been fixed during the last year. What prevents me to think I'm completely confused is that the is.na()- usage is proposed in: ?NA, S Programming, the R Language Definition manual, R's News file, but I cannot find it in the green book right now. Uwe Ligges __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] SIMCA algorithm implementation
Dear All Is there a SIMCA (Soft Independent Modelling Class Analogy) implementation on R or does anyone know if is it possible to replicate the SIMCA algorithm using existing R functions? Thanks Mike White __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] plotting results from leaps library
Hi In trying to fit a linear model , I use the leaps() function to determine wich predictors I should include in my model. I would like to plot the Mallow's Cp criteria against p with the indexes of selected model variates as points labels Is there already such a function? (I could not find it) Thanks Anne [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Generating automatic plots
Hello,
I have been trying to write a small program to generate automatic plots. What
I want is to draw boxplots for some variables in my data frame, contrasting
with a variable called 'missing' that has value 1 if some variable in a
concrete case has at least one missing value, in order to check if the cases
that don't enter in my analysis are biased.
The first attempt was to start generating a list with the variables of
contrast and then apply the list to the plots:
varlist - c(var1, var2, var3, var4, ...)
for (i in 1:length(varlist)) {
+ jpeg(varlist[i].jpg, width=480, heigth=480)
+ boxplot (varlist[i] ~ missing, xlab=missing values, ylab=varlist[i])
+ }
But I got 'Error in model.frame(formula, rownames, variables, varnames,
extranames: invalid variable type'
I suppose that is because I forget something related to the extraction of
values in vectors, but I can't find it on the R manual neither in other books
about R that I have checked.
Is there a way to draw plots using commands like the one above?
Thanks,
Xavier
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Note this behaviour:
a-a
a-NA
mode(a)
[1] logical
a-a
is.na(a) - T
mode(a)
[1] character
However after either way of assigning NA to a, is.na(a) is true,
and it prints as NA, so I can't see it's ever likely to matter. [Why
do I say these things? Expect usual flood of examples where it
does matter.]
Also if a is a character vector, a[2] - NA coerces the NA to
as.character(NA); again, just as one would hope/expect.
I have to echo Richard O'K's remark: if - NA can ever go wrong,
is that not a bug rather than a feature?
Simon Fear
Senior Statistician
Syne qua non Ltd
Tel: +44 (0) 1379 69
Fax: +44 (0) 1379 65
email: [EMAIL PROTECTED]
web: http://www.synequanon.com
Number of attachments included with this message: 0
This message (and any associated files) is confidential and\...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] plotting results from leaps library
Anne Piotet wrote: Hi In trying to fit a linear model , I use the leaps() function to determine wich predictors I should include in my model. I would like to plot the Mallow's Cp criteria against p with the indexes of selected model variates as points labels Is there already such a function? (I could not find it) Thanks Anne [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help try subsets() in John Fox's car package. You need to use regsubsets() in package leaps, not leaps(), but the plot.subsets() function in car would appear to produce what you want. HTH Gav -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [T] +44 (0)20 7679 5522 ENSIS Research Fellow [F] +44 (0)20 7679 7565 ENSIS Ltd. ECRC [E] [EMAIL PROTECTED] UCL Department of Geography [W] http://www.ucl.ac.uk/~ucfagls/cv/ 26 Bedford Way[W] http://www.ucl.ac.uk/~ucfagls/ London. WC1H 0AP. %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Generating automatic plots
Xavier Fernández i Marín wrote:
...
varlist - c(var1, var2, var3, var4, ...)
Instead of a character vector with the names, it'd make life easier if
you had a list of the vectors...
# make sure you use the naming - makes life easier later.
mylist - list(var1=var1, var2=var2, var3=var3, var4=var4)
Then...
for(i in seq(along=mylist)) { # seq(along=...) is safest.
jpeg(names(mylist)[i], width=...)
boxplot(mylist[[i]] ~ missing, xlab=...)
dev.off() # you didn't have this. it's required.
}
# but I don't see how missing gets correctly passed in your example,
# so I'm not sure what it's supposed to do.
I suppose that is because I forget something related to the extraction of
values in vectors, but I can't find it on the R manual neither in other books
about R that I have checked.
Nope. It's about passing a strings of text to plot, instead of actual
data. You'd have to actually parse the string to produce the value.
And that's getting tricky (at least, it's tricky to my brain). The
above is pretty simple by comparison.
Cheers
Jason
--
Indigo Industrial Controls Ltd.
http://www.indigoindustrial.co.nz
64-21-343-545
[EMAIL PROTECTED]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
On Wed, 8 Oct 2003, Simon Fear wrote: Note this behaviour: a-a a-NA mode(a) [1] logical a-a is.na(a) - T mode(a) [1] character However after either way of assigning NA to a, is.na(a) is true, and it prints as NA, so I can't see it's ever likely to matter. [Why do I say these things? Expect usual flood of examples where it does matter.] Also if a is a character vector, a[2] - NA coerces the NA to as.character(NA); again, just as one would hope/expect. I have to echo Richard O'K's remark: if - NA can ever go wrong, is that not a bug rather than a feature? I don't think it can ever `go wrong', but it can do things other than the user intends. The intention of is.na- is clearer, and so perhaps user error is less likely? That is the thinking behind the function, anyway. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Why does a[which(b == c[d])] not work?
Dear list, I can not understand why the expression in the subject does not work correct: dcrn[which(fn == inve[2])] numeric(0) inve[2] [1] 406.7 dcrn[which(fn == 406.7)] [1] 1.3994e-07 1.3988e-07 1.3953e-07 1.3966e-07 1.3953e-07 1.3968e-07 Is this a kick self problem or an bug? Thaks very much Thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
On Wednesday 08 October 2003 11:27, Thomas Bock wrote: Dear list, I can not understand why the expression in the subject does not work correct: dcrn[which(fn == inve[2])] numeric(0) inve[2] [1] 406.7 dcrn[which(fn == 406.7)] [1] 1.3994e-07 1.3988e-07 1.3953e-07 1.3966e-07 1.3953e-07 1.3968e-07 The reason is that you shouldn't compare numeric output like that, consider the following example: R x - 406.7 + 1e-5 R x [1] 406.7 R x == 406.7 [1] FALSE whereas R x - 406.7 + 1e-20 R x [1] 406.7 R x == 406.7 [1] TRUE that is 1.) `==' comparisons have a certain tolerance 2.) the print output is not necessarily precisely your number Instead of using `==' you should use a comparison with a certain tolerance you can specify... hth, Z Is this a kick self problem or an bug? Thaks very much Thomas __ [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: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Well, that's a convincing argument, but maybe
it's the name that's worrying some of us. Maybe it would be
more intuitive if called set.na (sorry, I mean setNA).
Also is.na- cannot be used to create a new variable of
NAs, so is not a universal method, which is a shame for its
advocates.
I note also that for a vector you can assign a new NA using
either TRUE or FALSE:
a - 1:3
is.na(a[4])-F
a
[1] 1 2 3 NA
For a list, assigning F leaves the new element set to NULL.
Mind you, I suspect this would be a particularly stupid thing
to do, so I'm not going to lose any sleep over R's reaction to it.
-Original Message-
From: Prof Brian Ripley [mailto:[EMAIL PROTECTED]
I don't think it can ever `go wrong', but it can do things other than
the
user intends. The intention of is.na- is clearer, and so
perhaps user
error is less likely? That is the thinking behind the
function, anyway.
Simon Fear
Senior Statistician
Syne qua non Ltd
Tel: +44 (0) 1379 69
Fax: +44 (0) 1379 65
email: [EMAIL PROTECTED]
web: http://www.synequanon.com
Number of attachments included with this message: 0
This message (and any associated files) is confidential and\...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
Thomas Bock [EMAIL PROTECTED] writes: Dear list, I can not understand why the expression in the subject does not work correct: dcrn[which(fn == inve[2])] numeric(0) inve[2] [1] 406.7 dcrn[which(fn == 406.7)] [1] 1.3994e-07 1.3988e-07 1.3953e-07 1.3966e-07 1.3953e-07 1.3968e-07 Is this a kick self problem or an bug? Kick self, most likely. What is inve[2] - 406.7 -- 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] fitdistr, mle's and gamma distribution
Thanks, your advice worked. I don't have much experience with maths, and
therefore tried to stay away from dealing with optimization, but going
down to this level opens a lot of possibilities. For the record, the
code I used, as you suggested:
###
shape - mean(data)^2/var(data)
scale - var(data)/mean(data)
gamma.param1 - shape
gamma.param2 - scale
log.gamma.param1 - log(gamma.param1)
log.gamma.param2 - log(gamma.param2)
gammaLoglik - function(params,
negative=TRUE){
lglk - sum(dgamma(data, shape=exp(params[1]), scale=exp(params[2]),
log=TRUE))
if(negative)
return(-lglk)
else
return(lglk)
}
optim.list - optim(c(log.gamma.param1, log.gamma.param2), gammaLoglik)
gamma.param1 - exp(optim.list$par[1])
gamma.param2 - exp(optim.list$par[2])
optim converges and the parameters provide a much better fit to the data
than the method of moments parameters do.
Armed with this new knowledge, I attempted to write a log(likelihood)
optimization method for the Pareto distribution. My ignorance of math
however shows, as the code does not work.
Here is what I did:
pareto.density - function(x, alpha, beta, log=FALSE){
# Pareto distribution not defined for x alpha
# Test for values x alpha, taking into account floating point error
test.vals - x-alpha
error.vals - test.vals[test.vals0]
if (length(error.vals)0)
stop(\nERROR: x alpha in pareto.distr(x, alpha, beta,
+ log=FALSE))
density - beta * (alpha^beta) * (x^(-beta-1))
if(log)
return(log(density))
else
return(density)
}
data.sorted - sort(data)
alpha.val - data.sorted[1]
beta.val - 1/((1/n) * sum(log(data.sorted/alpha.val)))
log.alpha.val - log(alpha.val)
log.beta.val - log(beta.val)
paretoLoglik - function(params, negative=TRUE){
lglk - sum(pareto.density(data.sorted, alpha=exp(params[1]),
+ beta=exp(params[2]), log=TRUE))
if(negative)
return(-lglk)
else
return(lglk)
}
optim.list - optim(c(log.alpha.val, log.beta.val), paretoLoglik,
+ method=L-BFGS-B, lower=c(log.alpha.val, 0), upper=c(log.alpha.val,
+ Inf))
pareto.param1 - exp(optim.list$par[1])
pareto.param2 - exp(optim.list$par[2])
#
I fixed the alpha parameter as my Pareto density function produces an
error if a datavalue alpha.
I get the following output:
source(browsesessoffplotfitted.R)
Error in optim(c(log.alpha.val, log.beta.val), paretoLoglik, method =
+ L-BFGS-B, :
non-finite finite-difference value [0]
Any ideas would be appreciated, otherwise its back to method of moments
for the Pareto distribution for me :)
Thanks
Lourens
On Tue, 2003-09-30 at 22:07, Spencer Graves wrote:
In my experience, the most likely cause of this problem is that
optim may try to test nonpositive values for shape or scale. I avoid
this situation by programming the log(likelihood) in terms of log(shape)
and log(scale) as follows:
gammaLoglik -
+ function(x, logShape, logScale, negative=TRUE){
+ lglk - sum(dgamma(x, shape=exp(logShape), scale=exp(logScale),
+ log=TRUE))
+ if(negative) return(-lglk) else return(lglk)
+ }
tst - rgamma(10, 1)
gammaLoglik(tst, 0, 0)
[1] 12.29849
Then I then call optim directly to minimize the negative of the
log(likelihood).
If I've guessed correctly, this should fix the problem. If not,
let us know.
hope this helps. spencer graves
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Contrast specified with C() - R vs S-Plus problem
Hi, For a n-level factor, I'd like to specify the first contrast and have the remaining n-2 constructed automatically so that the set is orthogonal. I then test the contrasts with summary.lm(anova-object). In S-Plus, the following works: y.anova - aov( y ~ C(CO2,c(1,0,-1)) ) summary.lm(y.anova) In R, it fails with the following error: levels(CO2) [1] A C E y.anova - aov(y + C(CO2,c(1,0,-1)) ) Error in contrasts-(*tmp*, value = contr) : wrong number of contrast matrix rows What is the way to do this in R? Thanks Pascal -- Dr. Pascal A. Niklaus Institute of Botany University of Basel Schönbeinstrasse 6 CH-4056 Basel / Switzerland __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Contrast specified with C() - R vs S-Plus problem
On Wed, 8 Oct 2003, Pascal A. Niklaus wrote: Hi, For a n-level factor, I'd like to specify the first contrast and have the remaining n-2 constructed automatically so that the set is orthogonal. I then test the contrasts with summary.lm(anova-object). In S-Plus, the following works: y.anova - aov( y ~ C(CO2,c(1,0,-1)) ) summary.lm(y.anova) I can't reproduce that in S-PLUS 6.1, and it is not as documented: contr what contrasts to use. May be one of four standard names ( helmert, poly, treatment, or sum), a function, or a matrix with as many rows as there are levels to the factor. You haven't given a matrix, and your vector gets converted to a 3x1 matrix when there are four levels. I suspect you have not got the same data in S-PLUS and R, but you haven't given us anything to check that. In R, it fails with the following error: levels(CO2) [1] A C E y.anova - aov(y + C(CO2,c(1,0,-1)) ) Error in contrasts-(*tmp*, value = contr) : wrong number of contrast matrix rows What is the way to do this in R? There are four levels, so CO2 - factor(c(, A, C, E)) attributes(C(CO2, as.matrix(1:4))) $levels [1] A C E $class [1] factor $contrasts [,1] [,2] [,3] 1 0.0236068 0.5472136 A2 -0.4393447 -0.7120227 C3 0.8078689 -0.2175955 E4 -0.3921311 0.3824045 does as you ask. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] fitdistr, mle's and gamma distribution
I'm sorry, but I don't have time to read all your code. However,
I saw that you tested for x alpha in your Pareto distribution
example. Have you considered reparameterizing to estimate log.del =
log(alpha-min(x))? Pass log.del as part of the vector of parameters to
estimate, then immediately inside the function, compute
alpha - (min(x)+exp(log.del))
I've fixed many convergence problems with simple tricks like this.
hope this helps. spencer graves
Lourens Olivier Walters wrote:
Thanks, your advice worked. I don't have much experience with maths, and
therefore tried to stay away from dealing with optimization, but going
down to this level opens a lot of possibilities. For the record, the
code I used, as you suggested:
###
shape - mean(data)^2/var(data)
scale - var(data)/mean(data)
gamma.param1 - shape
gamma.param2 - scale
log.gamma.param1 - log(gamma.param1)
log.gamma.param2 - log(gamma.param2)
gammaLoglik - function(params,
negative=TRUE){
lglk - sum(dgamma(data, shape=exp(params[1]), scale=exp(params[2]),
log=TRUE))
if(negative)
return(-lglk)
else
return(lglk)
}
optim.list - optim(c(log.gamma.param1, log.gamma.param2), gammaLoglik)
gamma.param1 - exp(optim.list$par[1])
gamma.param2 - exp(optim.list$par[2])
optim converges and the parameters provide a much better fit to the data
than the method of moments parameters do.
Armed with this new knowledge, I attempted to write a log(likelihood)
optimization method for the Pareto distribution. My ignorance of math
however shows, as the code does not work.
Here is what I did:
pareto.density - function(x, alpha, beta, log=FALSE){
# Pareto distribution not defined for x alpha
# Test for values x alpha, taking into account floating point error
test.vals - x-alpha
error.vals - test.vals[test.vals0]
if (length(error.vals)0)
stop(\nERROR: x alpha in pareto.distr(x, alpha, beta,
+ log=FALSE))
density - beta * (alpha^beta) * (x^(-beta-1))
if(log)
return(log(density))
else
return(density)
}
data.sorted - sort(data)
alpha.val - data.sorted[1]
beta.val - 1/((1/n) * sum(log(data.sorted/alpha.val)))
log.alpha.val - log(alpha.val)
log.beta.val - log(beta.val)
paretoLoglik - function(params, negative=TRUE){
lglk - sum(pareto.density(data.sorted, alpha=exp(params[1]),
+ beta=exp(params[2]), log=TRUE))
if(negative)
return(-lglk)
else
return(lglk)
}
optim.list - optim(c(log.alpha.val, log.beta.val), paretoLoglik,
+ method=L-BFGS-B, lower=c(log.alpha.val, 0), upper=c(log.alpha.val,
+ Inf))
pareto.param1 - exp(optim.list$par[1])
pareto.param2 - exp(optim.list$par[2])
#
I fixed the alpha parameter as my Pareto density function produces an
error if a datavalue alpha.
I get the following output:
source(browsesessoffplotfitted.R)
Error in optim(c(log.alpha.val, log.beta.val), paretoLoglik, method =
+ L-BFGS-B, :
non-finite finite-difference value [0]
Any ideas would be appreciated, otherwise its back to method of moments
for the Pareto distribution for me :)
Thanks
Lourens
On Tue, 2003-09-30 at 22:07, Spencer Graves wrote:
In my experience, the most likely cause of this problem is that
optim may try to test nonpositive values for shape or scale. I avoid
this situation by programming the log(likelihood) in terms of log(shape)
and log(scale) as follows:
gammaLoglik -
+ function(x, logShape, logScale, negative=TRUE){
+ lglk - sum(dgamma(x, shape=exp(logShape), scale=exp(logScale),
+ log=TRUE))
+ if(negative) return(-lglk) else return(lglk)
+ }
tst - rgamma(10, 1)
gammaLoglik(tst, 0, 0)
[1] 12.29849
Then I then call optim directly to minimize the negative of the
log(likelihood).
If I've guessed correctly, this should fix the problem. If not,
let us know.
hope this helps. spencer graves
__
[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] Why does a[which(b == c[d])] not work?
Your question has been answered by Achim and Peter Dalgaard (at least). Just a note: Using a[which(logic)] looks like a clumsy and inefficient way of writing a[ logic ] and I think you shouldn't propagate its use ... 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] New R - recompiling all packages
Hi All, I'm running R 1.7.1, and I've installed some additional packages such a Bioconductor. Do I've to re-install all the additional packages when ugrading to R 1.8.0 (i.e. are there compile in dependencies)? thanks for your help, Arne __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] R-1.8.0 is released
I've rolled up R-1.8.0.tgz a short while ago. This is a new version
with major changes (see below). Notably, the Macintosh version for OS
X has been substantially improved; the old Carbon interface is no
longer being supported.
Also notice that the underscore will no longer work as an assignment
operator.
There is also a bunch of new functions and an assortment of bugs have
been fixed.
You can get it from
http://cran.us.r-project.org/src/base/R-1.8.0.tgz
or wait for it to be mirrored at a CRAN site nearer to you. Binaries
for various platforms will appear in due course.
There is also a version split for floppies.
These are the md5sums for the freshly created files, in case you wish
to check that they are uncorrupted:
bb019d7e12e38ac8a8bc247f06cc6b42 R-1.8.0.tgz
80fd20fdd2b995ab69cfd8cde80267cf R-1.8.0.tgz-split.aa
074140872c9c015089543c9497a3be87 R-1.8.0.tgz-split.ab
f0e6005491839dd27636020f6475f4e2 R-1.8.0.tgz-split.ac
e7e817911d57e3c88959d6b86c061dd5 R-1.8.0.tgz-split.ad
52e9aff83357b65926f422abececba92 R-1.8.0.tgz-split.ae
6c6965e7f97b627280326ec6c436ab48 R-1.8.0.tgz-split.af
00ce1b1384f403ba76d8de302e92da1d R-1.8.0.tgz-split.ag
On behalf of the R Core Team,
Peter Dalgaard
Here's the relevant part of the NEWS file
CHANGES IN R VERSION 1.8.0
MACOS CHANGES
o As from this release there is only one R port for the
Macintosh, which runs only on MacOS X. (The `Carbon' port has
been discontinued, and the `Darwin' port is part of the new
version.) The current version can be run either as a
command-line application or as an `Aqua' console. There is a
`Quartz' device quartz(), and the download and installation of
both source and binary packages is supported from the Aqua
console. Those CRAN and BioC packages which build under MacOS
X have binary versions updated daily.
USER-VISIBLE CHANGES
o The defaults for glm.control(epsilon=1e-8, maxit=25) have been
tightened: this will produce more accurate results, slightly
slower.
o sub, gsub, grep, regexpr, chartr, tolower, toupper, substr,
substring, abbreviate and strsplit now handle missing values
differently from NA.
o Saving data containing name space references no longer warns
about name spaces possibly being unavailable on load.
o On Unix-like systems interrupt signals now set a flag that is
checked periodically rather than calling longjmp from the
signal handler. This is analogous to the behavior on Windows.
This reduces responsiveness to interrupts but prevents bugs
caused by interrupting computations in a way that leaves the
system in an inconsistent state. It also reduces the number
of system calls, which can speed up computations on some
platforms and make R more usable with systems like Mosix.
CHANGES TO THE LANGUAGE
o Error and warning handling has been modified to incorporate a
flexible condition handling mechanism. See the online
documentation of 'tryCatch' and 'signalCondition'. Code that
does not use these new facilities should remain unaffected.
o A triple colon operator can be used to access values of internal
variables in a name space (i.e. a:::b is the value of the internal
variable b in name space a).
o Non-syntactic variable names can now be specified by inclusion
between backticks `Like This`. The deparse() code has been
changed to output non-syntactical names with this convention,
when they occur as operands in expressions. This is controlled
by a `backtick' argument, which is by default TRUE for
composite expressions and FALSE for single symbols. This
should give minimal interference with existing code.
o Variables in formulae can be quoted by backticks, and such
formulae can be used in the common model-fitting functions.
terms.formula() will quote (by backticks) non-syntactic names
in its term.labels attribute. [Note that other code using
terms objects may expect syntactic names and/or not accept
quoted names: such code will still work if the new feature is
not used.]
NEW FEATURES
o New function bquote() does partial substitution like LISP backquote.
o capture.output() takes arbitrary connections for `file' argument.
o contr.poly() has a new `scores' argument to use as the base set
for the polynomials.
o cor() has a new argument `method = c(pearson,spearman,kendall)'
as cor.test() did forever. The two rank based measures do work with
all three missing value strategies.
o New utility function cov2cor() {Cov - Corr matrix}.
o cut.POSIXt() now allows `breaks' to be more general intervals
as allowed for the `by' argument to seq.POSIXt().
o data() now
Re: [R] Contrast specified with C() - R vs S-Plus problem
Thanks for the reply. Below is the solution and the S-Plus and R code that does the same (for documentation). I can't reproduce that in S-PLUS 6.1, and it is not as documented: In S-Plus 2000, C() complements the contrast matrix with orthogonal contrasts if only the first is given. CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum of Sq Mean Sq F Value Pr(F) C(CO2, c(1, 0, -1)) 2 16.0 8.00 7.259259 0.004013532 Residuals 21 23.14286 1.102041 summary.lm(aov(m ~ C(CO2,c(1,0,-1 Coefficients: Value Std. Error t value Pr(|t|) (Intercept) 2.5000 0.214311.6667 0. C(CO2, c(1, 0, -1))1 -1. 0.2624-3.8103 0.0010 C(CO2, c(1, 0, -1))2 0. 0.3712 0. 1. Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088 F-statistic: 7.259 on 2 and 21 degrees of freedom, the p-value is 0.004014 Correlation of Coefficients: (Intercept) C(CO2, c(1, 0, -1))1 C(CO2, c(1, 0, -1))1 0 C(CO2, c(1, 0, -1))2 0 0 In the meanwhile, I found that the same can be done in R like this: library(gregmisc) CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) contrastsCO2 - rbind(A vs E=c(1,0,-1)) a - aov(m ~ CO2, contrasts=list( CO2=make.contrasts(contrastsCO2))) summary(a) Df Sum Sq Mean Sq F value Pr(F) CO2 2 16.000 8.000 7.2593 0.004014 ** Residuals 21 23.143 1.102 summary.lm(a) Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.500e+00 2.143e-0111.67 1.22e-10 *** CO2A vs E -2.000e+00 5.249e-01-3.81 0.00102 ** CO2C29.774e-17 3.712e-01 2.63e-16 1.0 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088, Adjusted R-squared: 0.3525 F-statistic: 7.259 on 2 and 21 DF, p-value: 0.004014 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] New R - recompiling all packages
Most packages do not need to be re-installed, but R 1.8.0 has some nicer formatting of help pages so you may want to do so, and I am reinstalling all packages that make use of the methods package and so have saved images. For Bioconductor you need to ask on their list, but my understanding is that many of the packages have been updated for 1.8.0 and so you need to get new versions. Packages effects and tree have updated versions for 1.8.0 which update.packages() will pick up once your CRAN mirror gets updated (they are currently in 1.8.0/Other). On Wed, 8 Oct 2003 [EMAIL PROTECTED] wrote: Hi All, I'm running R 1.7.1, and I've installed some additional packages such a Bioconductor. Do I've to re-install all the additional packages when ugrading to R 1.8.0 (i.e. are there compile in dependencies)? thanks for your help, Arne -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] fitdistr, mle's and gamma distribution
Thanks for the help, the wrapper function was very useful. I managed to
solve the problem using Spencer Graves' suggestion. I am analyzing the
interarrival times between HTTP packets on a campus network. The dataset
actually has more than 14 Million entries! It represents the traffic
generated by approximately 3000 users browsing the web for 30 days. I
have to be careful to always remove unused objects from my workspace,
but otherwise I have so far managed to cope with 512Mb of memory on a
Pentium 600Mhz.
Lourens
On Tue, 2003-09-30 at 23:25, Ben Bolker wrote:
PS. 11 MILLION entries??
On Tue, 30 Sep 2003, Ben Bolker wrote:
Spencer Graves's suggestion of using shape and scale parameters on a log
scale is a good one.
To do specifically what you want (check values for which the objective
function is called and see what happens) you can do the following
(untested!), which makes a local copy of dgamma that you can mess with:
dgamma.old - dgamma
dgamma - function(x,shape,rate,...) {
d - dgamma.old(x,shape,rate,...)
cat(shape,rate,d,\n)
return(d)
}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Bootstrap Question
I have a question regarding bootstrap coverage. I am trying to understand the
benefits of using the bootstrap for small sample sets. To do this I created a normal
population and then picked 10 from the populations and applied both traditional
statistical methods and the Bootstrap (bcanon, 5000 bootstrap samples) to calculate a
95% confidence interval of on the mean. I saved the width of the confidence interval
and how many misses, repeated this 1000 times, and output the summary on the CI width
and the number of misses for each method (actual script below) . I had expected to see
about 5% of the CI to miss the actual mean. For the traditional method about 6%
missed but for the Bootstrap it was over 11%.
Traditional:
Min. 1st Qu. MedianMean 3rd Qu.Max.
0.4531 1.1460 1.3550 1.3690 1.5900 2.4330
[1] 0.062
Bootstrap:
Min. 1st Qu. MedianMean 3rd Qu.Max.
0.3661 0.9455 1.1290 1.1380 1.3220 2.0160
[1] 0.113
The bootstrap method consistently missed almost twice as many times as the traditional
method.
What am I missing? Is this the best that can be expected when working with a sample
set of only 10 pieces?
Although this experiment was done using a normal distribution, my real datasets would
be non-normal of an unknown distribution.
Thanks for any help,
Art Nuzzo
Motorola
***
library(bootstrap)
CImiss = function(M, lower, upper) {lower M || upper M }
CIr = function(lower, upper) {upper - lower}
C = c(); B = c() # CI Range
Ccov = 0; Bcov = 0 # Number of Ci Misses
cnt = 1000; # reps
x = rnorm(1) # create population
m = mean(x)
for (i in 1:cnt) {
s = sample(x,10,replace=F) # sample population
tresults = t.test(s)
attach(tresults)
C[i] = CIr(conf.int[1],conf.int[2])
if (CImiss(m,conf.int[1],conf.int[2])) Ccov = Ccov + 1
detach(tresults)
bcaresults - bcanon(s,5000,mean,alpha=c(.025,.975))
attach(bcaresults)
B[i] = CIr(confpoints[1,2],confpoints[2,2])
if (CImiss(m,confpoints[1,2],confpoints[2,2])) Bcov = Bcov + 1
detach(bcaresults)
}
print(summary (C))
print(Ccov/cnt)
print(summary (B))
print(Bcov/cnt)
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] updating via CRAN and http
Hello, thanks for the tips on updating packages for 1.8.0. The updating is a real problem for me, since I've to do it sort of manually using my web-browser or wget. I'm behind a firewall that requires http/ftp authentification (username and passwd) for every request it sends to a server outside our intranet. Therefore all the nice tools for automatic updating (cran, cpan ...) don't for me (I've tried). I understand that the non-paranoid rest of the world can't be bothered, but is there any intenstion to include such authentification into the update procedures of R? I think for ftp it's kind of tricky, but at least for http the authentification seems to be straight forward. kind regards, Arne __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Contrast specified with C() - R vs S-Plus problem
On Wed, 8 Oct 2003, Pascal A. Niklaus wrote: Thanks for the reply. Below is the solution and the S-Plus and R code that does the same (for documentation). I can't reproduce that in S-PLUS 6.1, and it is not as documented: In S-Plus 2000, C() complements the contrast matrix with orthogonal contrasts if only the first is given. And so does R. You just did not specify it correctly in R. The code you quote works in R as well, but your original example had four levels, not three. CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum Sq Mean Sq F value Pr(F) C(CO2, c(1, 0, -1)) 2 16.000 8.000 7.2593 0.004014 Residuals 21 23.143 1.102 Please do not make out your user errors to be S-R differences. Did you actually try your example in R? CO2 - factor(rep(c(A,C,E),each=8)) m - rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8) summary.aov(aov(m ~ C(CO2,c(1,0,-1 Df Sum of Sq Mean Sq F Value Pr(F) C(CO2, c(1, 0, -1)) 2 16.0 8.00 7.259259 0.004013532 Residuals 21 23.14286 1.102041 summary.lm(aov(m ~ C(CO2,c(1,0,-1 Coefficients: Value Std. Error t value Pr(|t|) (Intercept) 2.5000 0.214311.6667 0. C(CO2, c(1, 0, -1))1 -1. 0.2624-3.8103 0.0010 C(CO2, c(1, 0, -1))2 0. 0.3712 0. 1. Residual standard error: 1.05 on 21 degrees of freedom Multiple R-Squared: 0.4088 F-statistic: 7.259 on 2 and 21 degrees of freedom, the p-value is 0.004014 Correlation of Coefficients: (Intercept) C(CO2, c(1, 0, -1))1 C(CO2, c(1, 0, -1))1 0 C(CO2, c(1, 0, -1))2 0 0 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Rdinfo verbosity increased in 1.8.0?
Hi, my out-of-the-box installation of R-1.8.0 on Tru64 (OSFv5.1) results in an tremendously increased verbosity of Rdinfo. Is it really intended that Rdinfo complaints about minor mistakes like missing empty lines at the end given the very variable quality of Rd-pages? Is there a way to decrease Rdinfo's standard verbosity? Regards Michael -- Michael T. Mader Institute for Bioinformatics/MIPS, GSF Ingolstaedter Landstrasse 1 D-80937 Neuherberg 0049-89-3187-3576 XML can be a simplifying choice or a complicating one. There is a lot of hype surrounding it, but don't become a fashion victim by either adopting or rejecting it uncritically. Choose carefully and bear the KISS principle in mind. Eric Steven Raymond: The Art of Unix Programming, 2003 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Also, presumably is.na- could be redefined by the user for particular classes so if you got in the habit of setting NAs that way it would generalize better. --- Date: Wed, 8 Oct 2003 11:49:29 +0100 (BST) From: Prof Brian Ripley [EMAIL PROTECTED] I don't think it can ever `go wrong', but it can do things other than the user intends. The intention of is.na- is clearer, and so perhaps user error is less likely? That is the thinking behind the function, anyway. ___ No banners. No pop-ups. No kidding. Introducing My Way - http://www.myway.com __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
At Wednesday 03:06 PM 10/8/2003 +0200, Martin Maechler wrote: Your question has been answered by Achim and Peter Dalgaard (at least). Just a note: Using a[which(logic)] looks like a clumsy and inefficient way of writing a[ logic ] and I think you shouldn't propagate its use ... What then is the recommended way of treating an NA in the logical subset as a FALSE? (Or were you just talking about the given example, which didn't have this issue. However, you admonition seemed more general.) As in: x - 1:4 y - c(1,2,NA,4) x[y %% 2 == 0] [1] 2 NA 4 x[which(y %% 2 == 0)] [1] 2 4 Sometimes one might want the first result, but more usually, I want the second, and using which() seems a convenient way to get it. -- Tony Plate __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] fitdistr, mle's and gamma distribution
Are you interested in turning that into a monitor, processing each
day's data sequentially or even each entry as it arrived? If yes, you
may wish to evaluate the Foundations of Monitoring documents
downloadable from www.prodsyse.com. If you have any questions about
that, I might be able to help.
hope this helps. spencer graves
Lourens Olivier Walters wrote:
Thanks for the help, the wrapper function was very useful. I managed to
solve the problem using Spencer Graves' suggestion. I am analyzing the
interarrival times between HTTP packets on a campus network. The dataset
actually has more than 14 Million entries! It represents the traffic
generated by approximately 3000 users browsing the web for 30 days. I
have to be careful to always remove unused objects from my workspace,
but otherwise I have so far managed to cope with 512Mb of memory on a
Pentium 600Mhz.
Lourens
On Tue, 2003-09-30 at 23:25, Ben Bolker wrote:
PS. 11 MILLION entries??
On Tue, 30 Sep 2003, Ben Bolker wrote:
Spencer Graves's suggestion of using shape and scale parameters on a log
scale is a good one.
To do specifically what you want (check values for which the objective
function is called and see what happens) you can do the following
(untested!), which makes a local copy of dgamma that you can mess with:
dgamma.old - dgamma
dgamma - function(x,shape,rate,...) {
d - dgamma.old(x,shape,rate,...)
cat(shape,rate,d,\n)
return(d)
}
__
[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] updating via CRAN and http
Sorry, I didn' mean it the nasty way. I wouldn't have been surprised if the R-team had told me the authentification with the firewall is my problem (i.e. a special case that cannot be dealt with by th R-team). Yess, and off course I should have had a much closer lookk into the docu. Thanks again for the hint + please forgive! +regards, Arne -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: 08 October 2003 17:20 To: Muller, Arne PH/FR Cc: [EMAIL PROTECTED] Subject: Re: [R] updating via CRAN and http On Wed, 8 Oct 2003 [EMAIL PROTECTED] wrote: Hello, thanks for the tips on updating packages for 1.8.0. The updating is a real problem for me, since I've to do it sort of manually using my web-browser or wget. I'm behind a firewall that requires http/ftp authentification (username and passwd) for every request it sends to a server outside our intranet. Therefore all the nice tools for automatic updating (cran, cpan ...) don't for me (I've tried). I understand that the non-paranoid rest of the world can't be bothered, but is there any intenstion to include such authentification into the update procedures of R? I think for ftp it's kind of tricky, but at least for http the authentification seems to be straight forward. It's available for http: see ?download.file, and you can even configure that to use wget. Your comments are very much misplaced: we *have* bothered to provide the facilities for you. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] updating via CRAN and http
On Wed, 8 Oct 2003 [EMAIL PROTECTED] wrote: Hello, thanks for the tips on updating packages for 1.8.0. The updating is a real problem for me, since I've to do it sort of manually using my web-browser or wget. I'm behind a firewall that requires http/ftp authentification (username and passwd) for every request it sends to a server outside our intranet. Therefore all the nice tools for automatic updating (cran, cpan ...) don't for me (I've tried). I understand that the non-paranoid rest of the world can't be bothered, but is there any intenstion to include such authentification into the update procedures of R? I think for ftp it's kind of tricky, but at least for http the authentification seems to be straight forward. It's available for http: see ?download.file, and you can even configure that to use wget. Your comments are very much misplaced: we *have* bothered to provide the facilities for you. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] using split.screen() in Sweave
Dear R and sweave users A further problem, which I couldn't resolve, using the manual: In R I use the split.screen command to put e.g. two timecourses one above the other into one plot: split.screen(c(2,1)) screen(1) plot(stick,type='h', col=red,lwd=2) screen(2) plot(deconvolution.amplitude,type='h',col=blue,lwd=2) Is there a similar way, doing this in Sweave? many thanks Christoph -- Christoph Lehmann [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Lattice cloud() funtion bug in R1.8.0beta
Cloud() function does not display anything with R1.8.0beta in WindowsXP ... Does any one noticed this ? others functions from lattice seem working properly. does it work in the final 1.8.0 for windows ? __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] updating via CRAN and http
On Wed, Oct 08, 2003 at 04:25:13PM +0200, [EMAIL PROTECTED] wrote: Hello, thanks for the tips on updating packages for 1.8.0. The updating is a real problem for me, since I've to do it sort of manually using my web-browser or wget. I'm behind a firewall that requires http/ftp authentification (username and passwd) for every request it sends to a server outside our intranet. Therefore all the nice tools for automatic updating (cran, cpan ...) don't for me (I've tried). Please try to read the help pages, FAQ or mailing list archives. Trust me that you are not the only corporate user around here. This question comes up really frequently. I understand that the non-paranoid rest of the world can't be bothered, but is there any intenstion to include such authentification into the update procedures of R? I think for ftp it's kind of tricky, but at least for http the authentification seems to be straight forward. Years ago a small patch of mine was integrated which allowed for wget as a download methid. Set e.g. options(download.file.method=wget) options(CRAN=http://cran.us.r-project.org;) in ~/.Rprofile. Get yourself a copy of wget, set it up with a suitable ~/.wgetrc as e.g. per http_proxy=http://some.where.on.the.intra.net:8080 proxy=on proxy-user=abcdefgh proxy-passwd= and make sure you change the file to read-only. Test it with wget http://www.google.com and if that works you will be fine for downloads within R too, including package upgrades and normal http/ftp access as e.g. used by get.hist.quote() in the tseries package. This worked reliably for me in various environments, locations and operating system, including several Win* variants. Hth, Dirk -- Those are my principles, and if you don't like them... well, I have others. -- Groucho Marx __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
Here are some different ways of doing this. Don't know whether any could be considered superior to the others. # y[x==5] regarding NAs in x as not matching x - c(5, NA, 7, 5, NA, 3) y - c(1, 2, 3, 4, 5, 6) subset(y,x==5) y[x %in% 5] y[x %in% c(5)] y[which(x==5)] --- Date: Wed, 08 Oct 2003 08:58:32 -0600 From: Tony Plate [EMAIL PROTECTED] At Wednesday 03:06 PM 10/8/2003 +0200, Martin Maechler wrote: Your question has been answered by Achim and Peter Dalgaard (at least). Just a note: Using a[which(logic)] looks like a clumsy and inefficient way of writing a[ logic ] and I think you shouldn't propagate its use ... What then is the recommended way of treating an NA in the logical subset as a FALSE? (Or were you just talking about the given example, which didn't have this issue. However, you admonition seemed more general.) As in: x - 1:4 y - c(1,2,NA,4) x[y %% 2 == 0] [1] 2 NA 4 x[which(y %% 2 == 0)] [1] 2 4 Sometimes one might want the first result, but more usually, I want the second, and using which() seems a convenient way to get it. -- Tony Plate ___ No banners. No pop-ups. No kidding. Introducing My Way - http://www.myway.com __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Installing GLMMGibbs problems
Dear all; Installing the GLMMGibbs package to my Solaris Unix box, I got an compiling error: ars.c:497:10: missing terminating character ars.c: In function `dump_arse': ars.c:498: error: parse error before mylesj . The compiling error was reported to the list on Jul 3, 2003. According to Prof. Brain Ripley this is a known problem with the package and gcc 3.3, and the maintainer has been informed ... Because I was unable to reach the maintainer by email (my email to [EMAIL PROTECTED] on Oct. 2, bounced), I saved the GLMMGibbs package from the tmp directory, edited the offensive lines in src/ars.c and saved the file. Following a test compiling the edited ars.c without error, I re-tarred and gzipped the saved GLMMGibbs package, stored the new GLMMGibbs_0.5-1.tar.gz in a directory Rsources. I expected to install the package from the saved directory (similar to install package(s) from local zip files... on W2K) and used the install.packages(): install.packages(GLMMGibbs, lib=/PATH/Rsources/R-1.7.1/library, destdir = /PATH/Rsources) Instead of using the local zip file from the Rsources directory, R downloaded the GLMMGibbs_0.5-1.tar.gz from CRAN: trying URL `http://cran.r-project.org/src/contrib/GLMMGibbs_0.5-1.tar.gz' Content type `application/x-tar' length 346260 bytes opened URL . and produced the same compiling error. The question becomes how to install packages from local zip files on Solaris or Linux platforms? In the absence of the original maintainer, could someone in R team edit the source file ars.c in compliance with gcc 3.3? TIA, Richard __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] using split.screen() in Sweave
Christoph Lehmann wrote: Dear R and sweave users A further problem, which I couldn't resolve, using the manual: In R I use the split.screen command to put e.g. two timecourses one above the other into one plot: split.screen(c(2,1)) screen(1) plot(stick,type='h', col=red,lwd=2) screen(2) plot(deconvolution.amplitude,type='h',col=blue,lwd=2) Is there a similar way, doing this in Sweave? I've never used split.screen. Check out ?par, and read the entries for mfcol and mfrow. Also check out ?layout. HTH Jason -- Indigo Industrial Controls Ltd. http://www.indigoindustrial.co.nz 64-21-343-545 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
Achim Zeileis wrote: On Wednesday 08 October 2003 11:27, Thomas Bock wrote: ... I can not understand why the expression in the subject does not work correct: dcrn[which(fn == inve[2])] numeric(0) inve[2] [1] 406.7 ... 1.) `==' comparisons have a certain tolerance 2.) the print output is not necessarily precisely your number Instead of using `==' you should use a comparison with a certain tolerance you can specify... I usually specify... tol - sqrt(.Machine$double.eps) dcrn[(fn - inve[2]) tol] See ?.Machine for details. Cheers Jason -- Indigo Industrial Controls Ltd. http://www.indigoindustrial.co.nz 64-21-343-545 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] binomial glm warnings revisited
This seems to me to be a special case of the general problem of a parameter on a boundary. Another example is the case of a variance component that is zero. For this latter problem, Pinhiero and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer, sec. 2.4.1) present simulation results showing that a 50-50 mixture of chi-square(0) and chi-square(1), for example, provide an excellent approximation to the actual sampling distribution of the 2*log(likelihood ratio). Recent discussions of this and related questions on this list and elsewhere produced the following list of articles that may be helpful: Donald Andrews (2001) Testing When a Parameter In on the Boundary of the Maintained Hypothesis, Econometrica, 69: 683-734. Donald Andrews (2000) Inconsistency of the Bootstrap When a Parameter Is on the Boundary of the Parameter Space, Econometrica, 68: 388-405. Donald Andrews (1999) Estimation When a Parameter Is on a Boundary, Econometrica, 67: 1341-1383. Rousseeuw, P. J. and Christmann, A. (2003) Robustness against separations and outliers in logistic regression, Computational Statistics Data Analysis, Vol. 43, pp. 315-332 ### Unfortunately, I have not had time to review these, so I can't comment further. hope this helps. spencer graves Tord Snall wrote: Dear all, Last autumn there was some discussion on the list of the warning Warning message: fitted probabilities numerically 0 or 1 occurred in: (if (is.empty.model(mt)) glm.fit.null else glm.fit)(x = X, y = Y, when fitting binomial GLMs with many 0 and few 1. Parts of replies: You should be able to tell which coefficients are infinite -- the coefficients and their standard errors will be large. When this happens the standard errors and the p-values reported by summary.glm() for those variables are useless. My guess is that the deviances and coefficients are entirely ok. I'd expect that problems in the general area that Thomas mentions to reveal themselves as a failure to converge. I have this problem with my data. In a GLM, I have 269 zeroes and only 1 one: summary(dbh) Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 0.1659 3.8781 0.0430.966 dbh -0.5872 0.5320 -1.1040.270 drop1(dbh, test = Chisq) Single term deletions Model: MPext ~ dbh Df Deviance AIC LRT Pr(Chi) none 9.9168 13.9168 dbh 1 13.1931 15.1931 3.2763 0.07029 . I now wonder, is the drop1() function output 'reliable'? If so, is then the estimates from MASS confint() also 'reliable'? It gives the same warning. Waiting for profiling to be done... 2.5 % 97.5 % (Intercept) -6.503472 -0.77470556 abund -1.962549 -0.07496205 There were 20 warnings (use warnings() to see them) Thanks in advance for your reply. Sincerely, Tord --- Tord Snäll Avd. f växtekologi, Evolutionsbiologiskt centrum, Uppsala universitet Dept. of Plant Ecology, Evolutionary Biology Centre, Uppsala University Villavägen 14 SE-752 36 Uppsala, Sweden Tel: 018-471 28 82 (int +46 18 471 28 82) (work) Tel: 018-25 71 33 (int +46 18 25 71 33) (home) Fax: 018-55 34 19 (int +46 18 55 34 19) (work) E-mail: [EMAIL PROTECTED] Check this: http://www.vaxtbio.uu.se/resfold/snall.htm! __ [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] Why does a[which(b == c[d])] not work?
Whoops. Hit send too quickly. Jason Turner wrote: tol - sqrt(.Machine$double.eps) dcrn[(fn - inve[2]) tol] that should be dcrn[abs(fn - inve[2]) tol] -- Indigo Industrial Controls Ltd. http://www.indigoindustrial.co.nz 64-21-343-545 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] 2 questions regarding base-n and identifing digits
Dear listers, I have two questions: (1) Is there a way in R to change the base-n of the calculations. I wnat to run some calculations either in binary (base-2) or base-4. Is there a way to specify that in R - to chnage from the decimal? (2) I also want to extract the digits from a larger number and store them as a vector. I could do it through converting top string, parsing the string and converting back. Is there another way. It would look something like that: a-1234 a [1] 1234 b-foo(a) b [,1] [,2] [,3] [,4] [1,]1234 Thanks! Andrej _ Andrej Kveder, M.A. researcher Institute of Medical Sciences SRS SASA; Novi trg 2, SI-1000 Ljubljana, Slovenia phone: +386 1 47 06 440 fax: +386 1 42 61 493 [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] binomial glm warnings revisited
Spencer Graves [EMAIL PROTECTED] writes: This seems to me to be a special case of the general problem of a parameter on a boundary. Umm, no... I have this problem with my data. In a GLM, I have 269 zeroes and only 1 one: I don't think that necessarily gets you a parameter estimate on the boundary. Only if the single 1 is smaller or bigger than all the others should that happen. summary(dbh) Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 0.1659 3.8781 0.0430.966 dbh -0.5872 0.5320 -1.1040.270 drop1(dbh, test = Chisq) Single term deletions Model: MPext ~ dbh Df Deviance AIC LRT Pr(Chi) none 9.9168 13.9168 dbh 1 13.1931 15.1931 3.2763 0.07029 . I now wonder, is the drop1() function output 'reliable'? If so, is then the estimates from MASS confint() also 'reliable'? It gives the same warning. (Intercept) -6.503472 -0.77470556 abund -1.962549 -0.07496205 There were 20 warnings (use warnings() to see them) During profiling, you may be pushing one of the parameter near the extremes and get a model where the fitted p's are very close to 0/1. That's not necessarily a sign of unreliability -- the procedure is to set one parameter to a sequence of fixed values and optimize over the other, and it might just be the case that the optimizations have been wandering a bit far from the optimum. (I'd actually be more suspicious about the fact that the name of the predictor suddenly changed) However, if you have only one 1 you are effectively asking whether one observation has a different mean than the other 269, and you have to consider the sensitivity to the distribution of the predictor. As far as I can see, you end up with the test of the null hypothesis beta==0 being essentially equivalent to a two sample t test between the mean of the 0 group and that of the 1 group, so with only one observation in one of the groups, the normal approximation of the test hinges quite strongly on a normal distribution of the predictor itself. -- 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] binomial glm warnings revisited
Thanks, Peter: You are absolutely correct. Thanks again for the correction. Spencer Graves Peter Dalgaard BSA wrote: Spencer Graves [EMAIL PROTECTED] writes: This seems to me to be a special case of the general problem of a parameter on a boundary. Umm, no... I have this problem with my data. In a GLM, I have 269 zeroes and only 1 one: I don't think that necessarily gets you a parameter estimate on the boundary. Only if the single 1 is smaller or bigger than all the others should that happen. summary(dbh) Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) 0.1659 3.8781 0.0430.966 dbh -0.5872 0.5320 -1.1040.270 drop1(dbh, test = Chisq) Single term deletions Model: MPext ~ dbh Df Deviance AIC LRT Pr(Chi) none 9.9168 13.9168 dbh 1 13.1931 15.1931 3.2763 0.07029 . I now wonder, is the drop1() function output 'reliable'? If so, is then the estimates from MASS confint() also 'reliable'? It gives the same warning. (Intercept) -6.503472 -0.77470556 abund -1.962549 -0.07496205 There were 20 warnings (use warnings() to see them) During profiling, you may be pushing one of the parameter near the extremes and get a model where the fitted p's are very close to 0/1. That's not necessarily a sign of unreliability -- the procedure is to set one parameter to a sequence of fixed values and optimize over the other, and it might just be the case that the optimizations have been wandering a bit far from the optimum. (I'd actually be more suspicious about the fact that the name of the predictor suddenly changed) However, if you have only one 1 you are effectively asking whether one observation has a different mean than the other 269, and you have to consider the sensitivity to the distribution of the predictor. As far as I can see, you end up with the test of the null hypothesis beta==0 being essentially equivalent to a two sample t test between the mean of the 0 group and that of the 1 group, so with only one observation in one of the groups, the normal approximation of the test hinges quite strongly on a normal distribution of the predictor itself. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] winedt or another editor for linux
Have got anybody experience using winedt and R in Linux - perhaps with wine Or exist another editor with the ability to parse r code into R in Linux, because k-edit etc. using only syntax-highlithing ??? many thanks, christian [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] winedt or another editor for linux
On Wednesday 08 October 2003 22:02, Christian Schulz wrote: Have got anybody experience using winedt and R in Linux - perhaps with wine Or exist another editor with the ability to parse r code into R in Linux, because k-edit etc. using only syntax-highlithing ??? I guess you have looked at (X)Emacs in connection with ESS (Emacs Speaks Statistics)? That is of course much more than an editor and not only available for Linux... hth, Z many thanks, christian [[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] 2 questions regarding base-n and identifing digits
From: Andrej Kveder [mailto:[EMAIL PROTECTED] Dear listers, I have two questions: (1) Is there a way in R to change the base-n of the calculations. I wnat to run some calculations either in binary (base-2) or base-4. Is there a way to specify that in R - to chnage from the decimal? I don't have any good answers, but just this: Arithmetics in any base is the same. I believe you just need a way to convert decimal to other bases and printed out as such. E.g., 1 + 5 = 6 = 110 in binary. You just need something that returns 110 when given 6. Personally I'd write such a thing in C because I don't know any better. (2) I also want to extract the digits from a larger number and store them as a vector. I could do it through converting top string, parsing the string and converting back. Is there another way. It would look something like that: a-1234 a [1] 1234 b-foo(a) b [,1] [,2] [,3] [,4] [1,]1234 Would the following do what you want? a - c(1234, 5678, 90) foo - function(x) strsplit(as.character(x), ) foo(a) [[1]] [1] 1 2 3 4 [[2]] [1] 5 6 7 8 [[3]] [1] 9 0 HTH, Andy Thanks! Andrej _ Andrej Kveder, M.A. researcher Institute of Medical Sciences SRS SASA; Novi trg 2, SI-1000 Ljubljana, Slovenia phone: +386 1 47 06 440 fax: +386 1 42 61 493 [[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: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Simon Fear [EMAIL PROTECTED] suggested that a-a a-NA mode(a) [1] logical a-a is.na(a) - T mode(a) [1] character might be a relevant difference between assigning NA and using is.na. But the analogy is flawed: is.na(x) - operates on the _elements_ of x, while x - affects the variable x. When you assign NA to _elements_ of a vector, the mode does not change: a - a is.na(a) - TRUE mode(a) [1] character b - b b[TRUE] - NA mode(b) [1] character c - c c[1] - NA mode(c) [1] character __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Concerning x[i] - NA vs is.na(x[i]) - TRUE
Brian Ripley wrote:
I don't think it can ever `go wrong', but it can do things other
than the user intends.
If the user writes x[i] - NA, the user has clearly indicated his intention
that the i element(s) of x should become NA. There isn't any clearer way to
say that. The only way it could ever do something other than the user
intends is if the mode of x changes or the selected elements are set to
something other than NA.
The ?NA help page *hints* that this might be the case, but does not give
an example.
The question remains, *WHAT* can x[i]-NA do that might be other than
what the user intends? An example (especially one added to the ?NA help)
would be very useful.
The intention of is.na- is clearer,
I find this extremely puzzling. x[i] - NA is an extremely clear and
obvious way of saying I want the i element(s) of x to become NA.
is.na(x) - ... is not only an indirect way of doing this, it is a way
which is confusing and error-prone.
Bear in mind that one way of implementing something is is.na() would be
to associate a bit with each element of a vector; is.na() would test and
is.na-() would set that bit. It would be possible to have a language
exactly like R -except- that
x - 1
is.na(x) - TRUE
x
= NA
is.na(x) - FALSE
x
= 1
would work. The very existence of an is.na- which accepts a logical
vector containing FALSE as well as TRUE strongly suggests this. But it
doesn't work like that. As I've pointed out,
is.logical(m) length(m) == length(x) done{is.na(x) - m}
= identical(is.na(x), m)
is the kind of behaviour that has been associated with well-behaved
sinister function calls for several decades, and yet this is not a fact
about R.
and so perhaps user error is less likely?
I see no reason to believe this; the bad behaviour of is.na- surely
makes user error *more* likely rather than *less* likely.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] (no subject)
Good afternoon, We currently have R installed on our HP Superdome. However, we are getting ready to migrate from RISC to Itanium 2 chips running HP-UX (not Linux). Does the latest version of R run on HP-UX Itanium 2? Any information would be greatly appreciated. Michael ___ R. Michael Sheetz, PhD High Performance Computing Group 326 A McVey Hall University of Kentucky Lexington, KY 40506 Phone: (859) 257-2900 ext 289 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Why does a[which(b == c[d])] not work?
Achim Zeileis [EMAIL PROTECTED] wrote: R x - 406.7 + 1e-20 R x [1] 406.7 R x == 406.7 [1] TRUE that is 1.) `==' comparisons have a certain tolerance No, all.equal() supports tolerance, == does not. Consider .Machine$double.eps [1] 2.220446e-16 That is, the smallest relative difference that can be represented on my (fairly typical IEEE-conforming) machine is 2 in 10 to the 16th. 1e-20/406.7 is 2.458815e-23, which is a factor of 10 million too small a relative difference for the hardware to be able to represent. So in x - 406.7 + 1e-20 the value of x is identical to 406.7 in every bit. Instead of using `==' you should use a comparison with a certain tolerance you can specify... Such as ?all.equal __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Sorry! Previous subject was cell means model in LME
Hi everybody I want to specify the contrasts to build a cell means model on LME (each coefficient is the mean value of the independent variable for the specific category of a factor variable) when there are several factors as fixed effect in the model and also interactions between them. Can anybody give me a hint on how to do accomplish this? Also, how can I override the default Contr.Treatment for linear models? I tried removing the intercept but this gave me the mean values just for the first factor included in the model and then used contrast treatment for all the other factors. Thanks for your help! Francisco _ Instant message in style with MSN Messenger 6.0. Download it now FREE! http://msnmessenger-download.com __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: is.na(v)-b (was: Re: [R] Beginner's query - segmentation fault)
Tongue in cheek But surely is.na(x) - is.na(x) is clearer than x[is.na(x)] - NA (neither of which is a no-op). /Tongue in cheek __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] R-1.8.0 is released
I'd just like to thank the R-Core team for the new version (all OS's), and especially Stefano Iacus for his work on the old Carbon MacOS port and his continuing work on the OS X port. Cheers, Simon. Simon Blomberg, PhD Depression Anxiety Consumer Research Unit Centre for Mental Health Research Australian National University http://www.anu.edu.au/cmhr/ [EMAIL PROTECTED] +61 (2) 6125 3379 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Re: Mixed effects with nlme
Hi, R-users: Last week I send a request for help to this list. I have receive until now two kindly responses from Spencer Graves and Renauld Lancelot. They both point interesting things to fit an adequate model to my data but unfortunately it persists without a satisfactory solution. I propose again the same problem but using with a different dataset (Assay), taken from Pinheiro and Bates' book on page 163, that is relevant with crossed random effects. I have fitted the same model (p. 165) fmAssay - lme(logDens ~ sample * dilut, Assay, random=pdBlocked(list(, pdIdent(~1), pdIdent(~sample-1),pdIdent(~dilut-1))) ) and the results with anova function (p. 166) are numDF denDF F-value p-value (Intercept) 129 537.6294 .0001 sample 529 11.2103 .0001 dilut429 420.5458 .0001 sample:dilut2029 1.6072 0.1192 The problem is that with this approach one obtains correct F-values, but using a common residual term for DenDF that is a combination of (Block + Block:sample + Block:dilut). Then probability values are different to that obtained when we used the classical AOV funtion to fit the same model, because in this case each term is mapped with a error term (so sample uses Block:sample, dilut uses Block:dilut, and sample:dilut uses Block:sample:dilut): mod-aov(logDens ~ sample*dilut + Error(Block+Block/sample+Block/dilut), data=Assay) summary(mod) Error: Block DfSum Sq Mean Sq F value Pr(F) Residuals 1 0.0083115 0.0083115 Error: Block:sample Df Sum Sq Mean Sq F value Pr(F) sample 5 0.276153 0.055231 11.213 0.009522 Residuals 5 0.024627 0.004925 Error: Block:dilut Df Sum Sq Mean Sq F valuePr(F) dilut 4 3.7491 0.9373 420.79 1.684e-05 Residuals 4 0.0089 0.0022 Error: Within Df Sum Sq Mean Sq F value Pr(F) sample:dilut 20 0.055525 0.002776 1.6069 0.1486 Residuals20 0.034555 0.001728 Obviously, the interest on linear mixed effects is open with the possibility of modeling with correlated and/or heterocedastic errors, and this end cannot be pursued with AOV function. Summarizing, the problem is that I have not found a way to obtain with NLME the same results (DF, F-ratios and probabilities) that I get with AOV and multistratum errors. Is this an inconvenience of program?, probably due to the impossibility to use multiple nested arguments as random(~1|Block/sample+dilut) or random(~1|Block/sample*dilut) SAS MIXED can also fit these data and obtain correct results by means of a combination of random and repeated arguments: model = sample dilut sample*dilut; random = Block sample*Block dilut*Block; repeated /type=cs Sub=Block; Type 3 Tests of Fixed Effects Num Den Effect DF DFF ValuePr F sample 5 5 11.210.0095 dilut 4 4 420.79.0001 sample*dilut 20 20 1.610.1486 May be possible obtain the same results with NLME function? I would appreciate any kind of help. Best regards, Manuel Ato University of Murcia Spain e-mail: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] 2 questions regarding base-n and identifing digits
On 08-Oct-03 Liaw, Andy wrote: From: Andrej Kveder [mailto:[EMAIL PROTECTED] I have two questions: (1) Is there a way in R to change the base-n of the calculations. I wnat to run some calculations either in binary (base-2) or base-4. Is there a way to specify that in R - to chnage from the decimal? I don't have any good answers, but just this: Arithmetics in any base is the same. I believe you just need a way to convert decimal to other bases and printed out as such. E.g., 1 + 5 = 6 = 110 in binary. You just need something that returns 110 when given 6. Personally I'd write such a thing in C because I don't know any better. Indeed! A pity that R's sprintf (though described as a wrapper for C's sprintf) does not accept the C format conversion specifiers %X (for hexadecimal) or, better stil, %o (for octal), because then Andy's strsplit usage below could give you one of 16 or 8 characters which could be used to index a lookup of a bit pattern! Ted. Would the following do what you want? a - c(1234, 5678, 90) foo - function(x) strsplit(as.character(x), ) foo(a) [[1]] [1] 1 2 3 4 [[2]] [1] 5 6 7 8 [[3]] [1] 9 0 E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 167 1972 Date: 09-Oct-03 Time: 00:09:16 -- XFMail -- __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Scoping rules
Dear List members:
I'm using R1.7.1 (Windows 2000) and having difficulty with scoping.
I've studied the FAQ and on-line manuals and think I have identified
the
source of my difficulty, but cannot work out the solution.
For the purposes of illustration. I have three functions as defined
below:
fnA - function(my.x)
{
list(first=as.character(substitute(my.x)), second=sqrt(my.x))
}
fnB - function(AA)
{
tmp.x - get(AA$first)
tmp.x^2
}
fnC - function()
{
x - 1:2
y - fnA(x)
z - fnB(y)
c(x,y,z)
}
fnA() has a vector as an argument and returns the name of the vector
and the square root of its elements in a list. fn(B) takes the result
of fn(A) as its argument, gets the appropriate vector and computes the
square of its elements. These work fine when called at the command
line.
fnC() defines a local vector x and calls fnA() which operates on this
vector. Then fnB() is called, but it operates on a global vector x in
GlobalEnv (or returns an error is x doesn't exist there) - but I want
it
to operate on the local vector.
I think this is related to the enclosing environment of all three
functions being GlobalEnv (since they were created at the command
line),
but the parent environment of fnB() being different when invoked from
within fnC().
My questions:
1 Have I correctly understood the issue ?
2 How do I make fnB() operate on the local vector rather than the
global one ?
3 And, as an aside, I have used as.character(substitute(my.x)) to
pass
the name - but deparse(substitute(my.x)) also works. Is there any
reason to prefer one over the other?
Thank you ...
Peter Alspach
__
The contents of this e-mail are privileged and/or confidential to the
named recipient and are not to be used by any other person and/or
organisation. If you have received this e-mail in error, please notify
the sender and delete all material pertaining to this e-mail.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] 2 questions regarding base-n and identifing digits
Check out: http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/decodeencode/decodeencode.rev --- Date: Wed, 8 Oct 2003 21:39:50 +0200 From: Andrej Kveder [EMAIL PROTECTED] Dear listers, I have two questions: (1) Is there a way in R to change the base-n of the calculations. I wnat to run some calculations either in binary (base-2) or base-4. Is there a way to specify that in R - to chnage from the decimal? (2) I also want to extract the digits from a larger number and store them as a vector. I could do it through converting top string, parsing the string and converting back. Is there another way. It would look something like that: a-1234 a [1] 1234 b-foo(a) b [,1] [,2] [,3] [,4] [1,] 1 2 3 4 Thanks! Andrej _ Andrej Kveder, M.A. researcher Institute of Medical Sciences SRS SASA; Novi trg 2, SI-1000 Ljubljana, Slovenia phone: +386 1 47 06 440 fax: +386 1 42 61 493 ___ No banners. No pop-ups. No kidding. Introducing My Way - http://www.myway.com __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Scoping rules
It seems like you want in fnB
get(AA$first, envir = parent.frame(1))
but I'm entirely clear on why your original function doesn't work. My
understanding was that get() should search through the parent frames.
-roger
Peter Alspach wrote:
Dear List members:
I'm using R1.7.1 (Windows 2000) and having difficulty with scoping.
I've studied the FAQ and on-line manuals and think I have identified
the
source of my difficulty, but cannot work out the solution.
For the purposes of illustration. I have three functions as defined
below:
fnA - function(my.x)
{
list(first=as.character(substitute(my.x)), second=sqrt(my.x))
}
fnB - function(AA)
{
tmp.x - get(AA$first)
tmp.x^2
}
fnC - function()
{
x - 1:2
y - fnA(x)
z - fnB(y)
c(x,y,z)
}
fnA() has a vector as an argument and returns the name of the vector
and the square root of its elements in a list. fn(B) takes the result
of fn(A) as its argument, gets the appropriate vector and computes the
square of its elements. These work fine when called at the command
line.
fnC() defines a local vector x and calls fnA() which operates on this
vector. Then fnB() is called, but it operates on a global vector x in
GlobalEnv (or returns an error is x doesn't exist there) - but I want
it
to operate on the local vector.
I think this is related to the enclosing environment of all three
functions being GlobalEnv (since they were created at the command
line),
but the parent environment of fnB() being different when invoked from
within fnC().
My questions:
1 Have I correctly understood the issue ?
2 How do I make fnB() operate on the local vector rather than the
global one ?
3 And, as an aside, I have used as.character(substitute(my.x)) to
pass
the name - but deparse(substitute(my.x)) also works. Is there any
reason to prefer one over the other?
Thank you ...
Peter Alspach
__
The contents of this e-mail are privileged and/or confidential to the
named recipient and are not to be used by any other person and/or
organisation. If you have received this e-mail in error, please notify
the sender and delete all material pertaining to this e-mail.
__
[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] double precision
Hi, I am new to R, but have extensive experience in matlab. I have searched on the web and in Venabels Ripley book but I was unable to find the equivalent of the eps function in matlab. eps returns the Floating point relative accuracy. thanks __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] double precision
Check out ?.Machine -roger christoff pale wrote: Hi, I am new to R, but have extensive experience in matlab. I have searched on the web and in Venabels Ripley book but I was unable to find the equivalent of the eps function in matlab. eps returns the Floating point relative accuracy. thanks __ [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] Debian packages of 1.8.0 available
Debian packages of R 1.8.0 were uploaded earlier for i386. Binaries for alpha, ia64, hppa and powerpc are already in the package pool; the arm, mipsel and s390 architecture have their built completled and will be added to the pool shortly while the remaining architecures should follow over the next few days. Similar to prior practice, we will probably make a testing release (for i386) available at CRAN as well. Regards, Dirk -- Those are my principles, and if you don't like them... well, I have others. -- Groucho Marx __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Scoping rules
On Wed, 8 Oct 2003, Roger D. Peng wrote:
It seems like you want in fnB
get(AA$first, envir = parent.frame(1))
but I'm entirely clear on why your original function doesn't work. My
understanding was that get() should search through the parent frames.
No, get() searches through the enclosing frames, like ordinary variable
lookup. The only way to get dynamic scope is to specify it explicitly
-thomas
-roger
Peter Alspach wrote:
Dear List members:
I'm using R1.7.1 (Windows 2000) and having difficulty with scoping.
I've studied the FAQ and on-line manuals and think I have identified
the
source of my difficulty, but cannot work out the solution.
For the purposes of illustration. I have three functions as defined
below:
fnA - function(my.x)
{
list(first=as.character(substitute(my.x)), second=sqrt(my.x))
}
fnB - function(AA)
{
tmp.x - get(AA$first)
tmp.x^2
}
fnC - function()
{
x - 1:2
y - fnA(x)
z - fnB(y)
c(x,y,z)
}
fnA() has a vector as an argument and returns the name of the vector
and the square root of its elements in a list. fn(B) takes the result
of fn(A) as its argument, gets the appropriate vector and computes the
square of its elements. These work fine when called at the command
line.
fnC() defines a local vector x and calls fnA() which operates on this
vector. Then fnB() is called, but it operates on a global vector x in
GlobalEnv (or returns an error is x doesn't exist there) - but I want
it
to operate on the local vector.
I think this is related to the enclosing environment of all three
functions being GlobalEnv (since they were created at the command
line),
but the parent environment of fnB() being different when invoked from
within fnC().
My questions:
1 Have I correctly understood the issue ?
2 How do I make fnB() operate on the local vector rather than the
global one ?
3 And, as an aside, I have used as.character(substitute(my.x)) to
pass
the name - but deparse(substitute(my.x)) also works. Is there any
reason to prefer one over the other?
Thank you ...
Peter Alspach
__
The contents of this e-mail are privileged and/or confidential to the
named recipient and are not to be used by any other person and/or
organisation. If you have received this e-mail in error, please notify
the sender and delete all material pertaining to this e-mail.
__
[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
Thomas Lumley Assoc. Professor, Biostatistics
[EMAIL PROTECTED] University of Washington, Seattle
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Scoping rules
On Wed, 8 Oct 2003, Roger D. Peng wrote:
It seems like you want in fnB
get(AA$first, envir = parent.frame(1))
but I'm entirely clear on why your original function doesn't work. My
understanding was that get() should search through the parent frames.
Where did you get that idea from? Argument `inherits' to get() defaults
to TRUE and is defined as
inherits: should the enclosing frames of the environment be inspected?
Note `enclosing', not `parent'. So the normal R scope rules apply when
looking for an object from the frame of fnB, which as its environment (aka
enclosing frame) is the workspace (.GlobalEnv) is to look in the local
frame and then along the search path.
[This does seem a fairly common misconception, so if you do have any idea
in which document it arises it would be good to know.]
Peter Alspach wrote:
Dear List members:
I'm using R1.7.1 (Windows 2000) and having difficulty with scoping.
I've studied the FAQ and on-line manuals and think I have identified
the
source of my difficulty, but cannot work out the solution.
For the purposes of illustration. I have three functions as defined
below:
fnA - function(my.x)
{
list(first=as.character(substitute(my.x)), second=sqrt(my.x))
}
fnB - function(AA)
{
tmp.x - get(AA$first)
tmp.x^2
}
fnC - function()
{
x - 1:2
y - fnA(x)
z - fnB(y)
c(x,y,z)
}
fnA() has a vector as an argument and returns the name of the vector
and the square root of its elements in a list. fn(B) takes the result
of fn(A) as its argument, gets the appropriate vector and computes the
square of its elements. These work fine when called at the command
line.
fnC() defines a local vector x and calls fnA() which operates on this
vector. Then fnB() is called, but it operates on a global vector x in
GlobalEnv (or returns an error is x doesn't exist there) - but I want
it
to operate on the local vector.
I am not sure what you really want to do here, but R works best if you
pass functions an object to work on, and not the name of an object.
Normally this sort of thing is best avoided, but if it is really needed it
is normally simpler to find the object in the calling function (your fnC).
I think this is related to the enclosing environment of all three
functions being GlobalEnv (since they were created at the command
line),
but the parent environment of fnB() being different when invoked from
within fnC().
My questions:
1 Have I correctly understood the issue ?
2 How do I make fnB() operate on the local vector rather than the
global one ?
3 And, as an aside, I have used as.character(substitute(my.x)) to
pass
the name - but deparse(substitute(my.x)) also works. Is there any
reason to prefer one over the other?
deparse() is preferred. One subtle reason is what happens with very long
expressions (and note deparse may give more than one line of output),
but the main reason is what happens with expressions such as calls:
fn1 - function(x) as.character(substitute(x))
fn2 - function(x) deparse(substitute(x))
fn1(log(x))
[1] log x
fn2(log(x))
[1] log(x)
It is normally dangerous to use get() on the result of
deparse(substitute()), as the latter may be the character representation
of an expression and not a simple name.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Scoping rules
I think I misinterpreted this section of the help file for get():
If `inherits' is `FALSE', only the first frame of the specified
environment is inspected. If `inherits' is `TRUE', the search is
continued up through the parent frames until a bound value of the
right mode is found.
Should this read parent environments instead of parent frames? I'm
taking this from R 1.7.1.
-roger
Prof Brian Ripley wrote:
On Wed, 8 Oct 2003, Roger D. Peng wrote:
It seems like you want in fnB
get(AA$first, envir = parent.frame(1))
but I'm entirely clear on why your original function doesn't work. My
understanding was that get() should search through the parent frames.
Where did you get that idea from? Argument `inherits' to get() defaults
to TRUE and is defined as
inherits: should the enclosing frames of the environment be inspected?
Note `enclosing', not `parent'. So the normal R scope rules apply when
looking for an object from the frame of fnB, which as its environment (aka
enclosing frame) is the workspace (.GlobalEnv) is to look in the local
frame and then along the search path.
[This does seem a fairly common misconception, so if you do have any idea
in which document it arises it would be good to know.]
Peter Alspach wrote:
Dear List members:
I'm using R1.7.1 (Windows 2000) and having difficulty with scoping.
I've studied the FAQ and on-line manuals and think I have identified
the
source of my difficulty, but cannot work out the solution.
For the purposes of illustration. I have three functions as defined
below:
fnA - function(my.x)
{
list(first=as.character(substitute(my.x)), second=sqrt(my.x))
}
fnB - function(AA)
{
tmp.x - get(AA$first)
tmp.x^2
}
fnC - function()
{
x - 1:2
y - fnA(x)
z - fnB(y)
c(x,y,z)
}
fnA() has a vector as an argument and returns the name of the vector
and the square root of its elements in a list. fn(B) takes the result
of fn(A) as its argument, gets the appropriate vector and computes the
square of its elements. These work fine when called at the command
line.
fnC() defines a local vector x and calls fnA() which operates on this
vector. Then fnB() is called, but it operates on a global vector x in
GlobalEnv (or returns an error is x doesn't exist there) - but I want
it
to operate on the local vector.
I am not sure what you really want to do here, but R works best if you
pass functions an object to work on, and not the name of an object.
Normally this sort of thing is best avoided, but if it is really needed it
is normally simpler to find the object in the calling function (your fnC).
I think this is related to the enclosing environment of all three
functions being GlobalEnv (since they were created at the command
line),
but the parent environment of fnB() being different when invoked from
within fnC().
My questions:
1 Have I correctly understood the issue ?
2 How do I make fnB() operate on the local vector rather than the
global one ?
3 And, as an aside, I have used as.character(substitute(my.x)) to
pass
the name - but deparse(substitute(my.x)) also works. Is there any
reason to prefer one over the other?
deparse() is preferred. One subtle reason is what happens with very long
expressions (and note deparse may give more than one line of output),
but the main reason is what happens with expressions such as calls:
fn1 - function(x) as.character(substitute(x))
fn2 - function(x) deparse(substitute(x))
fn1(log(x))
[1] log x
fn2(log(x))
[1] log(x)
It is normally dangerous to use get() on the result of
deparse(substitute()), as the latter may be the character representation
of an expression and not a simple name.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] nlme lme for complex mixed ANOVAs
Dear List, I downloaded R for the first time yesterday, in the hopes that I might deal more effectively with a complex repeated measures experimental design involving inbred strains of laboratory mice. The design below, somewhat simplified, cannot be computed with standard ANOVA, because something called the X'X matrix is too large. The design has the following factors: Between-subject factors (levels): inbred mouse strain (20, twenty) sex (2) Animals: 10 per sex*strain combination (400 total) Within-subject factors: drug (3) trial set (4) stimulus characteristic 1 (2) stimulus characteristic 2 (2) My question for the R community is, does R have in nlme or lme the ability to compute this problem on a typical desktop PC? In Stata, for instance, which has a matrix size limit of 11,000, the problem above will not fit, using standard ANOVA. Matrix size can be determined by listing all the terms of the full model and multiplying the levels of each factor withing a term and summing all terms. I said simplified in the first paragraph above, because if I include the day of drug challenge, the model oversteps the number of within-subject factors allowed. So, surprisingly to me, Stata/SE, which I bought for large data sets, is too small! Not that I don't like Stata, but I am annoyed that I must find another tool to use. I understand that SAS Proc Mixed will compute the problem, because it may handle the covariance matrix in some kind of piecemeal fashion (perhaps by animal but I've no confirmation of this, except that it zips through a comparable data set on someone else's computer). However, I am running Apple OS X and don't have SAS on my machine. I don't really understand what's going on underneath these programs respective hoods, but my question here is whether R computes mixed models in such a way as to require a matrix size like Stata or like SAS, when mixed models of this size are presented? Sincerely, -Dave __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] residual lag scatter plot
Dear all,
I am looking for a function that can scatter plot the residuals obtained
from a longitudinal model i.e. plot e_{i,j} e_{i,k} for
all j k = 1,..n ( I have 7 observations for each subject). something
similar to the pairs() function.
I can do it the long way by constructing a residual matrix and use pairs.
However, I am wondering if there is a builty in function to do this.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] run R under unix
Dear I used to use R under windows. I always save my code in a txt file. When I want to run something, I just copy and paste to the windows. I recently need to do a big simulation work which will cost 5 or 6 days. I am afraid the memory of windows can not cope with this situation. I now want to run the code under unix. However, I do not know how to run this code in txt file under unix. I just log in to the unix and enter the R, what should I do next? One more question is: if I log off my computer when R is running under unix (i.e., disconnect my computer from server), will I get the result when I log in my computer next time? Thanks! _ Get 10mb of inbox space with MSN Hotmail Extra Storage __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] run R under unix
Zhen Pang wrote: ... I now want to run the code under unix. However, I do not know how to run this code in txt file under unix. I just log in to the unix and enter the R, what should I do next? One more question is: if I log off my computer when R is running under unix (i.e., disconnect my computer from server), will I get the result when I log in my computer next time? ... You'll lose it if you run R in the normal, interactive way. Running it in the background will allow you to log out and still have it running, but! 1) If you're not the only person using this machine, you learn the command nice before you begin. 2) I'm not certain you'll be able to produce jpeg or png graphics when backgrounded; your backgrounded task needs access to the windowing system for graphics rendering, and local security policy might prohibit this. 3) Save early, save often. You probably already know that, but it bears repeating. Here are some suggested steps to run your simulation in background mode. Unfortunately, the exact commands will depend on which version of Unix you're using, and what command shells are available. Consult your local expert. 1) transfer the text file of commands to the unix machine. FTP, using ASCII mode is safest. 2) log onto the Unix machine. 3) run sh, ksh, or bash. (the syntax for what follows is different for the C shell, and I don't know it). 4) using a stripped-down version of your script which will complete in a short time (say, a minute or two), just to check that things work, type nohup nice 15 R my.small.script my.output 21 (again, learn what nice means before you use it. This may not be suitable, and it's impossible for me to tell from here if it is). I know that's not the full answer, but only someone who knows the local setup can give you that answer. Cheers Jason -- Indigo Industrial Controls Ltd. http://www.indigoindustrial.co.nz 64-21-343-545 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
