CODEWITHSUNDEEP
rdistributionlmpolynomials
Zhubarb
Zhubarb

Reputation: 11905

R - creating a linear model with fixed poly() coefficients

In R, one can build an lm() or glm() object with fixed coefficients, using the offset parameter in a formula.

x=seq(1,100)
y=x^2+3*x+7
# Forcing to fit the polynomial: 2x^2 + 4x + 8 
fixed_model = lm(y ~ 0 + offset(8 + 4*x + 2*I(x^2) ))

Is it possible to do the same thing using poly()? I tried the code below but it doesn't seem to work. 

fixed_model_w_poly <- lm(y ~ 0 + offset(poly(x, order=2, raw = TRUE, coefs= c(8, 4, 2))))

Error : number of offsets is 200, should equal 100 (number of observations)

I want to use poly() as a convenient interface to run iterations with a high number of fixed coefficients or order values, rather than having to manually code: offset(8 + 4*x + 2*I(x^2) ) for each order/coefficient combination.

P.S: Further but not essential information: This is to go inside an MCMC routine. So an example usage would be to generate (and then compare) model_current to model_next in the below code:

library(MASS)
coeffs_current <- c(8, 4, 2)
model_current <- lm(y ~ 0 + offset(poly(x, order=2, raw = TRUE, coefs= coeffs_current )))
cov <- diag(rep(1,3))
coeffs_next <- mvrnorm(1, mu = as.numeric(coeffs_current ),
                       Sigma = cov ) 
model_next <- lm(y ~ 0 + offset(poly(x, order=2, raw = TRUE, coeffs_next ))

Upvotes: 1

Views: 604

Answers (1)

IRTFM
IRTFM

Reputation: 263421

This demonstrates what I suggested. (Not to use poly.)

library(MASS)
# coeffs_current <- c(8, 4, 2) Name change for compactness.
cc <- c(8, 4, 2)
form <- as.formula(bquote(y~x+offset(.(cc[1])+x*.(cc[2])+.(cc[3])*I(x^2) )))
model_current <- lm(form, data=dat))

I really have no idea what you intend to do with this next code. Looks like you want something based on the inputs to the prior function, but doesn't look like you want it based on the results.

cov <- diag(rep(1,3))
coeffs_next <- mvrnorm(1, mu = as.numeric(cc ),
                       Sigma = cov )

The code works (at least as I intended) with a simple test case. The bquote function substitutes values into expressions (well actually calls) and the as.formula function evaluates its argument and then dresses the result up as a proper formula-object.

dat <- data.frame(x=rnorm(20), y=rnorm(20) )
cc <- c(8, 4, 2)
form <- as.formula( bquote(y~x+offset(.(cc[1])+x*.(cc[2])+.(cc[3])*I(x^2) )))
model_current <- lm(form, data=dat)
#--------
> model_current

Call:
lm(formula = form, data = dat)

Coefficients:
(Intercept)            x  
     -9.372       -5.326    # Bizarre results due to the offset.
#--------
form
#y ~ x + offset(8 + x * 4 + 2 * I(x^2))

Upvotes: 1

Related Questions