[R] inconsistent lm results with fixed response variable

[tyler]

Hi,

I'm analyzing a large number of simulations using lm(), a sample of the
resulting data is pasted below. In some simulations, the response
variable doesn't vary, ie:

> tmp[[2]]$richness
 [1] 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40

When I analyze this using R version 2.8.0 (2008-10-20) on a linux
cluster, I get an appropriate result:


## begin R ##

summary(lm(richness ~ het, data = tmp[[2]]))

Call:
lm(formula = richness ~ het, data = tmp[[2]])

Residuals:
   Min     1Q Median     3Q    Max
     0      0      0      0      0

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)       40          0     Inf   <2e-16 ***
het                0          0      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0 on 23 degrees of freedom
Multiple R-squared:      ,      Adjusted R-squared:
F-statistic:       on 1 and 23 DF,  p-value: NA

## end R ##

This is good, as when I extract the Adjusted R-squared and slope I get
NaN and 0, which are easily identified in my aggregate analysis, so I
can deal with them appropriately. 

However, this isn't always the case:

## begin R ##

 tmp[[1]]$richness
 [1] 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
[26] 40 40 40 40 40 40 40 40 40 40 40

 summary(lm(richness ~ het, data = tmp[[1]]))

Call:
lm(formula = richness ~ het, data = tmp[[1]])

Residuals:
       Min         1Q     Median         3Q        Max
-8.265e-14  1.689e-15  2.384e-15  2.946e-15  4.022e-15

Coefficients:
             Estimate Std. Error   t value Pr(>|t|)
(Intercept) 4.000e+01  8.418e-15 4.752e+15   <2e-16 ***
het         1.495e-14  4.723e-14 3.160e-01    0.754
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.44e-14 on 34 degrees of freedom
Multiple R-squared: 0.5112,     Adjusted R-squared: 0.4968
F-statistic: 35.56 on 1 and 34 DF,  p-value: 9.609e-07

## end R ##

This is a problem, as when I plot the adj. R sq as part of an aggregate
analysis of a large number of simulations, it appears to be a very
strong regression. I wouldn't have caught this except it was
exceptionally high for the simulation parameters. It also differs by
more than rounding error from the results with R 2.8.1 running on my
laptop (Debian GNU/Linux), i.e., adj. R sq 0.5042 vs 0.4968.
Furthermore, on my laptop, none of the analyses produce a NaN adj. R sq,
even for data that do produce that result on the cluster.

Both my laptop and the linux cluster have na.action set to na.omit. Is
there something else I can do to ensure that lm() returns slope == 0
and adj.R.sq == NaN when the response variable is fixed? 

Thanks for any suggestions,

Tyler

Data follows:

`tmp` <-
list(structure(list(richness = c(40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40), range = c(0.655084651733024, 0.579667533660137, 0.433092220907644, 
0.62937198839679, 0.787891987978164, 0.623511540624239, 0.542744487102066, 
0.905937570175433, 0.806802881350753, 0.680413208666325, 0.873426339019084, 
0.699982832956593, 0.697716600618959, 0.952729864926405, 0.782938474636578, 
1.03899695305995, 0.715075858219333, 0.579749205792549, 1.20648999819246, 
0.648677938600964, 0.651883559714785, 0.997318331273967, 0.926368116052012, 
0.91001274146868, 1.20737951037620, 1.12006560586723, 1.09806272133903, 
0.9750792390176, 0.356496202035743, 0.612018080768747, 0.701905693862144, 
0.735857916053381, 0.991787489781244, 1.07247435214078, 0.60061903319766, 
0.699733090379818), het = c(0.154538307084452, 0.143186508136608, 
0.0690948358402777, 0.132337152911839, 0.169037344105692, 0.117783183361602, 
0.117524251767612, 0.221161206774407, 0.204574928003633, 0.170571000779693, 
0.204489357007294, 0.131749663515638, 0.154127894997213, 0.232672587431942, 
0.198610891796736, 0.260497696582693, 0.129028191256682, 0.128717975847452, 
0.254300896783617, 0.113546727236817, 0.142220347446853, 0.24828642688332, 
0.194340945175726, 0.190782985783610, 0.214676796387244, 0.252940213066992, 
0.22362832797347, 0.182423482989676, 0.0602332226418674, 0.145400861749859, 
0.141297315445974, 0.139798699247632, 0.222815139716421, 0.211971297234962, 
0.120813579628747, 0.150590744533818), n.rich = c(40, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40)), .Names = c("richness", "range", "het", "n.rich")), 
 structure(list(richness = c(40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40, 40, 40), range = c(0.753203162648624, 0.599708526308711, 
0.714477274087683, 0.892359682406808, 0.868440625159371, 0.753239521511417, 
1.20164969658467, 1.20462111558583, 1.13142122690491, 0.95241921975703, 
1.13214481653550, 0.827528954009827, 1.14827745443481, 0.936048043180592, 
0.874649332193952, 1.38844778296649, 0.985016220913809, 1.18166853164661, 
0.784679773255876, 0.94894149080785, 0.770312904574722, 1.10203660758219, 
1.15624067277321, 0.692776967548628, 0.79343712876973), 
het = c(0.170481207967181, 
0.108265674755723, 0.123316519598517, 0.220631611141464, 0.160460967122565, 
0.145032358811883, 0.293678286125082, 0.284769842125969, 0.258637372765782, 
0.18303781265474, 0.265304220319150, 0.194784967445680, 0.248055723803990, 
0.204658616507612, 0.167203828355069, 0.287030735881294, 0.247639113771915, 
0.269348295820692, 0.111409735752589, 0.209076579513581, 0.176890183224181, 
0.249378876987384, 0.260323833307383, 0.177061093736427, 0.172263958005774
), n.rich = c(40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40)), .Names = c("richness", 
"range", "het", "n.rich")))