[R] inconsistent lm results with fixed response variable

[Rolf Turner]

Oh for Pete's sake!

Computers use floating point arithmetic.  Your residual standard  
error in
case 2 (i.e. 1.44e-14) *is* 0, but floating point arithmetic can't  
quite see
that this is so. Put in a check for the RSE being 0, and ``over- 
ride'' the
adjusted R squared to be NA (or NaN, or whatever floats your boat) in  
such instances.
The all.equal() function might be useful to you:

 > x <- 1.44e-14
 > all.equal(x,0)
[1] TRUE

(Caution:  Trap for Young Players:  If x and y are ``really'' different,
then all.equal(x,y) doesn't return FALSE as you might expect, but rather
a description of the difference between x and y --- which may be  
complicated
if x and y are complicated objects.  The function isTRUE() is useful  
here.)

	cheers,

		Rolf Turner


On 21/01/2009, at 9:21 AM, tyler wrote:

> 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")))
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting- 
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

######################################################################
Attention: 
This e-mail message is privileged and confidential. If you are not the 
intended recipient please delete the message and notify the sender. 
Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal 
www.marshalsoftware.com
######################################################################