R Under development (unstable) (2024-08-17 r87027 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## recreating exact fit problem for categorical data > require(robustbase) Loading required package: robustbase > > ## some simple balanced dataset with one grouping variable > ngrp <- 10 > nrep <- 10 > set.seed(2) > data <- data.frame(y = rnorm(ngrp*nrep), grp=rep(letters[1:ngrp], each=nrep)) > ## this works fine > m1 <- lmrob(y ~ grp, data) > > ## now contaminate the dataset > data2 <- data > data2$y[1:48] <- 1e10 > try(m2 <- lmrob(y ~ grp, data2, trace.lev = 3)) lmrob_S(n = 100, nRes = 500): fast_s() [non-large n]: fast_s(*, s_y=4.8e+09, n=100, p=10, ipsi=1, ..) before INIT_WLS(): Subsampling 500 times to find candidate betas: refine_fast_s(s0=-1, convChk=FALSE): too many zeroes -> scale=0 & quit refinement Sample[ 0]: after refine_(*, conv=F): beta_ref : 1e+10 0 0 0 0 -1e+10 -1e+10 -1e+10 -1e+10 -1e+10 with ||beta_ref - beta_cand|| = 0, --> sc = 0 |s|=0: Have 62 (too many) exact zeroes -> leaving refinement! lmrob.S(): scale = 0; coeff.= [1] 1e+10 0e+00 0e+00 0e+00 0e+00 -1e+10 -1e+10 -1e+10 -1e+10 -1e+10 init *NOT* converged; init$scale = 0, init$coef: (Intercept) grpb grpc grpd grpe grpf 1e+10 0e+00 0e+00 0e+00 0e+00 -1e+10 grpg grph grpi grpj -1e+10 -1e+10 -1e+10 -1e+10 Warning messages: 1: In lmrob.S(x, y, control = control) : S-estimated scale == 0: Probably exact fit; check your data 2: In lmrob.fit(x, y, control, init = init) : initial estim. 'init' not converged -- will be return()ed basically unchanged > ## All observations of group "e" get rob. weight of 0: > weights <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1) ## from trace output > weights %*% m1$x (Intercept) grpb grpc grpd grpe grpf grpg grph grpi grpj [1,] 90 10 10 10 0 10 10 10 10 10 > > proc.time() user system elapsed 0.25 0.07 0.31