Re: [R] Help with a more flexible funtion for multiple comparisio n of means
Jose -
Before implementing SNK and Duncan's, you may want to be aware of some
criticisms of these methods:
From Hsu (1996),
Newman-Keuls multiple range test is not a confident inequalities method and
cannot be recommended.
Duncan's multiple range test is not a confident inequalities method and
cannot be recommended either. In the words of Tukey (1991), Duncan's
multiple range test was a 'distraction' in the history of multiple
comparisons, amounting to 'talking 5% while using more than 5%
simultaneous'
@Book{Hsu1996,
author = {Jason C. Hsu},
title ={Multiple Comparisons: Theory and Methods},
publisher ={Chapman Hall},
year = {1996}
}
-- Jim
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{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
Re: [R] Anova with Scheffe Tests
In good ol' Statview (now dearly departed) to complete a Scheffe
test you selected the independent variables and dependent variable
and it produced a table with the pairwise comparisons of the levels
of the factor. I'm looking for a system that is as basic, but can be
done using R and has documentation so I'm not guessing what I'm
doing. I'd rather not have to do plots in R and then run over to
dead software to do Scheffe's if possible.
I checked on google and there seems to be code for a couple of
functions out there, but I need something that has a manual.
Is there a Scheffe function out there that is reasonably well
documented, or should I consider some other method of dealing with
this data. We have been using Scheffe for this type of analysis as I
was under the impression it was very conservative. Tukey's HSD seems
to be conservative as well. Should I try this? Is there a different
approacch that is better and where can I read about it.
Thanks for any help you can provide.
Sam
It sounds like you are only interested in simple pairwise comparisons, in
which case you are better off using Tukey's HSD. This will still protect the
family-wise error rate, but will allow you to infer more differences than
you would using Scheffé's method. Scheffé's method is exact if you are truly
interested in all contrasts (a situation I have never come across), but it
is overly conservative when inferences are only made for pairwise
differences.
A geometric explanation of the difference between the two methods can be
found in Jason Hsu's book Multiple Comparisons: Theory and Methods.
As you can also read in that book, there is no good reason to precede
Tukey's HSD with an ANOVA F-test. Tukey's HSD controls the family-wise error
rate anyway, so an initial F-test is superfluous.
--Jim
James A. Rogers
Manager, Nonclinical Statistics
PGRD Groton Labs
Eastern Point Road (MS 8260-1331)
Groton, CT 06340
office: (860) 686-0786
fax: (860) 715-5445
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{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
Re: Re: [R] gnls or nlme : how to obtain confidence intervals of fitted values
I believe what Pierre means to be asking for are prediction intervals, or
possibly tolerance intervals, which are different from confidence
intervals. A great reference is:
Hahn G and Meeker WQ: Statistical Intervals. John Wiley and Sons, New York,
1991
For example, in a one-sample Normal problem, \bar{X} \pm 2
\hat{sigma}/\sqrt{n} is an approximate 95% confidence interval, while
\bar{X} \pm 2 \hat{sigma} is an approximate 95% prediction interval
(assuming \bar{X} and \hat{sigma} are high precision estimates; if they are
not you probably want to consider tolerance intervals, which are discussed
in the above reference).
Both intervals.lme{nlme} and estimable{gmodels} will give you confidence
intervals, but neither will give you prediction intervals. I am not aware of
any published R function that gives you prediction intervals or tolerance
intervals for lme models. It is not easy to write such a function for the
general case, but it may be relatively easy to write your own for special
cases of lme models.
Jim
Message: 9
Date: Mon, 4 Oct 2004 21:42:20 +1000
From: Andrew Robinson [EMAIL PROTECTED]
Subject: Re: [R] gnls or nlme : how to obtain confidence intervals of
fitted values
To: David Scott [EMAIL PROTECTED]
Cc: [EMAIL PROTECTED], Spencer Graves [EMAIL PROTECTED]
Message-ID: [EMAIL PROTECTED]
Content-Type: text/plain; charset=us-ascii
The function estimable() from the gmodels part of the gregmisc package
will do this, if applied appropriately.
It allows for the estimation of arbitrary linear combinations of the
parameter estimates of a model object, and calcualtes confidence
intervals.
I hope that this helps,
Andrew
On Tue, Oct 05, 2004 at 12:31:48AM +1300, David Scott wrote:
On Mon, 4 Oct 2004, Pierre MONTPIED wrote:
Thanks Spencer but the intervals function gives confidence intervals of
the parameters of the model not the predicted values. In the Soybean
example it would be the CI of predicted weight for a given time,
knowing
all the parameters (Asym, xmid, scal, variance function and residual)
and
their distributions.
And what I need to calculate is precisely this CI.
My question is therefore is there an analytical way to calculate such
CI,
whatever the model, or could I try some randomizing techniques such as
bootstrap or other ?
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] 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] anova and post hoc multicomparison tests
Guillaume,
Your comments are a compliment to R.
Undoubtedly other software is preferable if you want to do
Student-Newman-Keuls or Fisher's protected LSD (ANOVA F-test followed by
unadjusted T-tests). Perhaps the reason is that neither Student-Newman-Keuls
nor Fisher's protected LSD is a valid multiple comparison procedure.
Student-Newman-Keuls does not even control the probability of making at
least one false assertion of inequality (which is the almost the minimum one
could ask of a multiple comparison procedure). For details, including
examples of where these methods fail, see:
Hsu, J.C. (1996). Multiple Comparisons: Theory and Methods. Chapman Hall.
If you want to use R to perform valid multiple comparisons, such as
Dunnett's MCC or Tukey's HSD, see the function TukeyHSD and also the
multcomp package.
Jim Rogers
Message: 78
Date: Fri, 24 Sep 2004 10:54:55 +0200
From: BLANCHER Guillaume [EMAIL PROTECTED]
Subject: [R] anova and post hoc multicomparison tests
To: [EMAIL PROTECTED]
Message-ID: [EMAIL PROTECTED]
Content-Type: text/plain; charset=us-ascii
Hello everyone,
Like a lot of people, I have been looking for functions in R doing ANOVA
(ok) and performing multicomparisons (like Student-Newman-Keuls, etc.).
As I have been a little bit disappointed, I have bee looking through the
net for such open source softwares. I found one in:
http://www.statpages.org/miller/openstat/OS4.html
I have begun to use it, and it seems good and simple to understand (as
for a non-specialist like me).
Sorry for R, but I prefer OpenStat4 to R for ANOVAs and post hoc tests.
Guillaume
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] 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] confidence intervals
Robert,
I have this quick hack to obtain approximate Shewhart prediction intervals
for variance component models fit with lme (to nitpick slightly,
confidence intervals have the interpretation of containing parameters,
while prediction and tolerance intervals have the interpretation of
containing future observations or statistics).
Back of the envelope documentation: the only argument that probably needs
explaining is the reps argument to shewhart(). If your model is, e.g.
fixed = Y ~ 1, random = ~ 1 | Batch
then specify reps = c(1, 1) if you want to predict a single future
observation from a single future batch, reps = c(1, 2) if you want to
predict the mean of two future observations from a single future batch, reps
= c(2, 2) if you want to predict the mean of 4 observations spread evenly
over 2 future batches, ...
Leave mult.check = 1, unless you want to do a Bonferroni correction.
HTH,
Jim Rogers
valStats2 -
function (x, fixed, random, ...)
{
mod - lme(fixed = fixed, data = x, random = random, ...)
mn - fixef(mod)
vc - VarCorr(mod)
err - Expecting only random intercept terms and a single fixed
intercept.\n
if (length(mn) 1 || ncol(vc) 2)
stop(err)
rn - rownames(vc)
skip - regexpr(=, rn) 0
if (!any(skip))
vnms - attr(vc, title)
else vnms - grep(=, rn, value = TRUE)
vc - vc[!skip, ]
vnms - trim(sub(=.*, , vnms))
vnms - c(vnms, Residual)
vnms - paste(V, vnms, sep = .)
vars - as.numeric(vc[, Variance])
stats - c(mn, vars)
names(stats) - c(Intercept, vnms)
stats
}
shewhart -
function (x, meancol = Intercept, varcols = grep(^V\\., names(x), value
= TRUE), reps = c(1, 1), alpha = 0.02, mult.check = 1)
{
mn - x[[meancol]]
vr - as.matrix(x[varcols])
totvar - vr %*% (1/reps)
totsd - sqrt(totvar)
LL.mean - mn + qnorm(alpha/2/mult.check) * totsd
UL.mean - mn + qnorm(1 - alpha/2/mult.check) * totsd
out - data.frame(V.Total = totvar, LL.mean = LL.mean, UL.mean =
UL.mean)
out
}
### Example, where x is your data.frame:
foo - valStats2(x, fixed = Y ~ 1, random = ~ 1|Batch)
foo - as.data.frame(t(as.matrix(foo)))
data.frame(foo, shewhart(foo))
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of Spencer Graves
Sent: Friday, September 03, 2004 11:58 AM
To: [EMAIL PROTECTED]
Cc: Robert Waters; [EMAIL PROTECTED]
Subject: Re: [R] confidence intervals
Hi, Robert:
While it may be difficult to program this in general
(as suggested
by it's position on Doug's To Do list), all the pieces should be
available to support a special script for your specific application.
What fixed and random model(s) interest you most?
hope this helps. spencer graves
Douglas Bates wrote:
Robert Waters wrote:
Dear R users;
Im working with lme and Id like to have an idea of how
can I get CI for the predictions made with the model.
Im not a stats guy but, if Im not wrong, the CIs
should be different if Im predicting a new data point
or a new group. Ive been searching through the web and
in help-lists with no luck. I know this topic had been
asked before but without replies. Can anyone give an
idea of where can I found information about this or
how can I get it from R?
Thanks for any hint
That's not currently implemented in lme. It's on the To
Do list but
it is not very close to the top.
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] 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 nls to fit a four parameter logistic model
Shalini,
I think your hill equation is meant to just be an alternative
parameterization of the four parameter logistic (BTW, the hill
*coefficient* is a function of the slope parameter of the FPL, but I don't
believe hill equation is standard terminology). Note conc is the input
in this parameterization, not log(conc).
nls(log(il10)~A+(B-A)/(1+(conc/xmid )^scal),data=test,
+ start = list(A=3.5, B=15,
+ xmid=600,scal=1/2.5))
Nonlinear regression model
model: log(il10) ~ A + (B - A)/(1 + (conc/xmid)^scal)
data: test
A Bxmidscal
14.7051665 3.7964534 607.9822962 0.3987786
residual sum-of-squares: 0.1667462
To see the equivalence to the other parametrization that you used, note
1/2.507653
[1] 0.3987793
log(607.9822962)
[1] 6.410146
--Jim
Message: 17
Date: Mon, 16 Aug 2004 11:25:57 -0500
From: [EMAIL PROTECTED]
Subject: [R] using nls to fit a four parameter logistic model
To: [EMAIL PROTECTED]
Message-ID:
[EMAIL PROTECTED]
Content-Type: text/plain; charset=US-ASCII
I am working on what appears to be a fairly simple problem for the
following data
test=data.frame(cbind(conc=c(25000, 12500, 6250, 3125, 1513, 781, 391,
195, 97.7, 48.4, 24, 12, 6, 3, 1.5, 0.001),
il10=c(330269, 216875, 104613, 51372, 26842, 13256, 7255, 3049, 1849,
743,
480, 255, 241, 128, 103, 50)))
I am able to fit the above data to the equation
nls(log(il10)~A+(B-A)/(1+exp((xmid-log(conc))/scal)),data=test,
+ start = list(A=log(0.001), B=log(10),
+ xmid=log(6000),scal=0.8))
Nonlinear regression model
model: log(il10) ~ A + (B - A)/(1 + exp((xmid - log(conc))/scal))
data: test
A B xmid scal
3.796457 14.705159 6.410144 2.507653
residual sum-of-squares: 0.1667462
But in attempting to achieve a fit to what is commonly known as the hill
equation, which is a four parameter fit that is used widely in biological
data analysis
nls(log(il10)~A+(B-A)/(1+(log(conc)/xmid )^scal),data=test,
+ start = list(A=log(0.001), B=log(10), xmid=log(6000),scal=0.8))
Nonlinear regression model
model: log(il10) ~ A + (B - A)/(1 + (log(conc)/xmid )^scal)
Error in numericDeriv(form[[3]], names(ind), env) :
Missing value or an Infinity produced when evaluating the model
Please would someone offer a suggestion
Shalini
James A. Rogers
Manager, Nonclinical Statistics
PGRD Groton Labs
Eastern Point Road (MS 260-1331)
Groton, CT 06340
office: (860) 686-0786
fax: (860) 715-5445
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] 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] sporadic errors with nlrq() / optim()
Dear List,
Apologies if this is a known problem ... I wasn't able to find it on the bug
list, but it is a problem that does not seem to occur with a MAC build of R
2.0, so perhaps this problem has already been addressed for the future.
I am getting *sporadic* errors when refitting the same model to the same
data set, using nlrq() in the nlrq package. The algorithm is not stochastic,
so I would expect to get errors either every time, or never.
###
library(stats) # or library(nls) if using R 1.9.0
library(nlrq)
test - data.frame(x = c(7.60090245954208, 6.90775527898214,
6.21460809842219, 5.52146091786225, 4.60517018598809, 3.91202300542815,
3.2188758248682 , 2.52572864430826, 1.83258146374831, 7.60090245954208,
6.90775527898214, 6.21460809842219, 5.52146091786225, 4.60517018598809,
3.91202300542815, 3.2188758248682 , 2.52572864430826, 1.83258146374831,
6.21460809842219, 6.21460809842219),
y = c( 11.0161506644269, 9.84267541313937,
8.66146668057266, 7.48099216286952, 6.50578406012823, 6.24027584517077,
5.63121178182137, 5.71702770140622, 5.64190707093811, 10.8983676287705,
9.91857318995417, 8.74608021735751, 7.58120982619635, 6.361302477573 ,
5.91889385427315, 5.63835466933375 , 5.80211837537706, 5.64897423816121,
8.6195692580331 , 8.70367275835886)
)
i - 1
while(i 500) {# I usually hit an error within 50 iterations
cat(i, \n)
nlrq(y ~ SSfpl(x, A, B, xmid, scal), data = test)
i - i + 1
}
###
Errors occur with version 1.8.1 and 1.9.0 on both Windows (two different
machines) and UNIX, but not on version 2.0 on a MAC (these are the only R
version - OS permutations I was able to get reports on easily).
Anyone understand what is happening here?
Thanks,
Jim
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to generate a report with graphics and tables?
Just to expand on an earlier suggestion:
I have some data sets which change on a daily bases. So far I have
imported these sets into R, done all my evaluations resulting in a
couple of plots, charts and tables of numbers which I copypasted via
clipboard into Powerpoint.
The procedure is always the same and I wonder, whether there is no
easier way for doing so. Is there some type of report generator
available or some HowTo on this.
Can anybody give me a hint on where to look for?
Regards,
Olaf Bürger
...
One other possibility is to use the R2HTML package (and maybe the xtable
package, too) to write the `report' in HTML.
HTH,
Andy
R2HTML is not (yet) as far along as Sweave, but depending on how much
sympathy you have for your clients, html can be a very nice solution. An
important requirement for many reports is that computer-naive clients be
able to edit and copy from the reports to an MS application, preferably
without severe font and formatting problems and without needing to
install/understand foo2bar type applications. See comments to this
effect in the recent R2HTML article by the package author, Eric Lecoutre in
the most recent R News (vol 3/3). As Andy noted, xtable is nice to use in
conjunction with R2HTML (just write a little xtable method like the one
following this message):
In particular, a naive user viewing the html file with Internet Explorer (I
know v. 5.5 or greater will work) will have an option, under the File menu
to Edit with Microsoft Word for Windows. This assumes they have set MS
Word as their HTML editor within IE. It also assumes some minimal formatting
requirements for the HTML file (but you can worry about this, the client
won't have to).
I find this solution keeps everyone happy:
1. I can auto-generate text source for my documents;
2. clients without proprietary software can view the report;
3. clients with standard MS software but without a clue can easily edit and
copy.
Cheers,
Jim
HTML.xtable - function(x, file = .HTML.file, append = TRUE) {
sink(file = file, append = append)
print(x, type = html)
sink()
}
James A. Rogers
Manager, Nonclinical Statistics
PGRD Groton Labs
Eastern Point Road (MS 8260-1331)
Groton, CT 06340
office: (860) 686-0786
fax: (860) 715-5445
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] RE: Trellis graph and two colors display
From: Patrick Giraudoux [EMAIL PROTECTED]
Hello,
I would like to display two groups of dots in different colors
or style and additionnaly a linear regression for all the dots in
panel plots of a Trellis graph. Actually to get in each panel
the equivalent of:
plot(x[mois==3],y[mois==3],col=blue)
points(x[mois==9],y[mois==9],col=red)
abline(lm(y~x), col=green)
mois being a grouping variable and ID (see below) the
conditioning variable in a data.frame of 230 rows
After several really disatreous trials, I have tried the following:
xyplot(y~x|ID,panel=function(x,y){panel.superpose(x,y,subscript
s=1:230,groups=mois);panel.abline(lm(y~x),col=green)})
You don't need to mess around with the subscripts argument for this. I think
the problem is that you are not passing the necessary graphical parameters
to your panel function. You could do:
dat - expand.grid(ID = 1:4, mois = c(3, 9))
dat - dat[rep(seq(nrow(dat)), rep(10, nrow(dat))), ]
dat - data.frame(dat, x = rnorm(nrow(dat)), y = rnorm(nrow(dat)))
dat$ID - factor(dat$ID)
lset(list(superpose.symbol = list(col = c(blue, red), pch = c(1, 1
xyplot(y ~ x | ID, data = dat,
groups = mois,
panel = function(x, y, ...) {
panel.superpose(x, y, ...)
panel.abline(lm(y~x), col = green)
})
Cheers,
Jim
James A. Rogers
Manager, Biometrics
PGRD Groton Labs
Eastern Point Road (MS 8260-1331)
Groton, CT 06340
office: (860) 686-0786
fax: (860) 715-5445
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] minor documentation typo ?
Hi,
I believe the following is a very minor typo in the documentation for
merge() :
all logical; all=L is shorthand for all.x=L and all.y=L
^^ ^
Looks like 'L' should be 'FALSE'.
Someone please alert me if I am wrong. If I am not wrong, is it appropriate
to submit a bug for such minor detail?
Thanks,
Jim
James A. Rogers
Senior Coordinator, Biometrics
PGRD Groton Labs
Eastern Point Road (MS 8260-1331)
Groton, CT 06340
office: (860) 686-0786
fax: (860) 715-5445
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
