-1

I'm trying to extract OLS coefficient from a set of model that run under lapply:

The problem that not all submodel from my list have all levels and most of the time I end up with coefficients "out of bound". For example in the code below "No Reward" option is not present for Reward2 and result0$coef[3,1] will be out of bound as there is no third estimate reported. Question: is it possible to force "lm" to report all the coefficients specified in the model even if there is no estimate available?

I would like to apologize to the community for not presenting a reproducible code on my earlier attempt. Since my last attempt I solved the problem by checking within the function for presence of a particular estimate, but the question still remains and here is the code:

RewardList<-c("Reward1","Reward2")
set.seed(1234)
GL<-rnorm(10)
RewardName<-rep(c("Reward1","Reward2"),each=5)
NoReward<-c(0,1,0,1,0,  0,0,0,0,0)
Under1<-  c(1,0,0,0,1,  0,1,0,1,0)
Above1<-  c(0,0,1,0,0,  1,0,1,0,1)

tinput<-as.data.frame(RewardName)
tinput<-cbind(tinput,GL,NoReward,Under1,Above1)

regMF <- lapply(seq_along(RewardList), 
function (n) {
tinput <- tinput[tinput$RewardName==RewardList[n],]
result0 <- summary(lm(GL~NoReward+Under1+Above1-1,tinput))
result1 <- result0$coef[1,1] #no rebate
result2 <- result0$coef[2,1]
result3 <- result0$coef[3,1]
return(list(result1,result2,result3))})
Konstantin Mingoulin
  • 161
  • 1
  • 8

2 Answers2

1

To make your function more robust, you could use tryCatch or use dplyr::failwith. For instance, replacing 

result0 <- summary(lm(GL~`No Reward`+`Under 1%`-1,tinput))

with

require(dplyr)
result0 <- failwith(NULL, summary(lm(GL~`No Reward`+`Under 1%`-1,tinput)))

may work, though it is difficult to tell without a reproducable example. On how to produce a reproducable example, please have a look here.

Although not exactly your question, I would like to point out a potentially easier way to arrange your data using broom, which avoids the list of lists-output structure of your approach. With broom and dplyr, you can collect the output of your models in a dataframe for easier access. For instance, have a look at the output of 

 library(dplyr) 
 library(broom)
 mtcars %>% group_by(gear) %>% do(data.frame(tidy(lm(mpg~cyl, data=.), conf.int=T)))

Here, you can wrap the lm function around failwith as well.

Community
  • 1
  • 1
coffeinjunky
  • 11,254
  • 39
  • 57
1

As I said, use lmList:

library(nlme)
fits <- lmList(GL~NoReward+Under1+Above1-1 | RewardName, data = tinput)
coef(fits)

#         NoReward     Under1     Above1
#Reward1 -1.034134 -0.3889705  1.0844412
#Reward2        NA -0.5695960 -0.3102046
Roland
  • 127,288
  • 10
  • 191
  • 288