my.rlm scale, RFPE use of control

[harrysouthworth]

my.rlm returns a different estimate of scale to robustbase::lmrrob

Also, RFPE doesn't handle control properly. Calls to my.rlm within RFPE pass in rlm.control() rather than the user supplied object.

liver <- texmex::liver rbc <- robustbase::lmrob.control(tuning.psi=3.443689)

x <- model.matrix(~ log(ALT.B), data=liver) y <- log(liver$ALT.M)

robustbase::lmrob(log(ALT.M) ~ log(ALT.B), data=liver, control=rbc)$scale FRBmodelselection::my.rlm(x, y, control=rbc)$scale

FRBmodelselection::RFPE function (x, y, model, control = rlm.control(), sigma.full) { if (missing(sigma.full)) sigma.full <- my.rlm(x, y, control = rlm.control())$scale # <------ SHOULD BE control = control

These are part of a bigger issue that I'm trying figure out. RobStatTM::lmdetMM.RFPE, FBRmodelselection::RFPE and robust::lmRob.RFPE give different results that don't appear to be due to scaling factors (my own implementation of RFPE appears to give similar results to robust::lmRob.RFPE, but that doesn't mean that it's right).

Thanks, Harry

[msalibian]

Thanks Harry. The short answer is that my.rlm() is based on an old version of the original code behind robustbase::lmrob(). I will look at this carefully next week, but one possible difference between these different versions of the code may be on the definition of the residual M-scale estimator--sometimes it is defined as the solution \sigma to mean( rho( resid / sigma) ) = 1/2, where rho is the bi-square function with tuning constant "tuning.chi", and in some other cases it is defined as the solution to sum( rho( resid / sigma) ) / (n - p) = 1/2. This may not solve your problem, but I'd check it first, if you haven't yet.

And thanks for catching the bug in FRBmodelselection::RFPE, I'll fix it soon.

[msalibian]

Indeed, a quick check confirms the difference that I mention above between the residual scale M-estimator in robustbase::lmrob and the one in FRBmodelselection::my.rlm, both are asymptotically equivalent, but generally slightly different:

liver <- texmex::liver
rbc <- robustbase::lmrob.control(tuning.psi=3.443689)
x <- model.matrix(~ log(ALT.B), data=liver)
y <- log(liver$ALT.M)
(tmp0 <- robustbase::lmrob(log(ALT.M) ~ log(ALT.B), data=liver, control=rbc))$scale
(tmp <- FRBmodelselection::my.rlm(x, y, control=rbc))$scale

# do the residual M-scale estimators satisfy their "usual" definition?

# re-scaled rho() function so that sup_x rho(x) = 1
rho <- function(a, cc) Rho(x=a, cc=cc) * 6 / cc^2

# residuals of the FRBmodelselection::my.rlm fit
re <- as.vector(y -  x %*% tmp$coef.s )
mean( rho(re / tmp$scale, cc=rbc$tuning.chi) ) # * 6 / rbc$tuning.chi^2 ) 

# residuals of the robustbase::lmrob fit
re0 <- as.vector(y -  x %*% tmp0$init.S$coef )
# robustbase::lmrob uses a slightly different definition for
# the residual M-scale estimator, in particular:
mean( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) # =! 1/2 !!
# but
sum( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) / (n-2)

[harrysouthworth]

Thanks for looking into this. It's appreciated. Harry

On 19 October 2018, at 17:56, Matias Salibian Barrera notifications@github.com wrote:

Indeed, a quick check confirms the difference that I mention above between the residual scale M-estimator in robustbase::lmrob and the one in FRBmodelselection::my.rlm, both are asymptotically equivalent, but generally slightly different:

liver <- texmex::liver rbc <- robustbase::lmrob.control(tuning.psi=3.443689) x <- model.matrix(~ log(ALT.B), data=liver) y <- log(liver$ALT.M) (tmp0 <- robustbase::lmrob(log(ALT.M) ~ log(ALT.B), data=liver, control=rbc))$scale (tmp <- FRBmodelselection::my.rlm(x, y, control=rbc))$scale # do the residual M-scale estimators satisfy their "usual" definition? # re-scaled rho() function so that sup_x rho(x) = 1 rho <- function(a, cc) Rho(x=a, cc=cc) 6 / cc^2 # residuals of the FRBmodelselection::my.rlm fit re <- as.vector(y - x %% tmp$coef.s ) mean( rho(re / tmp$scale, cc=rbc$tuning.chi) ) # 6 / rbc$tuning.chi^2 ) # residuals of the robustbase::lmrob fit re0 <- as.vector(y - x %% tmp0$init.S$coef ) # robustbase::lmrob uses a slightly different definition for # the residual M-scale estimator, in particular: mean( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) # =! 1/2 !! # but sum( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) / (n-2)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/msalibian/FRBmodelselection","title":"msalibian/FRBmodelselection","subtitle":"GitHub repository","main_image_url":"https://assets-cdn.github.com/images/email/message_cards/header.png","avatar_image_url":"https://assets-cdn.github.com/images/email/message_cards/avatar.png","action":{"name":"Open in GitHub","url":"https://github.com/msalibian/FRBmodelselection"}},"updates":{"snippets":[{"icon":"PERSON","message":"@msalibian in #1: Indeed, a quick check confirms the difference that I mention above between the residual scale M-estimator in robustbase::lmrob and the one in FRBmodelselection::my.rlm, both are asymptotically equivalent, but generally slightly different: \r\n\r\n\r\nliver \u003c- texmex::liver\r\nrbc \u003c- robustbase::lmrob.control(tuning.psi=3.443689)\r\nx \u003c- model.matrix(~ log(ALT.B), data=liver)\r\ny \u003c- log(liver$ALT.M)\r\n(tmp0 \u003c- robustbase::lmrob(log(ALT.M) ~ log(ALT.B), data=liver, control=rbc))$scale\r\n(tmp \u003c- FRBmodelselection::my.rlm(x, y, control=rbc))$scale\r\n\r\n# do the residual M-scale estimators satisfy their \"usual\" definition?\r\n\r\n# re-scaled rho() function so that sup_x rho(x) = 1\r\nrho \u003c- function(a, cc) Rho(x=a, cc=cc) * 6 / cc^2\r\n\r\n# residuals of the FRBmodelselection::my.rlm fit\r\nre \u003c- as.vector(y - x %*% tmp$coef.s )\r\nmean( rho(re / tmp$scale, cc=rbc$tuning.chi) ) # * 6 / rbc$tuning.chi^2 ) \r\n\r\n# residuals of the robustbase::lmrob fit\r\nre0 \u003c- as.vector(y - x %*% tmp0$init.S$coef )\r\n# robustbase::lmrob uses a slightly different definition for\r\n# the residual M-scale estimator, in particular:\r\nmean( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) # =! 1/2 !!\r\n# but\r\nsum( rho(re0 / tmp0$scale, cc=rbc$tuning.chi) ) / (n-2)\r\n"}],"action":{"name":"View Issue","url":"https://github.com/msalibian/FRBmodelselection/issues/1#issuecomment-431428953"}}} [ { "@context": "http://schema.org", "@type": "EmailMessage", "potentialAction": { "@type": "ViewAction", "target": "https://github.com/msalibian/FRBmodelselection/issues/1#issuecomment-431428953", "url": "https://github.com/msalibian/FRBmodelselection/issues/1#issuecomment-431428953", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { "@type": "Organization", "name": "GitHub", "url": "https://github.com" } }, { "@type": "MessageCard", "@context": "http://schema.org/extensions", "hideOriginalBody": "false", "originator": "AF6C5A86-E920-430C-9C59-A73278B5EFEB", "title": "Re: [msalibian/FRBmodelselection] my.rlm scale, RFPE use of control (#1)", "sections": [ { "text": "", "activityTitle": "Matias Salibian Barrera", "activityImage": "https://assets-cdn.github.com/images/email/message_cards/avatar.png", "activitySubtitle": "@msalibian", "facts": [ ] } ], "potentialAction": [ { "name": "Add a comment", "@type": "ActionCard", "inputs": [ { "isMultiLine": true, "@type": "TextInput", "id": "IssueComment", "isRequired": false } ], "actions": [ { "name": "Comment", "@type": "HttpPOST", "target": "https://api.github.com", "body": "{\n\"commandName\": \"IssueComment\",\n\"repositoryFullName\": \"msalibian/FRBmodelselection\",\n\"issueId\": 1,\n\"IssueComment\": \"{{IssueComment.value}}\"\n}" } ] }, { "name": "Close issue", "@type": "HttpPOST", "target": "https://api.github.com", "body": "{\n\"commandName\": \"IssueClose\",\n\"repositoryFullName\": \"msalibian/FRBmodelselection\",\n\"issueId\": 1\n}" }, { "targets": [ { "os": "default", "uri": "https://github.com/msalibian/FRBmodelselection/issues/1#issuecomment-431428953" } ], "@type": "OpenUri", "name": "View on GitHub" }, { "name": "Unsubscribe", "@type": "HttpPOST", "target": "https://api.github.com", "body": "{\n\"commandName\": \"MuteNotification\",\n\"threadId\": 398424769\n}" } ], "themeColor": "26292E" } ]