[R] Full or non constrained/reparametrized model.matrix

Gregor Gorjanc gregor.gorjanc at bfro.uni-lj.si
Wed May 10 08:57:21 CEST 2006


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.




More information about the R-help mailing list