[R] confidence interval around Aalen-Johansen if no events

[Andreas Leha]

Hi all,

I'd like to calculate confidence intervals around Aalen-Johansen
estimates at given time points.  And later I'd like to be able to
compare the Aalen-Johansen estimates at given time points between
groups.
My problem is, that there are no events recorded.

So my question is:
Does anyone know how to get confidence intervals around Aalen-Johansen
estimates in the absence of (or prior to) events?
(And ideally also how to compare the Aalen-Johansen estimates at given
time points between groups in that situation.)

Note that I can do the same thing for Kaplan-Meier estimates (using
the bpcp package).

Many thanks in advance!

Best,
Andreas

PS:  Some code that illustrates my problem:

-
library("dplyr")
library("survival")

## read data
sdat <- structure(list(time = c(2, 189, 182, 2, 184, 179, 179, 159, 18, 
177, 177, 182, 184, 178, 179, 178, 179, 179, 177, 178, 175, 184, 
179, 179, 176, 176, 175, 193, 183, 182, 179, 175, 178, 90, 189, 
177, 185, 179, 184, 171, 182, 181, 179, 179, 175, 104, 185, 180, 
180, 177, 176, 184, 181, 178, 184, 178, 182, 176, 176, 180, 180, 
176, 183, 182, 175, 184, 182, 182, 177, 182, 64, 178, 179, 179, 
176, 93, 131, 184, 181), event = structure(c(2L, 1L, 1L, 2L, 
1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("censored", 
"competing_event", "event_of_interest"), class = "factor"), group = 
structure(c(2L, 
2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 
1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 
2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L), .Label = c("G1", 
"G2"), class = "factor")), row.names = c(NA, -79L), class = c("tbl_df", 
"tbl", "data.frame"))
sdat
## there are no events of interest:
sdat$event %>% table(exclude = NULL)

## fit works fine
fit <- survfit(Surv(time, event) ~ group, data = sdat)
## the estimate is at 0 (as expected)
fit$pstate[,3]
## but there is no confidence estimate
fit$lower[,3]
fit$upper[,3]
## cmprsk does not event return estimates for the event of interest
library("cmprsk")
fit <- cuminc(ftime = sdat$time,
  fstatus = sdat$event,
  cencode = "censored",
  group = sdat$group)
fit



## note that I can get confidence intervals from a Kaplan-Meier
## estimator:

## censor competing events:
sdat <-
  sdat %>%
  mutate(event2 = recode(event,
 competing_event = "censored")) %>%
  mutate(event2 = as.numeric(event2) - 1)
sdat$event2 %>% table(exclude = NULL)
## fit using the survival package
fit <- survfit(Surv(time, event2) ~ 1, data = sdat)
fit
## no (meaningful?) confidence intervals
fit$lower
fit$upper
## but through the bpcp package this is possible
library("bpcp")
sfit2 <- bpcpfit(Surv(time, event2) ~ 1, data = sdat)
sfit2
tidykmciLR(sfit2)
## can further compare between groups:
bpcp2samp(sdat$time, sdat$event2, sdat$group, testtime = 150,
  parmtype = "difference")

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for p1-p2 and plot the CI with bar chart

[AbouEl-Makarim Aboueissa]

Hi Jim:

Thank you very much for your help in this topic.

with many thanks
abou
__


*AbouEl-Makarim Aboueissa, PhD*

*Professor, Statistics and Data Science*
*Graduate Coordinator*

*Department of Mathematics and Statistics*
*University of Southern Maine*



On Sat, Nov 13, 2021 at 8:47 PM Jim Lemon  wrote:

> Hi Abou,
> Perhaps this will be helpful. Be aware that you will cop some flak for
> putting error bars on a bar plot.
>
> aadat<-data.frame(group=c(rep("Exp",50),rep("Con",50)),
>  v1=sample(0:1,100,TRUE),
>  v2=sample(0:1,100,TRUE),
>  v3=sample(0:1,100,TRUE),
>  v4=sample(0:1,100,TRUE),
>  v5=sample(0:1,100,TRUE))
> ggps<-function(x,group) {
>  gns<-as.vector(table(group))
>  return(by(x,group,sum)/gns)
> }
> testggps<-data.frame(
>  group=c("A","A","A","B","B","B","B","C","C","C","C","C"),
>  x=c(1,0,1,1,0,1,0,1,1,0,0,0))
> aaprop<-sapply(aadat[,2:6],ggps,aadat[,1])
> library(plotrix)
> barpos<-barp(aaprop,ylim=c(0,0.65),col=c(2,3),names.arg=colnames(aaprop))
> legend(2.5,0.65,c("Con","Exp"),fill=c(2,3))
> dispersion(barpos$x,barpos$y,ulim=aaprop/10)
>
> Jim
>
> On Sun, Nov 14, 2021 at 11:01 AM AbouEl-Makarim Aboueissa
>  wrote:
> >
> > Dear All:
> >
> >
> >
> > I do have a binary data set with multiple variables, event = 1 in all
> > variables. As an example, I attached a data set with 6 variables. The
> first
> > column is the grouping variable. Then the next 5 columns are the binary
> > data for 5 variables.
> >
> >
> >
> > - Can we compute the confidence interval for the difference between the
> two
> > proportions of the event = 1 in both groups (say: G1 – G2) for the 5
> > variables in one shut.
> >
> >
> >
> > - I also need to create the Bar plot of individual proportions (both
> groups
> > side-by-side) and add the confidence intervals bar for the 5 variables in
> > one graph.
> >
> >
> >
> >
> >
> > Example.Data <- read.table(file="F/Example_Data_for_R.csv", header=T,
> > sep=",")
> >
> > Example.Data
> >
> > attach(Example.Data)
> >
> >
> >
> >
> >
> >
> >
> >  For example, this is how I use the prop.test() function to get the
> CI
> > for p1-p2
> >
> >
> >
> > x12 <- c(x1, x2)
> >
> > n12 <- c(n1, n2)
> >
> > prop.test(x12, n12, conf.level = 0.95)$conf.int
> >
> >
> >
> > But, I am not sure how to use it for raw data, and for multiple pairs of
> > data in one shut if possible.
> >
> >
> >
> > With many thanks in advance
> >
> > Abou
> > __
> >
> >
> > *AbouEl-Makarim Aboueissa, PhD*
> >
> > *Professor, Statistics and Data Science*
> > *Graduate Coordinator*
> >
> > *Department of Mathematics and Statistics*
> > *University of Southern Maine*
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for p1-p2 and plot the CI with bar chart

[Jim Lemon]

Hi Abou,
Perhaps this will be helpful. Be aware that you will cop some flak for
putting error bars on a bar plot.

aadat<-data.frame(group=c(rep("Exp",50),rep("Con",50)),
 v1=sample(0:1,100,TRUE),
 v2=sample(0:1,100,TRUE),
 v3=sample(0:1,100,TRUE),
 v4=sample(0:1,100,TRUE),
 v5=sample(0:1,100,TRUE))
ggps<-function(x,group) {
 gns<-as.vector(table(group))
 return(by(x,group,sum)/gns)
}
testggps<-data.frame(
 group=c("A","A","A","B","B","B","B","C","C","C","C","C"),
 x=c(1,0,1,1,0,1,0,1,1,0,0,0))
aaprop<-sapply(aadat[,2:6],ggps,aadat[,1])
library(plotrix)
barpos<-barp(aaprop,ylim=c(0,0.65),col=c(2,3),names.arg=colnames(aaprop))
legend(2.5,0.65,c("Con","Exp"),fill=c(2,3))
dispersion(barpos$x,barpos$y,ulim=aaprop/10)

Jim

On Sun, Nov 14, 2021 at 11:01 AM AbouEl-Makarim Aboueissa
 wrote:
>
> Dear All:
>
>
>
> I do have a binary data set with multiple variables, event = 1 in all
> variables. As an example, I attached a data set with 6 variables. The first
> column is the grouping variable. Then the next 5 columns are the binary
> data for 5 variables.
>
>
>
> - Can we compute the confidence interval for the difference between the two
> proportions of the event = 1 in both groups (say: G1 – G2) for the 5
> variables in one shut.
>
>
>
> - I also need to create the Bar plot of individual proportions (both groups
> side-by-side) and add the confidence intervals bar for the 5 variables in
> one graph.
>
>
>
>
>
> Example.Data <- read.table(file="F/Example_Data_for_R.csv", header=T,
> sep=",")
>
> Example.Data
>
> attach(Example.Data)
>
>
>
>
>
>
>
>  For example, this is how I use the prop.test() function to get the CI
> for p1-p2
>
>
>
> x12 <- c(x1, x2)
>
> n12 <- c(n1, n2)
>
> prop.test(x12, n12, conf.level = 0.95)$conf.int
>
>
>
> But, I am not sure how to use it for raw data, and for multiple pairs of
> data in one shut if possible.
>
>
>
> With many thanks in advance
>
> Abou
> __
>
>
> *AbouEl-Makarim Aboueissa, PhD*
>
> *Professor, Statistics and Data Science*
> *Graduate Coordinator*
>
> *Department of Mathematics and Statistics*
> *University of Southern Maine*
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for p1-p2 and plot the CI with bar chart

[AbouEl-Makarim Aboueissa]

Dear All:



I do have a binary data set with multiple variables, event = 1 in all
variables. As an example, I attached a data set with 6 variables. The first
column is the grouping variable. Then the next 5 columns are the binary
data for 5 variables.



- Can we compute the confidence interval for the difference between the two
proportions of the event = 1 in both groups (say: G1 – G2) for the 5
variables in one shut.



- I also need to create the Bar plot of individual proportions (both groups
side-by-side) and add the confidence intervals bar for the 5 variables in
one graph.





Example.Data <- read.table(file="F/Example_Data_for_R.csv", header=T,
sep=",")

Example.Data

attach(Example.Data)







 For example, this is how I use the prop.test() function to get the CI
for p1-p2



x12 <- c(x1, x2)

n12 <- c(n1, n2)

prop.test(x12, n12, conf.level = 0.95)$conf.int



But, I am not sure how to use it for raw data, and for multiple pairs of
data in one shut if possible.



With many thanks in advance

Abou
__


*AbouEl-Makarim Aboueissa, PhD*

*Professor, Statistics and Data Science*
*Graduate Coordinator*

*Department of Mathematics and Statistics*
*University of Southern Maine*
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] confidence interval

[David Winsemius]

On 4/26/20 11:26 AM, Alex Serafim wrote:

there is a function called "confbands", which is no longer available in
software R. It is not "confband" or "confBands", but "confbands", I need to
use this object to add confidence intervals in my work. Why is this
function not available?
attached is an image that shows how the graph should look, in another image
is the error that appears in the software. "confbands" object not
found. Also attached is an image showing the function functioning normally.
In which even the functions of the object appear.



There was no attached image or images. That is because you have not 
carefully read the information provided here: 
https://stat.ethz.ch/mailman/listinfo/r-help and in the Posting Guide.


R is composed of packages. You need to know which package the 
"confbands" function came from before we can identify the cause of its 
absence.


It could be that you have not installed or loaded the package into your 
session. Or it could be that the author of the package in which it used 
to be made available has not maintained it on CRAN. Or perhaps it was 
never on CRAN and was instead on some blog or class assignment where it 
was defined upstream and you failed to notice it. All these situations 
have been seen on Rhelp posted by persons claiming that a function has 
gone missing.



--

David.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] confidence interval

[Alex Serafim]

there is a function called "confbands", which is no longer available in
software R. It is not "confband" or "confBands", but "confbands", I need to
use this object to add confidence intervals in my work. Why is this
function not available?
attached is an image that shows how the graph should look, in another image
is the error that appears in the software. "confbands" object not
found. Also attached is an image showing the function functioning normally.
In which even the functions of the object appear.
-- 



Livre
de vírus. www.avast.com
.
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence interval on quantile regression (quantreg)

[Cade, Brian]

Depending on the procedure used for estimating the CI, especially if the
default rankscore inversion method, then it is possible that legitimate end
points of the intervals for some quantiles with a given alpha (e.g., 0.05
for 95% CI) cannot be refined beyond plus or minus infinity.  Of course,
this typically happens for smaller sample sizes, more extreme taus, and
more complex models.  But unusual distributional characteristics of the
data distributions can also contribute to this issue.

Brian

Brian S. Cade, PhD

U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO  80526-8818

email:  ca...@usgs.gov 
tel:  970 226-9326


On Wed, Oct 19, 2016 at 11:05 AM, Frank Black 
wrote:

> Hi all,
>
> I am using the quantile regression package for estimating some models.
> However, in some cases the intervals' upper bounds are either abnormally
> high or low, with values such as -1.797693e+308 or 1.797693e+308. Actually,
> the number is the same in absolute terms.
>
> Does anyone know a reasonable explanation for this?
>
> Thanks.
>
> Kind regards,
> Frank
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence interval on quantile regression (quantreg)

[Frank Black]

Hi all,

I am using the quantile regression package for estimating some models.
However, in some cases the intervals' upper bounds are either abnormally
high or low, with values such as -1.797693e+308 or 1.797693e+308. Actually,
the number is the same in absolute terms.

Does anyone know a reasonable explanation for this?

Thanks.

Kind regards,
Frank

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for Ipred Bagging outputs ?

[David Winsemius]

> On Mar 19, 2016, at 12:36 AM, Majid Javanmard  
> wrote:
> 
> Hello everyone
> 
> here is the code that implements bagging using ipred package in R :
> 
> library(ipred)
> library(mlbench)
> data("BostonHousing")
> # Test set error (nbagg=25, trees pruned): 3.41 (Breiman, 1996a, Table 8)
> mod <- bagging(medv ~ ., data=BostonHousing, coob=TRUE)
> print(mod)
> pred <- predict(mod)
> pred<- as.data.frame(pred)
> 
> How can I have 95% Confidence interval for each predicted values !?

Perhaps you really mean prediction intervals, since none of those results 
really a parameters. I don't think it makes sense to talk about 95%CI's in the 
context of a bagging procedure because there really is no single model. In any 
case it has already been suggested that this is not really an R coding problem 
but rather a conceptual problem. You were advised to post further questions on 
stats.stackexchange.com (if you were unable to find an answered question), the 
first hit on a Google search for "confidence Interval randomForest":

 
http://stats.stackexchange.com/questions/56895/do-the-predictions-of-a-random-forest-model-have-a-prediction-interval

-- 

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for Ipred Bagging outputs ?

[Majid Javanmard]

Hello everyone

here is the code that implements bagging using ipred package in R :

library(ipred)
library(mlbench)
data("BostonHousing")
# Test set error (nbagg=25, trees pruned): 3.41 (Breiman, 1996a, Table 8)
mod <- bagging(medv ~ ., data=BostonHousing, coob=TRUE)
print(mod)
pred <- predict(mod)
pred<- as.data.frame(pred)

How can I have 95% Confidence interval for each predicted values !?

I appreciate if someone help me

Thanks

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for Ipred Bagging outputs (please help)

[Majid Javanmard]

Hello everyone

here is the code that implements bagging using ipred package :

library(ipred)
library(mlbench)
data("BostonHousing")
# Test set error (nbagg=25, trees pruned): 3.41 (Breiman, 1996a, Table 8)
mod <- bagging(medv ~ ., data=BostonHousing, coob=TRUE)
print(mod)
pred <- predict(mod)
pred<- as.data.frame(pred)

How can I have 95% Confidence interval for each predicted values !?

I appreciate if someone help me

Thanks

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for R-squared

[R Martinez]

To the list,

Does R have a function for computing the confidence interval of R-squared 
(coefficient of determination)?

I checked R's Help function with “confidence interval” as the query but none of 
the packages and functions that R found pertain to R-squared.

Thanks in advance for your help.

Raul Martinez
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for R-squared

[Sarah Goslee]

Hi Raul,

Searching for "confidence interval of R-squared" on rseek.org turns up
some packages that might be of use, including bootstrap and MBESS.



On Tue, Jan 26, 2016 at 3:11 PM, R Martinez  wrote:
> To the list,
>
> Does R have a function for computing the confidence interval of R-squared 
> (coefficient of determination)?
>
> I checked R's Help function with “confidence interval” as the query but none 
> of the packages and functions that R found pertain to R-squared.
>
> Thanks in advance for your help.
>
> Raul Martinez

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for R-squared

[Nordlund, Dan (DSHS/RDA)]

A Google search suggested the use of the boot package to bootstrap a confidence 
interval for R-squared.

http://www.statmethods.net/advstats/bootstrapping.html

Dan

Daniel Nordlund, PhD
Research and Data Analysis Division
Services & Enterprise Support Administration
Washington State Department of Social and Health Services

> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of R Martinez
> Sent: Tuesday, January 26, 2016 12:12 PM
> To: r-help@r-project.org
> Subject: [R] Confidence Interval for R-squared
> 
> To the list,
> 
> Does R have a function for computing the confidence interval of R-squared
> (coefficient of determination)?
> 
> I checked R's Help function with “confidence interval” as the query but none
> of the packages and functions that R found pertain to R-squared.
> 
> Thanks in advance for your help.
> 
> Raul Martinez
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence Interval for R-squared

[David L Carlson]

Sometimes Google just works better. My search on "confidence interval R2 R" 
turned up several hits on the first two pages. 

A Quick-R page on bootstrapping with an example using R-squared
http://www.statmethods.net/advstats/bootstrapping.html

Function ci.R2() in package MBESS
http://www.inside-r.org/packages/cran/MBESS/docs/ci.R2

Functions CI.Rsq() and CI.Rsqlm() in package psychometric
http://finzi.psych.upenn.edu/library/psychometric/html/CI.Rsq.html


-
David L Carlson
Department of Anthropology
Texas A University
College Station, TX 77840-4352


-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of R Martinez
Sent: Tuesday, January 26, 2016 2:12 PM
To: r-help@r-project.org
Subject: [R] Confidence Interval for R-squared

To the list,

Does R have a function for computing the confidence interval of R-squared 
(coefficient of determination)?

I checked R's Help function with “confidence interval” as the query but none of 
the packages and functions that R found pertain to R-squared.

Thanks in advance for your help.

Raul Martinez
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence interval for the mean in the random intercept and random slople model

[li li]

Hi all,
  I am fitting a random slope and random intercept model usign lme
fucntion as shown below. Type is factor with two levels. I would like
to to find a confidence interval for mean of this model. Note that the
variance we use in finding the confidence interval should include the
variariance component from random effect.  Any suggestions?


## using lme function
 mod_lme - lme(ti  ~ type*months, random=~ 1+months|lot, na.action=na.omit,
+ data=one, control = lmeControl(opt = optim))
 summary(mod_lme)
Linear mixed-effects model fit by REML
 Data: one
AIC   BIC   logLik
  -82.60042 -70.15763 49.30021

Random effects:
 Formula: ~1 + months | lot
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr
(Intercept) 8.907584e-03 (Intr)
months  6.039781e-05 -0.096
Residual4.471243e-02

Fixed effects: ti ~ type * months
 Value   Std.Error DF   t-value p-value
(Intercept) 0.25831245 0.016891587 31 15.292373  0.
type0.13502089 0.026676101  4  5.061493  0.0072
months  0.00804790 0.001218941 31  6.602368  0.
type:months -0.00693679 0.002981859 31 -2.326329  0.0267
 Correlation:
   (Intr) typPPQ months
type   -0.633
months -0.785  0.497
type:months  0.321 -0.762 -0.409

Standardized Within-Group Residuals:
  MinQ1   MedQ3   Max
-2.162856e+00 -1.962972e-01 -2.771184e-05  3.749035e-01  2.088392e+00

Number of Observations: 39
Number of Groups: 6

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] confidence interval for wilcox_test

[RomanGelzhaeuser]

Hello,
When using wilcox_test from coin I got a sample estimator which does not lie
in the confidence interval.(see below) How is that possible? (The
documentation said that both referre to some kind of pseudo median but it
seems to be the same for both the confidence interval and the samplle
estimator)

 wilcox_test(ErgebnisBefragung~V,daten,conf.int=T)

Asymptotic Wilcoxon Mann-Whitney Rank Sum Test

data:  ErgebnisBefragung by V (Burch, TVT)
Z = -2.0157, p-value = 0.04383
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 -2.056811e-05 -4.179582e-06
sample estimates:
difference in location 
 -3.251934e-05 

 



--
View this message in context: 
http://r.789695.n4.nabble.com/confidence-interval-for-wilcox-test-tp4702950.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Confidence interval, multiple imputation

[Mª Teresa Martinez Soriano]

Hi to
everyone, I have a big data set where rows are observations and columns are
variables. It contains a lot of missing values. I have used multiple imputation
with library mice and I get an “exact” prediction of each missing value. Now, I
would like to know the error I can commit or the confidence interval.

How can I
get this?

This is
part of my code

library(mice)

mod1-mice(dat,
method=c(,,rep(pmm,6)))

ro-round(cor(dat,
use = pair), 3)

 

predictor-quickpred(dat)#
esta matriz predictora se construye según las correlaciones



mod1-mice(dat,method=c(,,rep(pmm,6)),
pred=predictor)

imputados-complete(mod1,'long')

x.imp=split(imputados,
imputados$.imp)

acumula=x.imp[[1]][,-c(1,2)]

   for(j
in 2:length(x.imp))

   {
acumula=acumula+x.imp[[j]][,-c(1,2)]}

med.imp=acumula/5

 

 

Thanks in
advance 
  
[[alternative HTML version deleted]]

__
R-help@r-project.org 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] confidence interval in rpart

[puja makkar]

Hi,


I am using rpart with method=anova.Can we compute 95% confidence
intervals and prediction interval for the predicted mean at each node?

Thanks,
Puja

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval Calculation

[Tasnuva Tabassum]

Hello,
I have the following r-codes for solving a quasilikelihood estimating
equation:

library(geepack)
fit-geese(y~x1+x2+x3,jack=TRUE,scale.fix=TRUE,data=dat,mean.link =
logit, corstr=independence)


Now my question is how can I calculate the confidence interval of the
parameters of the above model fit?

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] confidence interval calculation for gee

[Tasnuva Tabassum]

Hello,
I have the following r-codes for solving a quasilikelihood estimating
equation:

library(geepack)
fit-geese(y~x1+x2+x3,jack=TRUE,id=id,scale.fix=TRUE,data=dat,mean.link =
logit, corstr=independence)


Now my question is how can I calculate the confidence interval of the
parameters of the above model fit?

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] confidence interval calculation for gee

[Duncan Mackay]

At 04:26 8/04/2013, you wrote:

Hello,
I have the following r-codes for solving a quasilikelihood estimating
equation:

library(geepack)
fit-geese(y~x1+x2+x3,jack=TRUE,id=id,scale.fix=TRUE,data=dat,mean.link =
logit, corstr=independence)


Now my question is how can I calculate the confidence interval of the
parameters of the above model fit?

[[alternative HTML version deleted]]

__
R-help@r-project.org 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.


It pays to look at the help page. usually with a regression there is 
a summary method


From ?geese and the index page of geepack there is summary.geese


names(summary(m2))
[1] meancorrelation 
scale   callmodel   control error 
clusz

 summary(m2)[[1]]
  estimatesan.se   wald  p
(Intercept)  1.3476092 0.1573571 73.3423807 0.
x0.1118360 0.1159304  0.9306116 0.33470405
trt -0.1068224 0.1936977  0.3041417 0.58129752
x:trt   -0.3023841 0.1710601  3.1247883 0.07710989

geepack uses sandwich standard errors (san.se)

summary(m2t)[[2]]
   estimate san.se waldp
alpha 0.5978407 0.08107156 54.37934 1.653122e-13

HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mac...@northnet.com.au

__
R-help@r-project.org 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] confidence interval for survfit

[Terry Therneau]

The first thing you are missing is the documentation -- try ?survfit.object.
 fit - survfit(Surv(time,status)~1,data)
fit$std.err will contain the standard error of the cumulative hazard or 
-log(survival)

The standard error of the survival curve is approximately S(t) * std(hazard), by the delta 
method.  This is what is printed by the summary function, because it is what user's 
expect, but it has very poor performance for computing confidence intervals.  A much 
better one is exp(-1* confidence interval for the cumulative hazard), which is the 
default.  In fact there are lots of better ones whose relative ranking depends on the 
details of your simulation study.  About the only really consistent result is that 
anything thoughtful beats S(t) +- 1.96 se(S), easily.  The default in R is the one that 
was best in the most recent paper I had read at the time I set the default.  If I were to 
rank them today using an average over all the comparison papers it would be second or 
third, but the good methods are so close that in practical terms it hardly matters.


Terry Therneau

On 03/15/2013 06:00 AM, r-help-requ...@r-project.org wrote:

Hi, I am wondering how the confidence interval for Kaplan-Meier estimator is 
calculated by survfit(). For example,?


  summary(survfit(Surv(time,status)~1,data),times=10)

Call: survfit(formula = Surv(rtime10, rstat10) ~ 1, data = mgi)

?time n.risk n.event survival std.err lower 95% CI upper 95% CI
?? 10 ?? 168? 55??? 0.761? 0.0282??? 0.707??? 0.818


I am trying to reproduce the upper and lower CI by using standard error. As far 
I understand, the default method for survfit() to calculate confidence interval 
is on the log survival scale, so:


upper CI = exp(log(0.761)+qnorm(0.975)*0.0282) = 0.804
lower CI = exp(log(0.761)-qnorm(0.975)*0.0282) = 0.720


they are not the same as the output from survfit().

Am I missing something?

Thanks

John


__
R-help@r-project.org 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] Confidence interval for Support Vector Regression

[João Oliveirinha]

Have you found a solution for this? I am also trying to find a way to 
retrieve the confidence intervals for the predictions.

Best regards,
João Oliveirinha

On Sunday, October 30, 2011 7:04:03 PM UTC, Muhammed Akbulut wrote:

 Hi,

 Is it possible to calculate confidence intervals for support vector
 regression?
 In the ksvm{kernlab} manual, it says that it supports confidence
 intervals for regression.
 However, i couldn't find any information about how to calculate confidence
 interval.
 Do you have any documents or examples about this or any other suggestions?

 Thanks in advance


 --

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org 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@r-project.org 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] Confidence interval for Whittle method

[Barun Saha]

Hello,

How could we get the confidence interval when using the whittleFit
method from fArma package?

--
Thanks,
Barun Saha

__
R-help@r-project.org 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] Confidence Interval from Moments?

[lambdatau]

Hi all,

I'm wondering whether it is possible to construct a confidence interval
using only the mean, variance, skewness and kurtosis, i.e. without any of
the population?

If anyone could help with this it'd be much appreciated (even if just a
confirmation of it being impossible!).

Thanks.

--
View this message in context: 
http://r.789695.n4.nabble.com/Confidence-Interval-from-Moments-tp4284937p4284937.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence Interval from Moments?

[R. Michael Weylandt]

Assuming a distribution defined solely by those moments it is possible
(e.g., z- or t-test confidence intervals) but this isn't really the
place to discuss such things since there's no R content to your
question: try stats.stackexchange.com

Michael

On Wed, Jan 11, 2012 at 4:56 AM, lambdatau lewis@rbs.com wrote:
 Hi all,

 I'm wondering whether it is possible to construct a confidence interval
 using only the mean, variance, skewness and kurtosis, i.e. without any of
 the population?

 If anyone could help with this it'd be much appreciated (even if just a
 confirmation of it being impossible!).

 Thanks.

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Confidence-Interval-from-Moments-tp4284937p4284937.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org 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@r-project.org 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] Confidence interval for Support Vector Regression

[Muhammed Akbulut]

Hi,

Is it possible to calculate confidence intervals for support vector
regression?
In the ksvm{kernlab} manual, it says that it supports confidence
intervals for regression.
However, i couldn't find any information about how to calculate confidence
interval.
Do you have any documents or examples about this or any other suggestions?

Thanks in advance






--

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval or p-value for difference in two c-statistics

[Bonnett, Laura]

Dear All,

Apologies if you have a seen a question like this from me before.  I am hoping 
that if I re-word my question more carefully someone may be able to offer more 
help than the last time I asked something similar.  I am using R 2.9.2 and 
Windows XP.

I am trying to determine if there is a statistically significant difference 
between two c-statistics (or equivalently D statistics).  In Stata I gather 
that this
can be done using the lincom command but I am not aware of anything similar in 
R.

I am doing a simulation study and have simulated two datasets (independent).  I 
can obtain c-statistics for each variable in each dataset using this code:

 rdev - rcorrcens(Surv(stimes1,eind1)~gendat1+neurodat1)
 rdev

Somers' Rank Correlation for Censored DataResponse
variable:Surv(stimes1, eind1)

  CDxy  aDxySDZ  Pn
gendat1   0.534  0.069 0.069 0.017 3.98 0.0001 1500
neurodat1 0.482 -0.036 0.036 0.011 3.18 0.0015 1500

 rval - rcorrcens(Surv(stimes2,eind2)~gendat2+neurodat2)
 rval

Somers' Rank Correlation for Censored DataResponse
variable:Surv(stimes2, eind2)

  CDxy  aDxySDZ Pn
gendat2   0.543  0.085 0.085 0.017 4.94 0e+00 1500
neurodat2 0.481 -0.038 0.038 0.011 3.44 6e-04 1500

Is there a way to determine if the c-statistics are significantly 
different/similar to one another across data sets (i.e. I want to compare the 
c-statistic for gendat1 with that for gendat2).  I was considering using the 
'tost' function within the equivalence package but that doesn't seem 
appropriate either.

Many thanks,
Laura



[[alternative HTML version deleted]]

__
R-help@r-project.org 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] confidence interval as shaded band (lme)

[bjmjarrett]

Hi all,

I’m trying to plot confidence intervals for the fitted values I get with my
lme model in R. 

Is there any way I can plot this in the form of a shaded band, like the
output of geom_smooth() in ggplot2 package. ggplot2 seems to use only lm,
glm, gam, loess and rlm as smoothing methods.

Any advice on the functions I should use to accomplish this will be very
helpful.

Thank you very much.

Ben

--
View this message in context: 
http://r.789695.n4.nabble.com/confidence-interval-as-shaded-band-lme-tp3727645p3727645.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] confidence interval as shaded band (lme)

[R. Michael Weylandt]

I don't know lme models very well, but if you have standard errors for your
values, this shouldn't be too hard (as a last resort) using polygon()

For example

x = 1:10
y = x^2
y.Err = 2*x
y.Up = y + y.Err; y.Dn =y-y.Err

# This graph is actually quite ugly so don't copy the formatting
plot(x,y,type=n)
polygon(c(x,rev(x)),c(y.Up,rev(y.Dn)),col=grey,border=red)
# Use the rev commands so the border moves logically around the shaded area
lines(x,y,type=b,lwd=3) # put the means back on top of the polygon

Still, this is a little brute force and I'd imagine that someone else will
shortly let you know how R can already do automatically.

Michael Weylandt

On Mon, Aug 8, 2011 at 1:07 PM, bjmjarrett bjmjarr...@gmail.com wrote:

 Hi all,

 I’m trying to plot confidence intervals for the fitted values I get with my
 lme model in R.

 Is there any way I can plot this in the form of a shaded band, like the
 output of geom_smooth() in ggplot2 package. ggplot2 seems to use only lm,
 glm, gam, loess and rlm as smoothing methods.

 Any advice on the functions I should use to accomplish this will be very
 helpful.

 Thank you very much.

 Ben

 --
 View this message in context:
 http://r.789695.n4.nabble.com/confidence-interval-as-shaded-band-lme-tp3727645p3727645.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org 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.


[[alternative HTML version deleted]]

__
R-help@r-project.org 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] confidence interval as shaded band (lme)

[Dennis Murphy]

Hi:

On Mon, Aug 8, 2011 at 10:07 AM, bjmjarrett bjmjarr...@gmail.com wrote:
 Hi all,

 I’m trying to plot confidence intervals for the fitted values I get with my
 lme model in R.

Which fitted values? The ones conditional on the random effects or the
ones averaged over the random effects? The standard errors of the two
sets of predictions are not the same.

 Is there any way I can plot this in the form of a shaded band, like the
 output of geom_smooth() in ggplot2 package. ggplot2 seems to use only lm,
 glm, gam, loess and rlm as smoothing methods.

You can fake it with geom_ribbon(), but it would be convenient to have
the endpoints of the CIs in a data frame you could input into ggplot2.

HTH,
Dennis

 Any advice on the functions I should use to accomplish this will be very
 helpful.

 Thank you very much.

 Ben

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/confidence-interval-as-shaded-band-lme-tp3727645p3727645.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org 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@r-project.org 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] confidence interval as shaded band (lme)

[Bert Gunter]

Would someone tell me how they propose to go from standard errors to
confidence intervals*. I suspect Doug Bates would probably like to
know, also, as he has expended a lot of effort on this over the years,
I believe.  :-)

-- Bert

* Note   +/- 2 std errors is almost certainly not the right answer --
how far wrong it is ... that's the question.


Cheers,
Bert

On Mon, Aug 8, 2011 at 1:07 PM, Dennis Murphy djmu...@gmail.com wrote:
 Hi:

 On Mon, Aug 8, 2011 at 10:07 AM, bjmjarrett bjmjarr...@gmail.com wrote:
 Hi all,

 I’m trying to plot confidence intervals for the fitted values I get with my
 lme model in R.

 Which fitted values? The ones conditional on the random effects or the
 ones averaged over the random effects? The standard errors of the two
 sets of predictions are not the same.

 Is there any way I can plot this in the form of a shaded band, like the
 output of geom_smooth() in ggplot2 package. ggplot2 seems to use only lm,
 glm, gam, loess and rlm as smoothing methods.

 You can fake it with geom_ribbon(), but it would be convenient to have
 the endpoints of the CIs in a data frame you could input into ggplot2.

 HTH,
 Dennis

 Any advice on the functions I should use to accomplish this will be very
 helpful.

 Thank you very much.

 Ben

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/confidence-interval-as-shaded-band-lme-tp3727645p3727645.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org 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@r-project.org 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.




-- 
Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions.

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics

__
R-help@r-project.org 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] Confidence interval from resampling

[Rolf Turner]

On 24/06/11 01:44, Adriana Bejarano wrote:

Dear R gurus,

I have the following code, but I still not know how to estimate and extract
confidence intervals (95%CI) from resampling.

SNIP

Some sound advice --- still sound after all these years I believe --- in 
respect

of  boot strap confidence intervals, is to be found in

Theoretical comparison of bootstrap confidence intervals,
Peter Hall, Annals of Statistics vol. 16 no. 3 1988, pp. 927 -- 953.


cheers,

Rolf Turner

__
R-help@r-project.org 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] Confidence interval from resampling

[Adriana Bejarano]

Dear R gurus,

I have the following code, but I still not know how to estimate and extract
confidence intervals (95%CI) from resampling.

Thanks!

~Adriana

#data
penta-c(770,729,640,486,450,410,400,340,306,283,278,260,253,242,240,229,201,198,190,186,180,170,168,151,150,148,147,125,117,110,107,104,85,83,80,74,70,66,54,46,45,43,40,38,10)
x-log(penta+1)
plot(ecdf(x), ylab=Probability, xlab=Concentration (Ln+1))

x.wei-fitdistr(x,weibull)
x.wei
 shapescale
  6.7291685   5.3769965
 (0.7807718) (0.1254696)
xwei.shape - x.wei$estimate[[1]]
xwei.scale -  x.wei$estimate[[2]]

x.wei-fitdistr(x,weibull)
x.wei
xwei.shape - x.wei$estimate[[1]]
xwei.scale -  x.wei$estimate[[2]]
curve(pweibull(x, shape=xwei.shape, scale = xwei.scale,lower.tail=TRUE,
log.p=FALSE), add=TRUE,col='green',lwd=3)

#draw random numbers from a weibull distribution 100 times with
shape=xwei.shape, scale = xwei.scale
draw - lapply(1:100, function(.x){
out-rweibull(x, shape=xwei.shape, scale = xwei.scale)
})
newx- data.frame(draw)

colnames(newx)-paste(x, 1:100, sep = )
newmat-data.matrix(newx)

# matrix of coefficients
rownum=2
colnum=100
ResultMat-matrix(NA, ncol=colnum, nrow=rownum)
rownum2=45
colnum2=100
ResultMat2-matrix(NA, ncol=colnum2, nrow=rownum2)

#loop through each column in the source matrix
for (i in 1:100)
{
sel_col-newmat[col(newmat)==i]
  {ResultMat[,i]-coef(fitdistr(sel_col,weibull))}
 xwei.shape- ResultMat[1,i]
   xwei.scale- ResultMat[2,i]
 curve(pweibull(x, shape=xwei.shape, scale=xwei.scale, lower.tail=TRUE,
log.p = FALSE), add=TRUE, col='grey',lwd=0.5)
 ResultMat2[,i]-pweibull(x, shape=xwei.shape, scale =
xwei.scale,lower.tail=TRUE, log.p=FALSE)
}

## convert dataframe to numeric
MatOut- as.matrix(ResultMat2)
mode(MatOut) - numeric

# initiate variables
mm-ml-mu-rep(0,length(MatOut[,1]))

# mean and upper/lower quantiles
for(i in 1:length(MatOut[,1])){
 mm[i]- mean(MatOut[i,])
 ml[i]- quantile(MatOut[i,], 0.025, na.rm=TRUE)
 mu[i]- quantile(MatOut[i,], 0.975, na.rm=TRUE)
}
#lines(x, mm, col=black)
lines(x, ml, col=blue, lwd=2)
lines(x, mu, col=blue, lwd=2)

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval from resampling

[David Winsemius]

On Jun 23, 2011, at 9:44 AM, Adriana Bejarano wrote:


Dear R gurus,

I have the following code, but I still not know how to estimate and  
extract

confidence intervals (95%CI) from resampling.



If you have a distribution of values, say resamp.stat, of a  
statistic from a properly performed resampling operation you can  
extract and display easily the 5th and 95th percentiles.


CI.stat - quantile(resamp.stat, c(0.05, 0.95) )
CI.stat

Note: I do not think that 100 replications would generally be  
sufficient for final work, although its probably acceptable for testing.


Your code as posted initially threw a bunch of errors since you did  
not include library(MASS), but fixing that fairly obvious problem  
shows that you can draw a nice plot. However, it remains unclear what  
statistic of what distribution you desire to assess. Mean, median, ...  
of what?


I do not think the error that arose on my machine from the wrapped  
text here:


#draw random numbers from a weibull distribution 100 times with
 ... shape=xwei.shape, scale = xwei.scale   - error

. was causing any problem, but there were a bunch of warnings that  
ought to be investigated:


Right after the loop I see ten of these:
Warning messages:
1: In dweibull(x, shape, scale, log) : NaNs produced

--
David


Thanks!

~Adriana

#data
penta- 
c 
(770,729,640,486,450,410,400,340,306,283,278,260,253,242,240,229,201,198,190,186,180,170,168,151,150,148,147,125,117,110,107,104,85,83,80,74,70,66,54,46,45,43,40,38,10 
)

x-log(penta+1)
plot(ecdf(x), ylab=Probability, xlab=Concentration (Ln+1))

x.wei-fitdistr(x,weibull)
x.wei
shapescale
 6.7291685   5.3769965
(0.7807718) (0.1254696)
xwei.shape - x.wei$estimate[[1]]
xwei.scale -  x.wei$estimate[[2]]

x.wei-fitdistr(x,weibull)
x.wei
xwei.shape - x.wei$estimate[[1]]
xwei.scale -  x.wei$estimate[[2]]
curve(pweibull(x, shape=xwei.shape, scale =  
xwei.scale,lower.tail=TRUE,

log.p=FALSE), add=TRUE,col='green',lwd=3)

#draw random numbers from a weibull distribution 100 times with
shape=xwei.shape, scale = xwei.scale
draw - lapply(1:100, function(.x){
out-rweibull(x, shape=xwei.shape, scale = xwei.scale)
})
newx- data.frame(draw)

colnames(newx)-paste(x, 1:100, sep = )
newmat-data.matrix(newx)

# matrix of coefficients
rownum=2
colnum=100
ResultMat-matrix(NA, ncol=colnum, nrow=rownum)
rownum2=45
colnum2=100
ResultMat2-matrix(NA, ncol=colnum2, nrow=rownum2)

#loop through each column in the source matrix
for (i in 1:100)
   {
   sel_col-newmat[col(newmat)==i]
 {ResultMat[,i]-coef(fitdistr(sel_col,weibull))}
xwei.shape- ResultMat[1,i]
  xwei.scale- ResultMat[2,i]
curve(pweibull(x, shape=xwei.shape, scale=xwei.scale, lower.tail=TRUE,
log.p = FALSE), add=TRUE, col='grey',lwd=0.5)
ResultMat2[,i]-pweibull(x, shape=xwei.shape, scale =
xwei.scale,lower.tail=TRUE, log.p=FALSE)
}

## convert dataframe to numeric
MatOut- as.matrix(ResultMat2)
mode(MatOut) - numeric

# initiate variables
mm-ml-mu-rep(0,length(MatOut[,1]))

# mean and upper/lower quantiles
for(i in 1:length(MatOut[,1])){
mm[i]- mean(MatOut[i,])
ml[i]- quantile(MatOut[i,], 0.025, na.rm=TRUE)
mu[i]- quantile(MatOut[i,], 0.975, na.rm=TRUE)
}
#lines(x, mm, col=black)
lines(x, ml, col=blue, lwd=2)
lines(x, mu, col=blue, lwd=2)

[[alternative HTML version deleted]]

__
R-help@r-project.org 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.


David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence interval from resampling

[Stephen Ellison]

Depending on how critical the problem is, you might also want to look at the 
literature on bootstrap CI's, perhaps starting from the references in boot.ci 
in the boot package. The simple quantiles are not necessarily the most 
appropriate. For example I seem to recall that BCa intervals were the intervals 
recommended for 'general use'  by Efron and Tibshirani (Introduction to the 
Bootstrap (1993) Chapman and Hall) with ABC intervals also getting an 
honourable mention. 1993 is a long time ago, though...

S Ellison



 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] On Behalf Of David Winsemius
 Sent: 23 June 2011 15:25
 To: Adriana Bejarano
 Cc: r-help@r-project.org
 Subject: Re: [R] Confidence interval from resampling
 
 
 On Jun 23, 2011, at 9:44 AM, Adriana Bejarano wrote:
 
  Dear R gurus,
 
  I have the following code, but I still not know how to estimate and 
  extract confidence intervals (95%CI) from resampling.
 
 
 If you have a distribution of values, say resamp.stat, of a 
 statistic from a properly performed resampling operation you 
 can extract and display easily the 5th and 95th percentiles.
 
 CI.stat - quantile(resamp.stat, c(0.05, 0.95) ) CI.stat
 
 Note: I do not think that 100 replications would generally be 
 sufficient for final work, although its probably acceptable 
 for testing.
 
 Your code as posted initially threw a bunch of errors since 
 you did not include library(MASS), but fixing that fairly 
 obvious problem shows that you can draw a nice plot. However, 
 it remains unclear what statistic of what distribution you 
 desire to assess. Mean, median, ...  
 of what?
 
 I do not think the error that arose on my machine from the 
 wrapped text here:
 
 #draw random numbers from a weibull distribution 100 times with
   ... shape=xwei.shape, scale = xwei.scale   - error
 
 . was causing any problem, but there were a bunch of 
 warnings that ought to be investigated:
 
 Right after the loop I see ten of these:
 Warning messages:
 1: In dweibull(x, shape, scale, log) : NaNs produced
 
 --
 David
 
  Thanks!
 
  ~Adriana
 
  #data
  penta-
  c
  
 (770,729,640,486,450,410,400,340,306,283,278,260,253,242,240,229,201,1
  
 98,190,186,180,170,168,151,150,148,147,125,117,110,107,104,85,83,80,74
  ,70,66,54,46,45,43,40,38,10
  )
  x-log(penta+1)
  plot(ecdf(x), ylab=Probability, xlab=Concentration (Ln+1))
 
  x.wei-fitdistr(x,weibull)
  x.wei
  shapescale
   6.7291685   5.3769965
  (0.7807718) (0.1254696)
  xwei.shape - x.wei$estimate[[1]]
  xwei.scale -  x.wei$estimate[[2]]
 
  x.wei-fitdistr(x,weibull)
  x.wei
  xwei.shape - x.wei$estimate[[1]]
  xwei.scale -  x.wei$estimate[[2]]
  curve(pweibull(x, shape=xwei.shape, scale = 
  xwei.scale,lower.tail=TRUE, log.p=FALSE), 
 add=TRUE,col='green',lwd=3)
 
  #draw random numbers from a weibull distribution 100 times with 
  shape=xwei.shape, scale = xwei.scale draw - lapply(1:100, 
  function(.x){ out-rweibull(x, shape=xwei.shape, scale = xwei.scale)
  })
  newx- data.frame(draw)
 
  colnames(newx)-paste(x, 1:100, sep = )
  newmat-data.matrix(newx)
 
  # matrix of coefficients
  rownum=2
  colnum=100
  ResultMat-matrix(NA, ncol=colnum, nrow=rownum)
  rownum2=45
  colnum2=100
  ResultMat2-matrix(NA, ncol=colnum2, nrow=rownum2)
 
  #loop through each column in the source matrix for (i in 1:100)
 {
 sel_col-newmat[col(newmat)==i]  
  {ResultMat[,i]-coef(fitdistr(sel_col,weibull))}
  xwei.shape- ResultMat[1,i]
xwei.scale- ResultMat[2,i]
  curve(pweibull(x, shape=xwei.shape, scale=xwei.scale, 
 lower.tail=TRUE, 
  log.p = FALSE), add=TRUE, col='grey',lwd=0.5) 
  ResultMat2[,i]-pweibull(x, shape=xwei.shape, scale = 
  xwei.scale,lower.tail=TRUE, log.p=FALSE) }
 
  ## convert dataframe to numeric
  MatOut- as.matrix(ResultMat2)
  mode(MatOut) - numeric
 
  # initiate variables
  mm-ml-mu-rep(0,length(MatOut[,1]))
 
  # mean and upper/lower quantiles
  for(i in 1:length(MatOut[,1])){
  mm[i]- mean(MatOut[i,])
  ml[i]- quantile(MatOut[i,], 0.025, na.rm=TRUE)
  mu[i]- quantile(MatOut[i,], 0.975, na.rm=TRUE) } #lines(x, mm, 
  col=black) lines(x, ml, col=blue, lwd=2) lines(x, mu, 
 col=blue, 
  lwd=2)
 
  [[alternative HTML version deleted]]
 
  __
  R-help@r-project.org 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.
 
 David Winsemius, MD
 West Hartford, CT
 
 __
 R-help@r-project.org 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.
 ***
This email and any

[R] confidence interval on mean in log-regression

[Alice Shelly]

Hello-

I am looking for R function that will give me some proper confidence
intervals on un-transformed mean prediction when performing a linear
regression on log-transformed data. I am referring to the UMVU estimate, the
el-shaarawi and viveros (1997) estimate, or the Wu, wong, and Wei (2005)
estimates, for example. Can't find anything, hoping not to have to program
it myself!

 

Thanks-

 

Alice Shelly

Environmental Statistician

TerraStat Consulting Group

5002 Lodge View Lane

Austin, TX 78731

(206) 362-3299

al...@blarg.net

 


[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval of Specral Density Plot

[Zablone Owiti]

Dear Users,

I wish to know at what confidence level is the  confidence interval 
provided in the  Spectrum function (plot.spec) plots.
The only information provided in the help regarding this is : a confidence 
interval will be plotted by plot.spec: this is asymmetric, and the width of 
the centre mark indicates the equivalent bandwidth. 

Attached is an example plot with the BLUE confidence interval.

Thanks


 ZABLONE OWITI

 GRADUATE STUDENT

College of Atmospheric Science

Nanjing University of Information, Science and Technology

Add: 219 Ning Liu Rd, Nanjing, Jiangsu, 21004, P.R. China

 Tel: +86-25-58731402

Fax: +86-25-58731456 
Mob. 15077895632 
Website:  www.nuist.edu.cn 



DO NOT PRINT THIS E-MAIL UNLESS NECESSARY. THE ENVIRONMENT CONCERNS US 
ALL.

 
__
R-help@r-project.org 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] Confidence Interval of Specral Density Plot

[David Winsemius]

On May 21, 2011, at 3:27 AM, Zablone Owiti wrote:


Dear Users,

I wish to know at what confidence level is the  confidence interval
provided in the  Spectrum function (plot.spec) plots.


The default level is clearly indicated in the arguments of the Usage  
section on the help page of plot.spec.


The only information provided in the help regarding this is : a  
confidence
interval will be plotted by plot.spec: this is asymmetric, and the  
width of

the centre mark indicates the equivalent bandwidth.


The authors of R consider the code to be the final arbiter of what  
documentation is provided. The code is readily available, since the  
plot.spec function is visible,  and it is not at all difficult to find  
the function that calculates the CI, since it is right at the top when  
you type: plot.spec




Attached is an example plot with the BLUE confidence interval.


Nothing attached, but in this case I don't think it would add  
anything. You probably need to read the Posting Guide before  
attempting to attach items. The mail-server is very picky about what  
it will accept.




Thanks


ZABLONE OWITI

GRADUATE STUDENT




David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence interval for difference in Harrell's c statistics (or equivalently Somers' D statistics)

[Laura Bonnett]

Dear All,

I am trying to calculate a 95% confidence interval for the difference in two
c statistics (or equivalently D statistics).  In Stata I gather that this
can be done using the lincom command.  Is there anything similar in R?

As you can see below I have two datasets (that are actually two independent
subsets of the same data) and the respective c statistics for the variables
in both cases.  What I would now like to do is to prove that there is no
statistically significant difference between the statistics (between the dev
and val datasets.)

Any help would be much appreciated.

 rdev - rcorrcens(Surv(stimes1,eind1)~gendat1+neurodat1)
 rdev

Somers' Rank Correlation for Censored DataResponse
variable:Surv(stimes1, eind1)

  CDxy  aDxySDZ  Pn
gendat1   0.534  0.069 0.069 0.017 3.98 0.0001 1500
neurodat1 0.482 -0.036 0.036 0.011 3.18 0.0015 1500

 rval - rcorrcens(Surv(stimes2,eind2)~gendat2+neurodat2)
 rval

Somers' Rank Correlation for Censored DataResponse
variable:Surv(stimes2, eind2)

  CDxy  aDxySDZ Pn
gendat2   0.543  0.085 0.085 0.017 4.94 0e+00 1500
neurodat2 0.481 -0.038 0.038 0.011 3.44 6e-04 1500

Many thanks,
Laura
P.S. I'm using Windows XP, R 2.9.2

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval for difference in Harrell's c statistics (or equivalently Somers' D statistics)

[David Winsemius]

On May 5, 2011, at 8:20 AM, Laura Bonnett wrote:


Dear All,

I am trying to calculate a 95% confidence interval for the  
difference in two
c statistics (or equivalently D statistics).  In Stata I gather that  
this

can be done using the lincom command.  Is there anything similar in R?


Have you looked at rcorrp.cens {Hmisc}?

snipped


[[alternative HTML version deleted]]



Please post in plain text.

--
David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence interval for difference in Harrell's c statistics (or equivalently Somers' D statistics)

[Laura Bonnett]

Dear David,

Thank you for your reply.  I have come across rcorrp.cens before.  However,
I'm not sure it does quite what I want it to.  It seems to compare whether
one predictor is more concordant than another within the same survival
function.  I want to see whether one predictor is more concordant than
another over two survival functions hence I fitted two rcorrcens functions.

E.g. if I have a development data set with a variable for age and a
validation data set for age then I want to know if the concordance is the
same over the development and validation data sets.

Thank you,
Laura

On 5 May 2011 16:09, David Winsemius dwinsem...@comcast.net wrote:


 On May 5, 2011, at 8:20 AM, Laura Bonnett wrote:

  Dear All,

 I am trying to calculate a 95% confidence interval for the difference in
 two
 c statistics (or equivalently D statistics).  In Stata I gather that this
 can be done using the lincom command.  Is there anything similar in R?


 Have you looked at rcorrp.cens {Hmisc}?

 snipped


[[alternative HTML version deleted]]



 Please post in plain text.

 --
 David Winsemius, MD
 West Hartford, CT



[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval for the difference between proportions - method used in prop.test()

[Stefanie Von Felten]

Hello,

Does anyone know which method from Newcombe (1998)* is implemented in prop.test 
for comparing two proportions?
I would guess it is the method based on the Wilson score (for single 
proportion), with and without continuity correction  for prop.test(..., 
correct=FALSE) and prop.test(..., correct=TRUE). These methods would correspond 
to no. 10 and 11 tested in Newcombe, respectively. Can someone confirm this? If 
not, which other methods are implemented by prop.test?

* Newcombe R.G. (1998) Two-Sided Confidence Intervals for the Single 
Proportion: Comparison of Seven Methods.  Statistics in Medicine  *17*, 857-872.

There is also the function ci.pd() from the R-package Epi, which should 
implement method no. 10 from Newcombe. However, prop.test(..., correct=FALSE) 
and ci.pd do not give the same result if I do the following:

successes - c(21, 41)
total - c(345, 345)
prop.test(successes, total, correct=FALSE)
library(Epi)
ci.pd(matrix(c(successes, total-successes),ncol=2, byrow=TRUE))

Can someone explain why?

Best wishes
Stefanie von Felten


Stefanie von Felten, PhD
Statistician
Clinical Trial Unit, CTU
University Hospital Basel
Schanzenstrasse 55
CH-4031 Basel

Phone: ++41(0)61 556 54 98

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval for the difference between proportions - method used in prop.test()

[Joshua Wiley]

Hi Stefanie,

Just to be clear, we are talking about differences in the third or
lower decimal place (at least with R version 2.13.0 alpha (2011-03-17
r54849), Epi_1.1.20).  This strikes me as small enough that both
functions may be implementing the same method, but maybe slightly
different ways of going about it?

If you are really concerned and need to know *exactly*, look at the
source code for both functions.  In case you did not know, if you type
the function name at the console with parentheses or any arguments,
just like:

 prop.test
 ci.pd

it will show the actual function code.  It looks to me like both of
them are implemented purely in R, and without even calling any other
complex functions (at least based on a quick glance through).  This
means if you have the Newscomb text, you should be able to sit down
and go through the code step by step comparing it.

Cheers,

Josh

FYI, you can use a matrix with prop.test, and then its transpose for ci.pd.
##
mymat - cbind(Successes = c(21, 41), Failures = c(345, 345) - c(21, 41))
require(Epi)
results - list(prop.test(mymat, correct=FALSE), ci.pd(t(mymat)))
results[[1]][[conf.int]] - results[[2]][6:7]

On Tue, Apr 5, 2011 at 3:38 AM, Stefanie Von Felten svonfel...@uhbs.ch wrote:
 Hello,

 Does anyone know which method from Newcombe (1998)* is implemented in 
 prop.test for comparing two proportions?
 I would guess it is the method based on the Wilson score (for single 
 proportion), with and without continuity correction  for prop.test(..., 
 correct=FALSE) and prop.test(..., correct=TRUE). These methods would 
 correspond to no. 10 and 11 tested in Newcombe, respectively. Can someone 
 confirm this? If not, which other methods are implemented by prop.test?

 * Newcombe R.G. (1998) Two-Sided Confidence Intervals for the Single 
 Proportion: Comparison of Seven Methods.  Statistics in Medicine  *17*, 
 857-872.

 There is also the function ci.pd() from the R-package Epi, which should 
 implement method no. 10 from Newcombe. However, prop.test(..., correct=FALSE) 
 and ci.pd do not give the same result if I do the following:

 successes - c(21, 41)
 total - c(345, 345)
 prop.test(successes, total, correct=FALSE)
 library(Epi)
 ci.pd(matrix(c(successes, total-successes),ncol=2, byrow=TRUE))

 Can someone explain why?

 Best wishes
 Stefanie von Felten


 Stefanie von Felten, PhD
 Statistician
 Clinical Trial Unit, CTU
 University Hospital Basel
 Schanzenstrasse 55
 CH-4031 Basel

 Phone: ++41(0)61 556 54 98

        [[alternative HTML version deleted]]

 __
 R-help@r-project.org 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.




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

__
R-help@r-project.org 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] Confidence interval for the difference between proportions - method used in prop.test()

[Joshua Wiley]

Dear Steffi,

On Tue, Apr 5, 2011 at 7:26 AM, Stefanie Von Felten svonfel...@uhbs.ch wrote:
 Dear Josh,

 Thanks for your help!

 Does your answer mean, that you agree the two methods should do the same,
 and what I was guessing, despite the small differences?

That would be my guess, but I have not actually read the reference in
discussion.  Still, the documentation for prop.test uses the same
Newcombe reference as ci.pd, so if method 10 is the clear winner, it
seems reasonable that prop.test is an implementation of it.


 What I prefer about ci.pd is, that the help clearly says which method is
 implemented, which is not the case for prop.test. But I do not know who has
 programmed the function.

Then for reporting, use ci.pd, and say it is method 10 from Newcombe.
You can always check the results with prop.test() which is part of R
core so you can be pretty confident whatever it does, it does it
correctly (and will be updated if necessary with future releases of
R).

Sincerely,

Josh

 Best wishes
 Steffi


 Stefanie von Felten, PhD
 Statistician
 Clinical Trial Unit, CTU
 University Hospital Basel
 Schanzenstrasse 55
 CH-4031 Basel

 Phone: ++41(0)61 556 54 98
 Joshua Wiley jwiley.ps...@gmail.com 05.04.11 15.59 Uhr 
 Hi Stefanie,

 Just to be clear, we are talking about differences in the third or
 lower decimal place (at least with R version 2.13.0 alpha (2011-03-17
 r54849), Epi_1.1.20). This strikes me as small enough that both
 functions may be implementing the same method, but maybe slightly
 different ways of going about it?

 If you are really concerned and need to know *exactly*, look at the
 source code for both functions. In case you did not know, if you type
 the function name at the console with parentheses or any arguments,
 just like:

 prop.test
 ci.pd

 it will show the actual function code. It looks to me like both of
 them are implemented purely in R, and without even calling any other
 complex functions (at least based on a quick glance through). This
 means if you have the Newscomb text, you should be able to sit down
 and go through the code step by step comparing it.

 Cheers,

 Josh

 FYI, you can use a matrix with prop.test, and then its transpose for ci.pd.
 ##
 mymat - cbind(Successes = c(21, 41), Failures = c(345, 345) - c(21, 41))
 require(Epi)
 results - list(prop.test(mymat, correct=FALSE), ci.pd(t(mymat)))
 results[[1]][[conf.int]] - results[[2]][6:7]

 On Tue, Apr 5, 2011 at 3:38 AM, Stefanie Von Felten svonfel...@uhbs.ch
 wrote:
 Hello,

 Does anyone know which method from Newcombe (1998)* is implemented in
 prop.test for comparing two proportions?
 I would guess it is the method based on the Wilson score (for single
 proportion), with and without continuity correction for prop.test(...,
 correct=FALSE) and prop.test(..., correct=TRUE). These methods would
 correspond to no. 10 and 11 tested in Newcombe, respectively. Can someone
 confirm this? If not, which other methods are implemented by prop.test?

 * Newcombe R.G. (1998) Two-Sided Confidence Intervals for the Single
 Proportion: Comparison of Seven Methods. Statistics in Medicine *17*,
 857-872.

 There is also the function ci.pd() from the R-package Epi, which should
 implement method no. 10 from Newcombe. However, prop.test(...,
 correct=FALSE) and ci.pd do not give the same result if I do the following:

 successes - c(21, 41)
 total - c(345, 345)
 prop.test(successes, total, correct=FALSE)
 library(Epi)
 ci.pd(matrix(c(successes, total-successes),ncol=2, byrow=TRUE))

 Can someone explain why?

 Best wishes
 Stefanie von Felten


 Stefanie von Felten, PhD
 Statistician
 Clinical Trial Unit, CTU
 University Hospital Basel
 Schanzenstrasse 55
 CH-4031 Basel

 Phone: ++41(0)61 556 54 98

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org 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.




 --
 Joshua Wiley
 Ph.D. Student, Health Psychology
 University of California, Los Angeles
 http://www.joshuawiley.com/




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

__
R-help@r-project.org 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] Confidence interval of Survival Curve of Weighted Cox Regression

[David Winsemius]

On Feb 15, 2011, at 9:05 PM, jeanneyue wrote:



Hi,

May I know how to obtain the confidence interval of the survival  
curve of

weighted Cox regression model?
I tried coxph, cph, and coxphw, but they did not work.
Any help would be much appreciated.


One possible reason that this question may have gone unanswered is  
that it is too vague. You don't indicate enough about your data  
situation or research goals to answer the question. And you are  
mentioning functions from different packages but not the name of the  
packages. I know for instance that cph is from rms and the it is not  
used to generate confidence intervals but that Predict would. Perhaps  
the same situation exists for coxph (presumably from survival) and  
coxphw whixh is not from either of hthose two packages.


You are asked to offered code and at the very least a description of  
the data.In that manner we can figure out what you mean by weights, a  
term that has at least two (and maybe more)  different interpretations  
in this context.


Since it appears you have not read the Posting Guide, please do so now  
and then post a more complete question.


--

David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence interval of Survival Curve of Weighted Cox Regression

[Terry Therneau]

--- begin included message ---
May I know how to obtain the confidence interval of the survival curve
of
weighted Cox regression model?
I tried coxph, cph, and coxphw, but they did not work.
Any help would be much appreciated. 

 end inclusion ---

Use the latest version of the survival pacakge from CRAN, coxph now
supports weighted survival curves.

Terry Therneau

__
R-help@r-project.org 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] Confidence interval of Survival Curve of Weighted Cox Regression

[jeanneyue]

Hi,

May I know how to obtain the confidence interval of the survival curve of
weighted Cox regression model?
I tried coxph, cph, and coxphw, but they did not work.
Any help would be much appreciated. 

Thanks,
Jeanne


-- 
View this message in context: 
http://r.789695.n4.nabble.com/Confidence-interval-of-Survival-Curve-of-Weighted-Cox-Regression-tp3308085p3308085.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence interval

[Peter Ehlers]

On 2011-02-13 18:31, Joshua Wiley wrote:

Hi,

The logical operators are actually vectorized, so I do not think you
need a loop.  Does this do what you want?

## Some data
set.seed(10)
dat- matrix(rnorm(500, sd = 3), nrow = 80)

## Hypothetical confidence interval
ci- c(-5, 5)

## Find the number of points outside interval
sum(dat  ci[1] | dat  ci[2], na.rm = TRUE)

Cheers,

Josh


Or you could use (no simpler) findInterval():

 fI - findInterval(dat, sort(ci))

## to see what's produced:
 table(fI)
#  0   1   2
# 28 512  20

## the 1s indicate values inside ci, etc,
## so we want

 sum(fI != 1)
# [1] 48

Peter Ehlers



On Sun, Feb 13, 2011 at 4:26 PM, Syd88jhea2...@uni.sydney.edu.au  wrote:


Hi,
I am trying to determine how many points fall ouside the confidence interval
range.

This is the code I have so far but it does not work. Any help would be
appreciated.

Count- vector ()
for (i in 1: nrow (dataname)){
if (dataname[i]l.ci.post[1]//
dataname[i]u.ci.post[i]){
count[i] -  1
}else
{count[i] -  0}
}



symbol
// = or - not sure if this is the right symbol though
--
View this message in context: 
http://r.789695.n4.nabble.com/Confidence-interval-tp3304258p3304258.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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@r-project.org 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] Confidence interval

[Syd88]

Hi,
I am trying to determine how many points fall ouside the confidence interval
range.

This is the code I have so far but it does not work. Any help would be
appreciated.

Count - vector ()
for (i in 1: nrow (dataname)){
if (dataname[i] l.ci.post[1]//
dataname[i] u.ci.post[i]){
count[i] - 1
}else
{count[i] - 0}
}



symbol
// = or - not sure if this is the right symbol though
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Confidence-interval-tp3304258p3304258.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence interval

[Joshua Wiley]

Hi,

The logical operators are actually vectorized, so I do not think you
need a loop.  Does this do what you want?

## Some data
set.seed(10)
dat - matrix(rnorm(500, sd = 3), nrow = 80)

## Hypothetical confidence interval
ci - c(-5, 5)

## Find the number of points outside interval
sum(dat  ci[1] | dat  ci[2], na.rm = TRUE)

Cheers,

Josh

On Sun, Feb 13, 2011 at 4:26 PM, Syd88 jhea2...@uni.sydney.edu.au wrote:

 Hi,
 I am trying to determine how many points fall ouside the confidence interval
 range.

 This is the code I have so far but it does not work. Any help would be
 appreciated.

 Count - vector ()
 for (i in 1: nrow (dataname)){
 if (dataname[i] l.ci.post[1]//
 dataname[i] u.ci.post[i]){
 count[i] - 1
 }else
 {count[i] - 0}
 }



 symbol
 // = or - not sure if this is the right symbol though
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Confidence-interval-tp3304258p3304258.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org 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.




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

__
R-help@r-project.org 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] Confidence interval based on MLE

[Jinsong Zhao]

On 2011-2-6 22:56, Ben Bolker wrote:

Jinsong Zhaojszhaoat  yeah.net  writes:



Hi there,

I have fitted a sample (with size 20) to a normal and/or logistic
distribution using fitdistr() in MASS or fitdist() in fitdistrplus
package. It's easy to get the parameter estimates. Now, I hope to report
the confidence interval for those parameter estimates. However, I don't
find a function that could give the confidence interval in R.

I hope to write a function, however, I don't find some detailed
information on the CI based on MLE. Would you please to give me some
hints on the CI calculation based on MLE?


Well, for the normal distribution I believe that the standard-error-
based confidence intervals are the same as those based on the MLE,
but in general I would suggest something along these lines:


library(bbmle)
z- rnorm(20)
m- mle2(z~dnorm(mean=mu,sd=sd),start=list(mu=0,sd=1),data=data.frame(z))

Warning message:
In dnorm(x, mean, sd, log) : NaNs produced

confint(m)

Profiling...
  2.5 %   97.5 %
mu -0.07880835 0.985382
sd  0.87314467 1.633600



Thank you very much for your kindly help and the way to get MLE through 
bbmle package. It works well.


I have a interval related question. I have a sample data set, with size 
20 or less. And I fit it to a three parameter distribution, e.g., a 
triangular distribution (oops, it cannot fitted by mle2 :-(). I get the 
quantile, q, for a given probability, p. Then, I hope to get the 
confidence (or prediction?) interval for the quantile, q. However, I 
don't know how to do.


I refer to some books on ecological data analysis. There's a explicit 
formula for CI to the normal distribution's q, based on delta method or 
Fieller's theorem. (And I think they should work for logistic 
distribution). But I don't find any thing that for other distribution.


BTW, is it possible to get a interval of p for a given q? Although, it's 
not a normal way in the view of statistics, it has a lot applications.


Any suggestions or comments will be really appreciated. Thanks in advance.

Regards,
Jinsong

__
R-help@r-project.org 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] Confidence interval based on MLE

[Jinsong Zhao]

Hi there,

I have fitted a sample (with size 20) to a normal and/or logistic 
distribution using fitdistr() in MASS or fitdist() in fitdistrplus 
package. It's easy to get the parameter estimates. Now, I hope to report 
the confidence interval for those parameter estimates. However, I don't 
find a function that could give the confidence interval in R.


I hope to write a function, however, I don't find some detailed 
information on the CI based on MLE. Would you please to give me some 
hints on the CI calculation based on MLE?


Any suggestions will be really appreciated. Thanks in advance.

Regards,
Jinsong

__
R-help@r-project.org 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] Confidence interval based on MLE

[Ben Bolker]

Jinsong Zhao jszhao at yeah.net writes:

 
 Hi there,
 
 I have fitted a sample (with size 20) to a normal and/or logistic 
 distribution using fitdistr() in MASS or fitdist() in fitdistrplus 
 package. It's easy to get the parameter estimates. Now, I hope to report 
 the confidence interval for those parameter estimates. However, I don't 
 find a function that could give the confidence interval in R.
 
 I hope to write a function, however, I don't find some detailed 
 information on the CI based on MLE. Would you please to give me some 
 hints on the CI calculation based on MLE?

   Well, for the normal distribution I believe that the standard-error-
based confidence intervals are the same as those based on the MLE,
but in general I would suggest something along these lines:

 library(bbmle)
 z - rnorm(20)
 m - mle2(z~dnorm(mean=mu,sd=sd),start=list(mu=0,sd=1),data=data.frame(z))
Warning message:
In dnorm(x, mean, sd, log) : NaNs produced
 confint(m)
Profiling...
 2.5 %   97.5 %
mu -0.07880835 0.985382
sd  0.87314467 1.633600

__
R-help@r-project.org 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] confidence interval

[Francesco Petrogalli]

Hi,
I have a circular shaped set of point on the plane (X,Y) centered in
zero. The distribution is more dense close to zero and less dense far
from zero.

I need to find the radius of a circle centered in zero that contains
65% of the points in the sample. Is there any R directive that can do
this?

I wanna start with 2D set of points, but the real case scenario is
with a 5D set of points.

Thanks,

Francesco

__
R-help@r-project.org 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] confidence interval

[David Winsemius]

On Jan 21, 2011, at 11:33 AM, Francesco Petrogalli wrote:


Hi,
I have a circular shaped set of point on the plane (X,Y) centered in
zero. The distribution is more dense close to zero and less dense far
from zero.

I need to find the radius of a circle centered in zero that contains
65% of the points in the sample. Is there any R directive that can do
this?

I wanna start with 2D set of points, but the real case scenario is
with a 5D set of points.


Something along the lines of

dxy= with(dfm, sqrt(x^2 +y^2))
quantile(dxy, probs=0.65)

The generalization to 5 dimensions appears trivial. Even the  
generalization to finding the radius around an arbitrary point seems  
trivial assuming an L2 norm.


--
David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence interval on quantile regression predictions

[Davey, Andrew]

I am using the quantreg package to build a quantile regression model and
wish to generate confidence intervals for the fitted values.

After fitting the model, I have tried running predict() and
predict.rq(), but in each case I obtain a vector of the fitted values
only.

For example:

library(quantreg)
y-rnorm(50,10,2)
x-seq(1,50,1)
model-rq(y~x,tau=0.5,method=br,model=TRUE)
model$fitted
predict(model,type=c(percentile),interval=c(confidence),level=0.95)
#returns same output as model$fitted

Am I missing something obvious? I have tried altering the settings for
the type and interval arguments but without success.

(I am running R v 2.10.1).

Thank you for your help.

Regards,


Andrew Davey
Senior Consultant

Tel:+ 44 (0) 1793 865023
Mob:+44 (0) 7747 863190
Fax:+ 44 (0) 1793 865001
E-mail: andrew.da...@wrcplc.co.uk
Website:www.wrcplc.co.uk
P Before printing, please think about the environment.

---
To read our latest newsletter visit 
http://www.wrcplc.co.uk/default.aspx?item=835 - Keeping our clients up-to-date 
with WRc's business developments
---
Visit our websites www.wrcplc.co.uk and www.wrcnsf.com, as well as 
www.waterportfolio.com for collaborative research projects.
---
The Information in this e-mail is confidential and is intended solely for the 
addressee. Access to this e-mail by any other party is unauthorised. If you are 
not the intended recipient, any disclosure, copying, distribution or any action 
taken in reliance on the information contained in this e-mail is prohibited and 
maybe unlawful. When addressed to WRc Clients, any opinions or advice contained 
in this e-mail are subject to the terms and conditions expressed in the 
governing WRc Client Agreement.
---
WRc plc is a company registered in England and Wales. Registered office 
address: Frankland Road, Blagrove, Swindon, Wiltshire SN5 8YF. Company 
registration number 2262098. VAT number 527 1804 53.
---
[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval on quantile regression predictions

[Roger Koenker]

You need to add explicitly newdata = list(x=x) and if you want percentile 
method you also
need se = boot.  


Roger Koenker
rkoen...@illinois.edu

On Jan 11, 2011, at 3:54 AM, Davey, Andrew wrote:

 I am using the quantreg package to build a quantile regression model and
 wish to generate confidence intervals for the fitted values.
 
 After fitting the model, I have tried running predict() and
 predict.rq(), but in each case I obtain a vector of the fitted values
 only.
 
 For example:
 
 library(quantreg)
 y-rnorm(50,10,2)
 x-seq(1,50,1)
 model-rq(y~x,tau=0.5,method=br,model=TRUE)
 model$fitted
 predict(model,type=c(percentile),interval=c(confidence),level=0.95)
 #returns same output as model$fitted
 
 Am I missing something obvious? I have tried altering the settings for
 the type and interval arguments but without success.
 
 (I am running R v 2.10.1).
 
 Thank you for your help.
 
 Regards,
 
 
 Andrew Davey
 Senior Consultant
 
 Tel:  + 44 (0) 1793 865023
 Mob:+44 (0) 7747 863190
 Fax:  + 44 (0) 1793 865001
 E-mail:   andrew.da...@wrcplc.co.uk
 Website:  www.wrcplc.co.uk
 P Before printing, please think about the environment.
 
 ---
 To read our latest newsletter visit 
 http://www.wrcplc.co.uk/default.aspx?item=835 - Keeping our clients 
 up-to-date with WRc's business developments
 ---
 Visit our websites www.wrcplc.co.uk and www.wrcnsf.com, as well as 
 www.waterportfolio.com for collaborative research projects.
 ---
 The Information in this e-mail is confidential and is intended solely for the 
 addressee. Access to this e-mail by any other party is unauthorised. If you 
 are not the intended recipient, any disclosure, copying, distribution or any 
 action taken in reliance on the information contained in this e-mail is 
 prohibited and maybe unlawful. When addressed to WRc Clients, any opinions or 
 advice contained in this e-mail are subject to the terms and conditions 
 expressed in the governing WRc Client Agreement.
 ---
 WRc plc is a company registered in England and Wales. Registered office 
 address: Frankland Road, Blagrove, Swindon, Wiltshire SN5 8YF. Company 
 registration number 2262098. VAT number 527 1804 53.
 ---
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org 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@r-project.org 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] confidence interval for logistic joinpoint regression from package ljr

[Vito Muggeo (UniPa)]

dear M.,
I do not know how to get the SE for the joinpoint (or breakpoint) from 
your ljr fit. However you can find useful the segmented package which 
works for any GLM (including the logistic one) and it returns 
(approximate) StErr (and Conf Int) also for the joinpoint (breakpoint in 
the segmented package). For some examples, see the R news paper

http://cran.r-project.org/doc/Rnews/Rnews_2008-1.pdf

Also, as regards to the APC, you could find interesting the following note

http://onlinelibrary.wiley.com/doi/10.1002/sim.3850/abstract

hope this helps you,
vito


Il 30/11/2010 23.54, moleps ha scritto:

I´m trying to run a logistic joinpoint regression utilising the ljr package. 
I´ve been using the forward selection technique to get the number of knots for 
the analysis, but I´m uncertain as to my results and the interpretation. The 
documentation is rather brief ( in the package and the stats in medicine 
article is quite technical) and without any good examples. At the moment I´m 
thinking
1)find the number of knots both using forward and backward techniques and see 
if they are close
2)present the annual percent change (APC) for each of the intervals, ie my 
present data (1950-2010 in 5 year intervals) is giving me

Variables  Coef
b0 Intercept -131.20404630
g0 t0.06146463
g1 max(t-tau1,0)   -0.51582466
g2 max(t-tau2,0)0.43429615

Joinpoints:

1 tau1= 1990.5
2 tau2= 1995.5

APC 1950-1990=exp(0.06)=1.06--6%
1990-1995=exp(0.06-0.51)=exp(-0.45)=0.63--  -37%
1995-2010=exp(0.06-0.51+0.43)---2%

3) Preferably a confidence interval for the APC should be given. However, this 
I havent figured out yet.

//M

__
R-help@r-project.org 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.



--

Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 23895240
fax: 091 485726/485612
http://dssm.unipa.it/vmuggeo

__
R-help@r-project.org 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] confidence interval for logistic joinpoint regression from package ljr

[moleps]

I´m trying to run a logistic joinpoint regression utilising the ljr package. 
I´ve been using the forward selection technique to get the number of knots for 
the analysis, but I´m uncertain as to my results and the interpretation. The 
documentation is rather brief ( in the package and the stats in medicine 
article is quite technical) and without any good examples. At the moment I´m 
thinking 
1)find the number of knots both using forward and backward techniques and see 
if they are close
2)present the annual percent change (APC) for each of the intervals, ie my 
present data (1950-2010 in 5 year intervals) is giving me

   Variables  Coef
b0 Intercept -131.20404630
g0 t0.06146463
g1 max(t-tau1,0)   -0.51582466
g2 max(t-tau2,0)0.43429615

Joinpoints:
  
1 tau1= 1990.5
2 tau2= 1995.5

APC 1950-1990=exp(0.06)=1.06--6%
1990-1995=exp(0.06-0.51)=exp(-0.45)=0.63-- -37%
1995-2010=exp(0.06-0.51+0.43)---2%

3) Preferably a confidence interval for the APC should be given. However, this 
I havent figured out yet.

//M

__
R-help@r-project.org 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] Confidence interval for response variable in mixed effects models

[Brian Willis]

HI,

I am using lmer() for a simple mixed effects model. The model is of the form
logit(y)~ x + (1|z), where x is an indicator variable and z a multi-level
factor.

I would like an estimate of the response variable (either y or logit y) with
an associated confidence interval for a given value of x.

There does not appear to be a predict function written for lmer().

The output for the fixed effects gives a standard error for the intercept,
the coefficient of x and the correlation. For n observations, I transform
the std errors to variances, and with the correlation, I think I can use the
formula Var(intercept + x) = Var(intercept) + Var(x) + 2Cov(x,intercept) to
get the variance of the fixed effects component of the response variable.

However, I would like to include the random effects component of the
variance, so I may derive a standard error and confidence interval for the
variance.

Does anyone know how to do this? Is there a ready made function like
predict() or does anyone know how to incorporate the variance of random
effects term to derive the std error of the response variable?

Regards,

Brian

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval for response variable in mixed effects models

[Ben Bolker]

Brian Willis brian.willis at manchester.ac.uk writes:

 I am using lmer() for a simple mixed effects model. The model is of the form
 logit(y)~ x + (1|z), where x is an indicator variable and z a multi-level
 factor.
 
 I would like an estimate of the response variable (either y or logit y) with
 an associated confidence interval for a given value of x.

[snippage: sorry to remove context, but I am posting via gmane, which
will complain if I have too much quoted context ...]

 Does anyone know how to do this? Is there a ready made function like
 predict() or does anyone know how to incorporate the variance of random
 effects term to derive the std error of the response variable?

  You should search the r-sig-mixed-models archive for answers, and
post there if you don't find what you need.  The problem is that it
can actually be a bit tricky to define these things properly for mixed
models, decide which random effects to include (or not) in the prediction
of the mean and include (or not) in the definition of the variance.
So far there has not been a confluence of people who want this,
people who know enough to construct a nice general solution, and
people who have time to do it ...

  Ben Bolker

__
R-help@r-project.org 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] Confidence interval for response variable in mixed effects models

[David Winsemius]

On Oct 30, 2010, at 11:33 AM, Brian Willis wrote:


HI,

I am using lmer() for a simple mixed effects model. The model is of  
the form
logit(y)~ x + (1|z), where x is an indicator variable and z a multi- 
level

factor.

I would like an estimate of the response variable (either y or logit  
y)


There is a fitted method for objects of class mer.

?mer-class


 with an associated confidence interval for a given value of x.


That is the sticking point as I understand it. See below.


There does not appear to be a predict function written for lmer().

The output for the fixed effects gives a standard error for the  
intercept,
the coefficient of x and the correlation. For n observations, I  
transform
the std errors to variances, and with the correlation, I think I can  
use the
formula Var(intercept + x) = Var(intercept) + Var(x) +  
2Cov(x,intercept) to
get the variance of the fixed effects component of the response  
variable.


However, I would like to include the random effects component of the
variance, so I may derive a standard error and confidence interval  
for the

variance.

Does anyone know how to do this? Is there a ready made function like
predict() or does anyone know how to incorporate the variance of  
random

effects term to derive the std error of the response variable?


I see Ben Bolker has offered a suggestion that you search the mixed  
models list.


An earlier question on this topic elicited this citation (from Bolker):

http://glmm.wikidot.com/faq

One of the available choices in Baron's search page is R-sig-mixed- 
models and here are a few links to threads (by people who know more  
than me) that may offer at least informed commentary and options if  
not a general solution (from  both R-help and the SIG-ME archives):


http://search.r-project.org/cgi-bin/namazu.cgi?query=+confidence+intervals+lmermax=100result=normalsort=scoreidxname=functionsidxname=Rhelp08idxname=Rhelp10idxname=R-sig-mixed-modelsidxname=Rhelp02

http://finzi.psych.upenn.edu/Rhelp10/2008-May/161516.html

http://rwiki.sciviews.org/doku.php?id=guides:lmer-tests

http://finzi.psych.upenn.edu/R/Rhelp02/archive/76742.html

http://finzi.psych.upenn.edu/R/Rhelp02/archive/81237.html

http://finzi.psych.upenn.edu/R/Rhelp02/archive/82567.html

--

David Winsemius, MD
West Hartford, CT

__
R-help@r-project.org 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] Confidence Interval in ARMA mean GARCH model

[Aditya Damani]

Hi,

How do I change Confidence Interval level (say from 95% to 80%) while
getting prediction intervals for ARMA mean models, and GARCH models in
R?

TIA
Aditya

__
R-help@r-project.org 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] Confidence interval calculation for intersection of two quadratic lines

[Leif Kirschenbaum]

How do I calculate the confidence interval for the value x given by the 
intersection of two quadratics (i.e. parabolas)?

I fit two quadratics of the form:

  y = C1 + B1*x + A1*x^2
  y = C2 + B2*x + A2*x^2

to two sets of points N1 and N2.

I test for whether they intersect, if they do then I calculate the roots of:
  0 = (C1 - C2) + (B1 - B2)*x + (A1 - A2)*x^2

to determine where they intersect and choose the proper root.
I have propagated errors through the quadratic equation so that I may calculate 
the variance of each root.

Naively I might take the variance of the root, find a Student-t quantile for 
degrees-of-freedom = (N1 + N2 - 6), and then calculate
C.I. = root +/- variance * t-quantile 

However I have several reservations about this approach:
(1) Let's say that one quadratic fit devolves to a linear fit (A1 ~ 0)*, then I 
have the intersection of a straight line and a quadratic.  (I code this 
situation separately with different error propagation formulas)  However 
intuitively if the intercept of the straight line varies by some dy up or down 
I would expect that the intersection with a parabola varies by a different 
amount if the straight line is moved up versus down, e.g. the confidence 
interval should be asymmetric.

(2) The variances for the coefficients A1, B1, C1 depend on N1 data points (and 
A2,B2,C2 on N2) and therefore intuitively it doesn't seem accurate to use a 
Student-t quantile derived from (N1 + N2).  Another way to think of it is if I 
constructed a linear model:
  x.1  - x  * (x   0)
  x.2  - x  * (x = 0)
  x2   - x^2
  x2.1 - x2 * (x   0)
  x2.2 - x2 * (x = 0)
  mylm - lm(y ~ x.1 + x.2 + x2.1 + x2.2) 
some coefficients are dependent on N1 data points (and their noise) and other 
coefficients depend on N2 data points (and their noise). (and I don't think lm 
would given me an aswer with such bifurcated 


Sincerely,


Leif S. Kirschenbaum, Ph.D., PMP, CRE


* In reality I will perform both quadratic and linear fits on N1 and N2, 
compare MSEs of the quadratic vs. linear on N1 to determine the best fit for N1 
points and similarly on N2.  The situation is that I am measuring many 
instances of an object (~16,000 devices per wafer) whose behavior is not well 
modelled at room temperature.  By inspection of ~4,000 plots of device data I 
see that there is sometimes a linear dependence and sometimes a dependence with 
curvature adequately represented by a quadratic fit.
The zero temperature device behavior is given by ( (y-yoffset)/yscale )^(2/3) + 
( (x-xcenter)/xscale )^(2/3) = 1 which is darn hard to fit with linear models 
even for perfect devices, which we do not have.  The room temperature device 
behavior is so poorly understood that it is only simulated with finite element 
models, which does not give me a neat analytic expression such as p^(2/3) + 
q^(2/3) = 1, and in any event the plots rarely look anything like x^(2/3) 
behavior.

Background:
I am developing a test system in LabVIEW, however I simulate data analysis 
offline in R, i.e. acquire raw data on the test system, import raw data into R 
and try analysis methods, then port the methods back to the test system.
Therefore use of any R functions more sophisticated than lm will be of limited 
use; LabVIEW does not have very sophisticated statistical functions (thank 
goodness it includes inverse-t).  I have already developed propagation of 
quadratic fit variances and written my own linear fit confidence interval 
calculations and now I am extending this to quadratics and to this intersection 
problem.

__
R-help@r-project.org 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] confidence interval for small non iid sample

[Kay Cichini]

hello,

can someone tell me if there is a way to estimate a confidence interval for
a small non iid sample. 
i computed a stratified boot.ci but i think it is not reasonable in the case
of such a small sample - are there any alternatives, can a conf.interval in
this case be estimated at all?

for example:
data (n = 10) nested within strata (1-3), like:

data = c(0.827189928, 0.728301515, 0.514646038, 
0.600658064, 0.558828789, 0.553728503, 0.755797662, 0.736486515, 
0.773758313, 0.722529142)

strata = as.factor(c(rep(1,2),rep(2,4),rep(3,4)))

thanks for any hints,
kay


-

Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck


-- 
View this message in context: 
http://r.789695.n4.nabble.com/confidence-interval-for-small-non-iid-sample-tp2132804p2132804.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence interval around a mean count (poisson based?)

[JiHO]

Hello all,

I am observing animals in a behavioural arena and recording their
distance from a specific point at regular time intervals (large enough
so that I can assume two successive positions are independent from
each other). Each animal provides a complete histogram of distances
which reflects its trajectory in the arena. I repeat those
observations with several animals in two scenari and I want to
describe the distribution of distances in each treatment.

I computed the mean histogram per treatment: per bin, I count the
number of distances falling in the bin for each animal and then
average this count over all animals, within treatment. Now I want to
represent the variability around this average count and compute a
confidence interval.

The data is counts so, unsurprisingly, it is not normal. I have less
than 30 animals in each treatment so I cannot assume that the mean
would be normally distributed. The means indeed look
Poisson-distributed, as the counts are. I tried to find ways to
compute confidence intervals for Poisson processes but everything I
come up with in R (poi.ci in NCStats, exactci in the package of the
same name) requires integers. This makes sense for raw counts but
makes me wonder if what I am trying to do is really right with those
means (which are floats).

I understand that this is more a statistics question than a R question
but I am a bit stumped and would appreciate any help I can get from
the experts on this list. Thank you very much in advance.

PS: what I did so far was just compute mean +/- SD. The result is here:
http://dl.dropbox.com/u/1047321/hist_mean-SD.pdf
Maybe the SD is already so large that it is not even worth trying to
pursue my goal above...

Sincerely,

JiHO
---
http://maururu.net

__
R-help@r-project.org 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] confidence interval for negatively skewed, leptokurtic sample

[pappnase]

Hello,

I´ve got a statistical problem that I hope you can help me with. It doesn´t
have to do directly with R, so if there´s another forum which would suit
better, please tell me! 


Now here´s the problem:

I want to derive confidence intervals for a variable X, which is - given the
descriptive statistics - obviously negatively skewed and leptokurtic (i.e.
peaked). My aim is to make a statement similar to this one: Given certain
values of the two explaining variables (see below), X will range from value
A to value B when drawing again from the same parent population.

The dataset that I´m using is pretty huge, it contains some 150,000 cases
and it can be seen as a sample out of the basic population which covers an
entire year (the sample).
The variable I´m interested in is the prediction error X of a given (daily
computed) wind power forecast which I compute as a difference of the
prediction value and the respective realisation value. To make things
clearer: the predictions yielding my data are calculated once a day, and
they cover three days so that there are three prediction values for each
realisation value.

Unfortunately, there is autocorrelation in the dataset because the there is
data for every quarter of an hour. That´s why I have to select some cases at
random (at least I think so). Second, and more important, I want to classify
the data in order to process the available information about a dependence of
X from the two explaining variables prediction horizon and prediction
level, i.e. the level of the predited power output in relation to the
maximum power output, the latter also called nominal power or rated power.
That´s why the sample I want to analyse is reduced down to about 300 cases.

As the mean of X is unsurprisingly always close to zero, I want to gather
information about the dispersion of X as a function of the explaining
variables. A regression however doesn´t seem appropriate to me because the
resulting confidence intervals of X subject to the explaining variables
would blur a lot of information hidden in the dataset (i.e. a stronger
dispersion for daytime predictions). That´s why I thought a classification
would meet my needs best.

My first aim is now to get some information about the standard deviation or
the variance of the parent population of X. I thought about bootstrapping:
drawing various samples from the same basic population would enable me to
calculate a confidence interval for the parameter of interest, i.e. the
standard deviation. Do you think that´s a suitable approach? I´m currently
using PASW (former SPSS) which is obviously not a very powerful software,
but I have access to Stata computers, too.

Assuming that I receive a confidence interval for, say, the standard
deviation, then the next problem arises: the distribution of X is still
negatively skewed and leptokurtic, so how can I anyhow derive a confidence
interval for X? Summing and subtracting the standard deviation multiplied by
1.96 would result in a symmetric confidence interval which is probably
wrong.


It would be great if someone could help me with this. I´m not making any
progress at the moment...

Best,
Andreas
-- 
View this message in context: 
http://n4.nabble.com/confidence-interval-for-negatively-skewed-leptokurtic-sample-tp1473062p1473062.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] confidence interval for return levels

[Madan Sigdel]

Dear List
 
From a sample size of 40 I have calculated retrun level values for 
2,5,10,25,50 and 100 years using a exponentail disrtibution as follows;
 Now how can I calculate confidence interval for 
each return levels?
seeking your help. Thank you
 






 yr
    return level

2
1.980421

5
4.598394

10
6.578815

25
9.196788

50
11.17721

100
13.15763
 

 


  
[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval on parameters from optim function

[Devred, Emmanuel]

Hi everyone,

I have two questions:

I would like to get confidence intervals on the coefficients derived
from the optim() function.
I apply optim() to a given function f
 res -
optim(c(0.08,0.04,1.),f,NULL,method=L-BFGS-B,lower=c(0.,0.,0.))
And I would like to get the p-value and confidence intervals associated
with 
 res$par

My second question deals with error message. I am doing a loop with the
optim() function in it, when I get an error message like below, the loop
is stopped, however I would like to change my initial values to avoid
this error message and keep the loop going, so if it crashes when I am
away the program can still run, any idea ?

Error in optim(c(0.08, 0.04, 1), f, NULL, method = L-BFGS-B, lower =
c(0,  : 
  L-BFGS-B needs finite values of 'fn'

Thank you for any information on these two problems.

  Emmanuel

---
Dr. Emmanuel Devred
Bedford Institute of Oceanography,
1 Challenger Drive,
Dartmouth, Nova Scotia, B2Y 4A2,
Canada

Ph:  (1) 902 426-4681
Fax: (1) 902 426-9388

devr...@mar.dfo-mpo.gc.ca

http://myweb.dal.ca/edevred/
---



[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval on parameters from optim function

[Ravi Varadhan]

Hi Emmanuel,

You can obtain standard error estimate from the Hessian matrix evaluated at
the optimum as: sqrt(diag(solve(ans$hessian))), where `ans' is the object
returned by `optim'. However, `optim' does not return the Hessian matrix by
default.  So, you need to specify `hessian = TRUE' when you call `optim'.
There are two issues that I would think about:

1.  You need to pay attention to the standard regularity conditions that are
required in order for the confidence interval to be valid.  One important
condition is that the solution be in the interior of the parameter space.
Since you are using box-constraints, I would check to make sure that the
solution is not on the boundary. 

2.  The Hessian returned by `optim' is a bit inaccurate.  Often, this is not
a big deal.  There are situations, however, where this matters, as you might
get an indefinite Hessian due to inaccuracy. So, I generally prefer to use
the `hessian' function in the numDeriv package.

Ravi.


-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Devred, Emmanuel
Sent: Wednesday, August 19, 2009 9:41 AM
To: r-help@r-project.org
Subject: [R] Confidence interval on parameters from optim function

Hi everyone,

I have two questions:

I would like to get confidence intervals on the coefficients derived
from the optim() function.
I apply optim() to a given function f
 res -
optim(c(0.08,0.04,1.),f,NULL,method=L-BFGS-B,lower=c(0.,0.,0.))
And I would like to get the p-value and confidence intervals associated
with 
 res$par

My second question deals with error message. I am doing a loop with the
optim() function in it, when I get an error message like below, the loop
is stopped, however I would like to change my initial values to avoid
this error message and keep the loop going, so if it crashes when I am
away the program can still run, any idea ?

Error in optim(c(0.08, 0.04, 1), f, NULL, method = L-BFGS-B, lower =
c(0,  : 
  L-BFGS-B needs finite values of 'fn'

Thank you for any information on these two problems.

  Emmanuel

---
Dr. Emmanuel Devred
Bedford Institute of Oceanography,
1 Challenger Drive,
Dartmouth, Nova Scotia, B2Y 4A2,
Canada

Ph:  (1) 902 426-4681
Fax: (1) 902 426-9388

devr...@mar.dfo-mpo.gc.ca

http://myweb.dal.ca/edevred/
---



[[alternative HTML version deleted]]

__
R-help@r-project.org 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@r-project.org 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] Confidence interval on parameters from optim function

[Ravi Varadhan]

Emmanuel,

I didn't answer your second question.

You can use the `try' function to capture errors and keep proceeding through
simulations without crashing out:  

?try

If `L-BFGS-B' does not work well, you could try the `spg' function in the
BB package.


Ravi.

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Devred, Emmanuel
Sent: Wednesday, August 19, 2009 9:41 AM
To: r-help@r-project.org
Subject: [R] Confidence interval on parameters from optim function

Hi everyone,

I have two questions:

I would like to get confidence intervals on the coefficients derived
from the optim() function.
I apply optim() to a given function f
 res -
optim(c(0.08,0.04,1.),f,NULL,method=L-BFGS-B,lower=c(0.,0.,0.))
And I would like to get the p-value and confidence intervals associated
with 
 res$par

My second question deals with error message. I am doing a loop with the
optim() function in it, when I get an error message like below, the loop
is stopped, however I would like to change my initial values to avoid
this error message and keep the loop going, so if it crashes when I am
away the program can still run, any idea ?

Error in optim(c(0.08, 0.04, 1), f, NULL, method = L-BFGS-B, lower =
c(0,  : 
  L-BFGS-B needs finite values of 'fn'

Thank you for any information on these two problems.

  Emmanuel

---
Dr. Emmanuel Devred
Bedford Institute of Oceanography,
1 Challenger Drive,
Dartmouth, Nova Scotia, B2Y 4A2,
Canada

Ph:  (1) 902 426-4681
Fax: (1) 902 426-9388

devr...@mar.dfo-mpo.gc.ca

http://myweb.dal.ca/edevred/
---



[[alternative HTML version deleted]]

__
R-help@r-project.org 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@r-project.org 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] Confidence interval on parameters from optim function

[Arun.stat]

Regrading your second question, I guess somehow you get undefined value like
logarithm of zero of your target function for some unfortunate parameter
values in the parameter space.



Devred, Emmanuel wrote:
 
 Hi everyone,
 
 I have two questions:
 
 I would like to get confidence intervals on the coefficients derived
 from the optim() function.
 I apply optim() to a given function f
 res -
 optim(c(0.08,0.04,1.),f,NULL,method=L-BFGS-B,lower=c(0.,0.,0.))
 And I would like to get the p-value and confidence intervals associated
 with 
 res$par
 
 My second question deals with error message. I am doing a loop with the
 optim() function in it, when I get an error message like below, the loop
 is stopped, however I would like to change my initial values to avoid
 this error message and keep the loop going, so if it crashes when I am
 away the program can still run, any idea ?
 
 Error in optim(c(0.08, 0.04, 1), f, NULL, method = L-BFGS-B, lower =
 c(0,  : 
   L-BFGS-B needs finite values of 'fn'
 
 Thank you for any information on these two problems.
 
   Emmanuel
 
 ---
 Dr. Emmanuel Devred
 Bedford Institute of Oceanography,
 1 Challenger Drive,
 Dartmouth, Nova Scotia, B2Y 4A2,
 Canada
 
 Ph:  (1) 902 426-4681
 Fax: (1) 902 426-9388
 
 devr...@mar.dfo-mpo.gc.ca
 
 http://myweb.dal.ca/edevred/
 ---
 
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/Confidence-interval-on-parameters-from-optim-function-tp25044388p25046772.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence interval on parameters from optim function

[Greg Snow]

If you f function is a likelihood (or better a log likelihood), then you can do 
profiling.  Take the parameter of interest and change it a small amount from 
the optimal value, rerun optim with this value fixed and let it optimize over 
everything else.  Repeat this for several values and see how the likelihood 
changes.  You can use this information to construct the confidence interval.  
One way is by using the chi-square approximation to log likelihood differences. 
 There are already some functions that use this procedure and you could either 
use those or look at the code to see how to do it yourself.  The mle, confint, 
and profile functions in the stats4 package are one set.

If your f function is not a likelihood or log likelihood, then we will need 
more information on what you are trying to accomplish. 

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Devred, Emmanuel
 Sent: Wednesday, August 19, 2009 7:41 AM
 To: r-help@r-project.org
 Subject: [R] Confidence interval on parameters from optim function
 
 Hi everyone,
 
 I have two questions:
 
 I would like to get confidence intervals on the coefficients derived
 from the optim() function.
 I apply optim() to a given function f
  res -
 optim(c(0.08,0.04,1.),f,NULL,method=L-BFGS-B,lower=c(0.,0.,0.))
 And I would like to get the p-value and confidence intervals associated
 with
  res$par
 
 My second question deals with error message. I am doing a loop with the
 optim() function in it, when I get an error message like below, the
 loop
 is stopped, however I would like to change my initial values to avoid
 this error message and keep the loop going, so if it crashes when I am
 away the program can still run, any idea ?
 
 Error in optim(c(0.08, 0.04, 1), f, NULL, method = L-BFGS-B, lower =
 c(0,  :
   L-BFGS-B needs finite values of 'fn'
 
 Thank you for any information on these two problems.
 
   Emmanuel
 
 ---
 Dr. Emmanuel Devred
 Bedford Institute of Oceanography,
 1 Challenger Drive,
 Dartmouth, Nova Scotia, B2Y 4A2,
 Canada
 
 Ph:  (1) 902 426-4681
 Fax: (1) 902 426-9388
 
 devr...@mar.dfo-mpo.gc.ca
 
 http://myweb.dal.ca/edevred/
 ---
 
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org 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@r-project.org 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] Confidence interval?

[Dieter Menne]

Thomas Seth Davis Thomas.Davis at nau.edu writes:

  I need help fitting/plotting a confidence interval to a frequency
distribution
 

In many medical journals, reviewers only want to see some error bars.
In 90% of the cases, these are wrong or misleading, but it is hopeless
to argue with medical professors playing statisticians. So I give them
the square-root of the counts, which is not much worse than anything else,
but I do not claim these to be confidence intervals.
Have some good excuse ready when some counts are low.

If you are working in a more serious field, you might consider giving
something like estimates of the parameters, but you cannot plot these
in your curves easily. Try fitdistr in package MASS for a starter.
 
  
  You must know your password to change your options (including changing
  the password, itself) or to unsubscribe.  It is:
  
davis234

This is really a great password, which needed 0.03 seconds on my 
computer to crack.

Dieter

__
R-help@r-project.org 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] Confidence interval?

[Thomas Seth Davis]

hi folks, 

 I need help fitting/plotting a confidence interval to a frequency 
distribution

 Can someone help with this? 

thanks, 

 tsd 



-Original Message-

 Date: Mon Apr 06 15:08:20 MST 2009
 From: r-help-requ...@r-project.org
 Subject: Welcome to the R-help mailing list
 To: t...@nau.edu

 Welcome to the R-help@r-project.org mailing list!
 
 To post to this list, send your email to:
 
   r-help@r-project.org
 
 General information about the mailing list is at:
 
   https://stat.ethz.ch/mailman/listinfo/r-help
 
 If you ever want to unsubscribe or change your options (eg, switch to
 or from digest mode, change your password, etc.), visit your
 subscription page at:
 
   https://stat.ethz.ch/mailman/options/r-help/tsd3%40nau.edu
 
 You can also make such adjustments via email by sending a message to:
 
   r-help-requ...@r-project.org
 
 with the word `help' in the subject or body (don't include the
 quotes), and you will get back a message with instructions.
 
 You must know your password to change your options (including changing
 the password, itself) or to unsubscribe.  It is:
 
   davis234
 
 Normally, Mailman will remind you of your r-project.org mailing list
 passwords once every month, although you can disable this if you
 prefer.  This reminder will also include instructions on how to
 unsubscribe or change your account options.  There is also a button on
 your options page that will email your current password to you.

__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear regression

[Thomas Lumley]

On Mon, 23 Mar 2009, Kevin J Emerson wrote:


 Now, it is simple enough to calculate the x-intercept itself ( - intercept / 
 slope ), but it is a whole separate process to generate the confidence 
 interval of it.  I can't figure out how to propagate the error of the slope 
 and intercept into the ratio of the two.  The option option I have tried to 
 figure out is to use the predict function to look for where the confidence 
 intervals cross the axis but this hasn't been too fruitful either.
 


A few people have written code to do nonlinear transformations automatically. 
There's one in my 'survey' package, and I think there's an example in Venables 
 Ripley's S Programming, and I'm sure there are others.

For example
 library(survey)
 m-lm(Dev.Rate~Temperature, data=t)
 svycontrast(m, quote(-`(Intercept)`/Temperature))
  nlcon SE
contrast 3.8061  0.1975

The `backquotes` are necessary because the coefficient name, (Intercept), isn't 
a legal variable name.

The svycontrast() function works on anything that has coef() and vcov() 
methods. The actual delta-method computations are in survey:::nlcon

-thomas

Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear regression

[Thomas Lumley]

I should probably also add the health warning that the delta-method approach to 
ratios of coefficients is REALLY BAD unless the denominator is well away from 
zero.

When the denominator is near zero nothing really works (Fieller's method is 
valid but only because 'valid' means less than you might like for a confidence 
interval method, and the bootstrap is useful but not necessarily valid)

 -thomas

Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear

[Ted Harding]

On 24-Mar-09 03:31:32, Kevin J Emerson wrote:
 Hello all,
 
 This is something that I am sure has a really suave solution in R,
 but I can't quite figure out the best (or even a basic) way to do it.
 
 I have a simple linear regression that is fit with lm for which I
 would like to estimate the x intercept with some measure of error
 around it (confidence interval). In biology, there is the concept
 of a developmental zero - a temperature under which development
 will not happen. This is often estimated by extrapolation of a
 curve of developmental rate as a function of temperature. This
 developmental zero is generally reported without error. I intend
 to change this! 
 There has to be some way to assign error to this term, I just have
 yet to figure it out.
 
 Now, it is simple enough to calculate the x-intercept itself ( -
 intercept / slope ), but it is a whole separate process to generate
 the confidence interval of it.  I can't figure out how to propagate
 the error of the slope and intercept into the ratio of the two.
 The option option I have tried to figure out is to use the predict
 function to look for where the confidence intervals cross the axis
 but this hasn't been too fruitful either.  
 
 I would greatly appreciate any insight you may be able to share.
 
 Cheers,
 Kevin
 
 Here is a small representative sample of some of my data where
 Dev.Rate ~ Temperature.

The difficulty with the situation you describe is that you are up
against what is known as the Calibration Problem in Statistics.
Essentially, the difficulty is that in the theory of the linear
regression there is sufficient probability for the slope to be
very close to zero -- and hence that the X-intercept maybe
very far out to the left or the right -- that the whole dconcept
of Standard error ceases to exist. Even the expected position of
the intercept does not exist. And this is true regardless of the
data (the one exception would be where it is absolutely certain
that the data will lie exactly on a straight line).

This has been addressed a number of times, and controversally,
in the statistical literature. One approach which I like (having
thought of it myself :)) is the use of likelihood-based confidence
intervals.

Given a linear regression

  Y = a + b*X + E

(where the error E has N(0, sig^2) distribution), suppose you
want to estimate the value X0 of X for which Y = Y0, a given value
(in your case Y0 = 0). For simplicity, measure X and Y from their
means sum(X)/N and sum(Y)/N.

The approach is based on comparing the likelihoods

[A] for unresticted ML estimation of a, b and sig
-- a.hat, b.hat, sig.hat
(where a.hat = 0 because of the origins of X and Y)
[B] for ML estimation assuming a particular value X1 for X0
-- a.tilde, b.tilde, sig.tilde

where
[A] a.hat = 0 (as above),
b.hat = (sum(X*Y))/(sum(X^2)
X0.hat= Y0/b.hat [ = -mean(Y)/b.hat in your case ]
sig.hat^2 = (sum(Y - b.hat*X)^2)/(N+1)

[B] a.tilde = (sum(X^2) - X0*sum(X*Y))/D
b.tilde = ((N+1)*sum(X*Y) + N*X1*Y0)/D
sig.hat.tilde^2
= (sum((Y - a.tilde - b.tilde*X)^2)
+ (Y0 - a.tilde - b.tilde*X1)^2)/(N+1)

where D = (N+1)*sum(X^2 + N*X1^2

Then let R(X1) = (sig.hat^2)/(sig.tilde^2)

A confidence interval for X0 is the set of values of X1 such
that R(X1)  R0 say, where Prob(R(X0)R0) = P, the desired
confidence level, when X0 is the true value.

It can be established that

  (N-2)*(1 - R(X0))/R(X0)

has the F distribution with 1 and (N-1) degrees of freedom.
Since this is independent of X0, the resulting set of values
of X1 constitute a confidence interval.

This was described in

  Harding, E.F. (1986). Modelling: the classical approach.
  The Statistician vol. 35, pp. 115-134
  (see pages 123-126)

and has been later referred to by P.J. Brown (thought I don't
have the references to hand just now).

When I have time for it (not today) I'll see if I can implement
this neatly in R. It's basically a question of solving

  (N-2)*(1 - R(X0))/R(X0) = qf(P,1,(N-1))

for X0 (two solutions, maybe one, if any exist). However, it is
quite possible that no solution exists (depending on P), so that
the confidence interval would then be (-inf,+inf); or, when only
one exists, it will be either (-inf,X0) or (X0,inf).
These possibilities of infinite intervals are directly related
to the possibility that, at significance level alpha = (1-P),
the estimated slope may not differ significantly from 0.

Maybe more later!
Ted.


E-Mail: (Ted Harding) ted.hard...@manchester.ac.uk
Fax-to-email: +44 (0)870 094 0861
Date: 24-Mar-09   Time: 10:04:30
-- XFMail --

__
R-help@r-project.org 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 

[R] CONFIDENCE INTERVAL FOR GLMER MODEL

[LiatL]

I've built a poisson regression model for multiple subjects by using the
GLMER function. I've also developed some curves for defining its limits but
I did not succeed in developing confidence interval for the model's curve
(confint or predict does not work - only for glm).
Does anyone know how can I produce confidence interva for a glmer model?
I'll appriciate any help...
Liat
-- 
View this message in context: 
http://www.nabble.com/CONFIDENCE-INTERVAL-FOR-GLMER-MODEL-tp22677543p22677543.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] CONFIDENCE INTERVAL FOR GLMER MODEL

[Douglas Bates]

On Tue, Mar 24, 2009 at 5:25 AM, LiatL liat.lam...@gmail.com wrote:

 I've built a poisson regression model for multiple subjects by using the
 GLMER function. I've also developed some curves for defining its limits but
 I did not succeed in developing confidence interval for the model's curve
 (confint or predict does not work - only for glm).
 Does anyone know how can I produce confidence interva for a glmer model?
 I'll appriciate any help...

You would need to be more explicit about what you mean by the model's
curve and what the confidence interval would represent, in the sense
of what components of the variability would be incorporated in the
confidence interval, before we could answer such a question.

May I suggest that the discussion be moved to the
r-sig-mixed-mod...@r-project.org mailing list, which I have cc:d on
this reply?

__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear

[Peter Dalgaard]

(Ted Harding) wrote:
 On 24-Mar-09 03:31:32, Kevin J Emerson wrote:
...
 When I have time for it (not today) I'll see if I can implement
 this neatly in R. It's basically a question of solving
 
   (N-2)*(1 - R(X0))/R(X0) = qf(P,1,(N-1))
 
 for X0 (two solutions, maybe one, if any exist). 
etc.

A quick and probably not-too-dirty way is to rewrite it as a nonlinear
model:

x - 1:10
Y - 2*x - 3 + rnorm(10, sd=.1)
cf - coef(lm(Y~x))
confint(nls(Y~beta*(x-x0), start=c(beta=cf[[2]],x0=-cf[[1]]/cf[[2]])))




-- 
   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
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear

[Kevin J Emerson]

I wanted to send out a quick thanks to all that replied to my query about 
estimating the confidence interval around the x-intercept of a linear 
regression.  The method I was able to implement in the most straightforward way 
was taken from Section 3.2 of Draper and Smith (1998). Applied Regression 
Analysis. The code follows - its not the most fancy code, but it gets the job 
done.  I will see if I can work out the other suggestions that were made and 
see how they all turn out.  

xInterceptCI - function(x, alpha = 0.05, ...) {
  intercept - coef(x)[1]  
  slope - coef(x)[2]  
  meanX - mean(x$model[,2])   
  n - length(x$model[,2]) 
  tstar - qt(alpha/2, n-2)
  sxx - sum(x$model[,2]^2) - sum(x$model[,2])^2 / n
  SSresidual - (1-cor(x$model[,1], x$model[,2])^2) * 
(sum(x$model[,1]^2)-sum(x$model[,1])^2/n)
  S - sqrt(SSresidual/(n-2))
  SEslope - S / sqrt(sxx)
  Xintercept - - intercept / slope
  y0 - 0
  g - (tstar / (slope/SEslope))^2
  left - (Xintercept - meanX) * g
  bottom - 1 - g
  Right - (tstar * S / slope) * sqrt( ((Xintercept - meanX)^2/sxx) + bottom/n)
  lower - Xintercept + (left + Right) / bottom
  upper - Xintercept + (left - Right) / bottom
  return(c(lower,upper))
}

On Tue, 24 Mar 2009 17:16:41 +0100, Peter Dalgaard p.dalga...@biostat.ku.dk 
wrote:
 (Ted Harding) wrote:
  On 24-Mar-09 03:31:32, Kevin J Emerson wrote:
 ...
  When I have time for it (not today) I'll see if I can implement
  this neatly in R. It's basically a question of solving
  
(N-2)*(1 - R(X0))/R(X0) = qf(P,1,(N-1))
  
  for X0 (two solutions, maybe one, if any exist). 
 etc.
 
 A quick and probably not-too-dirty way is to rewrite it as a nonlinear
 model:
 
 x - 1:10
 Y - 2*x - 3 + rnorm(10, sd=.1)
 cf - coef(lm(Y~x))
 confint(nls(Y~beta*(x-x0), start=c(beta=cf[[2]],x0=-cf[[1]]/cf[[2]])))
 
 
 
 
 -- 
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
 ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907
 


--
==
==
Kevin J Emerson
Bradshaw-Holzapfel Lab
Center for Ecology and Evolutionary Biology
1210 University of Oregon
Eugene, Oregon 97403
kemer...@uoregon.edu
__
R-help@r-project.org 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] confidence interval or error of x intercept of a linear regression

[Kevin J Emerson]

Hello all,

This is something that I am sure has a really suave solution in R, but I can't 
quite figure out the best (or even a basic) way to do it.

I have a simple linear regression that is fit with lm for which I would like to 
estimate the x intercept with some measure of error around it (confidence 
interval).  In biology, there is the concept of a developmental zero - a 
temperature under which development will not happen. This is often estimated by 
extrapolation of a curve of developmental rate as a function of temperature.  
This developmental zero is generally reported without error.  I intend to 
change this!  There has to be some way to assign error to this term, I just 
have yet to figure it out.

Now, it is simple enough to calculate the x-intercept itself ( - intercept / 
slope ), but it is a whole separate process to generate the confidence interval 
of it.  I can't figure out how to propagate the error of the slope and 
intercept into the ratio of the two.  The option option I have tried to figure 
out is to use the predict function to look for where the confidence intervals 
cross the axis but this hasn't been too fruitful either.  

I would greatly appreciate any insight you may be able to share.

Cheers,
Kevin

Here is a small representative sample of some of my data where Dev.Rate ~ 
Temperature.

t -
structure(list(Temperature = c(5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 9, 9, 10.5, 
10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 10.5, 
10.5, 10.5, 10.5, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
12, 12, 12, 12, 12, 12, 12, 14, 14, 14, 14, 16, 16, 16, 16, 16, 
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 18, 18, 
18, 18, 18, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 
20, 20, 22, 22), Dev.Rate = c(0.007518797, 0.007518797, 0.007518797, 
0.007194245, 0.007194245, 0.007194245, 0.007194245, 0.007194245, 
0.007194245, 0.006896552, 0.006896552, 0.012820513, 0.012820513, 
0.012195122, 0.012195122, 0.012195122, 0.012195122, 0.011363636, 
0.011363636, 0.011363636, 0.011363636, 0.011363636, 0.011363636, 
0.011363636, 0.010869565, 0.00952381, 0.00952381, 0.015151515, 
0.015151515, 0.022727273, 0.022727273, 0.022727273, 0.022727273, 
0.022727273, 0.022727273, 0.022727273, 0.022727273, 0.022727273, 
0.022727273, 0.022727273, 0.022727273, 0.02083, 0.02083, 
0.02083, 0.034482759, 0.029411765, 0.029411765, 0.029411765, 
0.029411765, 0.029411765, 0.029411765, 0.029411765, 0.029411765, 
0.029411765, 0.029411765, 0.029411765, 0.029411765, 0.029411765, 
0.029411765, 0.029411765, 0.027027027, 0.025, 0.038461538, 0.03030303, 
0.03030303, 0.03030303, 0.052631579, 0.052631579, 0.045454545, 
0.045454545, 0.045454545, 0.045454545, 0.045454545, 0.045454545, 
0.045454545, 0.045454545, 0.045454545, 0.038461538, 0.038461538, 
0.038461538, 0.038461538, 0.038461538, 0.038461538, 0.038461538, 
0.038461538, 0.047619048, 0.047619048, 0.047619048, 0.047619048, 
0.047619048, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 
0.0625, 0.0625, 0.0625, 0.0625, 0.052631579, 0.052631579, 0.052631579, 
0.052631579, 0.052631579, 0.076923077, 0.071428571)), .Names = c(Temperature, 
Dev.Rate), class = data.frame, row.names = c(NA, 107L))



--
==
==
Kevin J Emerson
Bradshaw-Holzapfel Lab
Center for Ecology and Evolutionary Biology
1210 University of Oregon
Eugene, Oregon 97403
kemer...@uoregon.edu

__
R-help@r-project.org 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] confidence interval for roc curves

[marc bernard]

Dear all,
 
I am looking for an R package that allows me to calculate and plot the 
confidence limits for the roc curve using for example some bootstrapping.
 I tried ROCR who seems doing such work but i couldn't find the right option 
to do it.
 
Many thanks
 
Bests
 
Marc
 
_


tarted! 

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval

[Yawar Amin]

rak1304 rkeyes87 at hotmail.com writes:
   I am new to R and Im 
some trouble with the following question... 

I'm starting to study stats and R again after almost a year, so I 
thought this is interesting. I think I have the answer. Here is how I 
arrived at it: 

  Generate 100 standard normal N(0,1) samples of size 100, 
X1(k),...,X100(k)  where k=1,...,100 (The k is and indicie in brackets) 
  Calculate the sample mean for each sample.   For each sample mean 
Xbark the 0.95-confidence interval for the mean mew=0  is given by... 
  Ik= ( Xbark plus or minus 1.96/10)   Find the number of intervals 
such that 0 does not belong to Ik. How many of  them do you expect to 
see?   Well so far I have come up with...   N-100; Nsamp-100  
A-matrix(rnorm(N*Nsamp,0,1),ncol=Nsamp)  means-apply(A,2,mean) 

This looks fine. Now you have to figure out I_k, the confidence 
interval. According to your formula for this, the lowest boundary in I_k 
must be XBar_k - 0.196, and the highest boundary XBar_k + 0.196. Now you 
have to figure out how many times 0 falls outside this interval. In R 
terms, XBar_k is your `means' variable. So you want to know how many 
times means - 0.196  0 or means + 0.196  0. 

A previous poster said that sum() counts the number of times something 
occurs if you're dealing with boolean data i.e. 1s and 0s, or TRUEs and 
FALSEs. R gives you these TRUE/FALSE values. E.g. if you have a variable 
w = 2, running w  0 will cause R to respond with TRUE. Same goes with 
vector data i.e. if w = c(0,1,2,3,4), then running w  0 will cause R to 
give you, this time, a vector or TRUEs and FALSEs: FALSE TRUE TRUE TRUE 
TRUE. Then, running sum(w0) will count only the TRUEs and give you 4. 
(Actually, R pretends it's a vector of 1s and 0s and sums it up: sum(0 1 
1 1 1) is 4.) 

  However I have no idea what I am doing and no idea if that even 
makes sense.  Any help would be greatly appreciated as I have no 
experience of statistical  software whatsoever. 

I hope this gives you some ideas as to how to turn your textbook 
statistical concepts into R language commands. Thanks to R's vector 
goodness, it's usually simpler than it seems--you don't need to 
overthink it. 

  Thanks in Advance.   Rachel 

Yawar

__
R-help@r-project.org 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] Confidence Interval

[rak1304]

I am new to R and Im some trouble with the following question...

Generate 100 standard normal N(0,1) samples of size 100, X1(k),...,X100(k) 
where k=1,...,100 (The k is and indicie in brackets)

Calculate the sample mean for each sample.

For each sample mean Xbark the 0.95-confidence interval for the mean mew=0
is given by...

Ik= ( Xbark plus or minus 1.96/10)

Find the number of intervals such that 0 does not belong to Ik.  How many of
them do you expect to see?

Well so far I have come up with...

N-100; Nsamp-100
A-matrix(rnorm(N*Nsamp,0,1),ncol=Nsamp)
means-apply(A,2,mean)

However I have no idea what I am doing and no idea if that even makes sense.
Any help would be greatly appreciated as I have no experience of statistical
software whatsoever. 

Thanks in Advance.

Rachel
-- 
View this message in context: 
http://www.nabble.com/Confidence-Interval-tp21521610p21521610.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] Confidence Interval

[Rachel Keyes]

I am new to R and Im some trouble with the following question... Generate 100 
standard normal N(0,1) samples of size 100, X1(k),...,X100(k)  where 
k=1,...,100 (The k is and indicie in brackets) Calculate the sample mean for 
each sample. For each sample mean Xbark the 0.95-confidence interval for the 
mean mew=0 is given by... Ik= ( Xbark plus or minus 1.96/10) Find the number of 
intervals such that 0 does not belong to Ik.  How many of them do you expect to 
see? Well so far I have come up with... N-100; Nsamp-100 
A-matrix(rnorm(N*Nsamp,0,1),ncol=Nsamp) means-apply(A,2,mean) However I have 
no idea what I am doing and no idea if that even makes sense. Any help would be 
greatly appreciated as I have no experience of statistical software whatsoever. 
Thanks in Advance. Rachel


_
Choose the perfect PC or mobile phone for you

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval

[rak1304]

Yes it is a homework problem, I included the whole question as I thought it
would make it easier to explain however I am unsure of how to do the
confidence interval part.  As far as I am aware I have set up a matrix with
my 100 samples of 100 and have calculated means.  Do I need to set up a new
matrix of these means and then do a confidence interval from there? Also
when I searched 'confint' in R I am confused as to how to adapt it to my
problem.

I am looking for the interval of each individual mean plus or minus 1.96/10

I wasnt aware of abs(x) and sum(means0) so thankyou for pointing them out
to me.




Ben Bolker wrote:
 
 Rachel Keyes rkeyes87 at hotmail.com writes:
 
 
 
 I am new to R and Im some trouble with the following question... Generate
 100
 standard normal N(0,1) samples
 of size 100, X1(k),...,X100(k)  where k=1,...,100 (The k is and indicie
 in
 brackets) Calculate the sample
 mean for each sample. For each sample mean Xbark the 0.95-confidence
 interval
 for the mean mew=0 is given
 by... Ik= ( Xbark plus or minus 1.96/10) Find the number of intervals
 such
 that 0 does not belong to Ik.  How
 many of them do you expect to see? Well so far I have come up with...
 N-100;
 Nsamp-100
 A-matrix(rnorm(N*Nsamp,0,1),ncol=Nsamp) means-apply(A,2,mean) However I
 have
 no idea what I am
 doing and no idea if that even makes sense. Any help would be greatly
 appreciated as I have no experience of
 statistical software whatsoever. Thanks in Advance. Rachel
 
 
   This sounds an awful lot like a homework problem.
 Can you convince us not by giving us a little more context?
 
 I will say that you seem to have made a pretty good start.
 A hint is that if you have a logical condition,
 sum(condition) will count the number of cases. For
 example, sum(means0) will count the number of positive
 means.  ?abs and ? may help too.
 
  By the way, I think that's supposed to be mu and not mew.
 
   Ben Bolker
 
 __
 R-help@r-project.org 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/Confidence-Interval-tp21523397p21523829.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org 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] confidence interval for glm

[Gavin Simpson]

On Sun, 2008-11-30 at 01:15 -0500, David Winsemius wrote:
 ?confint.glm   # ... in MASS

That provides confidence intervals on the parameters of the model, which
is not what the OP wanted. He wants confidence intervals on:

predict(mod, newdat)

One way to do this is to compute them in the normal way but do so on the
scale of the link function. Then transform those confidence intervals on
to the scale of the response by applying the inverse of the link
function used in the GLM. One way to do this is the following: [not
tested as no reproducible code to try it on]

## compute the predicted values on the scale of the link
preds - predict(mod, newdata = pred.data, se.fit = TRUE)

## where mod and pred.data are your glm model and 
## points at which you wish to evaluate the fitted function.
## leave out newdata = pred.data if you want the fitted values

## get inverse of link as a function
ilogit - family(mod)$linkinv

## compute standard confidence intervals (2 * se.fit)
ci.u - with(preds, fit + (2 * se.fit))
ci.l - with(preds, fit - (2 * se.fit))
## change 2 to be a quantile from the t distribution for your level
## of confidence (2 is approximately 95% for large n)
##
## something like for 95%, N is sample size;
## (hope I got the correction for df right here, no doubt someone will
## correct me if wrong), and use crit.t in place of 2 above.
crit.t - qt(0.975, df = N-1)

## compute predictions on scale of response by applying ilogit
preds - ilogit(preds$fit)

## transform conf int on to scale of response
ci.u - ilogit(ci.u)
ci.l - ilogit(ci.l)

Is that what you were looking for?

G

 
 On Nov 28, 2008, at 9:29 AM, Gerard M. Keogh wrote:
 
 
  Hi all,
 
  simple Q:
 
  how do I extract the upper and lower CI for predicted probabilities
  directly for a glm - I'm sure there's a one line to do it but I  
  can't find
  it.
  the predicted values I get with the predict (.. response)
 
  Thanks
 
  Gerard
 
 
  **
  The information transmitted is intended only for the person or  
  entity to which it is addressed and may contain confidential and/or  
  privileged material. Any review, retransmission, dissemination or  
  other use of, or taking of any action in reliance upon, this  
  information by persons or entities other than the intended recipient  
  is prohibited. If you received this in error, please contact the  
  sender and delete the material from any computer.  It is the policy  
  of the Department of Justice, Equality and Law Reform and the  
  Agencies and Offices using its IT services to disallow the sending  
  of offensive material.
  Should you consider that the material contained in this message is  
  offensive you should contact the sender immediately and also  
  mailminder[at]justice.ie.
 
  Is le haghaidh an duine nó an eintitis ar a bhfuil sí dírithe, agus  
  le haghaidh an duine nó an eintitis sin amháin, a bheartaítear an  
  fhaisnéis a tarchuireadh agus féadfaidh sé go bhfuil ábhar faoi rún  
  agus/nó faoi phribhléid inti. Toirmisctear aon athbhreithniú,  
  atarchur nó leathadh a dhéanamh ar an bhfaisnéis seo, aon úsáid eile  
  a bhaint aisti nó aon ghníomh a dhéanamh ar a hiontaoibh, ag daoine  
  nó ag eintitis seachas an faighteoir beartaithe. Má fuair tú é seo  
  trí dhearmad, téigh i dteagmháil leis an seoltóir, le do thoil, agus  
  scrios an t-ábhar as aon ríomhaire. Is é beartas na Roinne Dlí agus  
  Cirt, Comhionannais agus Athchóirithe Dlí, agus na nOifígí agus na  
  nGníomhaireachtaí a úsáideann seirbhísí TF na Roinne, seoladh ábhair  
  cholúil a dhícheadú.
  Más rud é go measann tú gur ábhar colúil atá san ábhar atá sa  
  teachtaireacht seo is ceart duit dul i dteagmháil leis an seoltóir 
  láithreach agus le mailminder[ag]justice.ie chomh maith.
  ***
 
 
 
  __
  R-help@r-project.org 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@r-project.org 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@r-project.org 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] confidence interval for glm

[David Winsemius]

?confint.glm   # ... in MASS

On Nov 28, 2008, at 9:29 AM, Gerard M. Keogh wrote:



Hi all,

simple Q:

how do I extract the upper and lower CI for predicted probabilities
directly for a glm - I'm sure there's a one line to do it but I  
can't find

it.
the predicted values I get with the predict (.. response)

Thanks

Gerard


**
The information transmitted is intended only for the person or  
entity to which it is addressed and may contain confidential and/or  
privileged material. Any review, retransmission, dissemination or  
other use of, or taking of any action in reliance upon, this  
information by persons or entities other than the intended recipient  
is prohibited. If you received this in error, please contact the  
sender and delete the material from any computer.  It is the policy  
of the Department of Justice, Equality and Law Reform and the  
Agencies and Offices using its IT services to disallow the sending  
of offensive material.
Should you consider that the material contained in this message is  
offensive you should contact the sender immediately and also  
mailminder[at]justice.ie.


Is le haghaidh an duine nó an eintitis ar a bhfuil sí dírithe, agus  
le haghaidh an duine nó an eintitis sin amháin, a bheartaítear an  
fhaisnéis a tarchuireadh agus féadfaidh sé go bhfuil ábhar faoi rún  
agus/nó faoi phribhléid inti. Toirmisctear aon athbhreithniú,  
atarchur nó leathadh a dhéanamh ar an bhfaisnéis seo, aon úsáid eile  
a bhaint aisti nó aon ghníomh a dhéanamh ar a hiontaoibh, ag daoine  
nó ag eintitis seachas an faighteoir beartaithe. Má fuair tú é seo  
trí dhearmad, téigh i dteagmháil leis an seoltóir, le do thoil, agus  
scrios an t-ábhar as aon ríomhaire. Is é beartas na Roinne Dlí agus  
Cirt, Comhionannais agus Athchóirithe Dlí, agus na nOifígí agus na  
nGníomhaireachtaí a úsáideann seirbhísí TF na Roinne, seoladh ábhair  
cholúil a dhícheadú.
Más rud é go measann tú gur ábhar colúil atá san ábhar atá sa  
teachtaireacht seo is ceart duit dul i dteagmháil leis an seoltóir 
láithreach agus le mailminder[ag]justice.ie chomh maith.

***



__
R-help@r-project.org 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@r-project.org 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] confidence interval for glm

[Gerard M. Keogh]

Hi all,

simple Q:

how do I extract the upper and lower CI for predicted probabilities
directly for a glm - I'm sure there's a one line to do it but I can't find
it.
the predicted values I get with the predict (.. response)

Thanks

Gerard


**
The information transmitted is intended only for the person or entity to which 
it is addressed and may contain confidential and/or privileged material. Any 
review, retransmission, dissemination or other use of, or taking of any action 
in reliance upon, this information by persons or entities other than the 
intended recipient is prohibited. If you received this in error, please contact 
the sender and delete the material from any computer.  It is the policy of the 
Department of Justice, Equality and Law Reform and the Agencies and Offices 
using its IT services to disallow the sending of offensive material.
Should you consider that the material contained in this message is offensive 
you should contact the sender immediately and also mailminder[at]justice.ie.

Is le haghaidh an duine nó an eintitis ar a bhfuil sí dírithe, agus le haghaidh 
an duine nó an eintitis sin amháin, a bheartaítear an fhaisnéis a tarchuireadh 
agus féadfaidh sé go bhfuil ábhar faoi rún agus/nó faoi phribhléid inti. 
Toirmisctear aon athbhreithniú, atarchur nó leathadh a dhéanamh ar an 
bhfaisnéis seo, aon úsáid eile a bhaint aisti nó aon ghníomh a dhéanamh ar a 
hiontaoibh, ag daoine nó ag eintitis seachas an faighteoir beartaithe. Má fuair 
tú é seo trí dhearmad, téigh i dteagmháil leis an seoltóir, le do thoil, agus 
scrios an t-ábhar as aon ríomhaire. Is é beartas na Roinne Dlí agus Cirt, 
Comhionannais agus Athchóirithe Dlí, agus na nOifígí agus na nGníomhaireachtaí 
a úsáideann seirbhísí TF na Roinne, seoladh ábhair cholúil a dhícheadú.
Más rud é go measann tú gur ábhar colúil atá san ábhar atá sa teachtaireacht 
seo is ceart duit dul i dteagmháil leis an seoltóir láithreach agus le 
mailminder[ag]justice.ie chomh maith. 
***



__
R-help@r-project.org 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] Confidence interval for binomial variance

[Robert A. LaBudde]

Based on simulations, I've come up with a simple function to compute 
the confidence interval for the variance of the binomial variance, 
where the true variance is


v   = rho*(1-rho)/n

where rho = true probability of success and n = # of trials.

For x = # successes observed in n trials, p = x / n as usual.

For p  0.25 or p  0.75, I use the proportion-based transformed 
confidence interval, as suggested by Moshe Olshansky. I used the 
Wilson interval, but some might prefer the Blaker interval. For p = 
0.25 or p = 0.75, I use the standard chi-square based confidence 
interval (which is very conservative for this problem in this range). 
The composite method gives reliable results over a wide range of rho.


This should suit the purpose until someone comes up with a more 
theoretically sound method. Bootstrap is not reliable for n  50, at least.


#09.26.08 02.50 binomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:

require('binGroup')

binomVarCI- function (n, x, conf=0.95) {
  p- x/n #proportion
  if (p0.25 | p0.75 | x==0 | x==n) {  #use proportion-based CI
pCI- binWilson(n, x, conf.level=conf)  #CI for proportion
vCI- sort(c(pCI[1]*(1-pCI[1])/(n-1), pCI[2]*(1-pCI[2])/(n-1)))
  } else {  #use chi-square-based CI
phiL- qchisq(0.025, n-1)/(n-1)
phiU- qchisq(0.975, n-1)/(n-1)
#vest- p*(1-p)/(n-1))  #variance estimate
vCI- c(vest/phiU, vest/phiL) #chi-square-based
  }
  return (vCI)
}

Here is a test program to measure coverage:

#09.26.08 03.10 tbinomVarCI.r
#copyright 2008 by Robert A LaBudde, all rights reserved
#test of CI for binomial sample variance
#created: 09.26.08 by r.a. labudde
#changes:

nReal - 1000

for (POD in c(0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.50)) {
  cat('\nPOD: ', sprintf('%8.4f', POD), '\n')
  for (nRepl in c(6, 12, 20, 50)) {
vtrue- POD*(1-POD)/nRepl
pcover- 0
for (iReal in 1:nReal) {
  x- rbinom(1, nRepl, POD)
  vCI- binomVarCI(nRepl, x)
  #vest- (x/nRepl)*(1 - x/nRepl)/(nRepl-1)
  if (vtrue = vCI[1]  vtrue= vCI[2]) pcover- pcover + 1
} #iReal
pcover- pcover/nReal
cat('n: ', sprintf('%4i', nRepl), ' Coverage: ', 
sprintf('%8.4f', pcover), '\n')

  } #nRepl
} #POD

Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: [EMAIL PROTECTED]
Least Cost Formulations, Ltd.URL: http://lcfltd.com/
824 Timberlake Drive Tel: 757-467-0954
Virginia Beach, VA 23464-3239Fax: 757-467-2947

Vere scire est per causas scire

__
R-help@r-project.org 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] confidence interval

[Sudheer Joseph]

I am a new user of R package and I needed a help
   regarding plotting confidence limits for FFT ( chisquare based) and
spectrum
   using the monte carlo simulations of background red noise . How do I do
it
   in R can any one help in this issue?
   sudheer

-- 
**
Dr. Sudheer Joseph
Scientist
Indian National Centre for Ocean Information Services (INCOIS)
Ocean Valley, Post Box No# 21,
IDA Jeedimetla P.O.
Hyderabad, Ranga Reddy District - 500 055
Andhra Pradesh, India.
TEl:+91-40-23044600(R),Tel:+91-9440832534(Mobile)
Tel:+91-40-23886047(O),Fax:+91-40-23895011(O)
E-mail:[EMAIL PROTECTED] [EMAIL PROTECTED];
[EMAIL PROTECTED]; [EMAIL PROTECTED]
Web- http://oppamthadathil.tripod.com
--* ---
The ultimate measure of a man is
not where he stands in moments of
comfort and convenience, but where
he stands at times of challenge and
controversy.
Martin Luther King, Jr.

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval

[Bernardo Rangel Tura]

Em Ter, 2008-08-19 às 23:25 -0300, Raphael Saldanha escreveu:
 Hi!
 
 With the following script, I'm trying to make a demonstration of a
 Confidence Interval, but I'm observing some differences on tails.

Raphael,

If you make demonstration of Confidence Interval why you don't use
ci.examp of TeachingDemos package?

It's a very good example


-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil

__
R-help@r-project.org 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] Confidence Interval

[Raphael Saldanha]

Thanks Professor and Bernardo,

What I'm trying to do is this: I have a macro for Minitab. His author says
it's a Monte Carlo simulation to estimate a confidence interval. But and
don't have Minitab, don't like to work with illegal licenses, and LOVE R.
So, I re-write the macro in a R script.

Bellow, is the original macro, for Minitab:

GMACRO

TesteMedia

do K1=1:1000
 random 200 C1;
 normal 100 2.
 mean C1 K2
 let C2(K1)=K2
enddo

sort C2 C3

let K3=C3(25)
let K4=C3(975)
print K3 K4

ENDMACRO

The macro's author think on this way: I have a y vector of sorted means with
1000 registers. So, my CI is between the 25º and 975º element.

So, I ask: Is this a Monte Carlo Simulation, and the nome is correct? The
way to isolate the inferior and superior means is correct?

About the graphics, I know the sample can't be reproduced because it's
random. But, re-running the script some times, I can see, some times,
differences in the right and left tail, under the normal curve. So, I wonder
to know where my code is wrong, but just for know, I agree with histogram
isn't the best way to illustrate confidence interval. I just show him on the
plot to illustrated the sample used. I will re-make the plot, just with the
curve and areas after the confidence interval, without the histogram.


Thanks for the help!






On Wed, Aug 20, 2008 at 12:16 AM, Prof Brian Ripley
[EMAIL PROTECTED]wrote:

 On Tue, 19 Aug 2008, Raphael Saldanha wrote:

  Hi!

 With the following script, I'm trying to make a demonstration of a
 Confidence Interval, but I'm observing some differences on tails.


 You need to tell us what you are trying to do.  You seem to be computing a
 parametric bootstrap interval, but incorrectly (you need to reflect
 percentile intervals).  See Davison  Hinkley (1997) 'The Bootstrap and its
 Application' for more details.

 In any case, your simulation is not repeatable (you set no seed), so we
 don't know what you saw and what 'differences' disturbed you.

 Your calculation of quantiles is not quite correct: use quantile().
 (The indices go from 1 to 1000, not 0 to 1000.)  E.g.

  quantile(1:1000, c(0.025, 0.975))

   2.5%   97.5%
  25.975 975.025

 and not 25, 975.

 As your results show, a histogram is not a good summary of the shape of a
 distribution on 1000 points.  We can do much better with an ecdf or a
 density estimate.


 # Teste de m?dia entre uma amostra e uma popula??o normal
 # Autor: Raphael de Freitas Saldanha
 # Agosto de 2008

 n- 200# Sample size
 xbar - 100# Sample mean
 s- 2  # Sample SD
 nc   - 0.95   # Confidence level (95% - 0.95)
 rep  - 1000   # Loops

 ###

 y - NULL# Vetor com as m?dias da amostra
 for (i in 1:rep){# Loop
 x - rnorm(n,xbar,s) # Gere uma amostra normal n elementos, xbar m?dia
 e
 s desvio-padr?o
 x - mean(x) # Calcule a m?dia (exata) dessa amostra
 y - c(y,x)  # Coloque essa m?dia em um registro em y
 }

 y - sort(y) # Ordene todas as m?dias geradas

 LI - y[((1-nc)/2)*rep] # Limite inferior,
 LS - y[rep-(((1-nc)/2)*rep)]   # Limite superior

 ### PARTE GR?FICA ###

 x - mean(y)

 xvals - seq(-LI, LI, length.out=5000)
 dvals - dnorm(xvals,mean(y), sd(y))[1:5000]

 xbvals - seq(LS, LS*2, length.out=5000)
 dbvals - dnorm(xbvals,mean(y), sd(y))[1:5000]

 ahist - hist(y, freq=FALSE, col=lightblue, main=Intervalo de
 confian?a)

 polygon(c(xvals,LI,LI), c(dvals,dnorm(LI,mean(y), sd(y)),min(dvals)),
 col=orange, lty=0)
 polygon(c(LS,LS,xbvals), c(min(dbvals),dnorm(LS,mean(y), sd(y)),dbvals),
 col=orange, lty=0)
 curve(dnorm(x,mean(y), sd(y)),add=TRUE, lty=1, col=red,lwd=2)

 ### Intervalo de Confian?a ###

 LI # Limite inferior
 LS # Limite superior

 --
 Atenciosamente,

 Raphael Saldanha
 [EMAIL PROTECTED]

[[alternative HTML version deleted]]



 --
 Brian D. Ripley,  [EMAIL PROTECTED]
 Professor of Applied Statistics,  
 http://www.stats.ox.ac.uk/~ripley/http://www.stats.ox.ac.uk/%7Eripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595




-- 
Atenciosamente,

Raphael Saldanha
[EMAIL PROTECTED]
Robert Orben  - To err is human - and to blame it on a computer is even
more so.

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence Interval

[Raphael Saldanha]

Hi!

Here is a better code, using the percentil value instead the position, and
some corrections on the graph code.

The reason in the difference on the tails is elementar, but I don't had see
it before. The right and left limits are calculated by simulation, and his
differences from the mean are not exactly equal, so the areas under the
normal curve! I'm feeling like dummie now...

Thanks for the help!

Raphael Saldanha
BRAZIL

On Wed, Aug 20, 2008 at 12:07 PM, Raphael Saldanha 
[EMAIL PROTECTED] wrote:

 Thanks Professor and Bernardo,

 What I'm trying to do is this: I have a macro for Minitab. His author says
 it's a Monte Carlo simulation to estimate a confidence interval. But and
 don't have Minitab, don't like to work with illegal licenses, and LOVE R.
 So, I re-write the macro in a R script.

 Bellow, is the original macro, for Minitab:

 GMACRO

 TesteMedia

 do K1=1:1000
  random 200 C1;
  normal 100 2.
  mean C1 K2
  let C2(K1)=K2
 enddo

 sort C2 C3

 let K3=C3(25)
 let K4=C3(975)
 print K3 K4

 ENDMACRO

 The macro's author think on this way: I have a y vector of sorted means
 with 1000 registers. So, my CI is between the 25º and 975º element.

 So, I ask: Is this a Monte Carlo Simulation, and the nome is correct? The
 way to isolate the inferior and superior means is correct?

 About the graphics, I know the sample can't be reproduced because it's
 random. But, re-running the script some times, I can see, some times,
 differences in the right and left tail, under the normal curve. So, I wonder
 to know where my code is wrong, but just for know, I agree with histogram
 isn't the best way to illustrate confidence interval. I just show him on the
 plot to illustrated the sample used. I will re-make the plot, just with the
 curve and areas after the confidence interval, without the histogram.


 Thanks for the help!






 On Wed, Aug 20, 2008 at 12:16 AM, Prof Brian Ripley [EMAIL PROTECTED]
  wrote:

 On Tue, 19 Aug 2008, Raphael Saldanha wrote:

  Hi!

 With the following script, I'm trying to make a demonstration of a
 Confidence Interval, but I'm observing some differences on tails.


 You need to tell us what you are trying to do.  You seem to be computing a
 parametric bootstrap interval, but incorrectly (you need to reflect
 percentile intervals).  See Davison  Hinkley (1997) 'The Bootstrap and its
 Application' for more details.

 In any case, your simulation is not repeatable (you set no seed), so we
 don't know what you saw and what 'differences' disturbed you.

 Your calculation of quantiles is not quite correct: use quantile().
 (The indices go from 1 to 1000, not 0 to 1000.)  E.g.

  quantile(1:1000, c(0.025, 0.975))

   2.5%   97.5%
  25.975 975.025

 and not 25, 975.

 As your results show, a histogram is not a good summary of the shape of a
 distribution on 1000 points.  We can do much better with an ecdf or a
 density estimate.


 # Teste de m?dia entre uma amostra e uma popula??o normal
 # Autor: Raphael de Freitas Saldanha
 # Agosto de 2008

 n- 200# Sample size
 xbar - 100# Sample mean
 s- 2  # Sample SD
 nc   - 0.95   # Confidence level (95% - 0.95)
 rep  - 1000   # Loops

 ###

 y - NULL# Vetor com as m?dias da amostra
 for (i in 1:rep){# Loop
 x - rnorm(n,xbar,s) # Gere uma amostra normal n elementos, xbar
 m?dia e
 s desvio-padr?o
 x - mean(x) # Calcule a m?dia (exata) dessa amostra
 y - c(y,x)  # Coloque essa m?dia em um registro em y
 }

 y - sort(y) # Ordene todas as m?dias geradas

 LI - y[((1-nc)/2)*rep] # Limite inferior,
 LS - y[rep-(((1-nc)/2)*rep)]   # Limite superior

 ### PARTE GR?FICA ###

 x - mean(y)

 xvals - seq(-LI, LI, length.out=5000)
 dvals - dnorm(xvals,mean(y), sd(y))[1:5000]

 xbvals - seq(LS, LS*2, length.out=5000)
 dbvals - dnorm(xbvals,mean(y), sd(y))[1:5000]

 ahist - hist(y, freq=FALSE, col=lightblue, main=Intervalo de
 confian?a)

 polygon(c(xvals,LI,LI), c(dvals,dnorm(LI,mean(y), sd(y)),min(dvals)),
 col=orange, lty=0)
 polygon(c(LS,LS,xbvals), c(min(dbvals),dnorm(LS,mean(y), sd(y)),dbvals),
 col=orange, lty=0)
 curve(dnorm(x,mean(y), sd(y)),add=TRUE, lty=1, col=red,lwd=2)

 ### Intervalo de Confian?a ###

 LI # Limite inferior
 LS # Limite superior

 --
 Atenciosamente,

 Raphael Saldanha
 [EMAIL PROTECTED]

[[alternative HTML version deleted]]



 --
 Brian D. Ripley,  [EMAIL PROTECTED]
 Professor of Applied Statistics,  
 http://www.stats.ox.ac.uk/~ripley/http://www.stats.ox.ac.uk/%7Eripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595




 --
 Atenciosamente,

 Raphael Saldanha
 [EMAIL PROTECTED]
 Robert Orben  - To err is human - and to blame it on a computer is even
 more so.




-- 
Atenciosamente,

Raphael Saldanha

[R] Confidence interval for the coefficient of variation

[Katrien Baert]

Dear,

We are trying to determine the (one-sided) CI for the coefficient of
variation in a small sample (say n = 10), with mean 100 and standard
deviation 21.
It appears though that the R-function ci.cv() and our simulation do not
agree.

The R-code:

library(MBESS)

n = 10
ci.cv(mean = 100, sd = 21, n = 10, conf.level = 0.9)
U10.95 - 0.3551754
ci.cv(mean = 100, sd = 21, n = 10, conf.level = 0.6)
U10.80 - 0.2742255

cv.x - c()
m10.95 - c()
m10.80 - c()

for(j in 1:10){
 for(i in 1:1){
  X - rnorm(n, 100, 21)
  cv.x[i] - sd(X)/mean(X)
 }
 m10.95[j] - mean(ifelse(cv.x  U10.95, 1, 0))
 m10.80[j] - mean(ifelse(cv.x  U10.80, 1, 0))
}
m10.95
m10.80
 m10.95
 [1] 0.0034 0.0054 0.0049 0.0045 0.0045 0.0050 0.0039 0.0043 0.0057 0.0042
 m10.80
 [1] 0.0935 0.0966 0.0976 0.0961 0.0943 0.0915 0.0911 0.0938 0.0968 0.0868

These probabilities are much lower than the expected 5% and 20%.
Does anyone know where our mistake is?

Thanks a lot.

[[alternative HTML version deleted]]

__
R-help@r-project.org 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] Confidence interval for the coefficient of variation

[Ben Bolker]

Katrien Baert katrien.baert at gmail.com writes:

 
 Dear,
 
 We are trying to determine the (one-sided) CI for the coefficient of
 variation in a small sample (say n = 10), with mean 100 and standard
 deviation 21.
 It appears though that the R-function ci.cv() and our simulation do not
 agree.
 

  You probably need to contact the package author.

  Ben Bolker

__
R-help@r-project.org 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.