1
$\begingroup$

I'm trying to run the lmrob() function in the R package robustbase, using data from a mixed factorial design study. My dependent variable is repeated (pre+ post) and I have two independent variables. I'm entering the post-measurement as the dependent variable, and controlling for the pre-measurement.

The function I'm trying to run is basically:

lmrob(Post ~ pre +IV1 + IV2, setting = "KS2014", data=mydataframe) 

When I run this, I get an error message saying:

S-estimated scale == 0:  Probably exact fit; check your data. 

When I don't control for the pre-measurement, I don't get this warning message. So I assume it has something to do with the repeated measures being correlated? Does anyone have any idea what might be causing the issue and how to resolve it?

This is the first time I'm posting here so please bear with me if I've not included enough detail or my question isn't clear enough, I will provide clarification if needed! I've copied in the output from the console, the session info, and the relevant extract from the dataset below.


Call:
lmrob(formula = Post_Intention ~ TPB_vs_no_TPB + Tailored_vs_Untailored + 
    Pre_Intention, data = dfa, setting = "KS2014")
 \--> method = "S"
Residuals:
   Min     1Q Median     3Q    Max 
   -80      0      0     10     90 

Exact fit detected

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)
(Intercept)                             0          0      NA       NA
TPB_vs_no_TPBNo TPB                     0          0      NA       NA
Tailored_vs_UntailoredUntailored        0          0      NA       NA
Pre_Intention                           1          0      NA       NA

Robustness weights: 
 88 observations c(1,5,10,13,16,19,21,25,27,28,29,31,32,34,36,40,42,50,53,56,57,61,66,68,75,79,83,86,87,93,98,104,108,112,116,119,124,125,128,134,135,136,138,139,144,145,148,153,156,161,164,165,171,173,175,176,177,179,182,184,190,193,197,198,199,203,213,215,216,217,218,219,220,226,232,233,235,236,239,241,242,243,245,246,247,250,251,252)
     are outliers with |weight| = 0 ( < 0.0004); 
 165 weights are ~= 1.
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4                bb       tuning.psi1       tuning.psi2 
       -5.000e-01         1.500e+00                NA         5.000e-01         5.000e-01        -5.000e-01         1.500e+00 
      tuning.psi3       tuning.psi4        refine.tol           rel.tol         solve.tol       eps.outlier             eps.x 
        9.500e-01                NA         1.000e-07         1.000e-07         1.000e-07         3.953e-04         1.819e-10 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max    maxit.scale      trace.lev            mts 
          1000            500             20              2           2000            200              0           1000 
    compute.rd fast.s.large.n 
             0           2000 
              setting                   psi           subsampling                   cov compute.outlier.stats 
             "KS2014"                 "lqq"         "nonsingular"             ".vcov.w"                "SMDM" 
seed : int(0) 
> source('~/R scripts/lmrob.R')
re-encoding from UTF-8
Warning message:
In lmrob.S(x, y, control = control, mf = mf) :
  S-estimated scale == 0:  Probably exact fit; check your data
> sessionInfo
function (package = NULL) 
{
    z <- list()
    z$R.version <- R.Version()
    z$platform <- z$R.version$platform
    if (nzchar(.Platform$r_arch)) 
        z$platform <- paste(z$platform, .Platform$r_arch, sep = "/")
    z$platform <- paste0(z$platform, " (", 8 * .Machine$sizeof.pointer, 
        "-bit)")
    z$locale <- Sys.getlocale()
    if (.Platform$OS.type == "windows") {
        z$running <- win.version()
    }
    else if (nzchar(Sys.which("uname"))) {
        uname <- system("uname -a", intern = TRUE)
        os <- sub(" .*", "", uname)
        z$running <- switch(os, Linux = if (file.exists("/etc/os-release")) {
            tmp <- readLines("/etc/os-release")
            t2 <- if (any(startsWith(tmp, "PRETTY_NAME="))) sub("^PRETTY_NAME=", 
                "", grep("^PRETTY_NAME=", tmp, value = TRUE)[1L]) else if (any(startsWith(tmp, 
                "NAME"))) sub("^NAME=", "", grep("^NAME=", tmp, 
                value = TRUE)[1L]) else "Linux (unknown distro)"
            sub("\\"(.*)\\"", "\\\\1", t2)
        } else if (file.exists("/etc/system-release")) {
            readLines("/etc/system-release")
        }, Darwin = {
            ver <- readLines("/System/Library/CoreServices/SystemVersion.plist")
            ind <- grep("ProductUserVisibleVersion", ver)
            ver <- ver[ind + 1L]
            ver <- sub(".*<string>", "", ver)
            ver <- sub("</string>$", "", ver)
            ver1 <- strsplit(ver, ".", fixed = TRUE)[[1L]][2L]
            sprintf("%s %s %s", ifelse(as.numeric(ver1) < 12, 
                "OS X", "macOS"), switch(ver1, `6` = "Snow Leopard", 
                `7` = "Lion", `8` = "Mountain Lion", `9` = "Mavericks", 
                `10` = "Yosemite", `11` = "El Capitan", `12` = "Sierra", 
                `13` = "High Sierra", ""), ver)
        }, SunOS = {
            ver <- system("uname -r", intern = TRUE)
            paste("Solaris", strsplit(ver, ".", fixed = TRUE)[[1L]][2L])
        }, uname)
    }
    if (is.null(package)) {
        package <- grep("^package:", search(), value = TRUE)
        keep <- vapply(package, function(x) x == "package:base" || 
            !is.null(attr(as.environment(x), "path")), NA)
        package <- .rmpkg(package[keep])
    }
    pkgDesc <- lapply(package, packageDescription, encoding = NA)
    if (length(package) == 0) 
        stop("no valid packages were specified")
    basePkgs <- sapply(pkgDesc, function(x) !is.null(x$Priority) && 
        x$Priority == "base")
    z$basePkgs <- package[basePkgs]
    if (any(!basePkgs)) {
        z$otherPkgs <- pkgDesc[!basePkgs]
        names(z$otherPkgs) <- package[!basePkgs]
    }
    loadedOnly <- loadedNamespaces()
    loadedOnly <- loadedOnly[!(loadedOnly %in% package)]
    if (length(loadedOnly)) {
        names(loadedOnly) <- loadedOnly
        pkgDesc <- c(pkgDesc, lapply(loadedOnly, packageDescription))
        z$loadedOnly <- pkgDesc[loadedOnly]
    }
    z$matprod <- as.character(options("matprod"))
    es <- extSoftVersion()
    z$BLAS <- as.character(es["BLAS"])
    z$LAPACK <- La_library()
    class(z) <- "sessionInfo"
    z
}
<bytecode: 0x0000000017538350>
<environment: namespace:utils>

Dataset

dput(dfa) structure(list(Pre_Intention = c(50, 10, 50, 100, 80, 100, 10, 0, 100, 50, 90, 100, 50, 100, 100, 50, 100, 100, 30, 0, 10, 0, 60, 10, 0, 0, 0, 0, 40, 10, 50, 0, 70, 50, 50, 20, 0, 50, 0, 70, 50, 60, 0, 100, 50, 100, 100, 100, 0, 50, 0, 10, 0, 60, 0, 50, 50, 20, 100, 100, 50, 100, 100, 0, 60, 20, 100, 10, 50, 100, 100, 100, 0, 50, 10, 0, 80, 50, 30, 100, 100, 80, 40, 100, 0, 10, 20, 100, 100, 100, 100, 0, 50, 100, 70, 100, 40, 20, 100, 100, 100, 100, 0, 50, 100, 100, 0, 50, 50, 90, 20, 50, 10, 40, 50, 20, 50, 0, 0, 100, 100, 10, 10, 80, 10, 100, 80, 50, 100, 10, 100, 100, 90, 50, 50, 50, 100, 50, 40, 100, 100, 40, 20, 10, 10, 20, 100, 20, 100, 0, 100, 0, 30, 0, 10, 10, 40, 0, 40, 0, 50, 0, 0, 30, 0, 0, 50, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 50, 20, 100, 50, 20, 90, 20, 100, 70, 20, 90, 10, 20, 80, 100, 90, 100, 100, 100, 10, 10, 20, 0, 60, 0, 10, 100, 100, 100, 10, 10, 0, 20, 0, 0, 40, 0, 30, 0, 40, 70, 0, 0, 30, 40, 0, 10, 20, 20, 50, 10, 10, 0, 0, 10, 0, 50, 0, 40, 0, 20, 0, 0, 20, 100, 0, 100, 30, 0, 50, 70, 10, 0, 0, 30, 10), Post_Intention = c(70, 10, 50, 100, 0, 100, 10, 0, 100, 90, 90, 100, 100, 100, 100, 80, 100, 100, 70, 0, 30, 0, 60, 10, 10, 0, 40, 20, 70, 10, 80, 40, 70, 40, 50, 30, 0, 50, 0, 100, 50, 100, 0, 100, 50, 100, 100, 100, 0, 60, 0, 10, 40, 60, 0, 100, 80, 20, 100, 100, 90, 100, 100, 0, 60, 30, 100, 0, 50, 100, 100, 100, 0, 50, 20, 0, 80, 50, 100, 100, 100, 80, 50, 100, 0, 0, 30, 100, 100, 100, 100, 0, 60, 100, 70, 100, 40, 30, 100, 100, 100, 100, 0, 100, 100, 100, 0, 100, 50, 90, 20, 70, 10, 40, 50, 50, 50, 0, 10, 100, 100, 10, 10, 20, 40, 100, 80, 80, 100, 10, 100, 100, 90, 80, 70, 70, 100, 80, 50, 100, 100, 40, 20, 100, 50, 20, 100, 30, 100, 0, 100, 0, 60, 0, 10, 20, 40, 0, 40, 0, 60, 0, 0, 10, 10, 0, 50, 0, 0, 0, 30, 10, 30, 0, 0, 30, 10, 50, 80, 100, 50, 50, 90, 0, 100, 70, 20, 90, 10, 30, 80, 100, 100, 100, 100, 100, 20, 20, 40, 0, 60, 0, 90, 100, 100, 100, 10, 10, 0, 20, 0, 0, 60, 0, 50, 20, 70, 80, 20, 20, 30, 40, 0, 10, 20, 40, 50, 10, 10, 0, 0, 20, 60, 50, 10, 60, 0, 20, 50, 0, 60, 70, 10, 100, 40, 30, 100, 70, 10, 30, 20, 40, 10), TPB_vs_no_TPB = structure(c(1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L), .Label = c("TPB", "No TPB"), class = "factor"), Tailored_vs_Untailored = structure(c(2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L), .Label = c("Tailored", "Untailored"), class = "factor")), .Names = c("Pre_Intention", "Post_Intention", "TPB_vs_no_TPB", "Tailored_vs_Untailored"), class = "data.frame", row.names = c(NA, -253L))

$\endgroup$
7
  • $\begingroup$ It looks like your pre and post measurements are perfectly correlated. Is that true? If not then more details are needed about your data-set. $\endgroup$
    – mdewey
    Commented Sep 25, 2017 at 17:17
  • $\begingroup$ No they are not perfectly correlated; the correlation is r=0.895. Does that help? If not let me know what further details you need to know about the dataset. $\endgroup$
    – J. Mueller
    Commented Sep 25, 2017 at 18:10
  • $\begingroup$ What happens if you fit the model just with the pre-score as the sole predictor? $\endgroup$
    – mdewey
    Commented Sep 26, 2017 at 10:09
  • $\begingroup$ Thanks for the response. I just ran the function again with only the pre-score as a predictor, and got the same error message... Is it possible R can't handle when the two scores are highly correlated (even if not quite perfect)? $\endgroup$
    – J. Mueller
    Commented Sep 27, 2017 at 17:59
  • 2
    $\begingroup$ There is of course a logical problem here as if the pre and post are so closely related your other variables have little to predict and what they do have may be mostly measurement error. $\endgroup$
    – mdewey
    Commented Sep 29, 2017 at 15:42

2 Answers 2

2
$\begingroup$

The initial S estimate is based on random sampling. It has itself quite a few tuning parameters, see help(lmrob.control) and probably ?lmrob.S . As @mdewey was also thinking, it could be that too many subsamples (of the random sample) gave perfect fits. I (as maintainer of robustbase) would be happy to investigate in more detail (when back at work), but for that I (or others helping) need to be able to get your data (or subset of your data which gives the same error). If it's small enough, you could use to dput(<dataframe>) and paste the result here.

Even before that you could try to use the currently most recommended setting="KS2014" and see if it helps (it already tries to use "better" S estimate tuning). Last but not least: Are you using the most current version of robustbase?

$\endgroup$
1
  • $\begingroup$ Thanks very much for your response. Regarding your queries/suggestions: I am using robustbase version 0.92-7, which I believe is the most up to date version. Also, I added the setting="KS2014" but this did not resolve the issue. I've copied in the output from dput() into my question above; if you are able to have a look at this once you are back at work, that would be great! I appreciate it. Many thanks. $\endgroup$
    – J. Mueller
    Commented Oct 4, 2017 at 7:52
1
$\begingroup$

Note you did not get an error, but a warning... and that's why you did get a result on which you could call summary(*) and got what you did show as output.

As a matter of fact that is already sufficient: You got 88 robustness weights == 0, and a perfect fit for the remaining observations (which had weights == 1: You can find that by using weight(fm, type = "robustness") where fm denotes the fitted model, i.e., fm <- lmrob(Post.. ~ Pre + .., dfa)`.

Your summary(fm) output indeed shows that the estimated coefficients are '1' for Pre and 0 for all other "predictors" (including the intercept). If you look more closely you can notice that from a total of 253 observations, in 165 cases, pre == post, identically, and hence a model where you reject 88 observations as outliers and use the remaining 165 ones gives a perfect prediction.

--- all clearly no error, and 'Ok' actually.

Now if you are more interested you can do the following,

fm1 <- lmrob(Post_Intention ~ ., setting = "KS2014", data = dfa)

Warning message:
In lmrob.S(x, y, control = control, mf = mf) :
S-estimated scale == 0:  Probably exact fit; check your data

Robustness weights, directly:

table(rw <- weights(fm1, type = "robustness"))

0   1 

88 165 #total: 253 observations

Model if we only look at the 88 cases that were weighed down to 0 (aka "thrown away"):

fms <- lmrob(Post_Intention ~ ., setting = "KS2014", data = dfa,
         subset = rw == 0)
summary(fms)

Call:
lmrob(formula = Post_Intention ~ ., data = dfa, subset = rw == 0, setting = "KS2014")
 \--> method = "SMDM"
Residuals:
Min      1Q  Median      3Q     Max 
-97.494  -9.791  -2.203  11.692  70.635 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      26.26067    4.98925   5.263 1.07e-06 ***
Pre_Intention                     0.96748    0.09624  10.053 4.56e-16 ***
TPB_vs_no_TPBNo TPB              -6.57097    4.45877  -1.474    0.144    
Tailored_vs_UntailoredUntailored  0.40615    4.52753   0.090    0.929    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 20.6 
Multiple R-squared:  0.5511,    Adjusted R-squared:  0.535 
Convergence in 12 IRWLS iterations

Robustness weights: 
 66 weights are ~= 1. The remaining 22 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0956  0.6930  0.8942  0.7596  0.9297  0.9981 
Algorithmic parameters: 
   .........
   .........
> 

This may (or may not) be interesting: Even for these cases, only the pre_intention has a significant effect.

$\endgroup$
1
  • $\begingroup$ Thank you for your helpful and detailed response. Could you help me understand why the model excludes the 88 responses? I understand that the model views them as outliers but in this case I don't think they should be excluded. A little background info: I was looking into ways to increase participants' intention to seek help from a doctor. It looks like for the majority (n=165), there was no change in intention. However for 88 participants, the intention changed: In only 10 cases it decreased, whereas in 78 cases it increased. I think this indicates that these weren't random outliers. $\endgroup$
    – J. Mueller
    Commented Oct 22, 2017 at 15:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.