Gregor Gorjanc
2006-May-09 21:43 UTC
[R] Full or non constrained/reparametrized model.matrix
Hello!
I have parameter estimates for a generalized linear model and would like
to produce fitted values i.e. fitted(). This can be easily done in R,
but my problem lies in fact that I have a vector of parameters from some
other software, that is is not constrained i.e. I have the following
estimates for model with one factor with 4 levels
beta = c(intercept group1 group2 group3 group4)
where group1:4 are estimated deviations from intercept i.e. sum to zero
contraint, but all parameter estimates are there! How can I create a
model matrix that will not contain any constraints since I would like to
compute
model.matrix("some_formula") %*% beta
I.e. I would like to have model.matrix of the form
1 1 0 0 0
1 0 1 0 0
1 0 0 1 0
1 0 0 0 1
and not of the following form with contr.treatment or any other contraints
1 0 0 0
1 1 0 0
1 0 1 0
1 0 0 1
I could remove group4 from beta and use sum to zero contraint for
contrast in fomula, but I would like to overcome this, as my model can
be richer in number or parameters. The following example, will show what
I would like to do:
## --- Setup ---
groupN <- 4
NPerGroup <- 5
min <- 1
max <- 5
g <- runif(n = groupN, min = min, max = max)
## --- Simulate ---
group <- factor(rep(paste("G", 1:groupN, sep = ""), each
= NPerGroup))
y <- vector(mode = "numeric", length = groupN * NPerGroup)
j <- 1
for (i in 1:groupN) {
y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i])
j <- (i * NPerGroup) + 1
}
## --- GLM ---
contrasts(group) <- contr.sum(groupN)
fit <- glm(y ~ group, family = "poisson")
coef(fit)
## Now this goes nicely
model.matrix(y ~ group) %*% coef(fit)
## But pretend I have the following vector of estimated parameters
beta <- c(coef(fit), 0 - sum(coef(fit)[-1]))
names(beta) <- c(names(beta)[1:4], "group4")
## I can not apply this as model matrix does not conform to beta
model.matrix(y ~ group) %*% beta
## Is there any general way of constructing full design/model matrix
## without any constraints/reparametrizations?
Thanks!
--
Lep pozdrav / With regards,
Gregor Gorjanc
----------------------------------------------------------------------
University of Ljubljana PhD student
Biotechnical Faculty
Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan
Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si
SI-1230 Domzale tel: +386 (0)1 72 17 861
Slovenia, Europe fax: +386 (0)1 72 17 888
----------------------------------------------------------------------
"One must learn by doing the thing; for though you think you know it,
you have no certainty until you try." Sophocles ~ 450 B.C.
Gabor Grothendieck
2006-May-10 00:01 UTC
[R] Full or non constrained/reparametrized model.matrix
On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote:> Hello! > > I have parameter estimates for a generalized linear model and would like > to produce fitted values i.e. fitted(). This can be easily done in R, > but my problem lies in fact that I have a vector of parameters from some > other software, that is is not constrained i.e. I have the following > estimates for model with one factor with 4 levels > > beta = c(intercept group1 group2 group3 group4) > > where group1:4 are estimated deviations from intercept i.e. sum to zero > contraint, but all parameter estimates are there! How can I create a > model matrix that will not contain any constraints since I would like to > compute > > model.matrix("some_formula") %*% beta > > I.e. I would like to have model.matrix of the form > > 1 1 0 0 0 > 1 0 1 0 0 > 1 0 0 1 0 > 1 0 0 0 1 > > and not of the following form with contr.treatment or any other contraints > > 1 0 0 0 > 1 1 0 0 > 1 0 1 0 > 1 0 0 1 > > I could remove group4 from beta and use sum to zero contraint for > contrast in fomula, but I would like to overcome this, as my model can > be richer in number or parameters. The following example, will show what > I would like to do: > > ## --- Setup --- > > groupN <- 4 > NPerGroup <- 5 > min <- 1 > max <- 5 > g <- runif(n = groupN, min = min, max = max) > > ## --- Simulate --- > > group <- factor(rep(paste("G", 1:groupN, sep = ""), each = NPerGroup)) > y <- vector(mode = "numeric", length = groupN * NPerGroup) > j <- 1 > for (i in 1:groupN) { > y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i]) > j <- (i * NPerGroup) + 1 > } > > ## --- GLM --- > > contrasts(group) <- contr.sum(groupN) > fit <- glm(y ~ group, family = "poisson") > coef(fit) > > ## Now this goes nicely > model.matrix(y ~ group) %*% coef(fit) > > ## But pretend I have the following vector of estimated parameters > beta <- c(coef(fit), 0 - sum(coef(fit)[-1])) > names(beta) <- c(names(beta)[1:4], "group4") > > ## I can not apply this as model matrix does not conform to beta > model.matrix(y ~ group) %*% betaTry this: model.matrix(y ~ group-1) %*% beta[-1] + beta[1]> > ## Is there any general way of constructing full design/model matrix > ## without any constraints/reparametrizations? > > Thanks! > > -- > Lep pozdrav / With regards, > Gregor Gorjanc > > ---------------------------------------------------------------------- > University of Ljubljana PhD student > Biotechnical Faculty > Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan > Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si > > SI-1230 Domzale tel: +386 (0)1 72 17 861 > Slovenia, Europe fax: +386 (0)1 72 17 888 > > ---------------------------------------------------------------------- > "One must learn by doing the thing; for though you think you know it, > you have no certainty until you try." Sophocles ~ 450 B.C. > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >
Gregor Gorjanc
2006-May-10 00:14 UTC
[R] Full or non constrained/reparametrized model.matrix
Hello! Thank you very much for the response. Please read within the lines Gabor Grothendieck wrote:> On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: > >> Hello! >> >> I have parameter estimates for a generalized linear model and would like >> to produce fitted values i.e. fitted(). This can be easily done in R, >> but my problem lies in fact that I have a vector of parameters from some >> other software, that is is not constrained i.e. I have the following >> estimates for model with one factor with 4 levels >> >> beta = c(intercept group1 group2 group3 group4) >> >> where group1:4 are estimated deviations from intercept i.e. sum to zero >> contraint, but all parameter estimates are there! How can I create a >> model matrix that will not contain any constraints since I would like to >> compute >> >> model.matrix("some_formula") %*% beta >> >> I.e. I would like to have model.matrix of the form >> >> 1 1 0 0 0 >> 1 0 1 0 0 >> 1 0 0 1 0 >> 1 0 0 0 1 >> >> and not of the following form with contr.treatment or any other >> contraints >> >> 1 0 0 0 >> 1 1 0 0 >> 1 0 1 0 >> 1 0 0 1 >> >> I could remove group4 from beta and use sum to zero contraint for >> contrast in fomula, but I would like to overcome this, as my model can >> be richer in number or parameters. The following example, will show what >> I would like to do: >> >> ## --- Setup --- >> >> groupN <- 4 >> NPerGroup <- 5 >> min <- 1 >> max <- 5 >> g <- runif(n = groupN, min = min, max = max) >> >> ## --- Simulate --- >> >> group <- factor(rep(paste("G", 1:groupN, sep = ""), each = NPerGroup)) >> y <- vector(mode = "numeric", length = groupN * NPerGroup) >> j <- 1 >> for (i in 1:groupN) { >> y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i]) >> j <- (i * NPerGroup) + 1 >> } >> >> ## --- GLM --- >> >> contrasts(group) <- contr.sum(groupN) >> fit <- glm(y ~ group, family = "poisson") >> coef(fit) >> >> ## Now this goes nicely >> model.matrix(y ~ group) %*% coef(fit) >> >> ## But pretend I have the following vector of estimated parameters >> beta <- c(coef(fit), 0 - sum(coef(fit)[-1])) >> names(beta) <- c(names(beta)[1:4], "group4") >> >> ## I can not apply this as model matrix does not conform to beta >> model.matrix(y ~ group) %*% beta > > > Try this: > > model.matrix(y ~ group-1) %*% beta[-1] + beta[1]This is a nice hack here, but does not work in general. Say I have another factor sex <- factor(rep(paste("S", 1:2, sep = ""), times=10)) model.matrix(y ~ group + sex - 1) produces a matrix of 5 columns and not 6 as I want to.>> >> ## Is there any general way of constructing full design/model matrix >> ## without any constraints/reparametrizations? >> >> Thanks! >> >> -- >> Lep pozdrav / With regards, >> Gregor Gorjanc >> >> ---------------------------------------------------------------------- >> University of Ljubljana PhD student >> Biotechnical Faculty >> Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan >> Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si >> >> SI-1230 Domzale tel: +386 (0)1 72 17 861 >> Slovenia, Europe fax: +386 (0)1 72 17 888 >> >> ---------------------------------------------------------------------- >> "One must learn by doing the thing; for though you think you know it, >> you have no certainty until you try." Sophocles ~ 450 B.C. >> >> ______________________________________________ >> R-help at stat.math.ethz.ch mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide! >> http://www.R-project.org/posting-guide.html >>-- Lep pozdrav / With regards, Gregor Gorjanc ---------------------------------------------------------------------- University of Ljubljana PhD student Biotechnical Faculty Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si SI-1230 Domzale tel: +386 (0)1 72 17 861 Slovenia, Europe fax: +386 (0)1 72 17 888 ---------------------------------------------------------------------- "One must learn by doing the thing; for though you think you know it, you have no certainty until you try." Sophocles ~ 450 B.C.
Gabor Grothendieck
2006-May-10 00:31 UTC
[R] Full or non constrained/reparametrized model.matrix
On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote:> Hello! > > Thank you very much for the response. Please read within the lines > > Gabor Grothendieck wrote: > > On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: > > > >> Hello! > >> > >> I have parameter estimates for a generalized linear model and would like > >> to produce fitted values i.e. fitted(). This can be easily done in R, > >> but my problem lies in fact that I have a vector of parameters from some > >> other software, that is is not constrained i.e. I have the following > >> estimates for model with one factor with 4 levels > >> > >> beta = c(intercept group1 group2 group3 group4) > >> > >> where group1:4 are estimated deviations from intercept i.e. sum to zero > >> contraint, but all parameter estimates are there! How can I create a > >> model matrix that will not contain any constraints since I would like to > >> compute > >> > >> model.matrix("some_formula") %*% beta > >> > >> I.e. I would like to have model.matrix of the form > >> > >> 1 1 0 0 0 > >> 1 0 1 0 0 > >> 1 0 0 1 0 > >> 1 0 0 0 1 > >> > >> and not of the following form with contr.treatment or any other > >> contraints > >> > >> 1 0 0 0 > >> 1 1 0 0 > >> 1 0 1 0 > >> 1 0 0 1 > >> > >> I could remove group4 from beta and use sum to zero contraint for > >> contrast in fomula, but I would like to overcome this, as my model can > >> be richer in number or parameters. The following example, will show what > >> I would like to do: > >> > >> ## --- Setup --- > >> > >> groupN <- 4 > >> NPerGroup <- 5 > >> min <- 1 > >> max <- 5 > >> g <- runif(n = groupN, min = min, max = max) > >> > >> ## --- Simulate --- > >> > >> group <- factor(rep(paste("G", 1:groupN, sep = ""), each = NPerGroup)) > >> y <- vector(mode = "numeric", length = groupN * NPerGroup) > >> j <- 1 > >> for (i in 1:groupN) { > >> y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i]) > >> j <- (i * NPerGroup) + 1 > >> } > >> > >> ## --- GLM --- > >> > >> contrasts(group) <- contr.sum(groupN) > >> fit <- glm(y ~ group, family = "poisson") > >> coef(fit) > >> > >> ## Now this goes nicely > >> model.matrix(y ~ group) %*% coef(fit) > >> > >> ## But pretend I have the following vector of estimated parameters > >> beta <- c(coef(fit), 0 - sum(coef(fit)[-1])) > >> names(beta) <- c(names(beta)[1:4], "group4") > >> > >> ## I can not apply this as model matrix does not conform to beta > >> model.matrix(y ~ group) %*% beta > > > > > > Try this: > > > > model.matrix(y ~ group-1) %*% beta[-1] + beta[1] > > This is a nice hack here, but does not work in general. Say I have > another factor > > sex <- factor(rep(paste("S", 1:2, sep = ""), times=10)) > > model.matrix(y ~ group + sex - 1) > > produces a matrix of 5 columns and not 6 as I want to.Try this: cbind(1, model.matrix(~group-1), model.matrix(~sex-1))> > >> > >> ## Is there any general way of constructing full design/model matrix > >> ## without any constraints/reparametrizations? > >> > >> Thanks! > >> > >> -- > >> Lep pozdrav / With regards, > >> Gregor Gorjanc > >> > >> ---------------------------------------------------------------------- > >> University of Ljubljana PhD student > >> Biotechnical Faculty > >> Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan > >> Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si > >> > >> SI-1230 Domzale tel: +386 (0)1 72 17 861 > >> Slovenia, Europe fax: +386 (0)1 72 17 888 > >> > >> ---------------------------------------------------------------------- > >> "One must learn by doing the thing; for though you think you know it, > >> you have no certainty until you try." Sophocles ~ 450 B.C. > >> > >> ______________________________________________ > >> R-help at stat.math.ethz.ch mailing list > >> https://stat.ethz.ch/mailman/listinfo/r-help > >> PLEASE do read the posting guide! > >> http://www.R-project.org/posting-guide.html > >> > > > -- > Lep pozdrav / With regards, > Gregor Gorjanc > > ---------------------------------------------------------------------- > University of Ljubljana PhD student > Biotechnical Faculty > Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan > Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si > > SI-1230 Domzale tel: +386 (0)1 72 17 861 > Slovenia, Europe fax: +386 (0)1 72 17 888 > > ---------------------------------------------------------------------- > "One must learn by doing the thing; for though you think you know it, > you have no certainty until you try." Sophocles ~ 450 B.C. > ---------------------------------------------------------------------- >
Gabor Grothendieck
2006-May-10 05:31 UTC
[R] Full or non constrained/reparametrized model.matrix
On 5/9/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: > > Hello! > > > > Thank you very much for the response. Please read within the lines > > > > Gabor Grothendieck wrote: > > > On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: > > > > > >> Hello! > > >> > > >> I have parameter estimates for a generalized linear model and would like > > >> to produce fitted values i.e. fitted(). This can be easily done in R, > > >> but my problem lies in fact that I have a vector of parameters from some > > >> other software, that is is not constrained i.e. I have the following > > >> estimates for model with one factor with 4 levels > > >> > > >> beta = c(intercept group1 group2 group3 group4) > > >> > > >> where group1:4 are estimated deviations from intercept i.e. sum to zero > > >> contraint, but all parameter estimates are there! How can I create a > > >> model matrix that will not contain any constraints since I would like to > > >> compute > > >> > > >> model.matrix("some_formula") %*% beta > > >> > > >> I.e. I would like to have model.matrix of the form > > >> > > >> 1 1 0 0 0 > > >> 1 0 1 0 0 > > >> 1 0 0 1 0 > > >> 1 0 0 0 1 > > >> > > >> and not of the following form with contr.treatment or any other > > >> contraints > > >> > > >> 1 0 0 0 > > >> 1 1 0 0 > > >> 1 0 1 0 > > >> 1 0 0 1 > > >> > > >> I could remove group4 from beta and use sum to zero contraint for > > >> contrast in fomula, but I would like to overcome this, as my model can > > >> be richer in number or parameters. The following example, will show what > > >> I would like to do: > > >> > > >> ## --- Setup --- > > >> > > >> groupN <- 4 > > >> NPerGroup <- 5 > > >> min <- 1 > > >> max <- 5 > > >> g <- runif(n = groupN, min = min, max = max) > > >> > > >> ## --- Simulate --- > > >> > > >> group <- factor(rep(paste("G", 1:groupN, sep = ""), each = NPerGroup)) > > >> y <- vector(mode = "numeric", length = groupN * NPerGroup) > > >> j <- 1 > > >> for (i in 1:groupN) { > > >> y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i]) > > >> j <- (i * NPerGroup) + 1 > > >> } > > >> > > >> ## --- GLM --- > > >> > > >> contrasts(group) <- contr.sum(groupN) > > >> fit <- glm(y ~ group, family = "poisson") > > >> coef(fit) > > >> > > >> ## Now this goes nicely > > >> model.matrix(y ~ group) %*% coef(fit) > > >> > > >> ## But pretend I have the following vector of estimated parameters > > >> beta <- c(coef(fit), 0 - sum(coef(fit)[-1])) > > >> names(beta) <- c(names(beta)[1:4], "group4") > > >> > > >> ## I can not apply this as model matrix does not conform to beta > > >> model.matrix(y ~ group) %*% beta > > > > > > > > > Try this: > > > > > > model.matrix(y ~ group-1) %*% beta[-1] + beta[1] > > > > This is a nice hack here, but does not work in general. Say I have > > another factor > > > > sex <- factor(rep(paste("S", 1:2, sep = ""), times=10)) > > > > model.matrix(y ~ group + sex - 1) > > > > produces a matrix of 5 columns and not 6 as I want to. > > Try this: > > cbind(1, model.matrix(~group-1), model.matrix(~sex-1)) >In thinking about this a bit more, try this: attr(group, "contrasts") <- diag(nlevels(group)) attr(sex, "contrasts") <- diag(nlevels(sex)) model.matrix(~ group + sex)
Gregor Gorjanc
2006-May-10 06:57 UTC
[R] Full or non constrained/reparametrized model.matrix
Gabor thanks! Last suggestion works like charm! Gabor Grothendieck wrote:> On 5/9/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote: > >> On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: >> > Hello! >> > >> > Thank you very much for the response. Please read within the lines >> > >> > Gabor Grothendieck wrote: >> > > On 5/9/06, Gregor Gorjanc <gregor.gorjanc at gmail.com> wrote: >> > > >> > >> Hello! >> > >> >> > >> I have parameter estimates for a generalized linear model and >> would like >> > >> to produce fitted values i.e. fitted(). This can be easily done >> in R, >> > >> but my problem lies in fact that I have a vector of parameters >> from some >> > >> other software, that is is not constrained i.e. I have the following >> > >> estimates for model with one factor with 4 levels >> > >> >> > >> beta = c(intercept group1 group2 group3 group4) >> > >> >> > >> where group1:4 are estimated deviations from intercept i.e. sum >> to zero >> > >> contraint, but all parameter estimates are there! How can I create a >> > >> model matrix that will not contain any constraints since I would >> like to >> > >> compute >> > >> >> > >> model.matrix("some_formula") %*% beta >> > >> >> > >> I.e. I would like to have model.matrix of the form >> > >> >> > >> 1 1 0 0 0 >> > >> 1 0 1 0 0 >> > >> 1 0 0 1 0 >> > >> 1 0 0 0 1 >> > >> >> > >> and not of the following form with contr.treatment or any other >> > >> contraints >> > >> >> > >> 1 0 0 0 >> > >> 1 1 0 0 >> > >> 1 0 1 0 >> > >> 1 0 0 1 >> > >> >> > >> I could remove group4 from beta and use sum to zero contraint for >> > >> contrast in fomula, but I would like to overcome this, as my >> model can >> > >> be richer in number or parameters. The following example, will >> show what >> > >> I would like to do: >> > >> >> > >> ## --- Setup --- >> > >> >> > >> groupN <- 4 >> > >> NPerGroup <- 5 >> > >> min <- 1 >> > >> max <- 5 >> > >> g <- runif(n = groupN, min = min, max = max) >> > >> >> > >> ## --- Simulate --- >> > >> >> > >> group <- factor(rep(paste("G", 1:groupN, sep = ""), each >> NPerGroup)) >> > >> y <- vector(mode = "numeric", length = groupN * NPerGroup) >> > >> j <- 1 >> > >> for (i in 1:groupN) { >> > >> y[j:(i * NPerGroup)] <- rpois(n = NPerGroup, lambda = g[i]) >> > >> j <- (i * NPerGroup) + 1 >> > >> } >> > >> >> > >> ## --- GLM --- >> > >> >> > >> contrasts(group) <- contr.sum(groupN) >> > >> fit <- glm(y ~ group, family = "poisson") >> > >> coef(fit) >> > >> >> > >> ## Now this goes nicely >> > >> model.matrix(y ~ group) %*% coef(fit) >> > >> >> > >> ## But pretend I have the following vector of estimated parameters >> > >> beta <- c(coef(fit), 0 - sum(coef(fit)[-1])) >> > >> names(beta) <- c(names(beta)[1:4], "group4") >> > >> >> > >> ## I can not apply this as model matrix does not conform to beta >> > >> model.matrix(y ~ group) %*% beta >> > > >> > > >> > > Try this: >> > > >> > > model.matrix(y ~ group-1) %*% beta[-1] + beta[1] >> > >> > This is a nice hack here, but does not work in general. Say I have >> > another factor >> > >> > sex <- factor(rep(paste("S", 1:2, sep = ""), times=10)) >> > >> > model.matrix(y ~ group + sex - 1) >> > >> > produces a matrix of 5 columns and not 6 as I want to. >> >> Try this: >> >> cbind(1, model.matrix(~group-1), model.matrix(~sex-1)) >> > > In thinking about this a bit more, try this: > > attr(group, "contrasts") <- diag(nlevels(group)) > attr(sex, "contrasts") <- diag(nlevels(sex)) > model.matrix(~ group + sex)-- Lep pozdrav / With regards, Gregor Gorjanc ---------------------------------------------------------------------- University of Ljubljana PhD student Biotechnical Faculty Zootechnical Department URI: http://www.bfro.uni-lj.si/MR/ggorjan Groblje 3 mail: gregor.gorjanc <at> bfro.uni-lj.si SI-1230 Domzale tel: +386 (0)1 72 17 861 Slovenia, Europe fax: +386 (0)1 72 17 888 ---------------------------------------------------------------------- "One must learn by doing the thing; for though you think you know it, you have no certainty until you try." Sophocles ~ 450 B.C.