Re: [R] gam confidence interval (package mgcv)

2011-06-28 Thread Simon Wood
not sure if I'm missing something here, but since you are using a log 
link, isn't the ratio you are looking for given by the `treatmentB' 
parameter in the summary (independent of X)


> summary(gfit)
[snip]

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.831830.03885  47.152 < 2e-16 ***
treatmentB   0.445130.05567   7.996  1.4e-09 ***
---
[snip]

Let mu = E(y), and b be a binary indicator for treatment B. You want
mu|b=1/mu|b=0

log(mu|b=1) = intercept + treatmentB + s(X)
log(mu|b=0) = intercept + s(X)

=> log(mu|b=1) - log(mu|b=0) = treatmentB

so mu|b=1/mu|b = exp(treatmentB)

So you can get the required interval by finding and interval for 
treatment B and exponentiating...


tB <- coef(gfit)[2]
se.tB <- sqrt(vcov(gfit)[2,2])
exp(c(tB - 2*se.tB,tB+2*se.tB))

On 06/28/2011 03:45 AM, Remko Duursma wrote:

Dear R-helpers,

I am trying to construct a confidence interval on a prediction of a
gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant
but I am not able to apply that to this, different, problem.

Any help is appreciated!

Basically I have a function Y = f(X) for two different treatments A
and B.  I am interested in the treatment ratios : Y(treatment = B) /
Y(treatment = A) as a function of X, including a confidence interval
for this treatment ratio (because we are testing this ratio against
some value, across the range of X).

The X values that Y is measured at differs between the treatments, but
the ranges are similar.


# Reproducible problem:
X1<- runif(20, 0.5, 4)
X2<- runif(20, 0.5, 4)

Y1<- 20*exp(-0.5*X1) + rnorm(20)
Y2<- 30*exp(-0.5*X2) + rnorm(20)

# Look at data:
plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5))
points(X2, Y2, pch=19, col="red")

# Full dataset
dfr<- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A",20),rep("B",20)))

# Fit gam
# I use a gamma family here although it is not necessary: in the real
problem it is, though.
gfit<- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log))

# I am interested in the relationship:
# Y(treatment =="B") / Y(treatment=="A") as a function of X, with a
confidence interval!

Do I just do a bootstrap here? Or is there a more appropriate method?

Thanks a lot for any help.

greetings,
Remko





-
Remko Duursma
Research Lecturer

Hawkesbury Institute for the Environment
University of Western Sydney
Hawkesbury Campus, Richmond

Mobile: +61 (0)422 096908
www.remkoduursma.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] gam confidence interval (package mgcv)

2011-06-27 Thread Remko Duursma
But that just gives me the prediction of Y for treatment A or B, not the ratio.

As I stated:

# I am interested in the relationship:
# Y(treatment =="B") / Y(treatment=="A") as a function of X, with a
confidence interval!


I can get the SE for either of them using predict.gam without a
problem, but I don't know how to get the CI for the ratio!

thanks,
remko

-
Remko Duursma
Research Lecturer

Hawkesbury Institute for the Environment
University of Western Sydney
Hawkesbury Campus, Richmond

Mobile: +61 (0)422 096908
www.remkoduursma.com



On Tue, Jun 28, 2011 at 2:09 PM, David Winsemius  wrote:
>
> On Jun 27, 2011, at 10:45 PM, Remko Duursma wrote:
>
>> Dear R-helpers,
>>
>> I am trying to construct a confidence interval on a prediction of a
>> gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant
>> but I am not able to apply that to this, different, problem.
>>
>> Any help is appreciated!
>>
>> Basically I have a function Y = f(X) for two different treatments A
>> and B.  I am interested in the treatment ratios : Y(treatment = B) /
>> Y(treatment = A) as a function of X, including a confidence interval
>> for this treatment ratio (because we are testing this ratio against
>> some value, across the range of X).
>>
>> The X values that Y is measured at differs between the treatments, but
>> the ranges are similar.
>>
>>
>> # Reproducible problem:
>> X1 <- runif(20, 0.5, 4)
>> X2 <- runif(20, 0.5, 4)
>>
>> Y1 <- 20*exp(-0.5*X1) + rnorm(20)
>> Y2 <- 30*exp(-0.5*X2) + rnorm(20)
>>
>> # Look at data:
>> plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5))
>> points(X2, Y2, pch=19, col="red")
>>
>> # Full dataset
>> dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2),
>> treatment=c(rep("A",20),rep("B",20)))
>>
>> # Fit gam
>> # I use a gamma family here although it is not necessary: in the real
>> problem it is, though.
>> gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log))
>>
>> # I am interested in the relationship:
>> # Y(treatment =="B") / Y(treatment=="A") as a function of X,
>
> Can't you use predict.gam?
>
>  plot(predict(gfit, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2),
>      treatment=c(rep("A",37),rep("B",37) ) ) )[1:37] )
> lines(predict(gfit3, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2),
>      treatment=c(rep("A",37),rep("B",37) ) ) )[-(1:37)])
>
>> with a confidence interval!
>
> There is an se.fit argument to predict.gam().
>
>
> --
>
> 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] gam confidence interval (package mgcv)

2011-06-27 Thread David Winsemius


On Jun 27, 2011, at 10:45 PM, Remko Duursma wrote:


Dear R-helpers,

I am trying to construct a confidence interval on a prediction of a
gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant
but I am not able to apply that to this, different, problem.

Any help is appreciated!

Basically I have a function Y = f(X) for two different treatments A
and B.  I am interested in the treatment ratios : Y(treatment = B) /
Y(treatment = A) as a function of X, including a confidence interval
for this treatment ratio (because we are testing this ratio against
some value, across the range of X).

The X values that Y is measured at differs between the treatments, but
the ranges are similar.


# Reproducible problem:
X1 <- runif(20, 0.5, 4)
X2 <- runif(20, 0.5, 4)

Y1 <- 20*exp(-0.5*X1) + rnorm(20)
Y2 <- 30*exp(-0.5*X2) + rnorm(20)

# Look at data:
plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5))
points(X2, Y2, pch=19, col="red")

# Full dataset
dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A", 
20),rep("B",20)))


# Fit gam
# I use a gamma family here although it is not necessary: in the real
problem it is, though.
gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log))

# I am interested in the relationship:
# Y(treatment =="B") / Y(treatment=="A") as a function of X,


Can't you use predict.gam?

 plot(predict(gfit, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2),
  treatment=c(rep("A",37),rep("B",37) ) ) )[1:37] )
lines(predict(gfit3, newdata=data.frame(X=rep(seq(0.4, 4, by=0.1), 2),
  treatment=c(rep("A",37),rep("B",37) ) ) )[-(1:37)])


with a confidence interval!


There is an se.fit argument to predict.gam().


--

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] gam confidence interval (package mgcv)

2011-06-27 Thread Remko Duursma
Dear R-helpers,

I am trying to construct a confidence interval on a prediction of a
gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant
but I am not able to apply that to this, different, problem.

Any help is appreciated!

Basically I have a function Y = f(X) for two different treatments A
and B.  I am interested in the treatment ratios : Y(treatment = B) /
Y(treatment = A) as a function of X, including a confidence interval
for this treatment ratio (because we are testing this ratio against
some value, across the range of X).

The X values that Y is measured at differs between the treatments, but
the ranges are similar.


# Reproducible problem:
X1 <- runif(20, 0.5, 4)
X2 <- runif(20, 0.5, 4)

Y1 <- 20*exp(-0.5*X1) + rnorm(20)
Y2 <- 30*exp(-0.5*X2) + rnorm(20)

# Look at data:
plot(X1, Y1, pch=19, col="blue", ylim=c(0,max(Y1,Y2)), xlim=c(0,5))
points(X2, Y2, pch=19, col="red")

# Full dataset
dfr <- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep("A",20),rep("B",20)))

# Fit gam
# I use a gamma family here although it is not necessary: in the real
problem it is, though.
gfit <- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log))

# I am interested in the relationship:
# Y(treatment =="B") / Y(treatment=="A") as a function of X, with a
confidence interval!

Do I just do a bootstrap here? Or is there a more appropriate method?

Thanks a lot for any help.

greetings,
Remko





-
Remko Duursma
Research Lecturer

Hawkesbury Institute for the Environment
University of Western Sydney
Hawkesbury Campus, Richmond

Mobile: +61 (0)422 096908
www.remkoduursma.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.