[R] how to create model matrix of second order terms

[Ben Bolker]

   I'm still not entirely clear what you want, in particular what you 
want to with factor variables ('squaring' them doesn't really seem to 
make sense).

Some combination of sprintf/paste/reformulate might do what you want:

## get all second-order terms (doesn't square numeric values)
quad_terms <- sprintf("(%s)^2", paste(b1, collapse = "+"))

## get squared terms
square_terms <- paste(sprintf("I(%s^2)", b1))

reformulate(c(quad_terms, square_terms))

The other alternative would be to construct a second-order polynomial 
across all of your terms (and then insert this into reformulate())

poly_term <- sprintf("poly(%s, degree = 2, raw = TRUE)",
              paste(b1, collapse = ", "))



On 2025-03-24 11:28 a.m., Stephen Bond via R-help wrote:
> Folks,
> 
> I appreciate your effort, but maybe I was not explicit enough, so let
> me try again.
> 
> The current setup for formulas does not allow for I(x^2) terms as
> explained in the MASS book at the end of Section 6.2 the x:x
> interaction is treated as x.
> 
> So I need to write my own code, which is clumsy unless you can refer me
> to a package that already exists or give me an idea how to improve my
> code. Also, writing out terms is not feasible when there are 1000
> variables, so the code needs to deal with taking a wide data frame or
> matrix with column names for convenience.
> 
> Let me know your ideas :-)
> 
> On Mon, 2025-03-24 at 02:43 -0700, Bert Gunter wrote:
>> Full disclosure: I did not attempt to decipher your code.
>>
>> But ~(A+B +C)^2 - (A + B + C)
>> gives all 2nd order interactions whether the terms are factors or
>> numeric.
>>
>> ~I(A^2) + I(B^2) gives quadratics in A and B, which must be numeric,
>> not factors, of course
>>
>> You can combine these as necessary to get a formula expression for
>> just 2nd order terms. Wrapping this in model.matrix() should then
>> give you the model matrix using "treatment" contrasts for the
>> contrasts involving factors (you can change the contrast types using
>> the 'contrasts.arg' argument of model.matrix())
>>
>> 1. Does this help?
>> 2. Do check this to make sure I'm correct
>>
>> Cheers,
>> Bert
>>
>> "An educated person is one who can entertain new ideas, entertain
>> others, and entertain herself."
>>
>>
>>
>> On Mon, Mar 24, 2025 at 12:22 AM Stephen Bond via R-help
>> <r-help using r-project.org> wrote:
>>> I am sending to this forum as stackoverflow has devolved into sth
>>> pretty bad.
>>> Below code shows how to get what I want in a clumsy way.
>>>
>>> cols <- letters[1:4]
>>> a1 <- outer(cols,cols,paste0)
>>> b1 <- a1[!lower.tri(a1)]
>>>
>>> X <- matrix(rnorm(80),ncol=4)
>>> colnames(X) <- cols
>>> X <- as.data.frame(X)
>>> XX <- matrix(0,nrow=nrow(X),ncol=length(b1))
>>> colnames(XX) <- b1
>>>
>>> for (k in 1:length(b1)){
>>>      XX[,k] <- X[,substr(b1[k],1,1)]*X[,substr(b1[k],2,2)]
>>> }
>>>
>>>
>>>
>>> Is there a way to get that using a formula or some neat trick? The
>>> above will not work for factors, so I will need to create the
>>> factor
>>> crossings using formula a*b*c and then cross with the numerics,
>>> which
>>> is even more clumsy.
>>> Thanks everybody
>>>
>>> ______________________________________________
>>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> https://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
 > E-mail is sent at my convenience; I don't expect replies outside of 
working hours.