

Hi,
I am trying to do BoxCox transformation, but I am not sure how to do it correctly. Here is an example showing what I am trying:
# example from MASS
require(MASS)
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
lambda = seq(0.05, 0.45, len = 20))
# Here is My attempt at getting the profile likelihood for the BoxCox parameter
lam < seq(0.05, 0.45, len = 20)
dev < rep(NA, length=20)
for (i in 1:20) {
a < lam[i]
ans < glm(((Days+1)^a1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data = quine)
dev[i] < ans$deviance
}
plot(lam, dev, type="b", xlab="lambda", ylab="deviance")
I am trying to create the profile likelihood for the BoxCox parameter, but obviously I am not getting it right. I am not sure that ans$deviance is the right thing to do.
I would appreciate any guidance.
Thanks & Best,
Ravi
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Hi Ravi,
Deviance is the SS in this case, but you need a normalizing constant
adjusted by the lambda to put them on the same scale. I modified your
example below to simplify slightly and use the normalization (see the
LL line).
Cheers,
Josh
######################################
require(MASS)
myp < function(y, lambda) (y^lambda1)/lambda
lambda < seq(0.05, 0.45, len = 20)
N < nrow(quine)
res < matrix(numeric(0), nrow = length(lambda), 2, dimnames =
list(NULL, c("Lambda", "LL")))
# scaling contant
C < exp(mean(log(quine$Days+1)))
for(i in seq_along(lambda)) {
r < resid(lm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
LL < ( (N/2) * log(sum((r/(C^lambda[i]))^2)))
res[i, ] < c(lambda[i], LL)
}
# box cox
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda)
# add our points on top to verify match
points(res[, 1], res[,2], pch = 16)
######################################
On Mon, Jul 7, 2014 at 11:33 AM, Ravi Varadhan < [hidden email]> wrote:
> Hi,
>
> I am trying to do BoxCox transformation, but I am not sure how to do it correctly. Here is an example showing what I am trying:
>
>
>
> # example from MASS
>
> require(MASS)
> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
> lambda = seq(0.05, 0.45, len = 20))
>
> # Here is My attempt at getting the profile likelihood for the BoxCox parameter
> lam < seq(0.05, 0.45, len = 20)
> dev < rep(NA, length=20)
>
> for (i in 1:20) {
> a < lam[i]
> ans < glm(((Days+1)^a1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data = quine)
> dev[i] < ans$deviance
> }
>
> plot(lam, dev, type="b", xlab="lambda", ylab="deviance")
>
> I am trying to create the profile likelihood for the BoxCox parameter, but obviously I am not getting it right. I am not sure that ans$deviance is the right thing to do.
>
> I would appreciate any guidance.
>
> Thanks & Best,
> Ravi
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.

Joshua F. Wiley
Ph.D. Student, UCLA Department of Psychology
http://joshuawiley.com/Senior Analyst, Elkhart Group Ltd.
http://elkhartgroup.comOffice: 260.673.5518
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Dear Ravi,
In my previous example, I used the residuals, so:
sum [ (r_i / scaling)^2 ]
If you want to use the deviance from glm, that gives you:
sum [ r_i^2 ]
and since the scaling factor is just a constant for any given lambda,
then the modification would be:
sum [ r_i^2 ] / ( scaling^2 )
and is given in the modified code below (posted back to Rhelp in case
any else has this question).
Hope this helps,
Josh
##########################################
require(MASS)
myp < function(y, lambda) (y^lambda1)/lambda
lambda < seq(0.05, 0.45, len = 20)
N < nrow(quine)
res < matrix(numeric(0), nrow = length(lambda), 2, dimnames =
list(NULL, c("Lambda", "LL")))
# scaling contant
C < exp(mean(log(quine$Days+1)))
for(i in seq_along(lambda)) {
SS < deviance(glm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
LL < ( (N/2) * log(SS/((C^lambda[i])^2)))
res[i, ] < c(lambda[i], LL)
}
# box cox
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda)
# add our points on top to verify match
points(res[, 1], res[,2], pch = 16)
##########################################
On Mon, Jul 7, 2014 at 11:57 PM, Ravi Varadhan < [hidden email]> wrote:
> Dear Josh,
> Thank you very much. I knew that the scaling had to be adjusted, but was not sure on how to do this.
>
> Can you please show me how to do this scaling with `glm'? In other words, how would I scale the deviance from glm?
>
> Thanks,
> Ravi
>
> Original Message
> From: Joshua Wiley [mailto: [hidden email]]
> Sent: Sunday, July 06, 2014 11:34 PM
> To: Ravi Varadhan
> Cc: [hidden email]
> Subject: Re: [R] Boxcox transformation
>
> Hi Ravi,
>
> Deviance is the SS in this case, but you need a normalizing constant adjusted by the lambda to put them on the same scale. I modified your example below to simplify slightly and use the normalization (see the LL line).
>
> Cheers,
>
> Josh
>
> ######################################
>
> require(MASS)
>
> myp < function(y, lambda) (y^lambda1)/lambda
>
>
> lambda < seq(0.05, 0.45, len = 20)
> N < nrow(quine)
> res < matrix(numeric(0), nrow = length(lambda), 2, dimnames = list(NULL, c("Lambda", "LL")))
>
> # scaling contant
> C < exp(mean(log(quine$Days+1)))
>
> for(i in seq_along(lambda)) {
> r < resid(lm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
> LL < ( (N/2) * log(sum((r/(C^lambda[i]))^2)))
> res[i, ] < c(lambda[i], LL)
> }
>
> # box cox
> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda) # add our points on top to verify match points(res[, 1], res[,2], pch = 16)
>
> ######################################
>
>
>
> On Mon, Jul 7, 2014 at 11:33 AM, Ravi Varadhan < [hidden email]> wrote:
>> Hi,
>>
>> I am trying to do BoxCox transformation, but I am not sure how to do it correctly. Here is an example showing what I am trying:
>>
>>
>>
>> # example from MASS
>>
>> require(MASS)
>> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
>> lambda = seq(0.05, 0.45, len = 20))
>>
>> # Here is My attempt at getting the profile likelihood for the BoxCox
>> parameter lam < seq(0.05, 0.45, len = 20) dev < rep(NA, length=20)
>>
>> for (i in 1:20) {
>> a < lam[i]
>> ans < glm(((Days+1)^a1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data =
>> quine) dev[i] < ans$deviance }
>>
>> plot(lam, dev, type="b", xlab="lambda", ylab="deviance")
>>
>> I am trying to create the profile likelihood for the BoxCox parameter, but obviously I am not getting it right. I am not sure that ans$deviance is the right thing to do.
>>
>> I would appreciate any guidance.
>>
>> Thanks & Best,
>> Ravi
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide
>> http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>
>
>
> 
> Joshua F. Wiley
> Ph.D. Student, UCLA Department of Psychology http://joshuawiley.com/ Senior Analyst, Elkhart Group Ltd.
> http://elkhartgroup.com> Office: 260.673.5518

Joshua F. Wiley
Ph.D. Student, UCLA Department of Psychology
http://joshuawiley.com/Senior Analyst, Elkhart Group Ltd.
http://elkhartgroup.comOffice: 260.673.5518
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Thank you. It is very helpful.
Ravi
Original Message
From: Joshua Wiley [mailto: [hidden email]]
Sent: Monday, July 07, 2014 4:15 PM
To: Ravi Varadhan
Cc: [hidden email]
Subject: Re: [R] Boxcox transformation
Dear Ravi,
In my previous example, I used the residuals, so:
sum [ (r_i / scaling)^2 ]
If you want to use the deviance from glm, that gives you:
sum [ r_i^2 ]
and since the scaling factor is just a constant for any given lambda, then the modification would be:
sum [ r_i^2 ] / ( scaling^2 )
and is given in the modified code below (posted back to Rhelp in case any else has this question).
Hope this helps,
Josh
##########################################
require(MASS)
myp < function(y, lambda) (y^lambda1)/lambda
lambda < seq(0.05, 0.45, len = 20)
N < nrow(quine)
res < matrix(numeric(0), nrow = length(lambda), 2, dimnames = list(NULL, c("Lambda", "LL")))
# scaling contant
C < exp(mean(log(quine$Days+1)))
for(i in seq_along(lambda)) {
SS < deviance(glm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
LL < ( (N/2) * log(SS/((C^lambda[i])^2)))
res[i, ] < c(lambda[i], LL)
}
# box cox
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda) # add our points on top to verify match points(res[, 1], res[,2], pch = 16)
##########################################
On Mon, Jul 7, 2014 at 11:57 PM, Ravi Varadhan < [hidden email]> wrote:
> Dear Josh,
> Thank you very much. I knew that the scaling had to be adjusted, but was not sure on how to do this.
>
> Can you please show me how to do this scaling with `glm'? In other words, how would I scale the deviance from glm?
>
> Thanks,
> Ravi
>
> Original Message
> From: Joshua Wiley [mailto: [hidden email]]
> Sent: Sunday, July 06, 2014 11:34 PM
> To: Ravi Varadhan
> Cc: [hidden email]
> Subject: Re: [R] Boxcox transformation
>
> Hi Ravi,
>
> Deviance is the SS in this case, but you need a normalizing constant adjusted by the lambda to put them on the same scale. I modified your example below to simplify slightly and use the normalization (see the LL line).
>
> Cheers,
>
> Josh
>
> ######################################
>
> require(MASS)
>
> myp < function(y, lambda) (y^lambda1)/lambda
>
>
> lambda < seq(0.05, 0.45, len = 20)
> N < nrow(quine)
> res < matrix(numeric(0), nrow = length(lambda), 2, dimnames =
> list(NULL, c("Lambda", "LL")))
>
> # scaling contant
> C < exp(mean(log(quine$Days+1)))
>
> for(i in seq_along(lambda)) {
> r < resid(lm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
> LL < ( (N/2) * log(sum((r/(C^lambda[i]))^2)))
> res[i, ] < c(lambda[i], LL)
> }
>
> # box cox
> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda) # add
> our points on top to verify match points(res[, 1], res[,2], pch = 16)
>
> ######################################
>
>
>
> On Mon, Jul 7, 2014 at 11:33 AM, Ravi Varadhan < [hidden email]> wrote:
>> Hi,
>>
>> I am trying to do BoxCox transformation, but I am not sure how to do it correctly. Here is an example showing what I am trying:
>>
>>
>>
>> # example from MASS
>>
>> require(MASS)
>> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
>> lambda = seq(0.05, 0.45, len = 20))
>>
>> # Here is My attempt at getting the profile likelihood for the
>> BoxCox parameter lam < seq(0.05, 0.45, len = 20) dev < rep(NA,
>> length=20)
>>
>> for (i in 1:20) {
>> a < lam[i]
>> ans < glm(((Days+1)^a1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data
>> =
>> quine) dev[i] < ans$deviance }
>>
>> plot(lam, dev, type="b", xlab="lambda", ylab="deviance")
>>
>> I am trying to create the profile likelihood for the BoxCox parameter, but obviously I am not getting it right. I am not sure that ans$deviance is the right thing to do.
>>
>> I would appreciate any guidance.
>>
>> Thanks & Best,
>> Ravi
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide
>> http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>
>
>
> 
> Joshua F. Wiley
> Ph.D. Student, UCLA Department of Psychology http://joshuawiley.com/ Senior Analyst, Elkhart Group Ltd.
> http://elkhartgroup.com> Office: 260.673.5518

Joshua F. Wiley
Ph.D. Student, UCLA Department of Psychology http://joshuawiley.com/ Senior Analyst, Elkhart Group Ltd.
http://elkhartgroup.comOffice: 260.673.5518
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

