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

Gregor Gorjanc gregor.gorjanc at gmail.com
Wed May 10 02:14:33 CEST 2006


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.




More information about the R-help mailing list