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@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.