[R] storing lm() results and other objects in a list

Idgarad idgarad at gmail.com
Wed Jul 15 21:17:57 CEST 2009


to clean up some code I would like to make a list of arbitrary length
to store various objects for use in a loop

sample code:


############ BEGIN SAMPLE ##############
# You can see the need for a loop already
linearModel1=lm(modelSource ~ .,mcReg)
linearModel2=step(linearModel1)
linearModel3=lm(modelSource ~ .-1,mcReg)
linearModel4=step(linearModel3)
#custom
linearModel5=lm(modelSource ~ . -ACF-MonthlyST1-MonthlyST2-MonthlyBLA,mcReg)

LinearModel1.res <- residuals(linearModel1)
LinearModel2.res <- residuals(linearModel2)
LinearModel3.res <- residuals(linearModel3)
LinearModel4.res <- residuals(linearModel4)
LinearModel5.res <- residuals(linearModel5)

#hmmm bolt on linearModel[x] as linearModel[x]$arma.fit?
arma1.fit <- auto.arima(LinearModel1.res)
arma2.fit <- auto.arima(LinearModel2.res)
arma3.fit <- auto.arima(LinearModel3.res)
arma4.fit <- auto.arima(LinearModel4.res)
arma5.fit <- auto.arima(LinearModel5.res,stepwise=T,trace=T)

#Ok what is left over after Regression and ARIMA that cannot
#be explained. Stupid outliers....
#AO's can be added to the cReg as a normal dummy variable
# but these are AOs from the model not the original data.
# is it better to handle AOs from the original data?

#linearModel[x]arma.ao?
arma1.ao <- detectAO(arma1.fit)
arma2.ao <- detectAO(arma2.fit)
arma3.ao <- detectAO(arma3.fit)
arma4.ao <- detectAO(arma4.fit)
arma5.ao <- detectAO(arma5.fit)

#What do I do with an innovative outlier? Transfer function or what?
#auto.arima doesn't handle the IO=c(...) stuff.... Umm...
#transfer functions, etc. are a deficency in the script at this point....

#linearModel[x]arma.io?
arma1.io <- detectIO(arma1.fit)
arma2.io <- detectIO(arma2.fit)
arma3.io <- detectIO(arma3.fit)
arma4.io <- detectIO(arma4.fit)
arma5.io <- detectIO(arma5.fit)

#Sample on how to auto-grab regressors from DetectAO and DetectIO and
#appened them to our regression array. You'd have to do this for each model
#as the residuals are where the outliers are coming from and diff models
#would have different residuals left over. IO is best left to arimax functions
#directly. I assume at this point that AO's can be added to Regression tables
#if that is the case then REM out the IO lines and pass the detectIO results

#into the arimax(x,y,z,IO=detectIO(blah))

#
# Need a better understanding of how to address the AO and IO's in
this script before implementing them
# (Repeat for each model, cReg1,cReg2,etc..)
#
#cReg1=cReg
#fReg1=fReg
#for(i in arma1.io$ind){ print(i);cReg1[,paste(sep="
","IO",i)]=1*(seq(cReg1[,2])==i)}
#for(i in arma1.ao$ind){ print(i);cReg1[,paste(sep="
","AO",i)]=1*(seq(cReg1[,2])==i)}
#for(i in arma1.io$ind){ print(i);fReg1[,paste(sep="
","IO",i)]=1*(seq(fReg1[,2]))}
#for(i in arma1.ao$ind){ print(i);fReg1[,paste(sep="
","AO",i)]=1*(seq(fReg1[,2]))}


#Get the pdq,PDQs into a variable so we can re-feed it if neccessary
#oh crap absorbing this into LinearModel[x] looks ugly for syntax....
arma1.fit$order=c(arma1.fit$arma[1],arma1.fit$arma[2],arma1.fit$arma[6])
arma2.fit$order=c(arma2.fit$arma[1],arma2.fit$arma[2],arma2.fit$arma[6])
arma3.fit$order=c(arma3.fit$arma[1],arma3.fit$arma[2],arma3.fit$arma[6])
arma4.fit$order=c(arma4.fit$arma[1],arma4.fit$arma[2],arma4.fit$arma[6])
arma5.fit$order=c(arma5.fit$arma[1],arma5.fit$arma[2],arma5.fit$arma[6])

arma1.fit$seasonal=c(arma1.fit$arma[3],arma1.fit$arma[4],arma1.fit$arma[7])
arma2.fit$seasonal=c(arma2.fit$arma[3],arma2.fit$arma[4],arma2.fit$arma[7])
arma3.fit$seasonal=c(arma3.fit$arma[3],arma3.fit$arma[4],arma3.fit$arma[7])
arma4.fit$seasonal=c(arma4.fit$arma[3],arma4.fit$arma[4],arma4.fit$arma[7])
arma5.fit$seasonal=c(arma5.fit$arma[3],arma5.fit$arma[4],arma5.fit$arma[7])

#these Two are used for linearModel2 and linearModel4, Get only the
#regressors that surived step removal.
newcReg=cReg[match(names(linearModel2$coeff[-1]),names(cReg))]
newfReg=fReg[match(names(linearModel2$coeff[-1]),names(fReg))]
newmcReg=mcReg[match(names(linearModel2$coeff[-1]),names(mcReg))]
newmfReg=mfReg[match(names(linearModel2$coeff[-1]),names(mfReg))]

#Scenario 1 - All Regressors Left In
newFit1.b <- Arima(modelSource,order=arma1.fit$order,seasonal=list(order=arma1.fit$seasonal),xreg=mcReg,include.drift=F)

#Scenario 2 - Step Removal of Regressors
newFit2.b <- Arima(modelSource,order=arma2.fit$order,seasonal=list(order=arma2.fit$seasonal),xreg=newmcReg,include.drift=F)

#Scenario 3 - All Regressors Left In with Intercept Removed
newFit3.b <- Arima(modelSource,order=arma3.fit$order,seasonal=list(order=arma3.fit$seasonal),xreg=mcReg,include.drift=F)

#Scenario 4 - Step Removal of Regressors with Intercept Removed (I
have a feeling this is identical to #2 in results
newFit4.b <- Arima(modelSource,order=arma4.fit$order,seasonal=list(order=arma4.fit$seasonal),xreg=newmcReg,include.drift=F)

#Scenario 5 - Robust1, For giggles and grins for now
newFit5.b <- Arima(modelSource,order=arma5.fit$order,seasonal=list(order=arma5.fit$seasonal),xreg=newmcReg,include.drift=F)

#All the drift terms ones are broke as the drift term doesn't match up
to the # of columns in fReg
#Still don't know how to fix at this time.
forecast1.b <-predict(newFit1.b,n.ahead=modLength-modMax,newxreg=mfReg)
forecast2.b <-predict(newFit2.b,n.ahead=modLength-modMax,newxreg=newmfReg)
forecast3.b <-predict(newFit3.b,n.ahead=modLength-modMax,newxreg=mfReg)
forecast4.b <-predict(newFit4.b,n.ahead=modLength-modMax,newxreg=newmfReg)
forecast5.b <-predict(newFit5.b,n.ahead=modLength-modMax,newxreg=newmfReg)

accuracy(forecast1.b$pred,verSource)
accuracy(forecast2.b$pred,verSource)
accuracy(forecast3.b$pred,verSource)
accuracy(forecast4.b$pred,verSource)
accuracy(forecast5.b$pred,verSource)

####### END SAMPLE #############

It's pretty clear I need to clean this up into a loop for each model
(as I like having an arbitrary number of models possible but R is not
very kind in trying to wrap my head around the construction of an
array of arrays (or list of lists)) and the subsequent matching.

For those curious this is to automatically forecast mainframe usage
based on a provided calendar of weekly events that are active (the
regressors, I have no control of the number of entries in the
calendar). It works fine in the existing form (sans the whole outlier
handling... haven't learned how to handle those yet, all self taught
so far...) but I need to clean this up into something useable in a
loop (ultimately I want to store models in a database and read them
into this).

Suggestions on how to handle the whole linearModel and the child data?




More information about the R-help mailing list