Re: [R] Fixed zeros in tables
Hi Andrew,
try the weights argument - apply zero weight to the structural zeros,
and 1 elsewhere.
Cheers
Andrew
On Sun, Nov 26, 2006 at 11:39:39AM +0700, A.R. Criswell wrote:
> Hello All R Users,
>
> Function loglm() in library MASS can be cajoled to accomodate
> structural zeros in a cross-classification table. An example from
> Fienberg demonstrates how this can be done.
>
> My question is: Can the function glm() perform the same task? Can
> glm() estimate a log-linear model with fixed zeros like loglm()?
>
> Thanks for your help,
> Andrew
>
> ## Fienberg, The Analysis of Cross-Classified Contingency Tables, 2nd
> ed., p.148.
> ## Results from survey of teenagers regarding their health concerns.
>
> health <- data.frame(expand.grid(CONCERNS = c("sex", "menstral",
> "healthy", "nothing"),
> AGE = c("12-15", "16-17"),
> GENDER = c("male", "female")),
> COUNT= c(4, 0, 42, 57, 2, 0, 7, 20,
> 9, 4, 19, 71, 7, 8, 10, 21))
>
> health <- xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = health)
>
> zeros <- data.frame(expand.grid(CONCERNS = c("sex", "menstral",
> "healthy", "nothing"),
> AGE = c("12-15", "16-17"),
> GENDER = c("male", "female")),
> COUNT= c(1, 0, 1, 1, 1, 0, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 1))
>
> zeros <- xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = zeros)
>
> library(MASS)
>
> fm.1 <- loglm(~ CONCERNS + AGE + GENDER,
> data = health, start = zeros, fitted = TRUE)
>
> fm.1; round(fm.1$fitted, 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
> and provide commented, minimal, self-contained, reproducible code.
--
Andrew Robinson
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
__
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
and provide commented, minimal, self-contained, reproducible code.
[R] plot p(Y=1) vs as
I am trying to fit a logistic regression model for this data set.
Firstly, I want to plot P(Y=1) vs As and P(Y=1) vs Aa.
Does any body know how to do these in R.
Thanks,
Aimin
> p5 <- read.csv("http://www.public.iastate.edu/~aiminy/data/p_5_2.csv";)
> str(p5)
'data.frame': 1030 obs. of 6 variables:
$ P : Factor w/ 5 levels "821p","8ABP",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Aa : Factor w/ 19 levels "ALA","ARG","ASN",..: 12 16 7 18 11 10 19 19
19 1 ...
$ As : num 126.9 64.1 82.7 7.6 42.0 ...
$ Ms : num 92.4 50.7 75.3 17.2 57.7 ...
$ Cur: num -0.1320 -0.0977 -0.0182 0.2368 0.1306 ...
$ Y : int 0 0 1 1 0 0 1 0 1 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
and provide commented, minimal, self-contained, reproducible code.
[R] Fixed zeros in tables
Hello All R Users,
Function loglm() in library MASS can be cajoled to accomodate
structural zeros in a cross-classification table. An example from
Fienberg demonstrates how this can be done.
My question is: Can the function glm() perform the same task? Can
glm() estimate a log-linear model with fixed zeros like loglm()?
Thanks for your help,
Andrew
## Fienberg, The Analysis of Cross-Classified Contingency Tables, 2nd
ed., p.148.
## Results from survey of teenagers regarding their health concerns.
health <- data.frame(expand.grid(CONCERNS = c("sex", "menstral",
"healthy", "nothing"),
AGE = c("12-15", "16-17"),
GENDER = c("male", "female")),
COUNT= c(4, 0, 42, 57, 2, 0, 7, 20,
9, 4, 19, 71, 7, 8, 10, 21))
health <- xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = health)
zeros <- data.frame(expand.grid(CONCERNS = c("sex", "menstral",
"healthy", "nothing"),
AGE = c("12-15", "16-17"),
GENDER = c("male", "female")),
COUNT= c(1, 0, 1, 1, 1, 0, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1))
zeros <- xtabs(COUNT ~ CONCERNS + AGE + GENDER, data = zeros)
library(MASS)
fm.1 <- loglm(~ CONCERNS + AGE + GENDER,
data = health, start = zeros, fitted = TRUE)
fm.1; round(fm.1$fitted, 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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem loading package Hmisc
Sorry about the question. I had installed the package locally and
hence the library command was not working as stated. library(Hmisc,
lib.loc="~/R/") did the trick.
On 11/25/06, Ritwik Sinha <[EMAIL PROTECTED]> wrote:
> Hi,
>
> I installed the package Hmisc with the command
> install.packages("Hmisc") without errors. When I try to load the
> library with command library(Hmisc) I get the error
>
> > library(Hmisc)
> Error in library(Hmisc) : there is no package called 'Hmisc'
>
>
> > version
>_
> platform i386-pc-linux-gnu
> arch i386
> os linux-gnu
> system i386, linux-gnu
> status
> major 2
> minor 4.0
> year 2006
> month 10
> day03
> svn rev39566
> language R
> version.string R version 2.4.0 (2006-10-03)
>
>
> Any help is appreciated. Unfortunately the RSiteSearch is not working
> now, so the question may not be well researched.
>
> Thanks
> Ritwik Sinha
> Graduate Student
> Epidemiology and Biostatistics
> Case Western Reserve University
> [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
>
--
Ritwik Sinha
Graduate Student
Epidemiology and Biostatistics
Case Western Reserve University
[EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
__
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
and provide commented, minimal, self-contained, reproducible code.
[R] problem loading package Hmisc
Hi,
I installed the package Hmisc with the command
install.packages("Hmisc") without errors. When I try to load the
library with command library(Hmisc) I get the error
> library(Hmisc)
Error in library(Hmisc) : there is no package called 'Hmisc'
> version
_
platform i386-pc-linux-gnu
arch i386
os linux-gnu
system i386, linux-gnu
status
major 2
minor 4.0
year 2006
month 10
day03
svn rev39566
language R
version.string R version 2.4.0 (2006-10-03)
Any help is appreciated. Unfortunately the RSiteSearch is not working
now, so the question may not be well researched.
Thanks
Ritwik Sinha
Graduate Student
Epidemiology and Biostatistics
Case Western Reserve University
[EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] My own correlation structure with nlme
I haven't seen a reply to this post, so I'll offer some comments,
even though I can't answer the question directly.
1. Given your question, I assume you have consulted Pinheiro and
Bates (2000) Mixed-Effects Models in S and S-Plus (Springer). If you
have not, I'm confident you will find material relevant to your question
there, especially in chapters 6-8.
2. Are you aware that the standard R installation includes a
subdirectory "~R\library\nlme\scripts", which contain copies of R
commands to create all the analyses in the book? In particular,
"ch06.R" and "ch08.R" might be particularly relevant to your question.
If you have not made local copies of these and walked through the code
line by line, I suggest you do so. I've learned a lot doing that.
3. Which versions of R and 'nlme' are you using? Some minor
enhancements to the help files were added a few months ago, and there
might be something helpful in the current examples that wasn't there
before.
4. What have you done to develop simple toy examples to help you
test your understanding of different aspects of the code? This
technique has helped me solve many problems.
5. Are you familiar with the use of the "methods" command to
trace logic flow? Consider for example the following:
> methods("corMatrix")
[1] corMatrix.corAR1* corMatrix.corARMA* corMatrix.corCAR1*
[4] corMatrix.corCompSymm* corMatrix.corIdent*corMatrix.corNatural*
[7] corMatrix.corSpatial* corMatrix.corStruct* corMatrix.corSymm*
[10] corMatrix.pdBlocked* corMatrix.pdCompSymm* corMatrix.pdDiag*
[13] corMatrix.pdIdent* corMatrix.pdMat* corMatrix.reStruct*
Non-visible functions are asterisked
6. Are you familiar with using "getAnywhere" to obtain code that
may be hidden with namespaces? For example, I just learned that
'corAR1' and 'corMatrix.corAR1' are two distinct functions. I found the
latter with this "methods" command, and I got the code to it using
"getAnywhere".
7. Are you familiar with the 'debug' command (in the 'base'
package, not the 'debug' package, which is probably much more powerful
but I haven't used the latter)? This allows a user to trace through
code line by line, examining the contents of objects -- and even
modifying them if you want. This is an extremely powerful way to learn
what a piece of code does.
8. If the above does not produce the answers you seek with a
reasonable additional effort, please submit another post. To increase
your chances of getting replies that are quicker and more helpful than
this one, please include commented, minimal, self-contained,
reproducible code. I'm confident that I could have helped you more with
less effort if your 'pairCorr' code had been included a sample call that
I could copy from your email into R and see if I get the same error
message you got. If I failed to get the same error as you saw, that
suggests a problem in your installation. If I got the same error
message, there is a good chance that I figure out how to get around that
and provide more helpful suggestions in less time than I've been
thinking about this.
Hope this helps.
Spencer Graves
Mohand wrote:
> Dear all,
>
> I am trying to define my own corStruct which is different from the
> classical one available in nlme. The structure of this correlation is
> given below.
> I am wondering to know how to continue with this structure by using
> specific functions (corMatrix, getCovariate, Initialize,...) in order to
> get a structure like corAR1, corSymm which will be working for my data.
>
> Thanks in advance.
>
> Regards
>
> M. Feddag
>
>
>
>
> *Correlation structure _
> _*
>
>
> pairCorr <- function(A, B, lowerTriangle = TRUE){
> n <- length(A)
> m <- length(B)
> if (n != m) stop("A and B must have same length")
> result <- matrix(0, n, n)
> same.A <- outer(A, A, "==")
> same.B <- outer(B, B, "==")
> A.matches.B <- outer(A, B, "==")
> B.matches.A <- outer(B, A, "==")
> result <- result + 0.5*same.A + 0.5*same.B -
> 0.5*A.matches.B - 0.5*B.matches.A
> rownames(result) <- colnames(result) <- paste(A, B, sep = "")
> if (lowerTriangle) result <- result[row(result) > col(result)]
> result
> }
>
> __
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] hausman Test
Hi, Stefan web.de> writes: > > Does anyone know how to do an Hausman test? > > I´ve estimate a modell (some alternatives) with clogit an wanted to test the > IIA test (Independence of Irrelevant > Alternatives) after estimating a multinomial logit model? My (privately) package includes Hausmann's IIA test. It is available from following URL: http://www.arumat.net/dc_0.1-3.zip (Windows Binary. If you need other binary or source format, email me) however you will easily do it: 1. Estimate 2 Logit models (full model, submodel), obtain coef (b.full, b.sub), covariance(vcov.full, vcov.sub), calculate their differences(db, dvcov) 2. calculate the statistic: t(db)%*%(dvcov)%*%(db) 3. df of this statistic is qr(dvcov)$rank 4. perform test by 1-pchisq(statistic, df) __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] cat not evaluated before file.choose
turn off the buffered write or use flush.console().
On 11/25/06, Kevin Middleton <[EMAIL PROTECTED]> wrote:
> As part of a larger function I have code similar to the reduced
> example below. The user is instructed to choose a file, which gets
> read using read.csv. In this example, I just have the name of the
> file print out.
>
> When I call this function with choosefile(), the file dialog box
> appears before the first cat line is printed. After I choose a file,
> both cat lines are printed.
>
> choosefile <- function (){
>cat("Choose the data file.\n")
>filename <- file.choose(new = FALSE)
>cat("You chose: ", filename, sep = "")
>}
>
> Is there a way to force the first cat line to print before the call
> to file.choose? I'm using R 2.4.0 Patched (2006-11-24 r39989) on OS
> X. Session info below.
>
> Any suggestions would be greatly appreciated.
>
> Kevin Middleton
>
> ---
>
> > sessionInfo()
> R version 2.4.0 Patched (2006-11-24 r39989)
> powerpc-apple-darwin8.8.0
>
> locale:
> en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] "stats" "graphics" "grDevices" "utils" "datasets"
> "methods" "base"
>
> __
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem you are trying to solve?
__
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
and provide commented, minimal, self-contained, reproducible code.
[R] cat not evaluated before file.choose
As part of a larger function I have code similar to the reduced
example below. The user is instructed to choose a file, which gets
read using read.csv. In this example, I just have the name of the
file print out.
When I call this function with choosefile(), the file dialog box
appears before the first cat line is printed. After I choose a file,
both cat lines are printed.
choosefile <- function (){
cat("Choose the data file.\n")
filename <- file.choose(new = FALSE)
cat("You chose: ", filename, sep = "")
}
Is there a way to force the first cat line to print before the call
to file.choose? I'm using R 2.4.0 Patched (2006-11-24 r39989) on OS
X. Session info below.
Any suggestions would be greatly appreciated.
Kevin Middleton
---
> sessionInfo()
R version 2.4.0 Patched (2006-11-24 r39989)
powerpc-apple-darwin8.8.0
locale:
en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] "stats" "graphics" "grDevices" "utils" "datasets"
"methods" "base"
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Fitting mixed-effects models with lme with fixed error term variances
On 11/22/06, Gregor Gorjanc <[EMAIL PROTECTED]> wrote:
> Hello!
>
> Douglas Bates wrote:
> > On 11/22/06, Gregor Gorjanc <[EMAIL PROTECTED]> wrote:
> >> Douglas Bates wrote:
> >> > On 11/21/06, Gregor Gorjanc <[EMAIL PROTECTED]> wrote:
> >> >> Douglas Bates stat.wisc.edu> writes:
> >> >> ...>
> >> >> > Can you be more specific about which parameters you want to fix and
> >> >> > which you want to vary in the optimization?
> >> >>
> >> >> It would be nice to have the ability to fix all variances i.e.
> >> >> variances of
> >> >> random effects.
> >> >
> ...
> >> > effects but allow the variance of a slope for the same grouping factor
> >> > to be estimated. Well, you could but only in the fortune("Yoda")
> >> > sense.
> >> >
> >>
> >> Yes, I agree here. Thank you for the detailed answer.
> >>
> >> > By the way, if you fix all the variances then what are you optimizing
> >> > over? The fixed effects? In that case the solution can be calculated
> >> > explicitly for a linear mixed model. The conditional estimates of the
> >> > fixed effects given the variance components are the solution to a
> >> > penalized linear least squares problem. (Yes, the solution can also
> >> > be expressed as a generalized linear least squares problem but there
> >> > are advantages to using the penalized least squares representation.
> >>
> >> Yup. It would really be great to be able to do that nicely in R, say use
> >> lmer() once and since this might take some time use estimates of
> >> variance components next time to get fixed and random effects. Is this
> >> possible with lmer or any related function - not in fortune("Yoda")
> >> sense ;)
> >
> > Not quite. There is now a capability in lmer to specify starting
> > estimates for the relative precision matrices which means that you can
> > use the estimates from one fit as starting estimates for another fit.
> > It looks like
> >
> > fm1 <- lmer(...)
> > fm2 <- lmer(y ~ x + (...), start = [EMAIL PROTECTED])
> >
> > I should say that in my experience this has not been as effective as I
> > had hoped it would be. What I see in the optimization is that the
> > first few iterations reduce the deviance quite quickly but the
> > majority of the time is spent refining the estimates near the optimum.
> > So, for example, if it took 40 iterations to converge from the rather
> > crude starting estimates calculated within lmer it might instead take
> > 35 iterations if you give it good starting estimates.
>
> Yes, I also have the same experience with other programs, say VCE[1].
> What I was hopping for was just solutions from linear mixed model i.e.
> either via GLS or PLS representations and no optimization if values for
> variance components are passed as input - there should be a way to
> distinguish that from just passing starting values..
The PLS representation in lmer is in terms of the relative
variance-covariance matrices of the random effects where "relative"
means relative to s^2. (In fact, the Omega slot contains the relative
precision matrices (inverses of the relative variance matrices) but
its the same idea. All variance components are expressed relative to
s^2.) If it is sufficient to fix these then you can easily do that.
Just follow the code in lmer until it gets to
.Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose)
LMEoptimize(mer) <- cv
return(new("lmer", mer,
and replace the first two lines by
.Call("mer_factor", mer, PACKAGE = "lme4")
A side-effect of performing the factorization is to evaluate the ML
and REML deviances and store the result in the deviance slot.
The "follow the code in lmer" part will get a bit tricky because of
the number of functions that are hidden in the lme4 namespace but I'm
sure you know how to get around that.
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] IIA tests
Hi Did you already know how to do an IIA test especial with Hausman? Thnak you __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple Conditional Tranformations
I have a program that is similar to your longer version, but I could
never get the syntax quite right. This will be a big help in
understanding how by works with functions.
Thanks,
Bob
-Original Message-
From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
Sent: Saturday, November 25, 2006 11:11 AM
To: Muenchen, Robert A (Bob)
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Multiple Conditional Tranformations
Here is a correction:
do.call(rbind, by(mydata, 1:nrow(mydata), function(x)
switch(as.character(x$gender),
m = transform(x, score1 = 3*q1+q2, score2 = 3.5*q1+q2),
f = transform(x, score1 = 2*q1+q2, score2 = 2.5*q1+q2),
transform(x, score1 = NA, score2 = NA))
))
On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> Here are some additional solutions. It appears that the SAS code is
performing
> the transformation row by row and for each row the code in your post
is
> specifying the transformation so if you want to do it that way we
> could use 'by'
> like this (where this time we have also added NA processing for the
gender):
>
>
> do.call(rbind, by(mydata, 1:nrow(mydata), function(x)
> switch(as.character(x$gender),
> m = transform(x, score1 = 3*q1+q2, score2 = 3.5*q1+q2),
> f = transform(x, score1 = 2*q1+q2, score2 = 2.5*q1+q2),
> NA)
> ))
>
> # or this somewhat longer version:
>
> do.call(rbind, by(mydata, 1:nrow(mydata), function(x) with(x, {
> if (is.na(gender)) {
> score1 <- score2 <- NA
> } else if (gender == "m") {
> score1 = 3 * q1 + q2
> score2 = 3.5 * q1 + q2
> } else if (gender == "f") {
> score1 = 2 * q1 + q2
> score2 = 2.5 * q1 + q2
> }
> cbind(x, score1, score2)
> })))
>
>
>
>
>
>
>
> On 11/25/06, Muehnchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > That's exactly what I'm looking for. Thanks so much for taking the
time
> > to do it that way.
> >
> > On the redundancy issue, I think SAS checks the "else if" condition
only
> > if the original "if" is false. The check for f when not m I put in
only
> > to exclude missing values for gender.
> >
> > Thanks!!
> > Bob
> >
> > -Original Message-
> > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> > Sent: Saturday, November 25, 2006 7:37 AM
> > To: Muenchen, Robert A (Bob)
> > Cc: r-help@stat.math.ethz.ch
> > Subject: Re: [R] Multiple Conditional Tranformations
> >
> > Firstly your outline does not check once, it checks twice. First it
> > check for "m" and then it redundantly checks for "f". On the other
> > hand the two variations in my post do check once.
> >
> > Although substantially longer than the solutions in my prior posts,
> > if you want the style shown in your post try this:
> >
> > mydata2 <- cbind(mydata, score1 = 0, score2 = 0)
> > is.m <- mydata$gender == "m"
> >
> > mydata2[is.m, ] <- transform(mydata[is.m, ],
> > score1 = 3 * q1 + q2,
> > score2 = 3.5 * q1 + q2
> > )
> >
> > mydata2[!is.m,] <- transform(mydata2[!is.m, ],
> > score1 = 2 * q1 + q2,
> > score2 = 2.5 * q1 + q2
> > )
> >
> > On 11/25/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > > Gabor,
> > >
> > > Those are handy variations! Perhaps my brain in still in SAS mode
on
> > > this. I'm expecting something like the code below that checks for
male
> > > only once, checks for female only when not male (skipping NAs) and
> > does
> > > all formulas under the appropriate conditions. The formulas I made
up
> > to
> > > keep the code short & may not be as easily modified to let the
logical
> > > 0/1 values fix them.
> > >
> > > if gender=="m" then do;
> > > Score1=...
> > > Score2=
> > > ...
> > > end;
> > > else if gender=="f" then do;
> > > Score1=...
> > > Score2=
> > > ...
> > > end;
> > >
> > > R may not have anything quite like that. R certainly has many
other
> > > features that SAS lacks.
> > >
> > > Thanks,
> > > Bob
> > >
> > > =
> > > Bob Muenchen (pronounced Min'-chen), Manager
> > > Statistical Consulting Center
> > > U of TN Office of Information Technology
> > > 200 Stokely Management Center, Knoxville, TN 37996-0520
> > > Voice: (865) 974-5230
> > > FAX: (865) 974-4810
> > > Email: [EMAIL PROTECTED]
> > > Web: http://oit.utk.edu/scc,
> > > News: http://listserv.utk.edu/archives/statnews.html
> > > =
> > >
> > >
> > > -Original Message-
> > > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> > > Sent: Saturday, November 25, 2006 12:39 AM
> > > To: Muenchen, Robert A (Bob)
> > > Cc: r-help@stat.math.ethz.ch
> > > Subject: Re: [R] Multiple Conditional Tranformations
> > >
> > > And here is a variation:
> > >
> > > transform(mydata,
> > > score1 = (2 + (gender == "m")) * q1 + q2,
> > > score2 = score1 + 0.5 * q1
> > > )
> > >
> > > or
> > >
> > > transform(
> > > transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2),
> > > score2 = score1 + 0.5 * q1
> > > )
> > >
> >
Re: [R] vector problem
thanks for your suggestions.
S Poetry looks really amazing, but also a lot of work for me.
I am pretty sure i´ll spent some hours over it.. (i should say weeks) :)
thanks for the zoo tips also..
but somehow i found out my problems are elsewhere.
i have to work with the following survey.
people were able to answer
in percentages that have to sum up to 100.
and answer is saved in one line of the answers.txt which i read in.
one answer from one participant can therefore content 1-3 rows.
e.g.:
participant1 says 10 and 20 and 70 and has therefore 3 rows
participant2 says 80 and nothing and 20 and has unfortunately 2 rows
what i wand to have is a matrix of all the answers with three
columns where every participant is a row.
should look like this:
10 20 70
80 0 20
I know data could have been saved better... but i do have this
dataset now and gotta rumble with it.
i will try again now, but i am still very glad about any help
because, this is my first real data apart from the practice stuff !
thanks in advance
Am 25.11.2006 um 11:24 schrieb Patrick Burns:
> Not an answer to your specific question, but S Poetry
> may be of use to you. In particular, it will say why it is
> not good practice to 'rbind' an object in a loop -- it is
> much better to create the matrix with its full size and then
> subscript into it.
>
> Patrick Burns
> [EMAIL PROTECTED]
> +44 (0)20 8525 0696
> http://www.burns-stat.com
> (home of S Poetry and "A Guide for the Unwilling S User")
>
> bunny , lautloscrew.com wrote:
>
>> Hello out there,
>>
>> i am not yet that experienced and trying to my best on a real
>> survey. but i am stuck with a little matrix / vector problem.
>> my vector of answers could have a length of 3 or only one. i want
>> to rbind all the answers into one matrix. (one vector for each
>> participant)
>> answers vectors for one participant could look like:
>>
>> p1: 100
>> p2: 20 80
>> p3: 40 10 50
>>
>> i have the following loop which should rbind them but i get no
>> proper matrix.
>> i´d like to have something like this
>>
>> 100 0 0
>> 20 80 0
>> 40 10 50
>>
>> Here´s my loop:
>>
>>
>> for(s in 1:length(fr))
>> {
>> ### ansmini is a complete survey of one participant
>> ansmini3=answers[relevant[s,],]
>>
>> for(x in 1:length(qidsb3))
>> {
>> # outputs all answer ids out of all data which belong to the
>> qids out of the question id vector
>> ansmin3=ansmini3[,1][ansmini3[,3]==qidsb3[x]]
>> newb3=rbind(newb3,ansmin3)
>> }
>> }
>>
>> this is doing the job for question with only one answer, so i
>> have only one answer id.
>> in this new case i have 1-3 possible answers.. that´s why i´am
>> not getting a nice newb3 matrix...
>>
>> i´d really be happy about some advice!
>> thanks in advance
>>
>> matthias
>>
>>
>>
>> [[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
>> and provide commented, minimal, self-contained, reproducible code.
>>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Nonlinear statistical modeling -- a comparison of R and AD Model Builder
On 11/24/06, dave fournier <[EMAIL PROTECTED]> wrote: > > >Dave > > Did you try supplying gradient information to nlminb? (I note that > nlminb is used for the optimization, but I don't see any gradient > information supplied to it.) I would suspect that supplying gradient > information would greatly speed up the computation (as you note in > comments at http://otter-rsch.ca/tresults.htm.) > Actually you should probably ask Norm Olsen these questions. > I am not proficient in R and am just using his code. Don't you find it somewhat disingenuous that you publish a comparison between the AD Model Builder software that you sell and R - a comparison that shows a tremendous advantage for your software - and then you write "I am not proficient in R"? Had you been proficient in R you might have known about the symbolic differentiation capabilities, specifically the deriv function, that have been part of the S language since the late 1980s. I believe that the 'AD' in "AD model builder" stands for automatic differentiation, which is actually something that John Chambers and I discussed at length when we were developing nonlinear modeling methods for S. In the end we went with symbolic differentiation rather than automatic differentiation because we felt that symbolic was more flexible. This is not to say that automatic differentiation isn't a perfectly legitimate technique. However, my recollection is that it would have required extensive revisions to the arithmetic expression evaluator, which is already very tricky code because of the "recycling rule" and the desire to shield users from knowledge of the internal representations and such details as whether you are using logical or integer or double precision operands or a combination. If you want to see these details you can, of course, examine the source code. I don't believe we would have the opportunity to examine how you implemented automatic differentiation. I must also agree with Spencer Graves that when I start reading a description of a nonlinear model with over 100 parameters, the example that you chose, I immediately start thinking of nonlinear mixed effects models. In my experience the only way in which a nonlinear model ends up with that number of parameters is through applying an underlying model with a low number of parameters to various groups within the data. Table 2 in the Schnute et al. paper to which you make reference states that the number of parameters in the model is T + A + 5 where T is the number of years of data and A is the number of age classes. To me that looks a lot like a nonlinear mixed effects model. Also, your choice of subject heading for your message seems deliberatively provocative. You seem to be implying that you are discussing a comparisons of AD Model Builder and R on all aspects of nonlinear statistical modeling but you are only discussing one comparison on simulated data using a model from the applications area for which you wrote AD Model Builder. Then you follow up by saying "I am not proficient in R" and your results for R are from applying code that someone else gave you. It seems that ADMB had a bit of a "home-field advantage" in this particular comparison. I view nonlinear statistical modeling differently. I have had a bit of experience in the area and I find that I want to examine the data carefully, usually through plots, before I embark on fitting complicated models. I like to have some assurance that the model makes sense in the context of the data. (In your example you don't need to worry about appropriateness of the model because the data were simulated.) I would never try to fit a nonlinear model with 100 parameters to data without carefully examining the data, and especially selected subsets of the data, first. For this the flexibility of the S language and tools like lattice graphics that were developed in this language are invaluable to me. The flexibility of data manipulation and graphics for interactive exploration of data is what attracted me to S in the first place. I realize that for many people the area of nonlinear statistical modeling is reduced to "Fit this model to these data and don't ask any questions. Just give me parameter estimates and p-values." If that is your situation then it would make sense to use software that gets you those estimates as quickly as possible with a minimum of effort. I'm just happy that I get to turn down people who ask me to do that. I like that fact that I can spend my time asking questions about the data and of the data. > However I can say that providing derivatives for such a model is a > highly nontrivial exercise. As I said in my posting, the R script and > data are available to anyone who feels that the exercise was not carried > out properly and would like to improve on it. Also one does not need > to provide derivatives to the AD Model Builder program. > > Finally suppose that you are very good at calculating derivatives
Re: [R] alphachannel on pdf-device
Thank you Paul Murrel for answering and thank you Brian Ripley for
correcting this behaviour in R-patched.
I am sorry, but I think there is just another 'misfeature' with
transparent figures since R-2.4.0 and it seems not to be corrected in
todays R-2.4.0 Patched (2006-11-24 r39989).
In earlier versions when producing a transparent figure without any
border (=NA), nothing was drawn. In newer versions there is something
like a white border line. You can see it in particular when drawing a
new dotted line above this border of the underlying figure. I should
demonstrate it with a small example:
pdf("test_alpha.pdf", version="1.4")
plot(1:8, 1:8, type="n", xlab="", ylab="")
lines(c(1,8), c(1,8), type="l", col="lightblue", lwd=15)
# rectangle with alphachannel=255 (no transparency)
goldenrod1 <- rgb(255,193,37,255, names=c("goldenrod1"), max=255)
rect(4,4,7,7, col=goldenrod1, border=NA)
lines(c(4,4,7,7,4), c(4,7,7,4,4), col="red", lty=2)
text(5,7, "no transparency")
# polygon with alphachannel=64 (transparency near 40%)
polygon(c(5,2,2,5,5), c(2,2,5,5,2), col=rgb(255,192,203,64, max = 255),
border=NA)
lines(c(5,2,2,5,5), c(2,2,5,5,2), col="red", lty=2)
text(4,2, "transparency 40%")
dev.off()
It is probably just a mistake. What do you think?
Rainer Hurling
Paul Murrell schrieb:
Hi
This appears to have been fixed in r-patched (by Brian Ripley)
Paul
Rainer Hurling wrote:
Dear R users,
since R-2.4.0 I am not able any more to draw figures without
transparency after drawing one figure with alphachannel set lower than 100%.
For better unterstanding here a small example of my problem:
---
pdf(paste("test_alpha_rgb.pdf",sep=""), version="1.4")
plot(1:8, 1:8, type="n", xlab="", ylab="")
lines(c(1,8), c(1,8), type="l", col="lightblue", lwd=15)
# polygon with alphachannel=64 (transparency 40%)
polygon(c(5,2,2,5,5), c(2,2,5,5,2), col=rgb(255,192,203,64, max = 255),
border=NA)
text(4,2, "transparency 40%")
# rectangle with alphachannel=255 (no transparency)
goldenrod1 <- rgb(255,193,37,255, names=c("goldenrod1"), max=255)
rect(4,4,7,7, col=goldenrod1, border=goldenrod1)
text(5,7, "no transparency")
dev.off()
---
In this example (see pdf in attachment) after plotting a line without
any transparency a polygon is drawn with alphachannel set to 40%. All
following figures and even text is shown with transparency.
The described behaviour seems to be the same on FreeBSD 7.0-CURRENT and
on WindowsXP.
Am I doing something wrong?
Many thanks for any help,
Rainer Hurling
pdf("test_alpha.pdf", version="1.4")
plot(1:8, 1:8, type="n", xlab="", ylab="")
lines(c(1,8), c(1,8), type="l", col="lightblue", lwd=15)
# rectangle with alphachannel=255 (no transparency)
goldenrod1 <- rgb(255,193,37,255, names=c("goldenrod1"), max=255)
rect(4,4,7,7, col=goldenrod1, border=NA)
lines(c(4,4,7,7,4), c(4,7,7,4,4), col="red", lty=2)
text(5,7, "no transparency")
# polygon with alphachannel=64 (transparency near 40%)
polygon(c(5,2,2,5,5), c(2,2,5,5,2), col=rgb(255,192,203,64, max = 255),
border=NA)
lines(c(5,2,2,5,5), c(2,2,5,5,2), col="red", lty=2)
text(4,2, "transparency 40%")
#for (x in 1:8)
# for (y in 1:8)
#rect(x-0.5,y-0.5, x+0.5, y+0.5, col=rgb(255,193,37, 255/(x*y), max=255))
#rect(5,5,7,7, col=rgb(255,0,0, 255, max=255), border=NA)
dev.off()
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Research Assistant Position
I posted this job opportunity on November 20 (see below). I would like to provide some background on the position for those who might be interested. It is in the University of California at the new Merced campus. This is the second academic year since the Merced campus opened, and the first in the all-new facilities. Merced is a 2-hour drive from the Bay Area, from Monterey, and from Yosemite National Park. Pay scales are competitive, and housing costs in Merced and neighboring communities are less than in most of California. There is considerable interest among the faculty in emphasizing the development, teaching and use of open source software as we build educational and research programs. Core courses for statistical data analysis for Engineering, Natural Sciences, and Cognitive and Social Sciences students using R are being taught or are in planning and development stages. This position is designed for a person experienced with programming and analysis using R. The advertised position includes duties supporting and participating in research and teaching development using R, including opportunities for publication or joint publication of research articles and R libraries. An example of relevant analyses and data visualization using R by the PI funding this position can be seen at http://www.sciencemag.org/cgi/rapidpdf/1128834.pdf . Anthony Westerling University of California Merced Merced, CA Programmer Analyst II/III (Research Assistant) Job Code SSNRI723A Open until filled. In the Sierra Nevada Research Institute at UC Merced, act in support of research in applied climatology and statistical modeling for wildfire, energy and water resource management applications and assist the Principle Investigator with the development of software, management of data sets and design, modification, and implementation of systems for modeling and analysis. Develop software libraries for the R statistical project for publication and use in the classroom environment. Requires experience and demonstrated expertise in programming and data visualization. Relevant programming experience includes R, Fortran and/or C. HTML is also desired; background in statistics, physics, climatology, hydrology, fire ecology or a similar field (Masters preferred); strong problem solving; demonstrated written communication and programming skills. A UC Merced job application, resume and cover letter are requested. For more information and to apply call 1-866-669-JOBS or visit http://jobs.ucmerced.edu/n/staff/position.jsf?positionId=723. EOE Anthony Westerling School of Engineering School of Social Sciences, Humanities, and Arts University of California, Merced http://tenaya.ucsd.edu/~westerli/westerling.html [EMAIL PROTECTED] (209) 228 4099 __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] low-variance warning in lmer
On 11/24/06, Ben Bolker <[EMAIL PROTECTED]> wrote:
> For block effects with small variance, lmer will sometimes
> estimate the variance as being very close to zero and issue
> a warning.
The "very close to zero" is merely a computational convenience. The
theoretical lower bound on the variance components is zero but
evaluating the profiled log-likelihood or log-restricted likelihood
when the variance-covariance matrix of the random effects is singular
would require a completely different algorithm. Evaluating the
gradient and the ECME update in this situation would be even more
complicated.
Instead of doing that I set the lower bounds to a value close to zero
but still large enough that the evaluation algorithm used elsewhere
works here too. The value I use for the relative variances is 5e-10.
The corresponding estimate of the variance component is 5e-10 * s^2.
So the algorithm hasn't really converged to 5e-10 * s^2. The estimate
of the variance component should be zero but I haven't allowed it to
go all the way to zero.
> I don't have a problem with this
I hope not. The maximum likelihood or REML estimates of variance
components can legitimately be zero.
> -- I've explored
> things a bit with some simulations (see below) and conclude that
> this is probably inevitable when trying to incorporate
> random effects with not very much data (the means and medians
> of estimates are plausibly close to the nominal values:
> the fraction of runs with warnings/near-zero estimates drops
> from about 50% when the between-block variance is 1% of
> the total (with 2 treatments, 12 blocks nested within treatment,
> 3 replicates per block), to 15% when between=30% of total,
> to near zero when between=50% of total.
>
> My question is: what should I suggest to students when this
> situation comes up?
If it is a variance component that is estimated as zero then just drop
that term from the model. Having the MLE or the REML estimate of a
variance component be zero is pretty strong evidence that it is not
significantly greater than zero (in the usual sense of failing to
reject the null hypothesis - it still could be the case that it is
greater than zero but we don't have strong enough evidence to reject
the case that it is zero). Hence we go with the simpler model that
eliminates this random effect.
The more difficult case is when the estimate of the
variance-covariance matrix of vector-valued random effects is
singular. Say you have a random effect for both the intercept and the
slope w.r.t. time for each subject. It is possible to converge to a
variance-covariance matrix for these random effects with a non-zero
variance for the intercept and a non-zero variance for the slope but
perfect correlation between these random effects. I'm not sure what
to suggest in that case as the reduced model is, as far as I can see,
not in the class of linear mixed models. I think I might just not
bother mentioning this in a classroom.
> Can anyone point me to appropriate references?
> (I haven't found anything relevant in Pinheiro and Bates, but
> I may not have looked in the right place ...)
>
> thanks
> Ben Bolker
>
>
> self-contained but unnecessarily complicated simulation
> code/demonstration:
> ---
> library(lme4)
> library(lattice)
>
> simfun <- function(reefeff,ntreat=2,nreef=12,
>nreefpertreat=3,
>t.eff=10,
>totvar=25,seed=NA) {
> if (!is.na(seed)) set.seed(seed)
> ntot = nreef*nreefpertreat
> npertreat=ntot/ntreat
> reef = gl(nreef,nreefpertreat)
> treat = gl(ntreat,npertreat)
> r.sd = sqrt(totvar*reefeff)
> e.sd = sqrt(totvar*(1-reefeff))
> y.det = ifelse(treat==1,0,t.eff)
> r.vals = rnorm(nreef,sd=r.sd)
> e.vals = rnorm(ntot,sd=e.sd)
> y <- y.det+r.vals[as.numeric(reef)]+e.vals
> data.frame(y,treat,reef)
> }
>
> getranvar <- function(x) {
> vc <- VarCorr(x)
> c(diag(vc[[1]]),attr(vc,"sc")^2)
> }
>
> estfun <- function(reefeff,...) {
> x <- simfun(reefeff=reefeff,...)
> ow <- options(warn=2)
> f1 <- try(lmer(y~treat+(1|reef),data=x))
> w <- (class(f1)=="try-error" && length(grep("effectively zero",f1))>0)
> options(ow)
> f2 <- lmer(y~treat+(1|reef),data=x)
> c(getranvar(f2),as.numeric(w))
> }
>
>
> rvec <- rep(c(0.01,0.05,0.1,0.15,0.2,0.3,0.5),each=100)
> X <- t(sapply(rvec,estfun))
> colnames(X) <- c("reefvar","resvar","warn")
> rfrac <- X[,"reefvar"]/(X[,"reefvar"]+X[,"resvar"])
> fracwarn <- tapply(X[,"warn"],rvec,mean)
> est.mean <- tapply(rfrac,rvec,mean)
>
>
> op <- par(mfrow=c(1,2))
> plot(rvec,rfrac,type="n",xlim=c(-0.02,0.55),axes=FALSE)
> axis(side=2)
> box()
> boxplot(rfrac~rvec,at=unique(rvec),add=TRUE,pars=list(boxwex=0.03),
> col="gray")
> points(jitter(rvec),rfrac,col=X[,"warn"]+1)
> lines(unique(rvec),fracwarn,col="blue",type="b",lwd=2)
> text(0.4,0.1,"fraction\nzero",col="blue")
> abline(a=0,b=1,lwd=2)
> points(unique(rvec),est.mean,cex=1.5,col=5)
> ##
> plot(rvec
Re: [R] trimming plotting area in a boxplot
Patrick Kuss wrote: > Hello all, > > I am trying to trim the plotting area in a boxplot, such that the space > between the plot margins (left and right) are of identical size as > between the boxes. > In the example below I want to get rid of the space outside of the > abline(). > > I appreciate any suggestions. > > Factor = LETTERS[1:5] > A = c(1:5); B = c(2:6); C = c(1:5); D = c(3:7); E = c(1:5); F = c(4:8) > df = data.frame(A,B,C,D,E,F) > boxplot(df) > abline(v=c(0.5,6.5)) par(xaxs="i") boxplot(df) Uwe Ligges > Cheers, Patrick > > __ > 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 > and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] OT: P(Z <= -1.46).
Duncan Murdoch <[EMAIL PROTECTED]> writes: > My copy of the CRC standard mathematical tables give 0.0721, without > citation. > > > Could two algorithms ``reasonably'' disagree in the 4th decimal > > place? > > One possible source for this error (if it is an error), would be someone > rounding to 5 places giving 0.07215, then rounding again to 4 places. > Is that reasonable? Wouldn't be surprised. I'm using an introductory textbook that has qnorm(.95) as alternatingly 1.64 and 1.65 in its tables, where the latter fairly clearly comes from rounding of 1.645 instead of 1.644854. -- 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple Conditional Tranformations
Here is a correction:
do.call(rbind, by(mydata, 1:nrow(mydata), function(x)
switch(as.character(x$gender),
m = transform(x, score1 = 3*q1+q2, score2 = 3.5*q1+q2),
f = transform(x, score1 = 2*q1+q2, score2 = 2.5*q1+q2),
transform(x, score1 = NA, score2 = NA))
))
On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> Here are some additional solutions. It appears that the SAS code is
> performing
> the transformation row by row and for each row the code in your post is
> specifying the transformation so if you want to do it that way we
> could use 'by'
> like this (where this time we have also added NA processing for the gender):
>
>
> do.call(rbind, by(mydata, 1:nrow(mydata), function(x)
> switch(as.character(x$gender),
> m = transform(x, score1 = 3*q1+q2, score2 = 3.5*q1+q2),
> f = transform(x, score1 = 2*q1+q2, score2 = 2.5*q1+q2),
> NA)
> ))
>
> # or this somewhat longer version:
>
> do.call(rbind, by(mydata, 1:nrow(mydata), function(x) with(x, {
> if (is.na(gender)) {
> score1 <- score2 <- NA
> } else if (gender == "m") {
> score1 = 3 * q1 + q2
> score2 = 3.5 * q1 + q2
> } else if (gender == "f") {
> score1 = 2 * q1 + q2
> score2 = 2.5 * q1 + q2
> }
> cbind(x, score1, score2)
> })))
>
>
>
>
>
>
>
> On 11/25/06, Muehnchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > That's exactly what I'm looking for. Thanks so much for taking the time
> > to do it that way.
> >
> > On the redundancy issue, I think SAS checks the "else if" condition only
> > if the original "if" is false. The check for f when not m I put in only
> > to exclude missing values for gender.
> >
> > Thanks!!
> > Bob
> >
> > -Original Message-
> > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> > Sent: Saturday, November 25, 2006 7:37 AM
> > To: Muenchen, Robert A (Bob)
> > Cc: r-help@stat.math.ethz.ch
> > Subject: Re: [R] Multiple Conditional Tranformations
> >
> > Firstly your outline does not check once, it checks twice. First it
> > check for "m" and then it redundantly checks for "f". On the other
> > hand the two variations in my post do check once.
> >
> > Although substantially longer than the solutions in my prior posts,
> > if you want the style shown in your post try this:
> >
> > mydata2 <- cbind(mydata, score1 = 0, score2 = 0)
> > is.m <- mydata$gender == "m"
> >
> > mydata2[is.m, ] <- transform(mydata[is.m, ],
> > score1 = 3 * q1 + q2,
> > score2 = 3.5 * q1 + q2
> > )
> >
> > mydata2[!is.m,] <- transform(mydata2[!is.m, ],
> > score1 = 2 * q1 + q2,
> > score2 = 2.5 * q1 + q2
> > )
> >
> > On 11/25/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > > Gabor,
> > >
> > > Those are handy variations! Perhaps my brain in still in SAS mode on
> > > this. I'm expecting something like the code below that checks for male
> > > only once, checks for female only when not male (skipping NAs) and
> > does
> > > all formulas under the appropriate conditions. The formulas I made up
> > to
> > > keep the code short & may not be as easily modified to let the logical
> > > 0/1 values fix them.
> > >
> > > if gender=="m" then do;
> > > Score1=...
> > > Score2=
> > > ...
> > > end;
> > > else if gender=="f" then do;
> > > Score1=...
> > > Score2=
> > > ...
> > > end;
> > >
> > > R may not have anything quite like that. R certainly has many other
> > > features that SAS lacks.
> > >
> > > Thanks,
> > > Bob
> > >
> > > =
> > > Bob Muenchen (pronounced Min'-chen), Manager
> > > Statistical Consulting Center
> > > U of TN Office of Information Technology
> > > 200 Stokely Management Center, Knoxville, TN 37996-0520
> > > Voice: (865) 974-5230
> > > FAX: (865) 974-4810
> > > Email: [EMAIL PROTECTED]
> > > Web: http://oit.utk.edu/scc,
> > > News: http://listserv.utk.edu/archives/statnews.html
> > > =
> > >
> > >
> > > -Original Message-
> > > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> > > Sent: Saturday, November 25, 2006 12:39 AM
> > > To: Muenchen, Robert A (Bob)
> > > Cc: r-help@stat.math.ethz.ch
> > > Subject: Re: [R] Multiple Conditional Tranformations
> > >
> > > And here is a variation:
> > >
> > > transform(mydata,
> > > score1 = (2 + (gender == "m")) * q1 + q2,
> > > score2 = score1 + 0.5 * q1
> > > )
> > >
> > > or
> > >
> > > transform(
> > > transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2),
> > > score2 = score1 + 0.5 * q1
> > > )
> > >
> > >
> > > On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> > > > Try this:
> > > >
> > > >
> > > > transform(mydata,
> > > > score1 = (2 + (gender == "m")) * q1 + q2,
> > > > score2 = (2.5 + (gender == "m")) * q1 + q2
> > > > )
> > > >
> > > >
> > > > On 11/24/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > > > > Mark,
> > > > >
> > > > > I finally got that approach to work by spreading
Re: [R] OT: P(Z <= -1.46).
Based on integration it appears that .0721 is correct. > integrate(function(x) exp(-x^2/2)/(2*pi)^.5, -Inf, -1.46) 0.07214504 with absolute error < 1.2e-07 On 11/25/06, [EMAIL PROTECTED] <[EMAIL PROTECTED]> wrote: > In checking over the solutions to some homework that I had assigned I > observed the fact that in R (version 2.4.0) pnorm(-1.46) gives > 0.07214504. The tables in the text book that I am using for the > course give the probability as 0.0722. > > Fascinated, I scanned through 5 or 6 other text books (amongst the > dozens of freebies from publishers that lurk on my shelf) and found > that some agree with R (giving P(Z <= -1.46) = 0.0721) and some agree > with the first text book, giving 0.0722. > > It is clearly of little-to-no practical import, but I'm curious as to > how such a discrepancy would arise in this era. Has anyone any > idea? Is there any possibility that the algorithm(s) used to > calculate this probability is/are not accurate to 4 decimal places? > > Could two algorithms ``reasonably'' disagree in the 4th decimal > place? >cheers, > >Rolf Turner >[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 > and provide commented, minimal, self-contained, reproducible code. > __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] OT: P(Z <= -1.46).
On 11/25/2006 10:21 AM, [EMAIL PROTECTED] wrote: > In checking over the solutions to some homework that I had assigned I > observed the fact that in R (version 2.4.0) pnorm(-1.46) gives > 0.07214504. The tables in the text book that I am using for the > course give the probability as 0.0722. > > Fascinated, I scanned through 5 or 6 other text books (amongst the > dozens of freebies from publishers that lurk on my shelf) and found > that some agree with R (giving P(Z <= -1.46) = 0.0721) and some agree > with the first text book, giving 0.0722. > > It is clearly of little-to-no practical import, but I'm curious as to > how such a discrepancy would arise in this era. Has anyone any > idea? Is there any possibility that the algorithm(s) used to > calculate this probability is/are not accurate to 4 decimal places? A text I've used gives the 0.0722 value, citing the 1962 edition of Lindgren's Statistical Theory. So it's not completely certain that this is "in this era". You can see parts of the 1993 version of Lindgren on books.google.com, and it repeats the 0.0722 value, but without citation (at least in the parts that are online). My copy of the CRC standard mathematical tables give 0.0721, without citation. > Could two algorithms ``reasonably'' disagree in the 4th decimal > place? One possible source for this error (if it is an error), would be someone rounding to 5 places giving 0.07215, then rounding again to 4 places. Is that reasonable? Duncan Murdoch > cheers, > > Rolf Turner > [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 > and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple Conditional Tranformations
Here are some additional solutions. It appears that the SAS code is performing
the transformation row by row and for each row the code in your post is
specifying the transformation so if you want to do it that way we
could use 'by'
like this (where this time we have also added NA processing for the gender):
do.call(rbind, by(mydata, 1:nrow(mydata), function(x)
switch(as.character(x$gender),
m = transform(x, score1 = 3*q1+q2, score2 = 3.5*q1+q2),
f = transform(x, score1 = 2*q1+q2, score2 = 2.5*q1+q2),
NA)
))
# or this somewhat longer version:
do.call(rbind, by(mydata, 1:nrow(mydata), function(x) with(x, {
if (is.na(gender)) {
score1 <- score2 <- NA
} else if (gender == "m") {
score1 = 3 * q1 + q2
score2 = 3.5 * q1 + q2
} else if (gender == "f") {
score1 = 2 * q1 + q2
score2 = 2.5 * q1 + q2
}
cbind(x, score1, score2)
})))
On 11/25/06, Muehnchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> That's exactly what I'm looking for. Thanks so much for taking the time
> to do it that way.
>
> On the redundancy issue, I think SAS checks the "else if" condition only
> if the original "if" is false. The check for f when not m I put in only
> to exclude missing values for gender.
>
> Thanks!!
> Bob
>
> -Original Message-
> From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> Sent: Saturday, November 25, 2006 7:37 AM
> To: Muenchen, Robert A (Bob)
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] Multiple Conditional Tranformations
>
> Firstly your outline does not check once, it checks twice. First it
> check for "m" and then it redundantly checks for "f". On the other
> hand the two variations in my post do check once.
>
> Although substantially longer than the solutions in my prior posts,
> if you want the style shown in your post try this:
>
> mydata2 <- cbind(mydata, score1 = 0, score2 = 0)
> is.m <- mydata$gender == "m"
>
> mydata2[is.m, ] <- transform(mydata[is.m, ],
> score1 = 3 * q1 + q2,
> score2 = 3.5 * q1 + q2
> )
>
> mydata2[!is.m,] <- transform(mydata2[!is.m, ],
> score1 = 2 * q1 + q2,
> score2 = 2.5 * q1 + q2
> )
>
> On 11/25/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > Gabor,
> >
> > Those are handy variations! Perhaps my brain in still in SAS mode on
> > this. I'm expecting something like the code below that checks for male
> > only once, checks for female only when not male (skipping NAs) and
> does
> > all formulas under the appropriate conditions. The formulas I made up
> to
> > keep the code short & may not be as easily modified to let the logical
> > 0/1 values fix them.
> >
> > if gender=="m" then do;
> > Score1=...
> > Score2=
> > ...
> > end;
> > else if gender=="f" then do;
> > Score1=...
> > Score2=
> > ...
> > end;
> >
> > R may not have anything quite like that. R certainly has many other
> > features that SAS lacks.
> >
> > Thanks,
> > Bob
> >
> > =
> > Bob Muenchen (pronounced Min'-chen), Manager
> > Statistical Consulting Center
> > U of TN Office of Information Technology
> > 200 Stokely Management Center, Knoxville, TN 37996-0520
> > Voice: (865) 974-5230
> > FAX: (865) 974-4810
> > Email: [EMAIL PROTECTED]
> > Web: http://oit.utk.edu/scc,
> > News: http://listserv.utk.edu/archives/statnews.html
> > =
> >
> >
> > -Original Message-
> > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
> > Sent: Saturday, November 25, 2006 12:39 AM
> > To: Muenchen, Robert A (Bob)
> > Cc: r-help@stat.math.ethz.ch
> > Subject: Re: [R] Multiple Conditional Tranformations
> >
> > And here is a variation:
> >
> > transform(mydata,
> > score1 = (2 + (gender == "m")) * q1 + q2,
> > score2 = score1 + 0.5 * q1
> > )
> >
> > or
> >
> > transform(
> > transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2),
> > score2 = score1 + 0.5 * q1
> > )
> >
> >
> > On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> > > Try this:
> > >
> > >
> > > transform(mydata,
> > > score1 = (2 + (gender == "m")) * q1 + q2,
> > > score2 = (2.5 + (gender == "m")) * q1 + q2
> > > )
> > >
> > >
> > > On 11/24/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > > > Mark,
> > > >
> > > > I finally got that approach to work by spreading the logical
> > condition
> > > > everywhere. That gets the lengths to match. Still, I can't help
> but
> > > > think there must be a way to specify the logic once per condition.
> > > >
> > > > Thanks,
> > > > Bob
> > > >
> > > > mydata$score1<-numeric(mydata$q1) #just initializing.
> > > > mydata$score2<-numeric(mydata$q1)
> > > > mydata$score1<-NA
> > > > mydata$score2<-NA
> > > > mydata
> > > >
> > > > mydata$score1[mydata$gender == "f"]<-
> > 2*mydata$q1[mydata$gender=="f"] +
> > > >
> > > > mydata$q2[mydata$gender=="f"]
> > > > mydata$score2[mydata$gender ==
> > "f"]<-2.5*my
[R] OT: P(Z <= -1.46).
In checking over the solutions to some homework that I had assigned I observed the fact that in R (version 2.4.0) pnorm(-1.46) gives 0.07214504. The tables in the text book that I am using for the course give the probability as 0.0722. Fascinated, I scanned through 5 or 6 other text books (amongst the dozens of freebies from publishers that lurk on my shelf) and found that some agree with R (giving P(Z <= -1.46) = 0.0721) and some agree with the first text book, giving 0.0722. It is clearly of little-to-no practical import, but I'm curious as to how such a discrepancy would arise in this era. Has anyone any idea? Is there any possibility that the algorithm(s) used to calculate this probability is/are not accurate to 4 decimal places? Could two algorithms ``reasonably'' disagree in the 4th decimal place? cheers, Rolf Turner [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 and provide commented, minimal, self-contained, reproducible code.
[R] hausman Test
Does anyone know how to do an Hausman test? I´ve estimate a modell (some alternatives) with clogit an wanted to test the IIA test (Independence of Irrelevant Alternatives) after estimating a multinomial logit model? Thank you __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Ryacas not working properly
On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote: > Ryacas starts up yacas with an initialization file R.ys which puts > the server in a mode which outputs XML (which is what Ryacas reads). > It does that by specifying --init /whatever/R.ys on the yacas command line > when it is run. > > You can run yacas without that init file or you can just enter into yacas > this to turn off XML output: > >PrettyPrinter() Yeah, the above works! Thanks, Gabor. Paul __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple Conditional Tranformations
That's exactly what I'm looking for. Thanks so much for taking the time to do it that way. On the redundancy issue, I think SAS checks the "else if" condition only if the original "if" is false. The check for f when not m I put in only to exclude missing values for gender. Thanks!! Bob -Original Message- From: Gabor Grothendieck [mailto:[EMAIL PROTECTED] Sent: Saturday, November 25, 2006 7:37 AM To: Muenchen, Robert A (Bob) Cc: r-help@stat.math.ethz.ch Subject: Re: [R] Multiple Conditional Tranformations Firstly your outline does not check once, it checks twice. First it check for "m" and then it redundantly checks for "f". On the other hand the two variations in my post do check once. Although substantially longer than the solutions in my prior posts, if you want the style shown in your post try this: mydata2 <- cbind(mydata, score1 = 0, score2 = 0) is.m <- mydata$gender == "m" mydata2[is.m, ] <- transform(mydata[is.m, ], score1 = 3 * q1 + q2, score2 = 3.5 * q1 + q2 ) mydata2[!is.m,] <- transform(mydata2[!is.m, ], score1 = 2 * q1 + q2, score2 = 2.5 * q1 + q2 ) On 11/25/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote: > Gabor, > > Those are handy variations! Perhaps my brain in still in SAS mode on > this. I'm expecting something like the code below that checks for male > only once, checks for female only when not male (skipping NAs) and does > all formulas under the appropriate conditions. The formulas I made up to > keep the code short & may not be as easily modified to let the logical > 0/1 values fix them. > > if gender=="m" then do; > Score1=... > Score2= > ... > end; > else if gender=="f" then do; > Score1=... > Score2= > ... > end; > > R may not have anything quite like that. R certainly has many other > features that SAS lacks. > > Thanks, > Bob > > = > Bob Muenchen (pronounced Min'-chen), Manager > Statistical Consulting Center > U of TN Office of Information Technology > 200 Stokely Management Center, Knoxville, TN 37996-0520 > Voice: (865) 974-5230 > FAX: (865) 974-4810 > Email: [EMAIL PROTECTED] > Web: http://oit.utk.edu/scc, > News: http://listserv.utk.edu/archives/statnews.html > = > > > -Original Message- > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED] > Sent: Saturday, November 25, 2006 12:39 AM > To: Muenchen, Robert A (Bob) > Cc: r-help@stat.math.ethz.ch > Subject: Re: [R] Multiple Conditional Tranformations > > And here is a variation: > > transform(mydata, > score1 = (2 + (gender == "m")) * q1 + q2, > score2 = score1 + 0.5 * q1 > ) > > or > > transform( > transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2), > score2 = score1 + 0.5 * q1 > ) > > > On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote: > > Try this: > > > > > > transform(mydata, > > score1 = (2 + (gender == "m")) * q1 + q2, > > score2 = (2.5 + (gender == "m")) * q1 + q2 > > ) > > > > > > On 11/24/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote: > > > Mark, > > > > > > I finally got that approach to work by spreading the logical > condition > > > everywhere. That gets the lengths to match. Still, I can't help but > > > think there must be a way to specify the logic once per condition. > > > > > > Thanks, > > > Bob > > > > > > mydata$score1<-numeric(mydata$q1) #just initializing. > > > mydata$score2<-numeric(mydata$q1) > > > mydata$score1<-NA > > > mydata$score2<-NA > > > mydata > > > > > > mydata$score1[mydata$gender == "f"]<- > 2*mydata$q1[mydata$gender=="f"] + > > > > > > mydata$q2[mydata$gender=="f"] > > > mydata$score2[mydata$gender == > "f"]<-2.5*mydata$q1[mydata$gender=="f"] + > > > > > > mydata$q2[mydata$gender=="f"] > > > mydata$score1[mydata$gender == "m"]<-3*mydata$q1[mydata$gender=="m"] > + > > > mydata$q2[mydata$gender=="m"] > > > mydata$score2[mydata$gender == > "m"]<-3.5*mydata$q1[mydata$gender=="m"] + > > > > > > mydata$q2[mydata$gender=="m"] > > > mydata > > > > > > = > > > Bob Muenchen (pronounced Min'-chen), Manager > > > Statistical Consulting Center > > > U of TN Office of Information Technology > > > 200 Stokely Management Center, Knoxville, TN 37996-0520 > > > Voice: (865) 974-5230 > > > FAX: (865) 974-4810 > > > Email: [EMAIL PROTECTED] > > > Web: http://oit.utk.edu/scc, > > > News: http://listserv.utk.edu/archives/statnews.html > > > = > > > > > > > > > -Original Message- > > > From: Leeds, Mark (IED) [mailto:[EMAIL PROTECTED] > > > Sent: Friday, November 24, 2006 8:45 PM > > > To: Muenchen, Robert A (Bob) > > > Subject: RE: [R] Multiple Conditional Tranformations > > > > > > I'm not sure if I understand your question but I don't think you > need > > > iflelse statements. > > > > > > myscore<-numeric(q1) ( because I'm not sure how to initialize a list > so > > > initialize a vec
Re: [R] Ryacas not working properly
Ryacas starts up yacas with an initialization file R.ys which puts
the server in a mode which outputs XML (which is what Ryacas reads).
It does that by specifying --init /whatever/R.ys on the yacas command line
when it is run.
You can run yacas without that init file or you can just enter into yacas
this to turn off XML output:
PrettyPrinter()
See the demo
demo("Ryacas-PrettyPrinter")
So if you are running yacas independently of R then
If you wish to run yacas outside of R just issue the command:
yacas
or maybe
yacas --init /dev/null
(untested)
On 11/25/06, Paul Smith <[EMAIL PROTECTED]> wrote:
> On 11/19/06, Marc Schwartz <[EMAIL PROTECTED]> wrote:
> > this might be a firewall/SELinux/"On Demand Services" problem, I have, I
> > think, nailed it down to two key things on FC6:
> >
> > 1. It requires the use of the '--enable-server' configure option for
> > yacas itself.
> >
> > 2. There is a directory permissions problem in the Ryacas package,
> > restricting a regular user's access to the files located in the 'yacdir'
> > subdirectory.
> >
> > Note that in the instructions below, I also install the GSL, which
> > presumably is not required, but I did it anyway after noting some
> > warnings during the yacas build process.
> >
> > One other thing, which is that when downloading the Ryacas package from
> > the Google site using Firefox, the file is downloaded as:
> >
> > Ryacas_0.2 3.tar.gz
> >
> > Note the missing '-' between the 2 and 3. Be sure to check for this
> > when saving the tarball to disk.
> >
> >
> > So here goes:
> >
> >
> > 1. Install the GSL (including the devel RPM) as root:
> >
> > yum install gsl*
> >
> >
> >
> >
> > 2. Install yacas from the source tarball using:
> >
> > ./ configure --enable-server
> > make
> >
> > then as root:
> >
> > make install
> >
> >
> >
> >
> > 4. Install the Ryacas package as root:
> >
> > R CMD INSTALL Ryacas_0.2-3.tar.gz
> >
> >
> > Be sure that you also have the 'XML' package from CRAN installed, which
> > is a dependency for Ryacas.
> >
> >
> >
> >
> > 3. Change the permissions for /usr/local/lib/R/library/Ryacas/yacdir:
> >
> > Note that the default permission for this directory after installation
> > is:
> >
> > drwxr--r-- 2 root root 4096 Nov 19 15:17 yacdir
> >
> >
> > Thus:
> >
> > su
> > chmod +x /usr/local/lib/R/library/Ryacas/yacdir
> >
> >
> > Now:
> >
> > drwxr-xr-x 2 root root 4096 Nov 19 15:17 yacdir
> >
> >
> >
> > Now in R as a regular user:
> >
> > > library(Ryacas)
> > Loading required package: XML
> >
> > > yacas("5/8 * 3/4")
> > [1] "Starting Yacas!"
> > Accepting requests from port 9734
> > expression(15/32)
> >
> > > yacas("3/7 * 5/9")
> > expression(5/21)
>
> There a side effect that I would like to remove: yacas compiled with
> the option '--enable-server' runs like this
>
> In> x := 2
>
> 2
>
> In> x
>
> 2
>
> In>
>
> i.e., with, apparently, xml commands. How can one run yacas normally,
> without those pieces of xml?
>
> Paul
>
> __
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problems with updating R-packages
Guten Tag Ich hab im R-Forum den Thread Re: [R] Problems with updating R-packages gefunden. Da ich das gleiche Problem (benütze ebenfalls die Linux Distribution Ubuntu) hatte und die Antworten darauf nicht hinreichend waren, habe ich selbst nach den Pakten umgesehen: Nach der Installation von den folgenden Paketen sollte der "lblas-3 Fehler" behoben sein. lapack3, lapack3-dev, lapack3-doc, refblas3, refblas3-dev, refblas3-doc Leider war ich nicht im Stande herauszufinden, wie man in diesem Forum Antworten schreiben kann, aus diesem Grunde melde ich mich direkt an euch und hoffe, dass ihr möglicherweise das Ins Forum schreiben könnted. Ich wünsche ein schönes Wochenende Gruss Reto Bürgin __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple Conditional Tranformations
Firstly your outline does not check once, it checks twice. First it check for "m" and then it redundantly checks for "f". On the other hand the two variations in my post do check once. Although substantially longer than the solutions in my prior posts, if you want the style shown in your post try this: mydata2 <- cbind(mydata, score1 = 0, score2 = 0) is.m <- mydata$gender == "m" mydata2[is.m, ] <- transform(mydata[is.m, ], score1 = 3 * q1 + q2, score2 = 3.5 * q1 + q2 ) mydata2[!is.m,] <- transform(mydata2[!is.m, ], score1 = 2 * q1 + q2, score2 = 2.5 * q1 + q2 ) On 11/25/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote: > Gabor, > > Those are handy variations! Perhaps my brain in still in SAS mode on > this. I'm expecting something like the code below that checks for male > only once, checks for female only when not male (skipping NAs) and does > all formulas under the appropriate conditions. The formulas I made up to > keep the code short & may not be as easily modified to let the logical > 0/1 values fix them. > > if gender=="m" then do; > Score1=... > Score2= > ... > end; > else if gender=="f" then do; > Score1=... > Score2= > ... > end; > > R may not have anything quite like that. R certainly has many other > features that SAS lacks. > > Thanks, > Bob > > = > Bob Muenchen (pronounced Min'-chen), Manager > Statistical Consulting Center > U of TN Office of Information Technology > 200 Stokely Management Center, Knoxville, TN 37996-0520 > Voice: (865) 974-5230 > FAX: (865) 974-4810 > Email: [EMAIL PROTECTED] > Web: http://oit.utk.edu/scc, > News: http://listserv.utk.edu/archives/statnews.html > = > > > -Original Message- > From: Gabor Grothendieck [mailto:[EMAIL PROTECTED] > Sent: Saturday, November 25, 2006 12:39 AM > To: Muenchen, Robert A (Bob) > Cc: r-help@stat.math.ethz.ch > Subject: Re: [R] Multiple Conditional Tranformations > > And here is a variation: > > transform(mydata, > score1 = (2 + (gender == "m")) * q1 + q2, > score2 = score1 + 0.5 * q1 > ) > > or > > transform( > transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2), > score2 = score1 + 0.5 * q1 > ) > > > On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote: > > Try this: > > > > > > transform(mydata, > > score1 = (2 + (gender == "m")) * q1 + q2, > > score2 = (2.5 + (gender == "m")) * q1 + q2 > > ) > > > > > > On 11/24/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote: > > > Mark, > > > > > > I finally got that approach to work by spreading the logical > condition > > > everywhere. That gets the lengths to match. Still, I can't help but > > > think there must be a way to specify the logic once per condition. > > > > > > Thanks, > > > Bob > > > > > > mydata$score1<-numeric(mydata$q1) #just initializing. > > > mydata$score2<-numeric(mydata$q1) > > > mydata$score1<-NA > > > mydata$score2<-NA > > > mydata > > > > > > mydata$score1[mydata$gender == "f"]<- > 2*mydata$q1[mydata$gender=="f"] + > > > > > > mydata$q2[mydata$gender=="f"] > > > mydata$score2[mydata$gender == > "f"]<-2.5*mydata$q1[mydata$gender=="f"] + > > > > > > mydata$q2[mydata$gender=="f"] > > > mydata$score1[mydata$gender == "m"]<-3*mydata$q1[mydata$gender=="m"] > + > > > mydata$q2[mydata$gender=="m"] > > > mydata$score2[mydata$gender == > "m"]<-3.5*mydata$q1[mydata$gender=="m"] + > > > > > > mydata$q2[mydata$gender=="m"] > > > mydata > > > > > > = > > > Bob Muenchen (pronounced Min'-chen), Manager > > > Statistical Consulting Center > > > U of TN Office of Information Technology > > > 200 Stokely Management Center, Knoxville, TN 37996-0520 > > > Voice: (865) 974-5230 > > > FAX: (865) 974-4810 > > > Email: [EMAIL PROTECTED] > > > Web: http://oit.utk.edu/scc, > > > News: http://listserv.utk.edu/archives/statnews.html > > > = > > > > > > > > > -Original Message- > > > From: Leeds, Mark (IED) [mailto:[EMAIL PROTECTED] > > > Sent: Friday, November 24, 2006 8:45 PM > > > To: Muenchen, Robert A (Bob) > > > Subject: RE: [R] Multiple Conditional Tranformations > > > > > > I'm not sure if I understand your question but I don't think you > need > > > iflelse statements. > > > > > > myscore<-numeric(q1) ( because I'm not sure how to initialize a list > so > > > initialize a vector with q1 elements ) > > > > > > myscore<-NA ( I think this should set all the values in myscore to > NA ) > > > myscore[mydata$gender == f]<-2*mydata$q1 + mydata$q2 > > > myscore[mydata$gender == m]<-3*mydata$q1 + mydata$q2 > > > > > > the above should do what you do in the first part of your code but I > > > don't know if that was your question ? > > > also, it does it making myscore a vector because I didn't know how > to > > > initialize a list. > > > Someone else may goive a better solution. I'm no expert. > > >
Re: [R] Multiple Conditional Tranformations
Gabor,
Those are handy variations! Perhaps my brain in still in SAS mode on
this. I'm expecting something like the code below that checks for male
only once, checks for female only when not male (skipping NAs) and does
all formulas under the appropriate conditions. The formulas I made up to
keep the code short & may not be as easily modified to let the logical
0/1 values fix them.
if gender=="m" then do;
Score1=...
Score2=
...
end;
else if gender=="f" then do;
Score1=...
Score2=
...
end;
R may not have anything quite like that. R certainly has many other
features that SAS lacks.
Thanks,
Bob
=
Bob Muenchen (pronounced Min'-chen), Manager
Statistical Consulting Center
U of TN Office of Information Technology
200 Stokely Management Center, Knoxville, TN 37996-0520
Voice: (865) 974-5230
FAX: (865) 974-4810
Email: [EMAIL PROTECTED]
Web: http://oit.utk.edu/scc,
News: http://listserv.utk.edu/archives/statnews.html
=
-Original Message-
From: Gabor Grothendieck [mailto:[EMAIL PROTECTED]
Sent: Saturday, November 25, 2006 12:39 AM
To: Muenchen, Robert A (Bob)
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Multiple Conditional Tranformations
And here is a variation:
transform(mydata,
score1 = (2 + (gender == "m")) * q1 + q2,
score2 = score1 + 0.5 * q1
)
or
transform(
transform(mydata, score1 = (2 + (gender == "m")) * q1 + q2),
score2 = score1 + 0.5 * q1
)
On 11/25/06, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> Try this:
>
>
> transform(mydata,
> score1 = (2 + (gender == "m")) * q1 + q2,
> score2 = (2.5 + (gender == "m")) * q1 + q2
> )
>
>
> On 11/24/06, Muenchen, Robert A (Bob) <[EMAIL PROTECTED]> wrote:
> > Mark,
> >
> > I finally got that approach to work by spreading the logical
condition
> > everywhere. That gets the lengths to match. Still, I can't help but
> > think there must be a way to specify the logic once per condition.
> >
> > Thanks,
> > Bob
> >
> > mydata$score1<-numeric(mydata$q1) #just initializing.
> > mydata$score2<-numeric(mydata$q1)
> > mydata$score1<-NA
> > mydata$score2<-NA
> > mydata
> >
> > mydata$score1[mydata$gender == "f"]<-
2*mydata$q1[mydata$gender=="f"] +
> >
> > mydata$q2[mydata$gender=="f"]
> > mydata$score2[mydata$gender ==
"f"]<-2.5*mydata$q1[mydata$gender=="f"] +
> >
> > mydata$q2[mydata$gender=="f"]
> > mydata$score1[mydata$gender == "m"]<-3*mydata$q1[mydata$gender=="m"]
+
> > mydata$q2[mydata$gender=="m"]
> > mydata$score2[mydata$gender ==
"m"]<-3.5*mydata$q1[mydata$gender=="m"] +
> >
> > mydata$q2[mydata$gender=="m"]
> > mydata
> >
> > =
> > Bob Muenchen (pronounced Min'-chen), Manager
> > Statistical Consulting Center
> > U of TN Office of Information Technology
> > 200 Stokely Management Center, Knoxville, TN 37996-0520
> > Voice: (865) 974-5230
> > FAX: (865) 974-4810
> > Email: [EMAIL PROTECTED]
> > Web: http://oit.utk.edu/scc,
> > News: http://listserv.utk.edu/archives/statnews.html
> > =
> >
> >
> > -Original Message-
> > From: Leeds, Mark (IED) [mailto:[EMAIL PROTECTED]
> > Sent: Friday, November 24, 2006 8:45 PM
> > To: Muenchen, Robert A (Bob)
> > Subject: RE: [R] Multiple Conditional Tranformations
> >
> > I'm not sure if I understand your question but I don't think you
need
> > iflelse statements.
> >
> > myscore<-numeric(q1) ( because I'm not sure how to initialize a list
so
> > initialize a vector with q1 elements )
> >
> > myscore<-NA ( I think this should set all the values in myscore to
NA )
> > myscore[mydata$gender == f]<-2*mydata$q1 + mydata$q2
> > myscore[mydata$gender == m]<-3*mydata$q1 + mydata$q2
> >
> > the above should do what you do in the first part of your code but I
> > don't know if that was your question ?
> > also, it does it making myscore a vector because I didn't know how
to
> > initialize a list.
> > Someone else may goive a better solution. I'm no expert.
> >
> >
> > -Original Message-
> > From: [EMAIL PROTECTED]
> > [mailto:[EMAIL PROTECTED] On Behalf Of Muenchen,
Robert
> > A (Bob)
> > Sent: Friday, November 24, 2006 8:27 PM
> > To: r-help@stat.math.ethz.ch
> > Subject: [R] Multiple Conditional Tranformations
> >
> > Greetings,
> >
> >
> >
> > I'm learning R and I'm stuck on a basic concept: how to specify a
> > logical condition once and then perform multiple transformations
under
> > that condition. The program below is simplified to demonstrate the
goal.
> > Its results are exactly what I want, but I would like to check the
> > logical state of gender only once and create both (or any number of)
> > scores at once.
> >
> >
> >
> > mystring<-
> >
> > ("id,group,gender,q1,q2,q3,q4
> >
> > 01,1,f,2,2,5,4
> >
> > 02,2,f,2,1,4,5
> >
> > 03,1,f,2,2,4,4
> >
> > 04,2,f,1,1,5,5
> >
> > 05,1,m,4,5,4,
> >
> > 06,2,m,5,4,5,5
> >
> > 07,1,m,3,3,4,5
> >
>
Re: [R] Diebold Mariano Test
Hello Graham
Pls find attached some old R code for the Diebold Mariano test. The code
is at least 6 years old and was not used in the meantime. Pls check
first before using it.
Best regards
Adrian
Dear List
Has anyone used R to distnguish between alternative forecasting models? In particular
is the Diebold Mariano test available for use within R.
Any assistance would be greatly appreciated.
Graham Leask
Lecturer in Strategy
Economics & Strategy Group
Aston University
Aston Triangle
Birmingham
B4 7ET
Tel: 0121 204 3150
E Mail: [EMAIL PROTECTED]
diebold.mariano.test <- function(x, alternative = c("two.sided", "less",
"greater"), k)
{
if (NCOL(x) > 1)
stop("x is not a vector or univariate time series")
if (any(is.na(x)))
stop("NAs in x")
alternative <- match.arg(alternative)
DNAME <- deparse(substitute(x))
n <- NROW(x)
cv <- acf(x, lag.max=k, type="covariance", plot=FALSE)$acf[,,1]
eps <- 1.0e-8
vr <- max(eps, sum(c(cv[1], 2*cv[-1])) / n)
STATISTIC <- mean(x) / sqrt(vr)
names(STATISTIC) <- "Standard Normal"
METHOD <- "Diebold-Mariano Test"
if (alternative == "two.sided")
PVAL <- 2 * pnorm(-abs(STATISTIC))
else if (alternative == "less")
PVAL <- pnorm(STATISTIC)
else if (alternative == "greater")
PVAL <- pnorm(STATISTIC, lower.tail = FALSE)
PARAMETER <- k
names(PARAMETER) <- "Truncation lag"
structure(list(statistic = STATISTIC, parameter = PARAMETER, alternative =
alternative,
p.value = PVAL, method = METHOD, data.name = DNAME),
class = "htest")
}
g <- function(x)
{
abs(x)
}
e1 <- rnorm(500)
e2 <- rnorm(500)
diebold.mariano.test(g(e1)-g(e2), k = 3)
e1 <- rnorm(500, sd=1)
e2 <- rnorm(500, sd=1.3)
diebold.mariano.test(g(e1)-g(e2), k = 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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Ryacas not working properly
On 11/19/06, Marc Schwartz <[EMAIL PROTECTED]> wrote:
> this might be a firewall/SELinux/"On Demand Services" problem, I have, I
> think, nailed it down to two key things on FC6:
>
> 1. It requires the use of the '--enable-server' configure option for
> yacas itself.
>
> 2. There is a directory permissions problem in the Ryacas package,
> restricting a regular user's access to the files located in the 'yacdir'
> subdirectory.
>
> Note that in the instructions below, I also install the GSL, which
> presumably is not required, but I did it anyway after noting some
> warnings during the yacas build process.
>
> One other thing, which is that when downloading the Ryacas package from
> the Google site using Firefox, the file is downloaded as:
>
> Ryacas_0.2 3.tar.gz
>
> Note the missing '-' between the 2 and 3. Be sure to check for this
> when saving the tarball to disk.
>
>
> So here goes:
>
>
> 1. Install the GSL (including the devel RPM) as root:
>
> yum install gsl*
>
>
>
>
> 2. Install yacas from the source tarball using:
>
> ./ configure --enable-server
> make
>
> then as root:
>
> make install
>
>
>
>
> 4. Install the Ryacas package as root:
>
> R CMD INSTALL Ryacas_0.2-3.tar.gz
>
>
> Be sure that you also have the 'XML' package from CRAN installed, which
> is a dependency for Ryacas.
>
>
>
>
> 3. Change the permissions for /usr/local/lib/R/library/Ryacas/yacdir:
>
> Note that the default permission for this directory after installation
> is:
>
> drwxr--r-- 2 root root 4096 Nov 19 15:17 yacdir
>
>
> Thus:
>
> su
> chmod +x /usr/local/lib/R/library/Ryacas/yacdir
>
>
> Now:
>
> drwxr-xr-x 2 root root 4096 Nov 19 15:17 yacdir
>
>
>
> Now in R as a regular user:
>
> > library(Ryacas)
> Loading required package: XML
>
> > yacas("5/8 * 3/4")
> [1] "Starting Yacas!"
> Accepting requests from port 9734
> expression(15/32)
>
> > yacas("3/7 * 5/9")
> expression(5/21)
There a side effect that I would like to remove: yacas compiled with
the option '--enable-server' runs like this
In> x := 2
2
In> x
2
In>
i.e., with, apparently, xml commands. How can one run yacas normally,
without those pieces of xml?
Paul
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using multiple objects in a for loop
You should use
assign("temp",get(on[j]))
or just
temp<-get(on[j])
as the on[j] IS just a character string. If you want to get the object
with that name, use get.
Ales Ziberna
Daniel Yanosky pravi:
> I have a situation where I want to perform the same manipulations to multiple
> R objects in series. I have constructed a vector () to serve as a list of the
> objects. However, my assign() assigns the object name (on[j]) as a character
> string to "temp" and not the object itself.
>
> ===
>
> fn<-c("FT1.1.01.RData","FT1.1.02.RData")
>
> on<-c("FT1101","FT1102")
>
> b<-1
> e<-2
>
> for (i in 1:2){
>
> for (j in b:e){
>
> load(fn[j])
> assign("temp",on[j])
>
> .
> .
> .
> rm(paste(stderr[i],".t",sep=""))
> rm(temp)
> rm(fn[i])
> }
>
> }
>
> ===
>
> Daniel
>
> __
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] predict and arima
Part of the fit is the Kalman filter state after running the model
forwards. Try reversing your series, fitting and then forecasting.
You might have more success in understanding arima0.
On Sat, 25 Nov 2006, Franck Arnaud wrote:
> Hi all,
>
> Forecasting from an arima model is easy with predict.
>
> But I can't manage to backcast : invent data from the model before the
> begining of the sample.
> The theory is easy : take your parameters, reverse your data, forecast, and
> then reverse the forecast
> I've tried to adapt the predict function to do that (i'm not sure that the
> statistical procedure is fine (with the residuals), but that's not my point
> right now) :
>
> mav.backcast.arima<-function(model,n.backcast,...)
> {
>if (class(model)[1]!="Arima") stop("argument must be an object
> of class 'Arima' (see ?arima)")
>
>model2<-model
>model$residuals<-rev(model$residuals)
>if (is.ts(model2$residuals))
> model$residuals<-ts(model$residuals,start=start(model2$residuals),
>frequency=frequency(model2$residuals))
>pred.before<-predict(model,n.ahead=n.backcast,...)
>
>freq<-frequency(model$residuals)
>startingdate<-per.sub(start(model2$residuals),n.backcast,freq=freq)
>
>pred<-ts(rev(pred.before$pred),start=startingdate,freq=freq)
>se<-ts(rev(pred.before$se),start=startingdate,freq=freq)
>
>return((list(pred = pred, se =se)))
> }
>
> This function does not work : it gives always the same result, it does not
> depend on the residuals (i've tried to insert
> a model$residuals<-rep(1,100) after the definition, to check that).
>
> Then i look at the code, with getS3method("predict","Arima"). And i get even
> more confused (!) :
> where does data play a role in the function ? residuals are loaded into rsd,
> but this variable is not used after...
> I looked at KalmanForecast and at the C code of KalmanFore, but it did not
> help me understand what was going on.
>
> thanks
>
> Franck A.
>
> btw, it has nothing to do with it, but i've done some stuff on time series
> (filtering with Hodrick prescott or Baxter King, for instance) that you can
> find on http://arnaud.ensae.net
>
> [[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
> and provide commented, minimal, self-contained, reproducible code.
>
--
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
and provide commented, minimal, self-contained, reproducible code.
