[R] Dropping terms from regression w/ poly()
Joshua Stults
joshua.stults at gmail.com
Thu Jun 4 16:07:00 CEST 2009
Hello r-help,
I'm fitting a model with lm() and using the orthogonal polynomials
from poly() as my basis:
dat <- read.csv("ConsolidatedData.csv", header=TRUE)
attach(dat)
nrows <- 1925
Rad <- poly(Radius, 2)
ntheta <- 14
Theta <- poly(T.Angle..deg., ntheta)
nbeta <- 4
Beta <- poly(B.Beta..deg., nbeta)
model.1 <- lm( Measurement ~ Block + Rad + Theta + Beta + Rad:Theta +
Rad:Beta + Theta:Beta)
Which works splendidly, and for my data set shows that the odd orders
in Theta and Beta are not significant (expected because it is a
symmetric measurement that I'm fitting). Now I want to simplify my
model by dropping all the odd-order terms. I couldn't figure out a
way to do that without defining new design matrices like this:
T2 <- matrix(0,nrows,ntheta)
for(i in 1:ntheta){
T2[,i] <- Theta[,i]
}
B2 <- matrix(0,nrows,nbeta)
for(i in 1:nbeta){
B2[,i] <- Beta[,i]
}
R2 <-matrix(0,nrows,2)
R2[,1] <- Rad[,1]
R2[,2] <- Rad[,2]
And then fitting a model using the individual columns of T2, R2 and
B2. That lets me drop all the odd terms, but it is very cumbersome:
model.2 <- lm( Measurement ~ Block + (R2[,1] + R2[,2]) + (T2[,1] +
T2[,2] + T2[,3] + T2[,4] + T2[,5] + T2[,6] + T2[,7] +
T2[,8] + T2[,9] + T2[,10] + T2[,11]+ T2[,12] + T2[,13] +
T2[,14]) + (B2[,1] + B2[,2] + B2[,3] + B2[,4]) + (R2[,1]
+ R2[,2]) : (T2[,1] + T2[,2] + T2[,3] + T2[,4] + T2[,5]
+ T2[,6] + T2[,7] + T2[,8] + T2[,9] + T2[,10] + T2[,11]+
T2[,12] + T2[,13] + T2[,14]) + (R2[,1] + R2[,2]) :
(B2[,1] + B2[,2] + B2[,3] + B2[,4]) + (T2[,1] + T2[,2] +
T2[,3] + T2[,4] + T2[,5] + T2[,6] + T2[,7] + T2[,8] +
T2[,9] + T2[,10] + T2[,11]+ T2[,12] + T2[,13] + T2[,14])
: (B2[,1] + B2[,2] + B2[,3] + B2[,4]))
because I have to call out each column explicitly. I'm new to R so I
have a suspicion that there is a much better way to accomplish this
(dropping odd terms from a model based on poly), but I haven't figured
it out. Any hints / suggestions? Thanks.
--
Joshua Stults
Website: http://j-stults.blogspot.com
More information about the R-help
mailing list