From: Gabor Grothendieck <ggrothendieck_at_gmail.com>

Date: Wed 10 May 2006 - 10:01:00 EST

R-help@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 Received on Wed May 10 10:06:41 2006

Date: Wed 10 May 2006 - 10:01:00 EST

On 5/9/06, Gregor Gorjanc <gregor.gorjanc@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]

*>
*

> ## 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@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
**>
*

R-help@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 Received on Wed May 10 10:06:41 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Wed 10 May 2006 - 12:10:13 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*