Hi,
I am trying to do Box-Cox 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 Box-Cox 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)^a-1)/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 Box-Cox 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/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, 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^lambda-1)/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 Box-Cox 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 Box-Cox 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)^a-1)/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 Box-Cox 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/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 F. Wiley Ph.D. Student, UCLA Department of Psychology http://joshuawiley.com/ Senior Analyst, Elkhart Group Ltd. http://elkhartgroup.com Office: 260.673.5518 ______________________________________________ [hidden email] 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. |
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 R-help in case any else has this question). Hope this helps, Josh ########################################## require(MASS) myp <- function(y, lambda) (y^lambda-1)/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] Box-cox 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^lambda-1)/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 Box-Cox 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 Box-Cox >> 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)^a-1)/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 Box-Cox 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/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 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.com Office: 260.673.5518 ______________________________________________ [hidden email] 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. |
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] Box-cox 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 R-help in case any else has this question). Hope this helps, Josh ########################################## require(MASS) myp <- function(y, lambda) (y^lambda-1)/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] Box-cox 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^lambda-1)/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 Box-Cox 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 >> Box-Cox 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)^a-1)/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 Box-Cox 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/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 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.com Office: 260.673.5518 ______________________________________________ [hidden email] 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. |
Free forum by Nabble | Edit this page |