I could be missing something -- but my multiple linear regression is seemingly working as expected:

> bigmodel <- lm(score ~ lean + gender + age, data = mydata)
> summary(bigmodel)

Call:
lm(formula = score ~ lean + gender + age, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.891  -4.354   0.892   6.240  18.537 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 70.96455    3.85275  18.419   <2e-16 ***
lean         0.62463    0.05938  10.518   <2e-16 ***
genderM     -2.24025    1.40362  -1.596   0.1121    
age          0.10783    0.06052   1.782   0.0764 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9 on 195 degrees of freedom
Multiple R-squared:  0.4188,    Adjusted R-squared:  0.4098 
F-statistic: 46.83 on 3 and 195 DF,  p-value: < 2.2e-16

> head(predict(bigmodel),20)
       1        2        3        4        5        6        7        8        9       10 
75.36711 74.43743 77.02533 78.76903 79.95515 79.09251 80.38647 81.65807 80.14846 78.96234 
      11       12       13       14       15       16       17       18       19       20 
82.39052 82.04468 81.05187 81.26753 84.50240 81.80667 80.92169 82.40895 81.76197 82.94809

But I can't wrap my head around the prediction after reading ?predict.lm. This output looks good to me for my original dataset -- but what if I want to run the prediction against a different dataset than the one I used to create bigmodel?

For example, if I import a .csv file into R called newmodel with 200 people complete with leans, gender, and age -- how can I use the regression formula from bigmodel to produce predictions for newmodel?

Thanks!