Re: [R] Logistic regression with more than two choices
>Dear all R-users, > >I am a new user of R and I am trying to build a discrete choice model (with >more than two alternatives A, B, C and D) using logistic regression. I have >data that describes the observed choice probabilities and some background >information. An example below describes the data: > >SexAge pr(A) pr(B) pr(C) pr(D) ... >1 11 0.5 0.5 0 0 >1 40 1 0 0 0 >0 34 0 0 0 1 >0 64 0.1 0.5 0.2 0.2 >... You can use multinom() Here is an exemple For example let this matrix to be analyzed: male female aborted factor 10 12 1 1.2 14 14 4 1.3 15 12 3 1.4 The data are to be entered in a text file like this: output factor n m 1.2 10 f 1.2 12 a 1.2 1 m 1.3 14 f 1.3 14 a 1.3 4 m 1.4 15 f 1.4 12 a 1.4 3 library(MASS) dt.plr <- multinom(output ~ factor, data=dt, weights=n, maxit=1000) dt.pr1<-predict(dt.plr, , type="probs") dt.pr1 >I have been able to model a case with only two alternatives "A" and "not A" >by using glm(). > >I do not know what functions are available to estimate such a model with >more than two alternatives. Multinom() is one possibility, but it only >allows the use of binary 0/1-data instead of observed probabilities. Did I >understand this correctly? > >Additionally, I am willing to use different independent variables for the >different alternatives in the model. Formally, I mean that: >Pr(A)=exp(uA)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) >Pr(B)=exp(uB)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) >... >where uA, uB, uC and uD are linear functions with different independent >variables, e.g. uA=alpha_A1*Age, uB=alpha_B1*Sex. > >Do you know how to estimate this type of models in R? I don't think it is possible... (at least simply, without writing all the script !) Note that I don't undrestand where the residual deviance from multinom() come from. I cant find the logic. Marc -- __ Marc Girondot, Pr Laboratoire Ecologie, Systématique et Evolution Equipe de Conservation des Populations et des Communautés CNRS, ENGREF et Université Paris-Sud 11 , UMR 8079 Bâtiment 362 91405 Orsay Cedex, France Tel: 33 1 (0)1.69.15.72.30 Fax: 33 1 (0)1 69 15 56 96 e-mail: [EMAIL PROTECTED] Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html Skype: girondot Fax in US: 1-425-732-6934 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Plotting rows (or columns) from a matrix in different graphs, not using "par"
On 6/14/05, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> On 6/14/05, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> > On 6/14/05, Camarda, Carlo Giovanni <[EMAIL PROTECTED]> wrote:
> > > Dear R-users,
> > > I would like to ask whether it's possible (for sure it would be), to
> > > plot each rows (or columns) in different graphs and in the same figure
> > > region without using the function "par" and then struggling around with
> > > "axes" and labels etc.
> > > Luckily, I would always have "rows + columns = even number" and the same
> > > "ylim".
> > >
> > > The next one could be a sort of example on what I would like to avoid
> > > plotting the rows of the matrix "mat":
> > >
> > > ### EXAMPLE ##
> > > dat <- sort(runif(16, 1, 1000))
> > > mat <- matrix(dat, ncol=4, byrow = T)
> > > y <- seq(1950, 1953)
> > > par(mfrow=c(2,2))
> > > plot(y, mat[1,], ylim=c(1,1000), axes=F, xlab="", ylab="simul.")
> > > box()
> > > axis(side=2, tick=T, labels=T)
> > > axis(side=3, tick=T, labels=T)
> > > plot(y, mat[2,], ylim=c(1,1000), axes=F, xlab="", ylab="")
> > > box()
> > > axis(side=3, tick=T, labels=T)
> > > plot(y, mat[3,], ylim=c(1,1000), axes=F, xlab="years", ylab="simul.")
> > > box()
> > > axis(side=1, tick=T, labels=T)
> > > axis(side=2, tick=T, labels=T)
> > > plot(y, mat[4,], ylim=c(1,1000), axes=F, xlab="years", ylab="")
> > > box()
> > > axis(side=1, tick=T, labels=T)
> > > ### END EXAMPLE
> > >
> > > Naturally something more compact would be even nicer.
> > >
> >
> > plot(ts(mat, start = 1950), nc = 2)
> >
>
> I just looked at this again and noticed that it did not have
> identical y axes. I tried with ylim= too but plot.ts seemed
> to ignore it.
>
> Here is another solution that uses the zoo library.
> Be sure to use zoo 1.0-1. It uses the screens= argument of
> plot.zoo which was not available in earlier versions of zoo.
>
> Instead of creating a ts series with 4 columns create a
> similar zoo series. Append to that 4 dummy series all
> with the same range. Use the screens= argument
> on plot.zoo to plot one real series and one dummy series
> in each graph. Use type = "l" to plot the real series
> and type = "n" for the dummy series. This causes the
> real series to be plotted with lines and the dummy
> series to be plotted invisibly yet still affect the
> y axis range.Note that we have used the fact that
> mat is square so if the real mat is not square be sure
> to modify this accordingly.
>
> library(zoo) #
>
> n <- nrow(mat)
> z <- zooreg(cbind(mat, max(mat) * diag(n)), start = 2000)
> plot(z, screens = c(1:4, 1:4), type = rep(c("l","n"), each = 4), nc = 2)
>
The problem with ylim= that seems to affect both plot.ts and plot.zoo
has been fixed in zoo so in the next version of zoo the above will
be reducable to just this:
plot(zooreg(mat, start = 2000), ylim = range(mat), nc = 2)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Puzzled in utilising summary.lm() to obtain Var(x)
On 6/14/05, Ajay Narottam Shah <[EMAIL PROTECTED]> wrote: > > > I have a program which is doing a few thousand runs of lm(). Suppose > > > it is a simple model > > > y = a + bx1 + cx2 + e > > > > > > I have the R object "d" where > > > d <- summary(lm(y ~ x1 + x2)) > > > > > > I would like to obtain Var(x2) out of "d". How might I do it? > > > > > > I can, of course, always do sd(x2). But it would be much more > > > convenient if I could snoop around the contents of summary.lm and > > > extract Var() out of it. I couldn't readily see how. Would you know > > > what would click? > > > > Is the question how to get the variance of a column of the > > model matrix for a model that is the sum of terms given only > > summary output and the column name but not the name of the > > data frame? If that is it then try this: > > > > d <- summary(lm(Sepal.Length ~ Sepal.Width, iris)) # test data > > var(model.matrix(eval(d$call))[,"Sepal.Width"]) > > Yes, this is indeed exactly what I was looking for :-) Thanks, > > The eval() pays the full cost of running d$call? Yes. It reruns it. If we can assume that the second arg to lm is data= then we could do this which simply grabs the indicated column from the data frame: f <- function(d, name) eval(substitute(with(eval(d$call[[3]]), name))) f(d, Sepal.Width) # same as iris$Sepal.Width # or f <- function(d, charname) eval(d$call[[3]])[[charname]] f(d, "Sepal.Width") __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Puzzled in utilising summary.lm() to obtain Var(x)
> > I have a program which is doing a few thousand runs of lm(). Suppose > > it is a simple model > > y = a + bx1 + cx2 + e > > > > I have the R object "d" where > > d <- summary(lm(y ~ x1 + x2)) > > > > I would like to obtain Var(x2) out of "d". How might I do it? > > > > I can, of course, always do sd(x2). But it would be much more > > convenient if I could snoop around the contents of summary.lm and > > extract Var() out of it. I couldn't readily see how. Would you know > > what would click? > > Is the question how to get the variance of a column of the > model matrix for a model that is the sum of terms given only > summary output and the column name but not the name of the > data frame? If that is it then try this: > > d <- summary(lm(Sepal.Length ~ Sepal.Width, iris)) # test data > var(model.matrix(eval(d$call))[,"Sepal.Width"]) Yes, this is indeed exactly what I was looking for :-) Thanks, The eval() pays the full cost of running d$call? -ans. -- Ajay Shah Consultant [EMAIL PROTECTED] Department of Economic Affairs http://www.mayin.org/ajayshah Ministry of Finance, New Delhi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] RGui crashes on wle call
Prof Brian Ripley wrote:
> On Mon, 13 Jun 2005, Chris Bergstresser wrote:
>>I'm seeing the following commands reliably produce a crash in RGui,
>> version 2.0.1, for both my home and office machine:
>>
>> > rm(list = ls(all = TRUE));
>> > load("dataset.R");
>> > library("wle");
>> > data.wle = wle.lm(abortion ~ year * lib.con + age + gender +
>> + urbanism + census + income + church.att + children + educ +
>> + religion.imp, data = data.set);
>>
> 1) Re-do the tests in the current version of R, preferably a beta of 2.1.1.
Yeah -- I upgraded to R 2.1.0, and it still reliably crashes.
> 2) Read the rw-FAQ, do the debugging reported there (with Dr MinGW or
> gdb) and find where it is crashing. (This is very likely to be in wle.)
> If it is in wle, send a report to the maintainer. If it is in R, send a
> report to R-bugs.
I'm a little loath to download and install a debugger, as I've never
done it before. I don't even know what to look for if I were to install it.
> In either case, supply enough data to reproduce the
> problem.
I can easily provide the datafile which seems to be causing it.
It's only 200k, so if anyone is interested in pursuing the matter I'd be
happy to send it to them. This is on Windows XP, btw.
-- Chris
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
Prof Brian Ripley wrote: > On Tue, 14 Jun 2005, Frank E Harrell Jr wrote: > >> Prof Brian Ripley wrote: >> >>> On Tue, 14 Jun 2005, James Salsman wrote: >>> >>> How can ordinary polynomial coefficients be calculated from an orthogonal polynomial fit? >>> >>> >>> >>> Why would you want to do that? predict() is perfectly happy with an >>> orthogonal polynomial fit and the `ordinary polynomial coefficients' >>> are rather badly determined in your example since the design matrix >>> has a very high condition number. >> >> >> Brian - I don't fully see the relevance of the high condition number >> nowadays unless the predictor has a really bad origin. Orthogonal >> polynomials are a mess for most people to deal with. > > > It means that if you write down the coeffs to a few places and then try > to reproduce the predictions you will do badly. The perturbation > analysis depends on the condition number, and so is saying that the > predictions are dependent on fine details of the coefficients. Right - I carry several digits of precision when I do this. > > Using (year-2000)/1000 or (year - 1970)/1000 would be a much better idea. > > Why do `people' need `to deal with' these, anyway. We have machines to > do that. The main application I think of is when we publish fitted models, but it wouldn't be that bad to restate fitted orthogonal polynomials in simpler notation. -Frank > >> >> Frank >> >>> >>> I'm trying to do something like find a,b,c,d from lm(billions ~ a+b*decade+c*decade^2+d*decade^3) but that gives: "Error in eval(expr, envir, enclos) : Object "a" not found" >>> >>> >>> >>> You could use >>> >>> lm(billions ~ decade + I(decade^2) + I(decade^3)) >>> >>> except that will be numerically inaccurate, since >>> >>> m <- model.matrix(~ decade + I(decade^2) + I(decade^3)) kappa(m) >>> >>> >>> [1] 3.506454e+16 >>> >>> >>> >>> > decade <- c(1950, 1960, 1970, 1980, 1990) > billions <- c(3.5, 5, 7.5, 13, 40) > # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg > > pm <- lm(billions ~ poly(decade, 3)) > > plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), main="average yearly inflation-adjusted dollar cost of extreme weather events worldwide") > curve(predict(pm, data.frame(decade=x)), add=TRUE) > # output: http://www.bovik.org/storms.gif > > summary(pm) Call: lm(formula = billions ~ poly(decade, 3)) Residuals: 1 2 3 4 5 0.2357 -0.9429 1.4143 -0.9429 0.2357 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept)13.800 0.882 15.647 0.0406 * poly(decade, 3)1 25.614 1.972 12.988 0.0489 * poly(decade, 3)2 14.432 1.972 7.318 0.0865 . poly(decade, 3)36.483 1.972 3.287 0.1880 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 1.972 on 1 degrees of freedom Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 > pm Call: lm(formula = billions ~ poly(decade, 3)) Coefficients: (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 13.80025.61414.432 6.483 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >>> >>> >> >> >> -- >> Frank E Harrell Jr Professor and Chair School of Medicine >> Department of Biostatistics Vanderbilt University >> >> > -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Calling C from Fortran
Dear Gilles,
You can try this.
Save the fortran file as 1.f
Save the C file as 2.c
Then
R CMD SHLIB 1.f 2.c
You should obtain a file named 1.so.
Then start you R program, then
dyn.load('/where/is/your/1.so')
.Fortran("testit")
You should obtain a random number which is normally distributed.
Good luck.
Jin
On 6/15/05, Gilles GUILLOT <[EMAIL PROTECTED]> wrote:
> I would like to call C routines from Fortran as suggested in section 5.6 of
> the "Writing R extensions" documentation.
>
> I'm familiar with Fortran but not with C.
> I understand the example provided in Fortran:
>
> subroutine testit()
> double precision normrnd, x
> call rndstart()
> x = normrnd()
> call dblepr("X was", 5, x, 1)
> call rndend()
> end
>
>
> but I don't understand the purpose of this C wrapper:
> #include
> void F77_SUB(rndstart)(void) { GetRNGstate(); }
> void F77_SUB(rndend)(void) { PutRNGstate(); }
> double F77_SUB(normrnd)(void) { return norm_rand(); }
>
> neither how I should compile it.
>
> Could anyone explain how I should compile and link
> the C and Fortran files above, and call the Fortran subroutine from R.
>
> Thanks,
>
> Gilles
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Plotting rows (or columns) from a matrix in different graphs, not using "par"
On 6/14/05, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> On 6/14/05, Camarda, Carlo Giovanni <[EMAIL PROTECTED]> wrote:
> > Dear R-users,
> > I would like to ask whether it's possible (for sure it would be), to
> > plot each rows (or columns) in different graphs and in the same figure
> > region without using the function "par" and then struggling around with
> > "axes" and labels etc.
> > Luckily, I would always have "rows + columns = even number" and the same
> > "ylim".
> >
> > The next one could be a sort of example on what I would like to avoid
> > plotting the rows of the matrix "mat":
> >
> > ### EXAMPLE ##
> > dat <- sort(runif(16, 1, 1000))
> > mat <- matrix(dat, ncol=4, byrow = T)
> > y <- seq(1950, 1953)
> > par(mfrow=c(2,2))
> > plot(y, mat[1,], ylim=c(1,1000), axes=F, xlab="", ylab="simul.")
> > box()
> > axis(side=2, tick=T, labels=T)
> > axis(side=3, tick=T, labels=T)
> > plot(y, mat[2,], ylim=c(1,1000), axes=F, xlab="", ylab="")
> > box()
> > axis(side=3, tick=T, labels=T)
> > plot(y, mat[3,], ylim=c(1,1000), axes=F, xlab="years", ylab="simul.")
> > box()
> > axis(side=1, tick=T, labels=T)
> > axis(side=2, tick=T, labels=T)
> > plot(y, mat[4,], ylim=c(1,1000), axes=F, xlab="years", ylab="")
> > box()
> > axis(side=1, tick=T, labels=T)
> > ### END EXAMPLE
> >
> > Naturally something more compact would be even nicer.
> >
>
> plot(ts(mat, start = 1950), nc = 2)
>
I just looked at this again and noticed that it did not have
identical y axes. I tried with ylim= too but plot.ts seemed
to ignore it.
Here is another solution that uses the zoo library.
Be sure to use zoo 1.0-1. It uses the screens= argument of
plot.zoo which was not available in earlier versions of zoo.
Instead of creating a ts series with 4 columns create a
similar zoo series. Append to that 4 dummy series all
with the same range. Use the screens= argument
on plot.zoo to plot one real series and one dummy series
in each graph. Use type = "l" to plot the real series
and type = "n" for the dummy series. This causes the
real series to be plotted with lines and the dummy
series to be plotted invisibly yet still affect the
y axis range.Note that we have used the fact that
mat is square so if the real mat is not square be sure
to modify this accordingly.
library(zoo) #
n <- nrow(mat)
z <- zooreg(cbind(mat, max(mat) * diag(n)), start = 2000)
plot(z, screens = c(1:4, 1:4), type = rep(c("l","n"), each = 4), nc = 2)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Logistic regression with more than two choices
Hi Koskinen For response variables with multiple categories, you may try polr() in MASS package, which implement a proportional odds model. And you may search the R archives, several threads discussed this problem before... Wuming On 6/15/05, Ville Koskinen <[EMAIL PROTECTED]> wrote: > Dear all R-users, > > I am a new user of R and I am trying to build a discrete choice model (with > more than two alternatives A, B, C and D) using logistic regression. I have > data that describes the observed choice probabilities and some background > information. An example below describes the data: > > Sex Age pr(A) pr(B) pr(C) pr(D) ... > 1 11 0.5 0.5 0 0 > 1 40 1 0 0 0 > 0 34 0 0 0 1 > 0 64 0.1 0.5 0.2 0.2 > ... > > I have been able to model a case with only two alternatives "A" and "not A" > by using glm(). > > I do not know what functions are available to estimate such a model with > more than two alternatives. Multinom() is one possibility, but it only > allows the use of binary 0/1-data instead of observed probabilities. Did I > understand this correctly? > > Additionally, I am willing to use different independent variables for the > different alternatives in the model. Formally, I mean that: > Pr(A)=exp(uA)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) > Pr(B)=exp(uB)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) > ... > where uA, uB, uC and uD are linear functions with different independent > variables, e.g. uA=alpha_A1*Age, uB=alpha_B1*Sex. > > Do you know how to estimate this type of models in R? > > Best regards, Ville Koskinen > > __ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Reducing the FPR (false positive rate)
Hello R-USERS, I think some people didn't understand my question. What I want is to use the training set to minimize the FALSE POSITIVE RATE. I think it is possible. I sacrifice ACCURACY to have less FALSE POSITIVES. I don't want a classifier result with 5% of FPR and, for example, 93% of ACCURACY. I want a 1% FPR sacrificing the 93% ACCURACY. Do you know how can I do this? I really need this because the requisite of the work I'm doing is only 1% of FPR. Thanks and best regards. * |Üñð뮧µñ| - Anderson de Rezende Rocha Bacharel em Ciência da Computação - UFLA Mestrando em Ciência da Computação - UNICAMP Esteganografia e esteganálise digital UNIVERSIDADE ESTADUAL DE CAMPINAS, SP - BRASIL < http://andersonrocha.cjb.net > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] anova.lme error
Hi,
I am working with R version 2.1.0, and I seem to have run into what looks
like a bug. I get the same error message when I run R on Windows as well as
when I run it on Linux.
When I call anova to do a LR test from inside a function, I get an error.
The same call works outside of a function. It appears to not find the right
environment when called from inside a function. I have provided the code
below.
Thanks,
Nisha Mulakken
myFunction <- function(myDataFrame) {
# Less restricted
fit1 <- gls(y ~ dose,
weights=varIdent(form=~1|dose),
data=myDataFrame)
# more restricted
fit2 <- gls(y ~ dose,
data=myDataFrame)
anova.results <- anova(fit1, fit2)
anova.results
}
df <- data.frame( y=c(12,3,45,1,53,6),
dose=c(0,10,200,0,10,200),
time=c("4.00 hrs", "4.00 hrs", "6.00 hrs", "6.00 hrs",
"8.00 hrs", "8.00 hrs"),
time.hours=c(4, 4, 6, 6, 8, 8),
rep=rep("a", 6)
)
## This leads to the following error:
## Error in anova.lme(object = fit1, fit2) : Object "fit2" not found
results <- myFunction(myDataFrame=df)
#
## The same thing outside of a function
# Less restricted
fit3 <- gls(y ~ dose,
weights=varIdent(form=~1|dose),
data=df)
# more restricted
fit4 <- gls(y ~ dose,
data=df)
## This works:
anova(fit3, fit4)
## The results:
## > anova(fit3, fit4)
## Model df AIC BIClogLik Test L.Ratio p-value
## fit3 1 5 57.98998 54.92145 -23.99499
## fit4 2 3 55.75284 53.91172 -24.87642 1 vs 2 1.76286 0.4142
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] within and between subject calculation
On 6/14/05, Kerry Bush <[EMAIL PROTECTED]> wrote:
> Dear helpers in this forum,
>
> I have the following question:
>
> Suppose I have the following data set:
>
> id x y
> 023 1 2
> 023 2 5
> 023 4 6
> 023 5 7
> 412 2 5
> 412 3 4
> 412 4 6
> 412 7 9
> 220 5 7
> 220 4 8
> 220 9 8
> ..
>
> and i want to calculate sum_{i=1}^k
> sum_{j=1}^{n_i}x_{ij}*y_{ij}
>
> is there a simple way to do this within and between
> subject summation in R?
You didn't make it clear what the indices i and j refer to. It seems
that part of the answer could use
> samp
id x y
1 23 1 2
2 23 2 5
3 23 4 6
4 23 5 7
5 412 2 5
6 412 3 4
7 412 4 6
8 412 7 9
9 220 5 7
10 220 4 8
11 220 9 8
> with(samp, tapply(x*y, id, sum))
23 220 412
71 139 109
but if you then sum this result you simply get the sum of x * y.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lattice, panel.grid, and scales=list(tick.number=XXX)
If you look at the code of panel.grid, you'll see why it doesn't work -- it
does not use any of the scale parameters. Moreover the Help page for
panel.grid explicitly warns that the h,v=-1 specification may not work, so
no promises have been broken.
I'm not sure how bwplot handles grids, since one of the axes is determined
by the number of groups and not the range of the data as in xyplot. You
might try specifying the at= argument for both the scales list and at the
top level call for panel.grid to pick up. i.e. bwplot( ..., scales = list
(at = where tic marks go),at= where tic marks go, ...).
Alternatively, wait for Deepayan to give a definitive answer.
-- Bert Gunter
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of M. K.
Sent: Tuesday, June 14, 2005 3:37 PM
To: R-help mailing list
Subject: [R] lattice, panel.grid, and scales=list(tick.number=XXX)
I have a Lattice plot in which I want to adjust the number of tick
marks used, and I want to have the drawn grid reflect that change.
Here is what I'm doing:
bwplot(var1 ~ var2, data=df, scales=list(tick.number=10),
panel=function(...) {
panel.grid(h=0,v=-1,...);
panel.stripplot(col="gray40", pch="|", cex=2, ...);
panel.bwplot(...);
})
Unfortunately this doesn't quite work. Although the bwplot's tick
marks are indeed increased as requested, the panel.grid produces the
same (3 line) grid as before, seemingly unaware of the changed # of
ticks.
Any suggestions on how to achieve what I want?
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Preparing timestamped data for fourier analysis
"Milos Zarkovic" <[EMAIL PROTECTED]> wrote in message news:[EMAIL PROTECTED] > I believe that FFT is not appropriate. However Lomb-Scargle periodogram > could be used. This may interest you: (Preprint of submitted paper) Detecting periodic patterns in unevenly spaced gene expression time series using Lomb-Scargle periodograms. http://research.stowers-institute.org/bioinfo/PDF/m2005_lomb-scargle_submitted.pdf R code here http://research.stowers-institute.org/efg/2005/LombScargle/ efg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Matrix stability problem
It looks to me like the x measurements are normalized to sum to 1. So
whatever measurement error there is gets "spread around", so to speak.
It would help if you could explain the setting a little more fully? Why is
it, for example, that you know the A values from experiment to experiment
(they do seem to vary) but there's no measurment error? Why the variation?
Do you know b, or is it an estimate from some measurements x and readings A?
I guess this is probably a (multivariate-response) regression problem, and
the question is where the error is and what its structure is. Imposing the
constraint that x sums to 1 would probably help. This makes an
overdetermined problem (two free parameters, three experiments) so you are
forced into regression.
Would you ever have more than three experiments? If so would that change the
formulation of the problem? More experiments + regression might be the
simplest way to get a more accurate solution.
Reid Huntsinger
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Lapointe, Pierre
Sent: Tuesday, June 14, 2005 5:45 PM
To: 'r-help@stat.math.ethz.ch'
Subject: [R] Matrix stability problem
Hello,
This is not a problem with R, the calculated results are mathematically
correct. This a matrix stability problem. Because of measuring errors, my
matrix solution is a bit off.
Here is what my equations look like:
A11 x11+A12 x12 +A13 x13 = b1
A21 x21+A22 x21 +A23 x23 = b2
A31 x31+A32 x31 +A33 x33 = b3
A is a reading, X is a measured weight, and b is total. The 3 experiments
give slightly different X values because of measurement errors.
For reproducibility, here's my A, x and b matrices and vectors
A <-matrix(
c(0.03,0.02,0.04,0.01,0.015,0.03,-0.01,-0.02,0.03),3,3,byrow=TRUE)
x <-matrix( c(0.2,0.3,0.5,0.205,0.305,0.49,0.19,0.29,0.52),3,3,byrow=TRUE)
b <-matrix( c(0.032,0.021325,0.0079),3,1)
As expected, rowSums(A*x) = b
Problem: Let's now assume I don't know x. I'd like to solve for x in Ax=b. I
am aware that my x is a matrix and solve(A,b) will give me a vector.
However, looking at the x matrix, one can easily see that the real x[,1]
(without measurement error) is close to 0.2, x[,2] is close to 0.3 and x[,3]
is close to 0.5
> x
[,1] [,2] [,3]
[1,] 0.200 0.300 0.50
[2,] 0.205 0.305 0.49
[3,] 0.190 0.290 0.52
However, solve(A,b) gives me a vector that is not close to the expected
solution: > solve(A,b)
[,1]
[1,] 0.214
[2,] 0.2612857 # Far from 0.2
[3,] 0.5088571
Do you know any function/package in R that could help me get a result closer
to the expected one?
Regards,
Pierre Lapointe
***
AVIS DE NON-RESPONSABILITE:\ Ce document transmis par courri...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lattice, panel.grid, and scales=list(tick.number=XXX)
I have a Lattice plot in which I want to adjust the number of tick
marks used, and I want to have the drawn grid reflect that change.
Here is what I'm doing:
bwplot(var1 ~ var2, data=df, scales=list(tick.number=10),
panel=function(...) {
panel.grid(h=0,v=-1,...);
panel.stripplot(col="gray40", pch="|", cex=2, ...);
panel.bwplot(...);
})
Unfortunately this doesn't quite work. Although the bwplot's tick
marks are indeed increased as requested, the panel.grid produces the
same (3 line) grid as before, seemingly unaware of the changed # of
ticks.
Any suggestions on how to achieve what I want?
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to fix false positve rates?
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Anderson de Rezende
Rocha
Sent: Tuesday, June 14, 2005 11:01 AM
To: r-help@stat.math.ethz.ch
Subject: [R] How to fix false positve rates?
Dear R-users,
I have a set of 12000 image samples. I can divide this set into two
categories: training and testing. I need to classify the test set into
a two qualitative outputs: true or false for some characteristic.
To do the classification I'm using the packages SVM in e1071 library
and LDA in the MASS library. However, I'm with a great number of FALSE
POSITIVE CASES in both classifiers, about 10/15%.
1) How can I use these classifiers with a fixed 1% FALSE POSITIVE RATE?
Surely you jest! Classification is not hypothesis testing.You do not get to
control either false positive or negative rates. You try to minimize (an
unbiased estimate of) mislclassification rate.
2) Could anynone indicate me a non-linear SVM classifier R-package?
3) Could anynone indicate me an ADA-BOOST classifier R-package?
Look!
RSiteSearch('boosting', restrict = 'functions')
Cheers,
Bert Gunter
Thanks for all. I'm anxiously wait the answers. Best regards.
*
|Üñð뮧µñ| - Anderson de Rezende Rocha
Bacharel em Ciência da Computação - UFLA
Mestrando em Ciência da Computação - UNICAMP
Esteganografia e esteganálise digital
UNIVERSIDADE ESTADUAL DE CAMPINAS, SP - BRASIL
< http://andersonrocha.cjb.net >
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] update.packages() - gregmisc
Is the code in your post intended to show what worked so others
will know what to do or is that code intended to show what you
did but did not work?
If its the latter, I successfully did it last week and don't
clearly remember my precise steps but I may have done this:
R CMD remove gdata
R CMD remove gmodels
R CMD remove gplots
R CMD remove gtools
R CMD remove gregmisc
I assume that using remove.packages, viz.
remove.packages("gregmisc")
etc.
would have given the same result.
On 6/14/05, Muhammad Subianto <[EMAIL PROTECTED]> wrote:
> Thanks.
> I do like this,
> > remove.packages("gregmisc", .libPaths()[1])
> > remove.packages("gtools", .libPaths()[1])
> > install.packages("gregmisc", .libPaths()[1])
> > update.packages()
> > update.packages()
> > install.packages("gtools", .libPaths()[1])
> > update.packages()
> > update.packages()
> > update.packages(ask='graphics')
>
> Regards,
> Muhammad Subianto
> R.2.1.0 on W2K
>
> On this day 6/14/2005 1:51 PM, Gabor Grothendieck wrote:
> > On 6/14/05, Muhammad Subianto <[EMAIL PROTECTED]> wrote:
> >
> >>Dear all,
> >>I have a problem to update package gregmisc.
> >>After I update,
> >> > update.packages(ask='graphics')
> >>trying URL
> >>'http://cran.at.r-project.org/bin/windows/contrib/2.1/gregmisc_2.0.8.zip'
> >>Content type 'application/zip' length 2465 bytes
> >>opened URL
> >>downloaded 2465 bytes
> >>
> >>package 'gregmisc' successfully unpacked and MD5 sums checked
> >>...
> >>
> >>then try to update again, still I must update package gregmisc, etc.
> >>I have tried 3,4,5, times with the same result.
> >>
> >
> >
> > This was discussed on r-devel recently. See:
> >
> > https://www.stat.math.ethz.ch/pipermail/r-devel/2005-June/033479.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
On 6/14/05, Richard Mott <[EMAIL PROTECTED]> wrote: > Does anyone know a way to do the following: > > Save a large number of R objects to a file (like load() does) but then > read back only a small named subset of them . As far as I can see, > load() reads back everything. > > The context is: > > I have an application which will generate a large number of large > matrices (approx 15000 matrices each of dimension 2000*30). I can > generate these matrices using an R-package I wrote, but it requires a > large amouint of memory and is slow so I want to do this only once. > However, I then want to do some subsequent processing, comprising a very > large number of runs in which small (~ 10) random selection of matrices > from the previously computed set are used for linear modeling. So I > need a way to load back named objects previously saved in a call to > save(). I can;t see anyway of doing this. Any ideas? Check out the g.data delayed data package on CRAN and the article in R News 2/3. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Puzzled in utilising summary.lm() to obtain Var(x)
On 6/14/05, Ajay Narottam Shah <[EMAIL PROTECTED]> wrote: > I have a program which is doing a few thousand runs of lm(). Suppose > it is a simple model > y = a + bx1 + cx2 + e > > I have the R object "d" where > d <- summary(lm(y ~ x1 + x2)) > > I would like to obtain Var(x2) out of "d". How might I do it? > > I can, of course, always do sd(x2). But it would be much more > convenient if I could snoop around the contents of summary.lm and > extract Var() out of it. I couldn't readily see how. Would you know > what would click? > Is the question how to get the variance of a column of the model matrix for a model that is the sum of terms given only summary output and the column name but not the name of the data frame? If that is it then try this: d <- summary(lm(Sepal.Length ~ Sepal.Width, iris)) # test data var(model.matrix(eval(d$call))[,"Sepal.Width"]) If that's not the question, try examining the contents of d using str(d) and examine the output of lm via str(eval(d$call)) and perhaps that will suggest the answer. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Plotting rows (or columns) from a matrix in different graphs, not using "par"
On 6/14/05, Camarda, Carlo Giovanni <[EMAIL PROTECTED]> wrote: > Dear R-users, > I would like to ask whether it's possible (for sure it would be), to > plot each rows (or columns) in different graphs and in the same figure > region without using the function "par" and then struggling around with > "axes" and labels etc. > Luckily, I would always have "rows + columns = even number" and the same > "ylim". > > The next one could be a sort of example on what I would like to avoid > plotting the rows of the matrix "mat": > > ### EXAMPLE ## > dat <- sort(runif(16, 1, 1000)) > mat <- matrix(dat, ncol=4, byrow = T) > y <- seq(1950, 1953) > par(mfrow=c(2,2)) > plot(y, mat[1,], ylim=c(1,1000), axes=F, xlab="", ylab="simul.") > box() > axis(side=2, tick=T, labels=T) > axis(side=3, tick=T, labels=T) > plot(y, mat[2,], ylim=c(1,1000), axes=F, xlab="", ylab="") > box() > axis(side=3, tick=T, labels=T) > plot(y, mat[3,], ylim=c(1,1000), axes=F, xlab="years", ylab="simul.") > box() > axis(side=1, tick=T, labels=T) > axis(side=2, tick=T, labels=T) > plot(y, mat[4,], ylim=c(1,1000), axes=F, xlab="years", ylab="") > box() > axis(side=1, tick=T, labels=T) > ### END EXAMPLE > > Naturally something more compact would be even nicer. > plot(ts(mat, start = 1950), nc = 2) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Dateticks
On 6/14/05, Bernard L. Dillard <[EMAIL PROTECTED]> wrote:
> Hello. I am having the worst time converting x-axis date ticks to real
> dates. I have tried several suggestions in online help tips and books to
> no avail.
>
> For example, the x-axis has 0, 50, 100, etc, and I want it to have
> "6/17/03", "8/6/03" etc. See attached (sample).
>
> Can anybody help me with this.
>
> Here's my code:
>
> ts.plot(date.attackmode.table[,1], type="l", col="blue", lty=2,ylab="IED
> Attacks", lwd=2,xlab="Attack Dates",main="Daily Summary of Attack
> Mode")
> grid()
>
The default format for dates in chron is the one needed here
so let us use that. (We could alternately use Date class
in defining tt and specify a format string as a second argument
to format in the last line.)
library(chron)
# test data. values are in y and times are in tt.
y <- rnorm(201)
tt <- chron("8/6/03") + seq(0, length = length(y))
# plot without axis
plot(tt, y, xaxt = "n")
# axis with ticks and labels every 50 days
idx50 <- seq(1, length(tt), 50)
axis(1, tt[idx50], format(tt[idx50]), cex.axis = 0.5)
or, to use Date class replace with this:
tt <- as.Date("2003-08-06") + seq(0, length = length(y))
plot(tt, y, xaxt = "n")
idx50 <- seq(1, length(tt), 50)
axis(1, tt[idx50], format(tt[idx50], "%m/%d/%y), cex.axis = 0.5)
More info on working with dates is in R News 4/1.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Manipulating dates
On 6/14/05, Richard Hillary <[EMAIL PROTECTED]> wrote: > Hello, > Given a vector of characters, or factors, denoting the date in > the following way: 28/03/2000, is there a method of > 1) Computing the earliest of these dates; > 2) Using this as a base, then converting all the other dates into merely > the number of days after this minimum date Convert dates to chron (whose default date format is the one needed here); convert that to numeric (which will be the number of days since some origin) and then do the indicated subtraction: library(chron) x <- as.character(x) # can omit if already character x.num <- as.numeric(chron(x)) x - x.num - min(x.num) Alternately convert it to Date and use julian (or convert x.Date to numeric and use subtraction, as with chron): x <- as.character(x) # can omit in latest R 2.1.0 patched x.Date <- as.Date(x, "%m/%d/%Y") julian(x.Date, min(x.Date)) The as.Date.factor method in the latest R 2.1.0 patched supports the format string (second argument of as.Date) for factors but older versions of R did not. If your x is a factor but you do not want to upgrade right now start off with the statement: x <- as.character(x) More info on dates is in R News 4/1. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Matrix stability problem
Hello,
This is not a problem with R, the calculated results are mathematically
correct. This a matrix stability problem. Because of measuring errors, my
matrix solution is a bit off.
Here is what my equations look like:
A11 x11+A12 x12 +A13 x13 = b1
A21 x21+A22 x21 +A23 x23 = b2
A31 x31+A32 x31 +A33 x33 = b3
A is a reading, X is a measured weight, and b is total. The 3 experiments
give slightly different X values because of measurement errors.
For reproducibility, here's my A, x and b matrices and vectors
A <-matrix(
c(0.03,0.02,0.04,0.01,0.015,0.03,-0.01,-0.02,0.03),3,3,byrow=TRUE)
x <-matrix( c(0.2,0.3,0.5,0.205,0.305,0.49,0.19,0.29,0.52),3,3,byrow=TRUE)
b <-matrix( c(0.032,0.021325,0.0079),3,1)
As expected, rowSums(A*x) = b
Problem: Let's now assume I don't know x. I'd like to solve for x in Ax=b. I
am aware that my x is a matrix and solve(A,b) will give me a vector.
However, looking at the x matrix, one can easily see that the real x[,1]
(without measurement error) is close to 0.2, x[,2] is close to 0.3 and x[,3]
is close to 0.5
> x
[,1] [,2] [,3]
[1,] 0.200 0.300 0.50
[2,] 0.205 0.305 0.49
[3,] 0.190 0.290 0.52
However, solve(A,b) gives me a vector that is not close to the expected
solution: > solve(A,b)
[,1]
[1,] 0.214
[2,] 0.2612857 # Far from 0.2
[3,] 0.5088571
Do you know any function/package in R that could help me get a result closer
to the expected one?
Regards,
Pierre Lapointe
***
AVIS DE NON-RESPONSABILITE:\ Ce document transmis par courri...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] within and between subject calculation
Dear helpers in this forum,
I have the following question:
Suppose I have the following data set:
id x y
023 1 2
023 2 5
023 4 6
023 5 7
412 2 5
412 3 4
412 4 6
412 7 9
220 5 7
220 4 8
220 9 8
..
and i want to calculate sum_{i=1}^k
sum_{j=1}^{n_i}x_{ij}*y_{ij}
is there a simple way to do this within and between
subject summation in R?
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem installing packages with compiled-from-source R.app on Mac OS X - Tiger
Thank you. I am compiling it now...
Costas
On 14 Ιουν 2005, at 10:13 ΜΜ, Prof Brian Ripley wrote:
> You need to use R-patched, nowadays 2.1.1 beta, not the released R
> 2.1.0.
> As far as I know it is fixed there.
>
> On Tue, 14 Jun 2005, Constantinos Antoniou wrote:
>
>
>> Hello all,
>>
>> This may be aimed for r-devel, but I encountered this as an R-user
>> and not an R-developer so I start here (having said that, please
>> direct me to R-devel if you think this is appropriate. I am not
>> cross-
>> posting, as I believe this is bad netiquette).
>>
>> I am a recent, but extremely happy R-user (especially after getting
>> my own copy of MASS 2002). My adventures started when I wanted to use
>> Rpy to use R also from within python. I compiled R (2.1.0) after
>> configuring it as follows:
>>
>> ./configure --enable-R-shlib --with-blas='-framework vecLib' --with-
>> lapack --with-aqua
>>
>> where --enable-R-shlib is required by Rpy.
>>
>> I then compiled Rpy-0.4.2.1 and I was able to call R from within
>> python. So far so good...
>>
>> I then -naively- tried to start R.app. After less than a bounce on
>> the dock... it would not start. [Sorry, I do not remember what
>> problem was showing up in the console...]
>>
>> So, I figured I had to also compile R.app. I got the code via:
>>
>> svn co https://svn.r-project.org/R-packages/trunk/Mac-GUI Mac-
>> GUI (executed this command yesterday)
>>
>> and compiled it as per the instructions... And R.app launches just
>> fine.
>>
>> Now, the issue comes when I try to install a(ny) package (via the GUI
>> package installer of R.app). I then get the following message:
>>
>> "Package installation failed. Package installation was not
>> successful. Please see the R console for details"
>>
>> And the R console says:
>>
>> "Error in install.packages(c("dyn"), lib = "/Library/Frameworks/
>> R.framework/Resources/library", :
>> couldn't find function ".install.macbinary""
>>
>> I see that this has been encountered before and resolved (at least
>> Prof. Ripley saw what was wrong) as per the following thread:
>>
>> https://stat.ethz.ch/pipermail/r-devel/2005-May/033115.html
>>
>> But I am not sure what needs to be done on my part so that I am not
>> affected from it.
>>
>> > str(.Platform)
>> List of 6
>> $ OS.type : chr "unix"
>> $ file.sep : chr "/"
>> $ dynlib.ext: chr ".so"
>> $ GUI : chr "X11"
>> $ endian: chr "big"
>> $ pkgType : chr "source"
>> > R.version.string
>> [1] "R version 2.1.0, 2005-04-18"
>>
>>
>> gcc version 4.0.0
>>
>> [Please let me know what other information may be relevant]
>>
>>
>> Thanking you in advance,
>>
>> Costas
>>
>>
>> --
>> Constantinos Antoniou, Ph.D.
>> Department of Transportation Planning and Engineering
>> National Technical University of Athens
>> 5, Iroon Polytechniou str. GR-15773, Athens, Greece
>>
>> __
>> R-help@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide! http://www.R-project.org/posting-
>> guide.html
>>
>>
>
> --
> Brian D. Ripley, [EMAIL PROTECTED]
> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel: +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax: +44 1865 272595
>
>
--
Constantinos Antoniou, Ph.D.
Department of Transportation Planning and Engineering
National Technical University of Athens
5, Iroon Polytechniou str. GR-15773, Athens, Greece
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Plotting rows (or columns) from a matrix in different graphs, not using "par"
On Tue, 2005-06-14 at 20:00 +0200, Camarda, Carlo Giovanni wrote:
> Dear R-users,
> I would like to ask whether it's possible (for sure it would be), to
> plot each rows (or columns) in different graphs and in the same figure
> region without using the function "par" and then struggling around with
> "axes" and labels etc.
> Luckily, I would always have "rows + columns = even number" and the same
> "ylim".
>
> The next one could be a sort of example on what I would like to avoid
> plotting the rows of the matrix "mat":
>
> ### EXAMPLE ##
> dat <- sort(runif(16, 1, 1000))
> mat <- matrix(dat, ncol=4, byrow = T)
> y <- seq(1950, 1953)
> par(mfrow=c(2,2))
> plot(y, mat[1,], ylim=c(1,1000), axes=F, xlab="", ylab="simul.")
> box()
> axis(side=2, tick=T, labels=T)
> axis(side=3, tick=T, labels=T)
> plot(y, mat[2,], ylim=c(1,1000), axes=F, xlab="", ylab="")
> box()
> axis(side=3, tick=T, labels=T)
> plot(y, mat[3,], ylim=c(1,1000), axes=F, xlab="years", ylab="simul.")
> box()
> axis(side=1, tick=T, labels=T)
> axis(side=2, tick=T, labels=T)
> plot(y, mat[4,], ylim=c(1,1000), axes=F, xlab="years", ylab="")
> box()
> axis(side=1, tick=T, labels=T)
> ### END EXAMPLE
>
> Naturally something more compact would be even nicer.
>
> Thanks in advance,
> Carlo Giovanni Camarda
If I am understanding what you are trying to do, then using xyplot()
from the lattice package (part of the standard R installation) may be
best:
Modify 'mat' so that all data is in 'df':
> df <- data.frame("Simul" = as.vector(mat),
"Years" = rep(1950:1953, 4),
"Group" = rep(LETTERS[1:4], each = 4))
> df
Simul Years Group
1 88.67956 1950 A
2 364.73579 1951 A
3 546.63525 1952 A
4 774.05869 1953 A
5 138.09586 1950 B
6 382.20986 1951 B
7 592.85512 1952 B
8 810.45569 1953 B
9 232.43912 1950 C
10 524.91085 1951 C
11 598.49688 1952 C
12 963.97328 1953 C
13 261.80598 1950 D
14 533.91143 1951 D
15 765.72522 1952 D
16 996.72192 1953 D
# Load the lattice package
> library(lattice)
# Draw the plot, note the use of standard R formulae in the syntax
# with the addition of the grouping vector after the '|'
> xyplot(Simul ~ Years | Group, data = df)
See:
library(lattice)
?xyplot
for more information.
HTH,
Marc Schwartz
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] trivial installation question
Hello, I try to install udunits and RNetCDF packages With the first one after R CMD INSTALL udunits_1.2.tar.gz the script looks for the library, does not find and says to edit udunits_1.0/udunits/src/Makevars.in with links to an include file and location of the library. When I put the correct path the script repeats the previous message still looking for the library in the default path. How do I communicate the location of libraries/include files to R CMD INSTALL? Thanks, Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bug in rpart?
Dear R-helpers,
Can you help me to see why "code 1" gives error
while "code 2" runs fine? The only difference in
the data is the distribution of age categories.
I am attaching the session after the code.
Many thanks.
XL
library(survival)
library(rpart)
# code 1
n <- 20
age <- rep(1:3, c(2, 3, 15))
eg<- data.frame(rexp(n), rbinom(n,1,prob=.3), age=age)
names(eg) <- c("surv", "status", "age")
rpart(Surv(surv, status)~age, data=eg)
# code 2
n <- 20
age <- rep(1:3, c(5, 5, 10))
eg<- data.frame(rexp(n), rbinom(n,1,prob=.3), age=age)
names(eg) <- c("surv", "status", "age")
rpart(Surv(surv, status)~age, data=eg)
# my session:
> library(rpart)
> # code 1
> n <- 20
> age <- rep(1:3, c(2, 3, 15))
> eg<- data.frame(rexp(n), rbinom(n,1,prob=.3),
age=age)
> names(eg) <- c("surv", "status", "age")
> rpart(Surv(surv, status)~age, data=eg)
Error in "$<-.data.frame"(`*tmp*`, "yval2", value =
c(1, 7)) :
replacement has 2 rows, data has 1
>
> # code 2
> n <- 20
> age <- rep(1:3, c(5, 5, 10))
> eg<- data.frame(rexp(n), rbinom(n,1,prob=.3),
age=age)
> names(eg) <- c("surv", "status", "age")
> rpart(Surv(surv, status)~age, data=eg)
n= 20
node), split, n, deviance, yval
* denotes terminal node
1) root 20 19.007310 1.000
2) age>=2.5 10 9.673372 0.8230355 *
3) age< 2.5 10 9.027225 1.1922660 *
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] some suggestion
Shame nobody has included this function in their package! Would be
useful to have as a standard function!
Sander.
Dimitris Rizopoulos wrote:
> take a look at this function by Kevin Wright
>
> RSiteSearch("sort.data.frame")
>
>
> Best,
> Dimitris
>
>
> Dimitris Rizopoulos
> Ph.D. Student
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
>
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/16/336899
> Fax: +32/16/337015
> Web: http://www.med.kuleuven.ac.be/biostat/
> http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
>
>
> - Original Message - From: "ronggui" <[EMAIL PROTECTED]>
> To:
> Sent: Tuesday, June 14, 2005 4:28 PM
> Subject: [R] some suggestion
>
>
>> it seems R has no function to sort the data.frame according to some
>> variable(s),though we can do these by order() and indexing.but why not
>> make sort() the a generic function,making it can sort vector and other
>> type of objects?
>>
>> maybe this is a silly suggestion,but i think it is quite usefull.
>>
>>
>>
>>
>> 2005-06-14
>>
>> --
>> Deparment of Sociology
>> Fudan University
>>
>> Blog:www.sociology.yculblog.com
>>
>>
>
>
>
>
>
>
>
>> __
>> R-help@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide!
>> http://www.R-project.org/posting-guide.html
>
>
>
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
--
Dr Sander P. Oom
Animal, Plant and Environmental Sciences,
University of the Witwatersrand
Private Bag 3, Wits 2050, South Africa
Tel (work) +27 (0)11 717 64 04
Tel (home) +27 (0)18 297 44 51
Fax +27 (0)18 299 24 64
Email [EMAIL PROTECTED]
Web www.oomvanlieshout.net/sander
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem installing packages with compiled-from-source R.app on Mac OS X - Tiger
You need to use R-patched, nowadays 2.1.1 beta, not the released R 2.1.0.
As far as I know it is fixed there.
On Tue, 14 Jun 2005, Constantinos Antoniou wrote:
> Hello all,
>
> This may be aimed for r-devel, but I encountered this as an R-user
> and not an R-developer so I start here (having said that, please
> direct me to R-devel if you think this is appropriate. I am not cross-
> posting, as I believe this is bad netiquette).
>
> I am a recent, but extremely happy R-user (especially after getting
> my own copy of MASS 2002). My adventures started when I wanted to use
> Rpy to use R also from within python. I compiled R (2.1.0) after
> configuring it as follows:
>
> ./configure --enable-R-shlib --with-blas='-framework vecLib' --with-
> lapack --with-aqua
>
> where --enable-R-shlib is required by Rpy.
>
> I then compiled Rpy-0.4.2.1 and I was able to call R from within
> python. So far so good...
>
> I then -naively- tried to start R.app. After less than a bounce on
> the dock... it would not start. [Sorry, I do not remember what
> problem was showing up in the console...]
>
> So, I figured I had to also compile R.app. I got the code via:
>
> svn co https://svn.r-project.org/R-packages/trunk/Mac-GUI Mac-
> GUI (executed this command yesterday)
>
> and compiled it as per the instructions... And R.app launches just fine.
>
> Now, the issue comes when I try to install a(ny) package (via the GUI
> package installer of R.app). I then get the following message:
>
> "Package installation failed. Package installation was not
> successful. Please see the R console for details"
>
> And the R console says:
>
> "Error in install.packages(c("dyn"), lib = "/Library/Frameworks/
> R.framework/Resources/library", :
> couldn't find function ".install.macbinary""
>
> I see that this has been encountered before and resolved (at least
> Prof. Ripley saw what was wrong) as per the following thread:
>
> https://stat.ethz.ch/pipermail/r-devel/2005-May/033115.html
>
> But I am not sure what needs to be done on my part so that I am not
> affected from it.
>
> > str(.Platform)
> List of 6
> $ OS.type : chr "unix"
> $ file.sep : chr "/"
> $ dynlib.ext: chr ".so"
> $ GUI : chr "X11"
> $ endian: chr "big"
> $ pkgType : chr "source"
> > R.version.string
> [1] "R version 2.1.0, 2005-04-18"
>
>
> gcc version 4.0.0
>
> [Please let me know what other information may be relevant]
>
>
> Thanking you in advance,
>
> Costas
>
>
> --
> Constantinos Antoniou, Ph.D.
> Department of Transportation Planning and Engineering
> National Technical University of Athens
> 5, Iroon Polytechniou str. GR-15773, Athens, Greece
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] using forecast() in dse2 with an ARMA model having a trend component
Scott This works for me: > arma.pred.without.trend$forecast[[1]][,1] [1] 106.0038 105.9789 105.9605 105.9396 105.9224 105.9052 105.8926 105.8849 [9] 105.8812 105.8880 105.9043 105.9240 105.9579 105.9878 105.9901 106.0095 [17] 106.0555 106.0782 106.0644 106.0427 106.0297 106.0072 106.0126 106.0125 > arma.pred.with.trend$forecast[[1]][,1] [1] 101.49566 97.46626 93.89069 90.70991 87.88602 85.37563 83.14843 [8] 81.17338 79.42193 77.87562 76.51150 75.30371 74.24589 73.30438 [15] 72.44303 71.69453 71.05658 70.47036 69.91559 69.41396 68.97527 [22] 68.57536 68.24555 67.94755 There was a bug that may be related to this fixed in the version dse_2005.4-1, which has been on CRAN for a month or so. The above result was actually with my working copy, but I don't recall any more recent changes that would affect this. If you actually have the 2005.4-1 version and are still having this problem then please let me know and I will check more carefully. (In that case, OS details would be helpful too.) Paul Gilbert Waichler, Scott R wrote: >(My apologies if this is a repeated posting. I couldn't find any trace >of my previous attempt in the archive.) > >I'm having trouble with forecast() in the dse2 package. It works fine >for me on a model without a trend, but gives me NaN output for the >forecast values when using a model with a trend. An example: > ># Set inputs and outputs for the ARMA model fit and test periods >arma.fit.input <- c(105.3332, 105.3573, 105.3113, 105.1493, 105.1209, >105.2111, 104.9161, >105.3654, 105.4682, 105.6789, 105.6297, 106.0155, >105.8454, 105.4322, >105.6062, 106.0739, 106.1109, 105.4470, 104.9739, >105.3427, 105.4305, >105.2563, 104.8501, 105.0358, 105.2827, 104.8977) > >arma.fit.output <- c(106.0376, 106.0514, 106.0716, 106.0570, 106.0442, >106.0414, 106.0375, > 106.0169, 106.0268, 106.0670, 106.1169, 106.1544, >106.1898, 106.2252, > 106.2605, 106.2959, 106.3324, 106.3974, 106.3460, >106.2357, 106.1897, > 106.1811, 106.1556, 106.1130, 106.0805, 106.0791) > >arma.pred.input <- c(104.9916, 104.8207, 104.8936, 104.8767, 104.9435, >104.8885, 104.9217, > 104.9029, 104.9508, 105.0065, 105.0557, 105.1982, >105.3392, 105.4007, > 105.6212, 105.5979, 105.2410, 105.4832, 105.8735, >105.5944, 105.1063, > 104.9809, 105.0821, 104.9362, 105.3037, 105.2322) >arma.pred.output <- c(106.0528, 106.0293, 106.0053, 105.9850, 105.9697, >105.9604, 105.9509, > 105.9430, 105.9357, 105.9314, 105.9333, 105.9420, >105.9640, 105.9994, > 106.0290, 106.0855, 106.1265, 106.1197, 106.1245, >106.1893, 106.2118, > 106.1503, 106.0883, 106.0511, 106.0194, 106.0221) > ># Set TSdata object >arma.fit.TSdata <- TSdata(input = arma.fit.input, output = >arma.fit.output) > ># Fit the model >arma.model.without.trend <- estVARXls(arma.fit.TSdata, max.lag=1, >trend=F) >arma.model.with.trend<- estVARXls(arma.fit.TSdata, max.lag=1, >trend=T) > ># Apply the model for the test period >arma.pred.TSdata <- TSdata(input = arma.pred.input, output = >arma.pred.output[1:2]) arma.pred.without.trend <- >forecast(TSmodel(arma.model.without.trend), arma.pred.TSdata) >arma.pred.with.trend<- forecast(TSmodel(arma.model.with.trend), >arma.pred.TSdata) > >The results: > > >>arma.pred.without.trend$forecast[[1]][,1] >> >> > [1] 106.0038 105.9789 105.9605 105.9396 105.9224 105.9052 105.8926 >105.8849 [9] 105.8812 105.8880 105.9043 105.9240 105.9579 105.9878 >105.9901 106.0095 [17] 106.0555 106.0782 106.0644 106.0427 106.0297 >106.0072 106.0126 106.0125 > > >>arma.pred.with.trend$forecast[[1]][,1] >> >> > [1] 5.76056e+228 NaN NaN NaN NaN > [6] NaN NaN NaN NaN NaN >[11] NaN NaN NaN NaN NaN >[16] NaN NaN NaN NaN NaN >[21] NaN NaN NaN NaN > >I read help on this function and the PDF manuals but can't see what I >might be missing. >Any ideas? > >Thanks, Sott Waichler >Pacific Northwest National Laboratory >[EMAIL PROTECTED] > >__ >R-help@stat.math.ethz.ch mailing list >https://stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Plotting quiver vector tensor arrows 2d field data
Hi All,
I'd like to plot something like
http://www.nawcwpns.navy.mil/~weather/mugu/mesodata/analysis.html
Looking through the galleries at
http://addictedtor.free.fr/graphiques/allgraph.php
http://r-spatial.sourceforge.net/gallery/
http://fawn.unibw-hamburg.de/cgi-bin/Rwiki.pl?GraphGallery
demo(graphics)
I did not find a function to plot a 2d field on a matrix. I did find
mention of a quiver function in the archives. Is this the best solution
or are there other tools I missed?
quiver<- function(u,v,scale=1,length=0.05)
# first stab at matlab's quiver in R
# from http://tolstoy.newcastle.edu.au/R/help/01c/2711.html
# Robin Hankin Tue 20 Nov 2001 - 13:10:28 EST
{
xpos <- col(u)
ypos <- max(row(u))-row(u)
speed <- sqrt(u*u+v*v)
maxspeed <- max(speed)
u <- u*scale/maxspeed
v <- v*scale/maxspeed
matplot(xpos,ypos,type="p",cex=0)
arrows(xpos,ypos,xpos+u,ypos+v,length=length*min(par.uin()))
}
par.uin <- function()
# determine scale of inches/userunits in x and y
# from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
# Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
{
u <- par("usr")
p <- par("pin")
c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
}
u <- matrix(rnorm(100),nrow=10)
v <- matrix(rnorm(100),nrow=10)
quiver(u,v)
I added these functions as an example to the Wiki:
http://fawn.unibw-hamburg.de/cgi-bin/Rwiki.pl?GraphGallery
http://fawn.unibw-hamburg.de/cgi-bin/Rwiki.pl?QuiverPlot
Thanks for your time,
Dave
--
Dr. David Forrest
[EMAIL PROTECTED](804)684-7900w
[EMAIL PROTECTED] (804)642-0662h
http://maplepark.com/~drf5n/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Survfit,newdata and continuous time varying covariates
Hi All, I am working with three time varying covariates in a coxph model. I cannot seem to figure out how to use survfit and the newdata argument to provide estimated survival curves for two scenarios of one covariate while holding the other two at the mean value. Is it possible to display how estimated survival depends on alternative scenarios/profiles of a time varying covariate? Alex Hanke Department of Fisheries and Oceans St. Andrews Biological Station 531 Brandy Cove Road St. Andrews, NB Canada E5B 2L9 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] ANN: BioC2005 Conference scheduled for August in Seattle
== BioC2005 Where Software and Biology Connect == http://bioconductor.org/meeting05/ About BioC2005 == This conference will highlight current developments within and beyond Bioconductor, a world-wide open source and open development software project for the analysis and comprehension of genomic data. Our goal is to provide a forum in which to discuss the use and design of software for analyzing data arising in biology with a focus on Bioconductor and genomic data. Where: Fred Hutchinson Cancer Research Center Seattle Wa. When: August 15 and 16, 2005 What: Morning Talks: 8:30-12:00 Afternoon Practicals: 2:00-5:00 Tuesday Evening 5:00-7:30 Posters and Wine & Cheese Fees: 250 USD for academic/research institute attendees 125 USD for enrolled full-time students DETAILS: http://bioconductor.org/meeting05/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to fix false positve rates?
Dear R-users, I have a set of 12000 image samples. I can divide this set into two categories: training and testing. I need to classify the test set into a two qualitative outputs: true or false for some characteristic. To do the classification I'm using the packages SVM in e1071 library and LDA in the MASS library. However, I'm with a great number of FALSE POSITIVE CASES in both classifiers, about 10/15%. 1) How can I use these classifiers with a fixed 1% FALSE POSITIVE RATE? 2) Could anynone indicate me a non-linear SVM classifier R-package? 3) Could anynone indicate me an ADA-BOOST classifier R-package? Thanks for all. I'm anxiously wait the answers. Best regards. * |Üñð뮧µñ| - Anderson de Rezende Rocha Bacharel em Ciência da Computação - UFLA Mestrando em Ciência da Computação - UNICAMP Esteganografia e esteganálise digital UNIVERSIDADE ESTADUAL DE CAMPINAS, SP - BRASIL < http://andersonrocha.cjb.net > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Plotting rows (or columns) from a matrix in different graphs, not using "par"
Dear R-users,
I would like to ask whether it's possible (for sure it would be), to
plot each rows (or columns) in different graphs and in the same figure
region without using the function "par" and then struggling around with
"axes" and labels etc.
Luckily, I would always have "rows + columns = even number" and the same
"ylim".
The next one could be a sort of example on what I would like to avoid
plotting the rows of the matrix "mat":
### EXAMPLE ##
dat <- sort(runif(16, 1, 1000))
mat <- matrix(dat, ncol=4, byrow = T)
y <- seq(1950, 1953)
par(mfrow=c(2,2))
plot(y, mat[1,], ylim=c(1,1000), axes=F, xlab="", ylab="simul.")
box()
axis(side=2, tick=T, labels=T)
axis(side=3, tick=T, labels=T)
plot(y, mat[2,], ylim=c(1,1000), axes=F, xlab="", ylab="")
box()
axis(side=3, tick=T, labels=T)
plot(y, mat[3,], ylim=c(1,1000), axes=F, xlab="years", ylab="simul.")
box()
axis(side=1, tick=T, labels=T)
axis(side=2, tick=T, labels=T)
plot(y, mat[4,], ylim=c(1,1000), axes=F, xlab="years", ylab="")
box()
axis(side=1, tick=T, labels=T)
### END EXAMPLE
Naturally something more compact would be even nicer.
Thanks in advance,
Carlo Giovanni Camarda
+
This mail has been sent through the MPI for Demographic Rese...{{dropped}}
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] any function to calculate the lamda coefficient?
i just write one for myself,but the result is different from the spss's.
i anyone can send me a copy of function calc the lamda coefficient?(i know i
should google first,but right now i can surf the internet.)
> lamda
function(table){
obs1.y<-sum(apply(table,1,max))
obs2.y<-max(apply(table,1,sum))
obs1.x<-sum(apply(table,2,max))
obs2.x<-max(apply(table,2,sum))
D.x<-(sum(table)-obs2.x)
D.y<-(sum(table)-obs2.y)
lamda.y<-(obs1.y-obs2.y)/D.y
lamda.x<-(obs1.x-obs2.x)/D.x
wt<-c(D.x,D.y)
lamda.xy<-weighted.mean(c(lamda.x,lamda.y),wt/sum(wt))
cat(
"the lamda.y coefficient is :",lamda.y,";","the lamda.x coefficient is
:",lamda.x,";",
"\n",
"the mean lamda coefficient is :",lamda.xy,".",fill=T)
}
> a<-matrix(c(44,20,54,8),2)
> a
[,1] [,2]
[1,] 44 54
[2,] 208
> lamda(a)
the lamda.y coefficient is : -0.8571429 ; the lamda.x coefficient is :
0.5483871 ;
the mean lamda coefficient is : 0.111 .
2005-06-15
--
Deparment of Sociology
Fudan University
Blog:www.sociology.yculblog.com
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] using forecast() in dse2 with an ARMA model having a trend component
(My apologies if this is a repeated posting. I couldn't find any trace of my previous attempt in the archive.) I'm having trouble with forecast() in the dse2 package. It works fine for me on a model without a trend, but gives me NaN output for the forecast values when using a model with a trend. An example: # Set inputs and outputs for the ARMA model fit and test periods arma.fit.input <- c(105.3332, 105.3573, 105.3113, 105.1493, 105.1209, 105.2111, 104.9161, 105.3654, 105.4682, 105.6789, 105.6297, 106.0155, 105.8454, 105.4322, 105.6062, 106.0739, 106.1109, 105.4470, 104.9739, 105.3427, 105.4305, 105.2563, 104.8501, 105.0358, 105.2827, 104.8977) arma.fit.output <- c(106.0376, 106.0514, 106.0716, 106.0570, 106.0442, 106.0414, 106.0375, 106.0169, 106.0268, 106.0670, 106.1169, 106.1544, 106.1898, 106.2252, 106.2605, 106.2959, 106.3324, 106.3974, 106.3460, 106.2357, 106.1897, 106.1811, 106.1556, 106.1130, 106.0805, 106.0791) arma.pred.input <- c(104.9916, 104.8207, 104.8936, 104.8767, 104.9435, 104.8885, 104.9217, 104.9029, 104.9508, 105.0065, 105.0557, 105.1982, 105.3392, 105.4007, 105.6212, 105.5979, 105.2410, 105.4832, 105.8735, 105.5944, 105.1063, 104.9809, 105.0821, 104.9362, 105.3037, 105.2322) arma.pred.output <- c(106.0528, 106.0293, 106.0053, 105.9850, 105.9697, 105.9604, 105.9509, 105.9430, 105.9357, 105.9314, 105.9333, 105.9420, 105.9640, 105.9994, 106.0290, 106.0855, 106.1265, 106.1197, 106.1245, 106.1893, 106.2118, 106.1503, 106.0883, 106.0511, 106.0194, 106.0221) # Set TSdata object arma.fit.TSdata <- TSdata(input = arma.fit.input, output = arma.fit.output) # Fit the model arma.model.without.trend <- estVARXls(arma.fit.TSdata, max.lag=1, trend=F) arma.model.with.trend<- estVARXls(arma.fit.TSdata, max.lag=1, trend=T) # Apply the model for the test period arma.pred.TSdata <- TSdata(input = arma.pred.input, output = arma.pred.output[1:2]) arma.pred.without.trend <- forecast(TSmodel(arma.model.without.trend), arma.pred.TSdata) arma.pred.with.trend<- forecast(TSmodel(arma.model.with.trend), arma.pred.TSdata) The results: > arma.pred.without.trend$forecast[[1]][,1] [1] 106.0038 105.9789 105.9605 105.9396 105.9224 105.9052 105.8926 105.8849 [9] 105.8812 105.8880 105.9043 105.9240 105.9579 105.9878 105.9901 106.0095 [17] 106.0555 106.0782 106.0644 106.0427 106.0297 106.0072 106.0126 106.0125 > arma.pred.with.trend$forecast[[1]][,1] [1] 5.76056e+228 NaN NaN NaN NaN [6] NaN NaN NaN NaN NaN [11] NaN NaN NaN NaN NaN [16] NaN NaN NaN NaN NaN [21] NaN NaN NaN NaN I read help on this function and the PDF manuals but can't see what I might be missing. Any ideas? Thanks, Sott Waichler Pacific Northwest National Laboratory [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Dateticks
try this example:
> x.1 <- strptime("6/17/03",'%m/%d/%y')
> x.1
[1] "2003-06-17"
> plot(0:250, xaxt='n')
> dates <- x.1 + c(0,50,100,150,200,250) * 86400 # 'dates' is in seconds,
so add the appropriate number of days
> dates
[1] "2003-06-17 00:00:00 EDT" "2003-08-06 00:00:00 EDT" "2003-09-25
00:00:00 EDT"
[4] "2003-11-13 23:00:00 EST" "2004-01-02 23:00:00 EST" "2004-02-21
23:00:00 EST"
> axis(1, at=c(0,50,100,150,200,250), labels=format(dates,"%m/%d/%y")) #
format the output
>
Jim
__
James Holtman"What is the problem you are trying to solve?"
Executive Technical Consultant -- Convergys Labs
[EMAIL PROTECTED]
+1 (513) 723-2929
"Bernard L. Dillard"
<[EMAIL PROTECTED]> To:
r-help@stat.math.ethz.ch
Sent by: cc:
[EMAIL PROTECTED]Subject: [R] Dateticks
ath.ethz.ch
06/14/2005 12:27
Hello. I am having the worst time converting x-axis date ticks to real
dates. I have tried several suggestions in online help tips and books to
no avail.
For example, the x-axis has 0, 50, 100, etc, and I want it to have
"6/17/03", "8/6/03" etc. See attached (sample).
Can anybody help me with this.
Here's my code:
ts.plot(date.attackmode.table[,1], type="l", col="blue", lty=2,ylab="IED
Attacks", lwd=2,xlab="Attack Dates",main="Daily Summary of Attack
Mode")
grid()
Thanks for your help if possible.
(See attached file: sample.pdf)
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
sample.pdf
Description: Adobe PDF document
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] problem installing packages with compiled-from-source R.app on Mac OS X - Tiger
Hello all,
This may be aimed for r-devel, but I encountered this as an R-user
and not an R-developer so I start here (having said that, please
direct me to R-devel if you think this is appropriate. I am not cross-
posting, as I believe this is bad netiquette).
I am a recent, but extremely happy R-user (especially after getting
my own copy of MASS 2002). My adventures started when I wanted to use
Rpy to use R also from within python. I compiled R (2.1.0) after
configuring it as follows:
./configure --enable-R-shlib --with-blas='-framework vecLib' --with-
lapack --with-aqua
where --enable-R-shlib is required by Rpy.
I then compiled Rpy-0.4.2.1 and I was able to call R from within
python. So far so good...
I then -naively- tried to start R.app. After less than a bounce on
the dock... it would not start. [Sorry, I do not remember what
problem was showing up in the console...]
So, I figured I had to also compile R.app. I got the code via:
svn co https://svn.r-project.org/R-packages/trunk/Mac-GUI Mac-
GUI (executed this command yesterday)
and compiled it as per the instructions... And R.app launches just fine.
Now, the issue comes when I try to install a(ny) package (via the GUI
package installer of R.app). I then get the following message:
"Package installation failed. Package installation was not
successful. Please see the R console for details"
And the R console says:
"Error in install.packages(c("dyn"), lib = "/Library/Frameworks/
R.framework/Resources/library", :
couldn't find function ".install.macbinary""
I see that this has been encountered before and resolved (at least
Prof. Ripley saw what was wrong) as per the following thread:
https://stat.ethz.ch/pipermail/r-devel/2005-May/033115.html
But I am not sure what needs to be done on my part so that I am not
affected from it.
> str(.Platform)
List of 6
$ OS.type : chr "unix"
$ file.sep : chr "/"
$ dynlib.ext: chr ".so"
$ GUI : chr "X11"
$ endian: chr "big"
$ pkgType : chr "source"
> R.version.string
[1] "R version 2.1.0, 2005-04-18"
gcc version 4.0.0
[Please let me know what other information may be relevant]
Thanking you in advance,
Costas
--
Constantinos Antoniou, Ph.D.
Department of Transportation Planning and Engineering
National Technical University of Athens
5, Iroon Polytechniou str. GR-15773, Athens, Greece
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] KMEANS output...
help(kmeans) describes the return value: a list with a) the complete classification b) the cluster centers c) the within-cluster sums of squares b) the cluster sizes. Thankfully, you don't need any parsing. Reid Huntsinger -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Omer Bakkalbasi Sent: Tuesday, June 14, 2005 12:40 PM To: r-help@stat.math.ethz.ch Subject: [R] KMEANS output... Using R 2.1.0 on Windows 2 questions: 1. Is there a way to parse the output from kmeans within R? 2. If the answer to 1. is convoluted or impossible, how do you save the output from kmeans in a plain text file for further processing outside R? Example: > ktx<-kmeans(x,12, nstart = 200) I would like to parse ktx within R to extract cluster sizes, sum-of-squares values, etc., OR save ktx in a plain text file in Windows to post-process it with a parser I would have to write. Thanks! Omer Cell: (914) 671-7447 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] KMEANS output...
Using R 2.1.0 on Windows 2 questions: 1. Is there a way to parse the output from kmeans within R? 2. If the answer to 1. is convoluted or impossible, how do you save the output from kmeans in a plain text file for further processing outside R? Example: > ktx<-kmeans(x,12, nstart = 200) I would like to parse ktx within R to extract cluster sizes, sum-of-squares values, etc., OR save ktx in a plain text file in Windows to post-process it with a parser I would have to write. Thanks! Omer Cell: (914) 671-7447 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] plots
Dear all, Is it possible to change the levels in a mosaic plot, the appearance of the level or the levels size? For instance: A C E B D Thanks for your help Adrián __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
> On Tue, 14 Jun 2005, Prof Brian Ripley wrote:
> If your file system does not like 15000 files you can always
> save in a DBMS.
Or, switch to a better/more appropriate file system:
http://en.wikipedia.org/wiki/Comparison_of_file_systems
ReiserFS would allow you to store up to about 1.2 million files in a directory.
> -Original Message-
> From: Prof Brian Ripley [mailto:[EMAIL PROTECTED]
> Sent: Tuesday, June 14, 2005 10:41 AM
> To: Barry Rowlingson
> Cc: r-help@stat.math.ethz.ch; Richard Mott
> Subject: Re: [R] load ing and saving R objects
>
>
> On Tue, 14 Jun 2005, Barry Rowlingson wrote:
>
> > Richard Mott wrote:
> >> Does anyone know a way to do the following:
> >>
> >> Save a large number of R objects to a file (like load()
> does) but then
> >> read back only a small named subset of them . As far as I can see,
> >> load() reads back everything.
> >
> > Save them to individual files when you generate them?
> >
> > for(i in 1:15000){
> >
> > m=generateBigMatrix(i)
> >
> > filename=paste("BigMatrix-",i,".Rdata",sep='')
> > save(m,file=filename)
> > }
> >
> > Note that load will always overwrite 'm', so to load a
> sample of them in
> > you'll need to do something like this:
> >
> > bigSamples=list()
> >
> > for(i in sample(15000,N)){
> >filename=paste("BigMatrix-",i,".Rdata",sep='')
> >load(filename)
> >bigSamples[[i]]=m
> > }
> >
> > But there may be a more efficient way to string up a big list like
> > that, I can never remember - get it working, then worry
> about optimisation.
>
> (Yes, use bigSamples <- vector("list", 15000) first.)
>
> > I hope your filesystem is happy with 15000 objects in it. I would
> > dedicate a folder or directory for just these objects'
> files, since it
> > then becomes near impossible to see anything other than the
> big matrix
> > files...
>
> .readRDS/.saveRDS might be a better way to do this, and avoids always
> restoring to "m".
>
> If your file system does not like 15000 files you can always
> save in a
> DBMS.
>
> I did once look into restoring just some of the objects in a save()ed
> file, but it is not really possible to do so efficiently due
> to sharing
> between objects.
>
> --
> Brian D. Ripley, [EMAIL PROTECTED]
> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel: +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax: +44 1865 272595
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Dateticks
Hello. I am having the worst time converting x-axis date ticks to real dates. I have tried several suggestions in online help tips and books to no avail. For example, the x-axis has 0, 50, 100, etc, and I want it to have "6/17/03", "8/6/03" etc. See attached (sample). Can anybody help me with this. Here's my code: ts.plot(date.attackmode.table[,1], type="l", col="blue", lty=2,ylab="IED Attacks", lwd=2,xlab="Attack Dates",main="Daily Summary of Attack Mode") grid() Thanks for your help if possible. sample.pdf Description: Adobe PDF document __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Calling C from Fortran
I would like to call C routines from Fortran as suggested in section 5.6 of
the "Writing R extensions" documentation.
I'm familiar with Fortran but not with C.
I understand the example provided in Fortran:
subroutine testit()
double precision normrnd, x
call rndstart()
x = normrnd()
call dblepr("X was", 5, x, 1)
call rndend()
end
but I don't understand the purpose of this C wrapper:
#include
void F77_SUB(rndstart)(void) { GetRNGstate(); }
void F77_SUB(rndend)(void) { PutRNGstate(); }
double F77_SUB(normrnd)(void) { return norm_rand(); }
neither how I should compile it.
Could anyone explain how I should compile and link
the C and Fortran files above, and call the Fortran subroutine from R.
Thanks,
Gilles
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Logistic regression with more than two choices
Dear all R-users, I am a new user of R and I am trying to build a discrete choice model (with more than two alternatives A, B, C and D) using logistic regression. I have data that describes the observed choice probabilities and some background information. An example below describes the data: Sex Age pr(A) pr(B) pr(C) pr(D) ... 1 11 0.5 0.5 0 0 1 40 1 0 0 0 0 34 0 0 0 1 0 64 0.1 0.5 0.2 0.2 ... I have been able to model a case with only two alternatives "A" and "not A" by using glm(). I do not know what functions are available to estimate such a model with more than two alternatives. Multinom() is one possibility, but it only allows the use of binary 0/1-data instead of observed probabilities. Did I understand this correctly? Additionally, I am willing to use different independent variables for the different alternatives in the model. Formally, I mean that: Pr(A)=exp(uA)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) Pr(B)=exp(uB)/(exp(uA)+exp(uB)+exp(uC)+exp(uD) ... where uA, uB, uC and uD are linear functions with different independent variables, e.g. uA=alpha_A1*Age, uB=alpha_B1*Sex. Do you know how to estimate this type of models in R? Best regards, Ville Koskinen __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] extracting components of a list
Dimitris, wouldn't this be more precise --- > sapply(jj,function(x) which(x$b[1]==4)) [[1]] [1] 1 [[2]] numeric(0) [[3]] [1] 1 John Dimitris wrote --- maybe something like this: jj <- list(list(a = 1, b = 4:7), list(a = 5, b = 3:6), list(a = 10, b = 4:5)) ### jj[sapply(jj, function(x) x$b[1] == 4)] I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.ac.be/biostat/ http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm - Original Message - From: "Robin Hankin" <[EMAIL PROTECTED]> To: "r-help" Sent: Monday, June 13, 2005 4:23 PM Subject: [R] extracting components of a list > Hi > > how do I extract those components of a list that satisfy a certain > requirement? If > > jj <- list(list(a=1,b=4:7),list(a=5,b=3:6),list(a=10,b=4:5)) > > > I want just the components of jj that have b[1] ==4 which in this > case > would be the first and > third of jj, vizlist (jj[[1]],jj[[3]]). > > How to do this efficiently? > > My only idea was to loop through jj, and set unwanted components to > NULL, but > FAQ 7.1 warns against this. > > > > > -- > Robin Hankin > Uncertainty Analyst > National Oceanography Centre, Southampton > European Way, Southampton SO14 3ZH, UK > tel 023-8059-7743 > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] t.test using RSPerl
This has nothing to do with RSPerl, instead it has to do what kind of
object you obtain and how these are print():ed. Typing the name of an
object, say, 'res', at R prompt and pressing ENTER;
> res
is equivalent as typing
> print(res)
This is for convenience to the user. Basically, this is why you can do
> 1+1
[1] 2
without having to do
> print(1+1)
[1] 2
When you "transfer" the object from R to Perl you will, as expect, only
receive the value of 'res', not the output from 'print(res)'. Here is an
example illustrating the behavior in R:
> res <- t.test(1:10,y=c(7:20))
> length(res)
[1] 9
> str(res)
List of 9
$ statistic : Named num -5.43
..- attr(*, "names")= chr "t"
$ parameter : Named num 22
..- attr(*, "names")= chr "df"
$ p.value: num 1.86e-05
$ conf.int : atomic [1:2] -11.05 -4.95
..- attr(*, "conf.level")= num 0.95
$ estimate : Named num [1:2] 5.5 13.5
..- attr(*, "names")= chr [1:2] "mean of x" "mean of y"
$ null.value : Named num 0
..- attr(*, "names")= chr "difference in means"
$ alternative: chr "two.sided"
$ method : chr "Welch Two Sample t-test"
$ data.name : chr "1:10 and c(7:20)"
- attr(*, "class")= chr "htest"
> print(res)
Welch Two Sample t-test
data: 1:10 and c(7:20)
t = -5.4349, df = 21.982, p-value = 1.855e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-11.052802 -4.947198
sample estimates:
mean of x mean of y
5.5 13.5
In your case, to capture what print(res) is outputting, you may want to
look at ?capture.out, that is
> output <- capture.output(res)
> output # ...that is, print(output)
[1] ""
[2] "\tWelch Two Sample t-test"
[3] ""
[4] "data: 1:10 and c(7:20) "
[5] "t = -5.4349, df = 21.982, p-value = 1.855e-05"
[6] "alternative hypothesis: true difference in means is not equal to 0 "
[7] "95 percent confidence interval:"
[8] " -11.052802 -4.947198 "
[9] "sample estimates:"
[10] "mean of x mean of y "
[11] " 5.5 13.5 "
[12] ""
and transfer 'output' to Perl.
Cheers
Henrik
Wagle, Mugdha wrote:
> Hi,
>
> I've just started using R and RSPerl. I have some code as follows:
>
> &R::initR("--no-save");
> &R::call("t.test", ([EMAIL PROTECTED], [EMAIL PROTECTED]));
>
> where @array1 and @array2 are both 1-dimensional arrays in Perl having 54675
> elements each. On execution the output is as follows:
>
> Calling R function name `t.test', # arguments: 3
> 1) Arg type 3
> Got a reference to a value 10
> Here now!2) Arg type 3
> Got a reference to a value 10
> Here now!Calling R
> t.test(c(0, 6.24280675278087, 6.35175793656943, 5.76925805661511,
> 7.0789316246711, 7.4636498661157, 8.13730810691084, 8.78203131644273,
> 9.64502765609435, 9.95631242346133, 5.83129579495516, 6.8798700754926,
> 7.31814159140937...(REST OF THE ARRAY ELEMENTS).
> 4.91632461462501, 3.38099467434464,
> 3.91800507710569, 3.23867845216438, 3.38439026334577, 4.64918707140487,
> 3.23474917402449, 3.62966009445396, 3.36729582998647, 3.91999117507732
> ))
> Performed the call, result has length 9
>
> My question is : with other functions such as sum and log10, the actual
> values of the result are displayed. Here the call seems to have worked but
> the output is not what you get when running t.test directly on the R command
> prompt..
>
> data: data4[2] and data4[3]
> t = 0.2186, df = 109.847, p-value = 0.8274
> alternative hypothesis: true difference in means is not equal to 0
> 95 percent confidence interval:
> -3722.830 4645.723
> sample estimates:
> mean of x mean of y
> 6185.139 5723.693
>
> which is what I had expected, after seeing the outputs in the case of simpler
> functions like sum. Could anyone please tell me how I can obtain the output I
> expect(i.e. the same as the command line outputgiving values of t,
> p-value and the means)?
>
> Thank you very much for the help!!
>
> Sincerely,
> Mugdha Wagle
> Hartwell Center for Bioinformatics and Biotechnology,
> St.Jude Children's Research Hospital, Memphis TN 38105
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] t.test using RSPerl
Hi,
I've just started using R and RSPerl. I have some code as follows:
&R::initR("--no-save");
&R::call("t.test", ([EMAIL PROTECTED], [EMAIL PROTECTED]));
where @array1 and @array2 are both 1-dimensional arrays in Perl having 54675
elements each. On execution the output is as follows:
Calling R function name `t.test', # arguments: 3
1) Arg type 3
Got a reference to a value 10
Here now!2) Arg type 3
Got a reference to a value 10
Here now!Calling R
t.test(c(0, 6.24280675278087, 6.35175793656943, 5.76925805661511,
7.0789316246711, 7.4636498661157, 8.13730810691084, 8.78203131644273,
9.64502765609435, 9.95631242346133, 5.83129579495516, 6.8798700754926,
7.31814159140937...(REST OF THE ARRAY ELEMENTS).
4.91632461462501, 3.38099467434464,
3.91800507710569, 3.23867845216438, 3.38439026334577, 4.64918707140487,
3.23474917402449, 3.62966009445396, 3.36729582998647, 3.91999117507732
))
Performed the call, result has length 9
My question is : with other functions such as sum and log10, the actual values
of the result are displayed. Here the call seems to have worked but the output
is not what you get when running t.test directly on the R command prompt..
data: data4[2] and data4[3]
t = 0.2186, df = 109.847, p-value = 0.8274
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3722.830 4645.723
sample estimates:
mean of x mean of y
6185.139 5723.693
which is what I had expected, after seeing the outputs in the case of simpler
functions like sum. Could anyone please tell me how I can obtain the output I
expect(i.e. the same as the command line outputgiving values of t, p-value
and the means)?
Thank you very much for the help!!
Sincerely,
Mugdha Wagle
Hartwell Center for Bioinformatics and Biotechnology,
St.Jude Children's Research Hospital, Memphis TN 38105
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
Thanks everyone for help on this. It looks like the solution is a zillion files. Richard -- Richard Mott | Wellcome Trust Centre tel 01865 287588 | for Human Genetics fax 01865 287697 | Roosevelt Drive, Oxford OX3 7BN __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] plotting confidence intervals for polynomial models? (was Re: ordinary polynomial coefficients from orthogonal polynomials?)
What I really want now are the 95% confidence intervals that I
mistakenly thought that linesHyperb.lm() would provide:
> pm <- lm(billions ~ I(decade) + I(decade^2) + I(decade^3))
> library(sfsmisc)
> linesHyperb.lm(pm)
Error in "names<-.default"(`*tmp*`, value = c("I(decade)",
"I(decade^2)", :
names attribute [3] must be the same length as the vector [1]
> pm <- lm(billions ~ poly(decade, 3))
> linesHyperb.lm(pm)
Warning message:
'newdata' had 100 rows but variable(s) found have 5 rows
Shouldn't curve(predict(...), add=TRUE) be able to plot confidence
interval bands?
Prof Brian Ripley wrote:
> Why do `people' need `to deal with' these, anyway. We have machines to
> do that.
Getting a 0.98 adjusted R^2 on the first try, compels me to
try to publish the fitted formula.
> decade <- c(1950, 1960, 1970, 1980, 1990)
> billions <- c(3.5, 5, 7.5, 13, 40)
> # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg
>
> pm <- lm(billions ~ poly(decade, 3))
>
> plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000),
main="average yearly inflation-adjusted dollar cost of extreme weather
events worldwide")
> curve(predict(pm, data.frame(decade=x)), add=TRUE)
> # output: http://www.bovik.org/storms.gif
>
> summary(pm)
Call:
lm(formula = billions ~ poly(decade, 3))
Residuals:
1 2 3 4 5
0.2357 -0.9429 1.4143 -0.9429 0.2357
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)13.800 0.882 15.647 0.0406 *
poly(decade, 3)1 25.614 1.972 12.988 0.0489 *
poly(decade, 3)2 14.432 1.972 7.318 0.0865 .
poly(decade, 3)36.483 1.972 3.287 0.1880
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.972 on 1 degrees of freedom
Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829
F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] use of gam
On Mon, 13 Jun 2005, Kerry Bush wrote: > > Suppose I fit the following model: > >> library(gam) > >> fit <- gam(y~x1+x2+s(x3),family=binomial) > > and then I use > >> fitf$coef > x1x2 s(x3) > 4.1947460 2.7967200 0.0788252 > > are the coefficients for x1 and x2 the estimated > coefficients? Yes. >what is the meaning of s(x3)? since this > is a non-parametric component, it may not belong here > as a coefficient, am I right? I believe that it is coefficient of the smooth term if the smooth term is standardized to a standard deviation of 1. If you consider the shape of the smooth term as fixed, the model looks like a glm(). This representation can be used to compute approximate standard errors for the linear terms (the approximation isn't very good if s(x3) has too many degrees of freedom, and Trevor Hastie et al have recently worked out a consistent estimator of the standard errors). -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] c(recursive=TRUE)
Hi R users,
I am currently using c(...,recursive=TRUE) to handle list-structured
objects.
This allows to represent something like:
> l1 = list(level1=1,level2=list(sub1=1,sub2=2))
as:
> (l1names = names(c(l1,recursive=TRUE)))
[1] "level1" "level2.sub1" "level2.sub2"
Then, one can use:
> (l1names = sapply(l1names,FUN=function(element)
strsplit(element,"\\.")))
$level1
[1] "level1"
$level2.sub1
[1] "level2" "sub1"
$level2.sub2
[1] "level2" "sub2"
Which allows to access to individuals terminal nodes of list using
indexing, as
l1[[c("level2","sub2")]] nicely works in R.
For my very specific uses, it unfortunatelt happens that some elements
of my list do have names that contains a dot (for example:
decimal.mark). This is very frustating as a consequence is that my
function does not work in that case...
One workaround would be to add an argument to c(), as delimitor='.' by
default, allowing to put a less usual delimitor (say "§"). But this
would mean changing in C (feature request?). Same for "unlist".
Is there any other easy way to get the structure of a list?
Eric
Eric Lecoutre
UCL / Institut de Statistique
Voie du Roman Pays, 20
1348 Louvain-la-Neuve
Belgium
tel: (+32)(0)10473050
[EMAIL PROTECTED]
http://www.stat.ucl.ac.be/ISpersonnel/lecoutre
If the statistics are boring, then you've got the wrong numbers. -Edward
Tufte
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] New Family object for GLM models...
On Tue, 14 Jun 2005, Abatih E. wrote:
> Dear R-Users, I wish to create a new family object based on the Binomial
> family. The only difference will be with the link function. Thus instead
> if using the 'logit(u)' link function, i plan to use '-log(i-u)'. So
> far, i have tried to write the function following that of the Binomial
> and Negative Binomial families. The major problem i have here is with
> the definition of the initial values. The definitions i have attempted
> so far don't yield convergence when implemented in a model.
Are you sure the problem is with initial values? Your link function
is capable of producing mu<0 from large enough negative eta. If the MLE
is on the boundary of the parameter space it will not solve the likelihood
equations and so won't be found by iterative reweighted least squares. It
might be found by step-halving in glm, but that isn't guaranteed.
If the mle is in the interior of the parameter space then in my experience
the starting values aren't terribly important (though my experience is
with the log-binomial rather than -log(1-mu) link).
-thomas
> I still
> can't figure out how the initial values are defined. I don't know if it
> is based on some property of the fitted model or the domain and /or
> range of the link function. I wish someone could help me provide a
> solution to this problem. I have appended the function i wrote.
>
>
> Add.haz<-function () {
>
>
>
>env <- new.env(parent = .GlobalEnv)
>
>assign(".nziek", nziek, envir = env)
>
>famname<-"Addhaza"
>
>link="addlink"
>
>linkfun<-function(mu) -log(1-mu)
>
>linkinv<-function(eta) 1-exp(-eta)
>
>variance<-function(mu) mu*(1-mu)
>
>validmu<-function(mu) all(mu > 0) && all(mu < 1)
>
>mu.eta<-function(mu) 1/(1-mu)
>
>dev.resids <- function(y, mu, wt) 2 * wt * (y * log(ifelse(y ==
>
>0, 1, y/mu)) + (1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 -
>
>mu
>
>aic <- function(y, n, mu, wt, dev) {
>
>m <- if (any(n > 1))
>
>n
>
>else wt
>
>-2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
>
>y), round(m), mu, log = TRUE))
>
>}
>
>
>
>initialize<- expression({
>
>n<-rep(1,nobs)
>
>if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
>
> mustart <- (weights * y + 0.5)/(weights + 1)
>
>m <- weights * y
>
>if (any(abs(m - round(m)) > 0.001)) warning("non-integer
> #successes in a binomial glm!")
>
>
>
>
>
> })
>
>
>
>environment(variance) <- environment(validmu) <- environment(dev.resids)
> <- environment(aic) <- env
>
>
>
>structure(list(family = famname, link = link, linkfun = linkfun,
>
>linkinv = linkinv, variance = variance, dev.resids = dev.resids,
>
>aic = aic, mu.eta =mu.eta, initialize = initialize,
>
>validmu = validmu), class = "family")
> }
>
> Thank you for your kind attention.
> Emmanuel
>
>
> EMMANUEL NJI ABATIH
> Rues des deux Eglises 140
> 1210 St Josse-Ten-Node, Bruxelles, BELGIUM
> MOBLIE: 0032486958988(anytime)
> Fix: 0032 2642 5038(8am to 6pm)
> http://www.iph.fgov.be/epidemio/epien
>
> -
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
Thomas Lumley Assoc. Professor, Biostatistics
[EMAIL PROTECTED] University of Washington, Seattle
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
On Tue, 14 Jun 2005, Barry Rowlingson wrote:
> Richard Mott wrote:
>> Does anyone know a way to do the following:
>>
>> Save a large number of R objects to a file (like load() does) but then
>> read back only a small named subset of them . As far as I can see,
>> load() reads back everything.
>
> Save them to individual files when you generate them?
>
> for(i in 1:15000){
>
> m=generateBigMatrix(i)
>
> filename=paste("BigMatrix-",i,".Rdata",sep='')
> save(m,file=filename)
> }
>
> Note that load will always overwrite 'm', so to load a sample of them in
> you'll need to do something like this:
>
> bigSamples=list()
>
> for(i in sample(15000,N)){
>filename=paste("BigMatrix-",i,".Rdata",sep='')
>load(filename)
>bigSamples[[i]]=m
> }
>
> But there may be a more efficient way to string up a big list like
> that, I can never remember - get it working, then worry about optimisation.
(Yes, use bigSamples <- vector("list", 15000) first.)
> I hope your filesystem is happy with 15000 objects in it. I would
> dedicate a folder or directory for just these objects' files, since it
> then becomes near impossible to see anything other than the big matrix
> files...
.readRDS/.saveRDS might be a better way to do this, and avoids always
restoring to "m".
If your file system does not like 15000 files you can always save in a
DBMS.
I did once look into restoring just some of the objects in a save()ed
file, but it is not really possible to do so efficiently due to sharing
between objects.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] problems with the fitted values of the function gls
Dear helpers, I am a beginner user of R and probably I am not using properly the function gls() of the library “nlme”. I estimated a model with AR(1) residuals but apparently the fitted and the predicted values aren’t constructed considering the autocorrelation of the residuals. It seems that instead of having: yt hat = (b hat)’ xt + Phi (ut-1 hat) the program gives just: yt hat = (b hat)’ xt To understand better this issue I tried to fit a model for a simulated AR(1) process with the following code: library(tseries) library(nlme) set.seed(2) ar1sim <- arima.sim ( list (order = c (1,0,0), ar = 0.5), n = 200 ) reg <- gls ( ar1sim ~ 1, correlation = corAR1() ) plot(ar1sim) points(predict(reg), type = “l”, col = “red”) - And I obtained the plot of the fitted values (just the intercept) compared with the actual ones. Probably I am doing something wrong, because I don’t think that the function should behave in this way, but I don’t understand what. Thank you so much for your help, Carlo Fezzi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] some suggestion
take a look at this function by Kevin Wright
RSiteSearch("sort.data.frame")
Best,
Dimitris
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
- Original Message -
From: "ronggui" <[EMAIL PROTECTED]>
To:
Sent: Tuesday, June 14, 2005 4:28 PM
Subject: [R] some suggestion
it seems R has no function to sort the data.frame according to some
variable(s),though we can do these by order() and indexing.but why
not make sort() the a generic function,making it can sort vector and
other type of objects?
maybe this is a silly suggestion,but i think it is quite usefull.
2005-06-14
--
Deparment of Sociology
Fudan University
Blog:www.sociology.yculblog.com
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
On Tue, 14 Jun 2005, Frank E Harrell Jr wrote: > Prof Brian Ripley wrote: >> On Tue, 14 Jun 2005, James Salsman wrote: >> >> >>> How can ordinary polynomial coefficients be calculated >>> from an orthogonal polynomial fit? >> >> >> Why would you want to do that? predict() is perfectly happy with an >> orthogonal polynomial fit and the `ordinary polynomial coefficients' are >> rather badly determined in your example since the design matrix has a very >> high condition number. > > Brian - I don't fully see the relevance of the high condition number nowadays > unless the predictor has a really bad origin. Orthogonal polynomials are a > mess for most people to deal with. It means that if you write down the coeffs to a few places and then try to reproduce the predictions you will do badly. The perturbation analysis depends on the condition number, and so is saying that the predictions are dependent on fine details of the coefficients. Using (year-2000)/1000 or (year - 1970)/1000 would be a much better idea. Why do `people' need `to deal with' these, anyway. We have machines to do that. > > Frank > >> >> >>> I'm trying to do something like find a,b,c,d from >>> lm(billions ~ a+b*decade+c*decade^2+d*decade^3) >>> but that gives: "Error in eval(expr, envir, enclos) : >>> Object "a" not found" >> >> >> You could use >> >> lm(billions ~ decade + I(decade^2) + I(decade^3)) >> >> except that will be numerically inaccurate, since >> >> >>> m <- model.matrix(~ decade + I(decade^2) + I(decade^3)) >>> kappa(m) >> >> [1] 3.506454e+16 >> >> >> >> decade <- c(1950, 1960, 1970, 1980, 1990) billions <- c(3.5, 5, 7.5, 13, 40) # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg pm <- lm(billions ~ poly(decade, 3)) plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), >>> >>> main="average yearly inflation-adjusted dollar cost of extreme weather >>> events worldwide") >>> curve(predict(pm, data.frame(decade=x)), add=TRUE) # output: http://www.bovik.org/storms.gif summary(pm) >>> >>> Call: >>> lm(formula = billions ~ poly(decade, 3)) >>> >>> Residuals: >>> 1 2 3 4 5 >>> 0.2357 -0.9429 1.4143 -0.9429 0.2357 >>> >>> Coefficients: >>> Estimate Std. Error t value Pr(>|t|) >>> (Intercept)13.800 0.882 15.647 0.0406 * >>> poly(decade, 3)1 25.614 1.972 12.988 0.0489 * >>> poly(decade, 3)2 14.432 1.972 7.318 0.0865 . >>> poly(decade, 3)36.483 1.972 3.287 0.1880 >>> --- >>> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 >>> >>> Residual standard error: 1.972 on 1 degrees of freedom >>> Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 >>> F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 >>> >>> pm >>> >>> Call: >>> lm(formula = billions ~ poly(decade, 3)) >>> >>> Coefficients: >>> (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 >>> 13.80025.61414.432 6.483 >>> >>> __ >>> R-help@stat.math.ethz.ch mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide! >>> http://www.R-project.org/posting-guide.html >>> >> >> > > > -- > Frank E Harrell Jr Professor and Chair School of Medicine > Department of Biostatistics Vanderbilt University > > -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] some suggestion
it seems R has no function to sort the data.frame according to some variable(s),though we can do these by order() and indexing.but why not make sort() the a generic function,making it can sort vector and other type of objects? maybe this is a silly suggestion,but i think it is quite usefull. 2005-06-14 -- Deparment of Sociology Fudan University Blog:www.sociology.yculblog.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Manipulating dates
try this:
x <- c("28/03/2000", "27/08/2001", "29/05/2002", "15/12/2003")
y <- as.Date(x, "%d/%m/%Y")
y - min(y)
I hope it helps.
Best,
Dimitris
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
- Original Message -
From: "Richard Hillary" <[EMAIL PROTECTED]>
To:
Sent: Tuesday, June 14, 2005 4:07 PM
Subject: [R] Manipulating dates
> Hello,
> Given a vector of characters, or factors, denoting the date
> in
> the following way: 28/03/2000, is there a method of
> 1) Computing the earliest of these dates;
> 2) Using this as a base, then converting all the other dates into
> merely
> the number of days after this minimum date
> Many thanks
> Richard Hillary
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Manipulating dates
Use POSIX. To convert: my.dates <- strptime(your.characters, format='%d/%m/%Y') once you have that, you can use 'min' to find the minimum. 'difftime' will give you the differences. Jim __ James Holtman"What is the problem you are trying to solve?" Executive Technical Consultant -- Convergys Labs [EMAIL PROTECTED] +1 (513) 723-2929 Richard Hillary <[EMAIL PROTECTED]To: r-help@stat.math.ethz.ch c.uk>cc: Sent by: Subject: [R] Manipulating dates [EMAIL PROTECTED] ath.ethz.ch 06/14/2005 10:07 Please respond to r.hillary Hello, Given a vector of characters, or factors, denoting the date in the following way: 28/03/2000, is there a method of 1) Computing the earliest of these dates; 2) Using this as a base, then converting all the other dates into merely the number of days after this minimum date Many thanks Richard Hillary __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Manipulating dates
Hello, Given a vector of characters, or factors, denoting the date in the following way: 28/03/2000, is there a method of 1) Computing the earliest of these dates; 2) Using this as a base, then converting all the other dates into merely the number of days after this minimum date Many thanks Richard Hillary __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
I think this is covered in the FAQ: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Models -roger James Salsman wrote: > How can ordinary polynomial coefficients be calculated > from an orthogonal polynomial fit? > > I'm trying to do something like find a,b,c,d from > lm(billions ~ a+b*decade+c*decade^2+d*decade^3) > but that gives: "Error in eval(expr, envir, enclos) : > Object "a" not found" > > > decade <- c(1950, 1960, 1970, 1980, 1990) > > billions <- c(3.5, 5, 7.5, 13, 40) > > # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg > > > > pm <- lm(billions ~ poly(decade, 3)) > > > > plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), > main="average yearly inflation-adjusted dollar cost of extreme weather > events worldwide") > > curve(predict(pm, data.frame(decade=x)), add=TRUE) > > # output: http://www.bovik.org/storms.gif > > > > summary(pm) > > Call: > lm(formula = billions ~ poly(decade, 3)) > > Residuals: >1 2 3 4 5 > 0.2357 -0.9429 1.4143 -0.9429 0.2357 > > Coefficients: > Estimate Std. Error t value Pr(>|t|) > (Intercept)13.800 0.882 15.647 0.0406 * > poly(decade, 3)1 25.614 1.972 12.988 0.0489 * > poly(decade, 3)2 14.432 1.972 7.318 0.0865 . > poly(decade, 3)36.483 1.972 3.287 0.1880 > --- > Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 > > Residual standard error: 1.972 on 1 degrees of freedom > Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 > F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 > > > pm > > Call: > lm(formula = billions ~ poly(decade, 3)) > > Coefficients: > (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 >13.80025.61414.432 6.483 > > __ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > -- Roger D. Peng http://www.biostat.jhsph.edu/~rpeng/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
Prof Brian Ripley wrote: > On Tue, 14 Jun 2005, James Salsman wrote: > > >>How can ordinary polynomial coefficients be calculated >>from an orthogonal polynomial fit? > > > Why would you want to do that? predict() is perfectly happy with an > orthogonal polynomial fit and the `ordinary polynomial coefficients' are > rather badly determined in your example since the design matrix has a very > high condition number. Brian - I don't fully see the relevance of the high condition number nowadays unless the predictor has a really bad origin. Orthogonal polynomials are a mess for most people to deal with. Frank > > >>I'm trying to do something like find a,b,c,d from >> lm(billions ~ a+b*decade+c*decade^2+d*decade^3) >>but that gives: "Error in eval(expr, envir, enclos) : >>Object "a" not found" > > > You could use > > lm(billions ~ decade + I(decade^2) + I(decade^3)) > > except that will be numerically inaccurate, since > > >>m <- model.matrix(~ decade + I(decade^2) + I(decade^3)) >>kappa(m) > > [1] 3.506454e+16 > > > > >>>decade <- c(1950, 1960, 1970, 1980, 1990) >>>billions <- c(3.5, 5, 7.5, 13, 40) >>># source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg >>> >>>pm <- lm(billions ~ poly(decade, 3)) >>> >>>plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), >> >>main="average yearly inflation-adjusted dollar cost of extreme weather >>events worldwide") >> >>>curve(predict(pm, data.frame(decade=x)), add=TRUE) >>># output: http://www.bovik.org/storms.gif >>> >>>summary(pm) >> >>Call: >>lm(formula = billions ~ poly(decade, 3)) >> >>Residuals: >> 1 2 3 4 5 >> 0.2357 -0.9429 1.4143 -0.9429 0.2357 >> >>Coefficients: >> Estimate Std. Error t value Pr(>|t|) >>(Intercept)13.800 0.882 15.647 0.0406 * >>poly(decade, 3)1 25.614 1.972 12.988 0.0489 * >>poly(decade, 3)2 14.432 1.972 7.318 0.0865 . >>poly(decade, 3)36.483 1.972 3.287 0.1880 >>--- >>Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 >> >>Residual standard error: 1.972 on 1 degrees of freedom >>Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 >>F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 >> >> >>>pm >> >>Call: >>lm(formula = billions ~ poly(decade, 3)) >> >>Coefficients: >> (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 >> 13.80025.61414.432 6.483 >> >>__ >>R-help@stat.math.ethz.ch mailing list >>https://stat.ethz.ch/mailman/listinfo/r-help >>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >> > > -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
I would suggest saving each object to an individual file with some sort of systematic file name. That way, you can implement a rudimentary key-value database and load only the objects you want. You might be interested in the 'serialize()' and 'unserialize()' functions for this purpose. If having ~15000 files is not desirable, then you need a database like GDBM. If you can live with something simpler, you might take a look at my 'filehash' package at http://sandybox.typepad.com/software/. It hasn't been tested much but it may suit your needs. -roger Richard Mott wrote: > Does anyone know a way to do the following: > > Save a large number of R objects to a file (like load() does) but then > read back only a small named subset of them . As far as I can see, > load() reads back everything. > > The context is: > > I have an application which will generate a large number of large > matrices (approx 15000 matrices each of dimension 2000*30). I can > generate these matrices using an R-package I wrote, but it requires a > large amouint of memory and is slow so I want to do this only once. > However, I then want to do some subsequent processing, comprising a very > large number of runs in which small (~ 10) random selection of matrices > from the previously computed set are used for linear modeling. So I > need a way to load back named objects previously saved in a call to > save(). I can;t see anyway of doing this. Any ideas? > > Thanks > > Richard Mott > > -- Roger D. Peng http://www.biostat.jhsph.edu/~rpeng/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
Richard Mott wrote:
> Does anyone know a way to do the following:
>
> Save a large number of R objects to a file (like load() does) but then
> read back only a small named subset of them . As far as I can see,
> load() reads back everything.
Save them to individual files when you generate them?
for(i in 1:15000){
m=generateBigMatrix(i)
filename=paste("BigMatrix-",i,".Rdata",sep='')
save(m,file=filename)
}
Note that load will always overwrite 'm', so to load a sample of them in
you'll need to do something like this:
bigSamples=list()
for(i in sample(15000,N)){
filename=paste("BigMatrix-",i,".Rdata",sep='')
load(filename)
bigSamples[[i]]=m
}
But there may be a more efficient way to string up a big list like
that, I can never remember - get it working, then worry about optimisation.
I hope your filesystem is happy with 15000 objects in it. I would
dedicate a folder or directory for just these objects' files, since it
then becomes near impossible to see anything other than the big matrix
files...
Baz
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] RES: Learning, "if" and "else"
Thank you very much! All methods are very useful.
Paulo M. Brando
Instituto de Pesquisa Ambiental da Amazonia (IPAM)
Santarem, PA, Brasil.
Av. Rui Barbosa, 136.
Fone: + 55 93 522 55 38
www.ipam.org.br
E-mail: [EMAIL PROTECTED]
-Mensagem original-
De: Adaikalavan Ramasamy [mailto:[EMAIL PROTECTED]
Enviada em: Tuesday, June 14, 2005 3:17 AM
Para: Paulo Brando
Cc: R-help@stat.math.ethz.ch
Assunto: Re: [R] Learning, "if" and "else"
You can also use switch() instead of ifelse(), which makes the code a
bit easier to read. The downside to this is that switch does not take
vectorised input and thus you need a loop. For example
# data
dbh <- c(30,29,15,14,30,29)
form <- factor(c("tree", "tree", "liana", "liana", "palm", "palm"))
df <- data.frame(dbh, form)
out <- numeric( nrow(df) )
for(i in 1:nrow(df)){
x <- as.numeric( df[i, 1] )
y <- as.character( df[i, 2] )
out[i] <- switch( y,
"tree" = exp( -4.898 + 4.512*log(x) - 0.319*(log(x))^2 ),
"liana" = 10^(0.07 + 2.17 * log10(x)),
NA
)
}
Or slightly more efficient solution is
out <- apply( df, 1, function(z){
x <- as.numeric(z[1]); y <- as.character(z[2])
switch( y,
"tree" = exp( -4.898 + 4.512*log(x) - 0.319*(log(x))^2 ),
"liana" = 10^(0.07 + 2.17 * log10(x)),
NA
)
})
Regards, Adai
On Mon, 2005-06-13 at 15:32 -0700, Paulo Brando wrote:
> Dear Rs,
>
> I have tried to use conditional expressions to calculate biomass for
> different life forms (trees, lianas, and palms).
>
> Here is an example:
>
> > lifeform
>
> dbh form
> 1 30 tree
> 2 29 tree
> 3 28 tree
> 4 27 tree
> 5 26 tree
> 6 25 tree
> 7 24 tree
> 8 23 tree
> 9 22 tree
> 10 21 tree
> 11 20 tree
> 12 15 liana
> 13 14 liana
> 14 13 liana
> 15 12 liana
> 16 11 liana
> 17 10 liana
> 18 9 liana
> 19 8 liana
> 20 7 liana
> 21 6 liana
> 22 5 liana
> 23 30 palm
> 24 29 palm
> 25 28 palm
> 26 27 palm
> 27 26 palm
> 28 25 palm
> 29 24 palm
> 30 23 palm
> 31 22 palm
> 32 21 palm
> 33 20 palm
>
> ### I want to include biomass
>
> lifeform$biomass <-
>
> {
> if(lifeform$form=="tree")
>exp(-4.898+4.512*log(dbh)-0.319*(log(dbh))^2)
>else{
> if (lifeform$form=="liana")
>10^(0.07 + 2.17 * log10 (dbh))
>else ("NA")
> }
> Warning message:
> the condition has length > 1 and only the first element will be used
in:
> if (lifeform$form == "tree") exp(-4.898 + 4.512 * log(dbh) -
>
>
> ### But I always get the message warning message above.
>
>
>
> I looked for similar examples in R mail list archive, but they did not
> help a lot.
>
> I am quite new to 'R'. Any material that covers this theme?
>
> Thank you very much!
>
> Paulo
>
> PS. Sorry about the last e-mail. I did not change the message title.
>
> Paulo M. Brando
> Instituto de Pesquisa Ambiental da Amazonia (IPAM)
> Santarem, PA, Brasil.
> Av. Rui Barbosa, 136.
> Fone: + 55 93 522 55 38
> www.ipam.org.br
> E-mail: [EMAIL PROTECTED]
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] load ing and saving R objects
This may not be quite the answer you're looking for, but I sometimes save each such object in its own file (usually .RData). Then, if you know which objects you're looking for, you know their names, and can load the individual files. Hope this helps, Matt Wiener -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Richard Mott Sent: Tuesday, June 14, 2005 9:03 AM To: r-help@stat.math.ethz.ch Subject: [R] load ing and saving R objects Does anyone know a way to do the following: Save a large number of R objects to a file (like load() does) but then read back only a small named subset of them . As far as I can see, load() reads back everything. The context is: I have an application which will generate a large number of large matrices (approx 15000 matrices each of dimension 2000*30). I can generate these matrices using an R-package I wrote, but it requires a large amouint of memory and is slow so I want to do this only once. However, I then want to do some subsequent processing, comprising a very large number of runs in which small (~ 10) random selection of matrices from the previously computed set are used for linear modeling. So I need a way to load back named objects previously saved in a call to save(). I can;t see anyway of doing this. Any ideas? Thanks Richard Mott -- Richard Mott | Wellcome Trust Centre tel 01865 287588 | for Human Genetics fax 01865 287697 | Roosevelt Drive, Oxford OX3 7BN __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Puzzled in utilising summary.lm() to obtain Var(x)
I have a program which is doing a few thousand runs of lm(). Suppose it is a simple model y = a + bx1 + cx2 + e I have the R object "d" where d <- summary(lm(y ~ x1 + x2)) I would like to obtain Var(x2) out of "d". How might I do it? I can, of course, always do sd(x2). But it would be much more convenient if I could snoop around the contents of summary.lm and extract Var() out of it. I couldn't readily see how. Would you know what would click? -- Ajay Shah Consultant [EMAIL PROTECTED] Department of Economic Affairs http://www.mayin.org/ajayshah Ministry of Finance, New Delhi __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] load ing and saving R objects
Does anyone know a way to do the following: Save a large number of R objects to a file (like load() does) but then read back only a small named subset of them . As far as I can see, load() reads back everything. The context is: I have an application which will generate a large number of large matrices (approx 15000 matrices each of dimension 2000*30). I can generate these matrices using an R-package I wrote, but it requires a large amouint of memory and is slow so I want to do this only once. However, I then want to do some subsequent processing, comprising a very large number of runs in which small (~ 10) random selection of matrices from the previously computed set are used for linear modeling. So I need a way to load back named objects previously saved in a call to save(). I can;t see anyway of doing this. Any ideas? Thanks Richard Mott -- Richard Mott | Wellcome Trust Centre tel 01865 287588 | for Human Genetics fax 01865 287697 | Roosevelt Drive, Oxford OX3 7BN __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] why does the unsubscribe not work
You can test this suggestion by asking for a password reminder from the https site. If you do get an email that says what your password is, then you are still subscribed. I actually asked the website to send me a password reminder to an account that is not subscribed. Although it says that "A reminder of your password has been emailed to you" on the website, the account in question did not receive any emails. Regards, Adai > What are the timeframes involved? There is always some time delay. In > particular, there is no way of stopping mails that have already been > sent by the list handling software. And if your mailbox is clogged, > things may stay in the queue for up to five days before the relevant > mail relay gives up. > > If your password has stopped working, it suggests to me that you have > in fact already been unsubscribed. > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] why does the unsubscribe not work
> "PD" Peter Dalgaard <[EMAIL PROTECTED]>
> on 14 Jun 2005 12:06:53 +0200 writes:
PD> [EMAIL PROTECTED] writes:
>> Hi I have tried to unsubscribe from this list now three
>> times, I desperately need to unsubscribe from this list
>> from this address because the list is choking this
>> mailbox.
>> Below are the various suggested ways of unsubscribing
>>
..
.. (whining whining ...)
..
PD> What are the timeframes involved? There is always some time delay. In
PD> particular, there is no way of stopping mails that have already been
PD> sent by the list handling software. And if your mailbox is clogged,
PD> things may stay in the queue for up to five days before the relevant
PD> mail relay gives up.
PD> If your password has stopped working, it suggests to me that you have
PD> in fact already been unsubscribed.
yes, and that's exactly true. He has been unsubscribed when he
wanted it -- and IIRC the confirmation e-mail even tells you
that you may continue to receive mail for a short while after
unsubscribing.
Yes, I ("r-help-owner") have seen your e-mail -- without any
politeness usually seen in letters including no name of sender --
among many other things that I must do or want to complete.
If "*-owner" was requested to deal with subscription and unsubscriptions
``manually'' each time someone unsubcscribes or subscribes,
this would amount to deal with about 20 such e-mails per day.
I'm usually inclined to reply to polite e-mail letters, but even
then only when more urgent duties have been fulfilled first.
It seems that 99% of all R-help subscribers can deal with the
system themselves.
Regards,
Martin Maechler, ETH Zurich
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] update.packages() - gregmisc
Thanks.
I do like this,
> remove.packages("gregmisc", .libPaths()[1])
> remove.packages("gtools", .libPaths()[1])
> install.packages("gregmisc", .libPaths()[1])
> update.packages()
> update.packages()
> install.packages("gtools", .libPaths()[1])
> update.packages()
> update.packages()
> update.packages(ask='graphics')
Regards,
Muhammad Subianto
R.2.1.0 on W2K
On this day 6/14/2005 1:51 PM, Gabor Grothendieck wrote:
> On 6/14/05, Muhammad Subianto <[EMAIL PROTECTED]> wrote:
>
>>Dear all,
>>I have a problem to update package gregmisc.
>>After I update,
>> > update.packages(ask='graphics')
>>trying URL
>>'http://cran.at.r-project.org/bin/windows/contrib/2.1/gregmisc_2.0.8.zip'
>>Content type 'application/zip' length 2465 bytes
>>opened URL
>>downloaded 2465 bytes
>>
>>package 'gregmisc' successfully unpacked and MD5 sums checked
>>...
>>
>>then try to update again, still I must update package gregmisc, etc.
>>I have tried 3,4,5, times with the same result.
>>
>
>
> This was discussed on r-devel recently. See:
>
> https://www.stat.math.ethz.ch/pipermail/r-devel/2005-June/033479.html
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bs() function of the splines package with intercept=FALSE
Hello, I'm implementing a function using non uniform B-Splines. Looking at the code of the bs() function, I realized that if the intercept was set to FALSE, the behavior of the function was the following (df is the number of degrees of freedom that I believe can be interpreted as the number of control points): - Compute df- ord + 1 _internal_ knots (ord is the order of the spline) - Add ord times the boundaries values at each extrem of the knots vector. - Compute the design matrix on this vector. This results in a design matrix with df+1 columns. The bs function the _removes_ the first column of the matrix (resulting as expected in a matrix with df columns). I'm a bit confused, does anyone see the mathematical justification of this manipulation? Thanks a lot, Laurent http://www.imagine-msn.com/hotmail/default.aspx?locale=fr-FR __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] update.packages() - gregmisc
On 6/14/05, Muhammad Subianto <[EMAIL PROTECTED]> wrote: > Dear all, > I have a problem to update package gregmisc. > After I update, > > update.packages(ask='graphics') > trying URL > 'http://cran.at.r-project.org/bin/windows/contrib/2.1/gregmisc_2.0.8.zip' > Content type 'application/zip' length 2465 bytes > opened URL > downloaded 2465 bytes > > package 'gregmisc' successfully unpacked and MD5 sums checked > ... > > then try to update again, still I must update package gregmisc, etc. > I have tried 3,4,5, times with the same result. > This was discussed on r-devel recently. See: https://www.stat.math.ethz.ch/pipermail/r-devel/2005-June/033479.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] update.packages() - gregmisc
Dear all, I have a problem to update package gregmisc. After I update, > update.packages(ask='graphics') trying URL 'http://cran.at.r-project.org/bin/windows/contrib/2.1/gregmisc_2.0.8.zip' Content type 'application/zip' length 2465 bytes opened URL downloaded 2465 bytes package 'gregmisc' successfully unpacked and MD5 sums checked ... then try to update again, still I must update package gregmisc, etc. I have tried 3,4,5, times with the same result. Best, Muhammad Subianto > version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor1.0 year 2005 month04 day 18 language R > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superscript in figures - basic question
On Tue, 14 Jun 2005, Martin Maechler wrote:
>> "Peter" == Peter Alspach <[EMAIL PROTECTED]>
>> on Tue, 14 Jun 2005 14:11:47 +1200 writes:
>
>Peter> Ben
>
>Peter> Others have pointed out plotmath. However, for some
>Peter> superscripts (including 2) it may be easier to use
>Peter> the appropriate escape sequence (at in Windows):
>
>Peter> ylab = 'BA (m\262/ha)'
>
> but please refrain from doing that way.
> You should write R code that is portable, and ``plotmath''
> has been designed to be portable. Windows escape codes are not,
> and may not even work in future (or current?) versions of
> Windows with `unusual' locale settings {fonts, etc}.
Indeed, it only works for superscript 2 and 3 and only in Western European
locales (I just happened to be testing Korean). Also, the spacing is not
done as carefully as plotmath.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
On Tue, 14 Jun 2005, James Salsman wrote: > How can ordinary polynomial coefficients be calculated > from an orthogonal polynomial fit? Why would you want to do that? predict() is perfectly happy with an orthogonal polynomial fit and the `ordinary polynomial coefficients' are rather badly determined in your example since the design matrix has a very high condition number. > I'm trying to do something like find a,b,c,d from > lm(billions ~ a+b*decade+c*decade^2+d*decade^3) > but that gives: "Error in eval(expr, envir, enclos) : > Object "a" not found" You could use lm(billions ~ decade + I(decade^2) + I(decade^3)) except that will be numerically inaccurate, since > m <- model.matrix(~ decade + I(decade^2) + I(decade^3)) > kappa(m) [1] 3.506454e+16 > > decade <- c(1950, 1960, 1970, 1980, 1990) > > billions <- c(3.5, 5, 7.5, 13, 40) > > # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg > > > > pm <- lm(billions ~ poly(decade, 3)) > > > > plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), > main="average yearly inflation-adjusted dollar cost of extreme weather > events worldwide") > > curve(predict(pm, data.frame(decade=x)), add=TRUE) > > # output: http://www.bovik.org/storms.gif > > > > summary(pm) > > Call: > lm(formula = billions ~ poly(decade, 3)) > > Residuals: > 1 2 3 4 5 > 0.2357 -0.9429 1.4143 -0.9429 0.2357 > > Coefficients: > Estimate Std. Error t value Pr(>|t|) > (Intercept)13.800 0.882 15.647 0.0406 * > poly(decade, 3)1 25.614 1.972 12.988 0.0489 * > poly(decade, 3)2 14.432 1.972 7.318 0.0865 . > poly(decade, 3)36.483 1.972 3.287 0.1880 > --- > Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 > > Residual standard error: 1.972 on 1 degrees of freedom > Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 > F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 > > > pm > > Call: > lm(formula = billions ~ poly(decade, 3)) > > Coefficients: > (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 > 13.80025.61414.432 6.483 > > __ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Learning, "if" and "else"
You can also use switch() instead of ifelse(), which makes the code a
bit easier to read. The downside to this is that switch does not take
vectorised input and thus you need a loop. For example
# data
dbh <- c(30,29,15,14,30,29)
form <- factor(c("tree", "tree", "liana", "liana", "palm", "palm"))
df <- data.frame(dbh, form)
out <- numeric( nrow(df) )
for(i in 1:nrow(df)){
x <- as.numeric( df[i, 1] )
y <- as.character( df[i, 2] )
out[i] <- switch( y,
"tree" = exp( -4.898 + 4.512*log(x) - 0.319*(log(x))^2 ),
"liana" = 10^(0.07 + 2.17 * log10(x)),
NA
)
}
Or slightly more efficient solution is
out <- apply( df, 1, function(z){
x <- as.numeric(z[1]); y <- as.character(z[2])
switch( y,
"tree" = exp( -4.898 + 4.512*log(x) - 0.319*(log(x))^2 ),
"liana" = 10^(0.07 + 2.17 * log10(x)),
NA
)
})
Regards, Adai
On Mon, 2005-06-13 at 15:32 -0700, Paulo Brando wrote:
> Dear Rs,
>
> I have tried to use conditional expressions to calculate biomass for
> different life forms (trees, lianas, and palms).
>
> Here is an example:
>
> > lifeform
>
> dbh form
> 1 30 tree
> 2 29 tree
> 3 28 tree
> 4 27 tree
> 5 26 tree
> 6 25 tree
> 7 24 tree
> 8 23 tree
> 9 22 tree
> 10 21 tree
> 11 20 tree
> 12 15 liana
> 13 14 liana
> 14 13 liana
> 15 12 liana
> 16 11 liana
> 17 10 liana
> 18 9 liana
> 19 8 liana
> 20 7 liana
> 21 6 liana
> 22 5 liana
> 23 30 palm
> 24 29 palm
> 25 28 palm
> 26 27 palm
> 27 26 palm
> 28 25 palm
> 29 24 palm
> 30 23 palm
> 31 22 palm
> 32 21 palm
> 33 20 palm
>
> ### I want to include biomass
>
> lifeform$biomass <-
>
> {
> if(lifeform$form=="tree")
>exp(-4.898+4.512*log(dbh)-0.319*(log(dbh))^2)
>else{
> if (lifeform$form=="liana")
>10^(0.07 + 2.17 * log10 (dbh))
>else ("NA")
> }
> Warning message:
> the condition has length > 1 and only the first element will be used in:
> if (lifeform$form == "tree") exp(-4.898 + 4.512 * log(dbh) -
>
>
> ### But I always get the message warning message above.
>
>
>
> I looked for similar examples in R mail list archive, but they did not
> help a lot.
>
> I am quite new to 'R'. Any material that covers this theme?
>
> Thank you very much!
>
> Paulo
>
> PS. Sorry about the last e-mail. I did not change the message title.
>
> Paulo M. Brando
> Instituto de Pesquisa Ambiental da Amazonia (IPAM)
> Santarem, PA, Brasil.
> Av. Rui Barbosa, 136.
> Fone: + 55 93 522 55 38
> www.ipam.org.br
> E-mail: [EMAIL PROTECTED]
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bs() function of the splines package
Laurent 14 juin 12:02 afficher les options De : "Laurent" <[EMAIL PROTECTED]> - Rechercher les messages de cet auteur Date : Tue, 14 Jun 2005 03:02:37 -0700 Local : Mar 14 juin 2005 12:02 Objet : bs() function of the splines package Répondre | Répondre à l'auteur | Transférer | Imprimer | Message individuel | Afficher l'original | Retirer | Signaler un cas d'utilisation abusive Hello, I'm implementing a function using non uniform B-Splines. Looking at the code of the bs() function, I realized that if the intercept was set to TRUE, the behavior of the function was the following (df is the number of degrees of freedom that I believe can be interpreted as the number of control points): - Compute df- ord + 1 _internal_ knots (ord is the order of the spline) - Add ord times the boundaries values at each extrem of the knots vector. - Compute the design matrix on this vector. This results in a design matrix with df+1 columns. The bs function then _removes_ the first column of the matrix (resulting as expected in a matrix with df columns). I'm a bit confused, does anyone see a mathematical justification to this manipulation? In this case, is it licit tu use df- ord + 2 internal knots and remove the last columns too? Thanks a lot, Laurent __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] why does the unsubscribe not work
[EMAIL PROTECTED] writes: > Hi I have tried to unsubscribe from this list now three times, I desperately > need to unsubscribe from this list from this address because the list is > choking > this mailbox. > Below are the various suggested ways of unsubscribing > > > To subscribe or unsubscribe via the World Wide Web, visit > > https://stat.ethz.ch/mailman/listinfo/r-help > > or, via email, send a message with subject or body 'help' to > > [EMAIL PROTECTED] > > I did the second the first time, I received an email suggesting that I send an > email to that address with unsubscribe in the subject nothing in the body, > which > I did, I received a confirmation request from that address and I confirmed, > after confirmation I received a response that I had been unsubscribed. I > continued to receive email from this list. > So I went ahead and tried to unsubscribe via the https address, unfortunately > it > said that it did not recognize my password, as I write down all my passwords I > doubt I've forgotten it and I am assuming that somewhere in the system I have > been unsubscribed but unfortunately not through the whole system. Anyway I > attempted again via the r-help-request, got the same response and of course > the > end result is I am still receiving the emails that are choking my system. > Yesterday in order to deal with this problem I tried to contact > [EMAIL PROTECTED] > and I have not received a response from him. This might very well be because > this mailbox is getting so damned choken by mail from this list that a number > of > emails from other lists etc. are being refused. I am cc'ing him on this > email. I > would like someone to get me unsubscribed from this list as I don't think it > very fair that I should have to abandon a mailbox that I have set to receiving > mails from several interesting lists and that I have been using for more than > a > year now. What are the timeframes involved? There is always some time delay. In particular, there is no way of stopping mails that have already been sent by the list handling software. And if your mailbox is clogged, things may stay in the queue for up to five days before the relevant mail relay gives up. If your password has stopped working, it suggests to me that you have in fact already been unsubscribed. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] ordinary polynomial coefficients from orthogonal polynomials?
How can ordinary polynomial coefficients be calculated from an orthogonal polynomial fit? I'm trying to do something like find a,b,c,d from lm(billions ~ a+b*decade+c*decade^2+d*decade^3) but that gives: "Error in eval(expr, envir, enclos) : Object "a" not found" > decade <- c(1950, 1960, 1970, 1980, 1990) > billions <- c(3.5, 5, 7.5, 13, 40) > # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg > > pm <- lm(billions ~ poly(decade, 3)) > > plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000), main="average yearly inflation-adjusted dollar cost of extreme weather events worldwide") > curve(predict(pm, data.frame(decade=x)), add=TRUE) > # output: http://www.bovik.org/storms.gif > > summary(pm) Call: lm(formula = billions ~ poly(decade, 3)) Residuals: 1 2 3 4 5 0.2357 -0.9429 1.4143 -0.9429 0.2357 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept)13.800 0.882 15.647 0.0406 * poly(decade, 3)1 25.614 1.972 12.988 0.0489 * poly(decade, 3)2 14.432 1.972 7.318 0.0865 . poly(decade, 3)36.483 1.972 3.287 0.1880 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 1.972 on 1 degrees of freedom Multiple R-Squared: 0.9957, Adjusted R-squared: 0.9829 F-statistic: 77.68 on 3 and 1 DF, p-value: 0.08317 > pm Call: lm(formula = billions ~ poly(decade, 3)) Coefficients: (Intercept) poly(decade, 3)1 poly(decade, 3)2 poly(decade, 3)3 13.80025.61414.432 6.483 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] why does the unsubscribe not work
Hi I have tried to unsubscribe from this list now three times, I desperately need to unsubscribe from this list from this address because the list is choking this mailbox. Below are the various suggested ways of unsubscribing > To subscribe or unsubscribe via the World Wide Web, visit > https://stat.ethz.ch/mailman/listinfo/r-help > or, via email, send a message with subject or body 'help' to > [EMAIL PROTECTED] I did the second the first time, I received an email suggesting that I send an email to that address with unsubscribe in the subject nothing in the body, which I did, I received a confirmation request from that address and I confirmed, after confirmation I received a response that I had been unsubscribed. I continued to receive email from this list. So I went ahead and tried to unsubscribe via the https address, unfortunately it said that it did not recognize my password, as I write down all my passwords I doubt I've forgotten it and I am assuming that somewhere in the system I have been unsubscribed but unfortunately not through the whole system. Anyway I attempted again via the r-help-request, got the same response and of course the end result is I am still receiving the emails that are choking my system. Yesterday in order to deal with this problem I tried to contact [EMAIL PROTECTED] and I have not received a response from him. This might very well be because this mailbox is getting so damned choken by mail from this list that a number of emails from other lists etc. are being refused. I am cc'ing him on this email. I would like someone to get me unsubscribed from this list as I don't think it very fair that I should have to abandon a mailbox that I have set to receiving mails from several interesting lists and that I have been using for more than a year now. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] protection stack overflow??
On Tue, 14 Jun 2005, Hu Chen wrote:
> Hi dear Rers,
> I am using SSOAP package to access SOAP service at NCBI.
> I followed the example code in SSOAP but failed.
>
>> z <-
>> .SOAP("http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/soap_adapter.cgi";,
>> method="run_eInfo", db="pubmed", action = I("einfo"))
> Error: protect(): protection stack overflow
>
> what's wrong?
Your code is overflowing the protection stack. Most likely due to heavy
recursion: see ?Startup as to how to increase the stack size.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] New Family object for GLM models...
Dear R-Users,
I wish to create a new family object based on the Binomial family. The only
difference will be with the link function. Thus instead if using the 'logit(u)'
link function, i plan to use '-log(i-u)'.
So far, i have tried to write the function following that of the Binomial and
Negative Binomial families.
The major problem i have here is with the definition of the initial values. The
definitions i have attempted so far don't yield convergence when implemented in
a model. I still can't figure out how the initial values are defined. I don't
know if it is based on some property of the fitted model or the domain and /or
range of the link function. I wish someone could help me provide a solution to
this problem. I have appended the function i wrote.
Add.haz<-function () {
env <- new.env(parent = .GlobalEnv)
assign(".nziek", nziek, envir = env)
famname<-"Addhaza"
link="addlink"
linkfun<-function(mu) -log(1-mu)
linkinv<-function(eta) 1-exp(-eta)
variance<-function(mu) mu*(1-mu)
validmu<-function(mu) all(mu > 0) && all(mu < 1)
mu.eta<-function(mu) 1/(1-mu)
dev.resids <- function(y, mu, wt) 2 * wt * (y * log(ifelse(y ==
0, 1, y/mu)) + (1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 -
mu
aic <- function(y, n, mu, wt, dev) {
m <- if (any(n > 1))
n
else wt
-2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
y), round(m), mu, log = TRUE))
}
initialize<- expression({
n<-rep(1,nobs)
if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
mustart <- (weights * y + 0.5)/(weights + 1)
m <- weights * y
if (any(abs(m - round(m)) > 0.001)) warning("non-integer #successes
in a binomial glm!")
})
environment(variance) <- environment(validmu) <- environment(dev.resids) <-
environment(aic) <- env
structure(list(family = famname, link = link, linkfun = linkfun,
linkinv = linkinv, variance = variance, dev.resids = dev.resids,
aic = aic, mu.eta =mu.eta, initialize = initialize,
validmu = validmu), class = "family")
}
Thank you for your kind attention.
Emmanuel
EMMANUEL NJI ABATIH
Rues des deux Eglises 140
1210 St Josse-Ten-Node, Bruxelles, BELGIUM
MOBLIE: 0032486958988(anytime)
Fix: 0032 2642 5038(8am to 6pm)
http://www.iph.fgov.be/epidemio/epien
-
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superscript in figures - basic question
Chuck Cleland <[EMAIL PROTECTED]> writes:
> See ?plotmath.
>
> plot(1:10, 1:10, ylab=expression(BA (m^{2}/ha)))
Ick... The parser will think that BA is a function and m^2/ha its
argument. This probably has consequences for the spacing.
Try expression("BA" * ("m"^2/"ha"))
--
O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superscript in figures - basic question
> "Peter" == Peter Alspach <[EMAIL PROTECTED]>
> on Tue, 14 Jun 2005 14:11:47 +1200 writes:
Peter> Ben
Peter> Others have pointed out plotmath. However, for some
Peter> superscripts (including 2) it may be easier to use
Peter> the appropriate escape sequence (at in Windows):
Peter> ylab = 'BA (m\262/ha)'
but please refrain from doing that way.
You should write R code that is portable, and ``plotmath''
has been designed to be portable. Windows escape codes are not,
and may not even work in future (or current?) versions of
Windows with `unusual' locale settings {fonts, etc}.
Martin Maechler
Peter> Cheers
Peter> Peter Alspach
"Benjamin M. Osborne" <[EMAIL PROTECTED]>
14/06/05 13:42:53 >>>
Peter> Although I see similar, but more complex, questions
Peter> addressed in the help archive, I'm having trouble
Peter> adding superscripted text to the y-axis labels of
Peter> some figures, and I can't find anything in the R
Peter> documentation on this.
Peter> I have:
Peter> ylab="BA (m2/ha)"
Peter> but I want the "2" to be superscripted. Thanks in
Peter> advence for the help, or for pointing out the
Peter> appropriate help file.
Peter> -Ben Osborne
Peter> __
Peter> R-help@stat.math.ethz.ch mailing list
Peter> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
Peter> do read the posting guide!
Peter> http://www.R-project.org/posting-guide.html
Peter> __
Peter> The contents of this e-mail are privileged and/or
Peter> confidenti...{{dropped}}
Peter> __
Peter> R-help@stat.math.ethz.ch mailing list
Peter> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
Peter> do read the posting guide!
Peter> http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] protection stack overflow??
Hi dear Rers,
I am using SSOAP package to access SOAP service at NCBI.
I followed the example code in SSOAP but failed.
> z <- .SOAP("http://www.ncbi.nlm.nih.gov/entrez/eutils/soap/soap_adapter.cgi";,
> method="run_eInfo", db="pubmed", action = I("einfo"))
Error: protect(): protection stack overflow
what's wrong?
Thanks very much.
Regards
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
