Thanks for your help. It will take until next week to look into this
and get back.
Joh
On Fri, Nov 28, 2008 at 11:34 AM, Matthias Kohl
<[EMAIL PROTECTED]> wrote:
> Dear Joh,
>
> in case of our packages ROptEst and RobLox we recently implemented the
> functions roptest and roblox, which enable the computation of (in our sense)
> optimally robust estimators via one single function.
>
> In case of our regression packages we do not have such functions so far.
> However, we plan to implement these functions in the near future also using
> the formula interface of R. For the time being you have to do the necessary
> steps by hand.
>
> After the installation of these packages, you should find a subfolder
> "scripts" in the package folders of ROptEst, ROptRegTS and RobRex, as well
> as a folder "tests" in case of package RobLox. These folders contain several
> examples how to use these packages.
>
> Let us for instance consider a regression of the form y ~ x.
>
> Some points you should think about are:
> What error distribution do you assume? If normal, go for package RobRex. If
> something different, go for package ROptRegTS.
> May x and y both be contaminated or only y?
>
> Two examples using RobRex ...
>
> ###############################################################################
> ## Example 1
> ###############################################################################
> require(MASS)
> require(robustbase)
> library(ISwR)
> data(thuesen)
> attach(thuesen)
>
> ## LS-estimator
> fit.LS <- lm(short.velocity ~ blood.glucose)
> ## M-estimator
> fit.M <- rlm(short.velocity ~ blood.glucose)
> ## MM-estimator
> fit.MM <- lmrob(short.velocity ~ blood.glucose)
> ## LTS-estimator
> fit.LTS <- ltsReg(short.velocity ~ blood.glucose)
>
> ## AL-estimators
> ## regressor distribution: design measure
> require(RobRex)
> K <- DiscreteMVDistribution(cbind(1,blood.glucose))
> ## ALc-estimator: conditional neighborhood; i.e., only y is contaminated
> ## about 10 sec. on my system
> system.time(IC1 <- rgsOptIC.ALc(r = 0.5, K = K, theta = fit.M$coeff, scale =
> fit.M$s))
> ALc1 <- oneStepEstimator(cbind(1,as.matrix(thuesen)), IC1, c(fit.M$coeff,
> fit.M$s))
> ## AL-estimator: unconditional neighborhood; i.e., x and y are contaminated
> ## about 85 sec. on my system
> system.time(IC2 <- rgsOptIC.AL(r = 0.5, K = K, theta = fit.MM$coeff, scale =
> fit.MM$s))
> AL1 <- oneStepEstimator(cbind(1,as.matrix(thuesen)), IC2, c(fit.MM$coeff,
> fit.MM$s))
>
> ## Plot
> require(RColorBrewer)
> myCol <- brewer.pal(6, "Set1")
> plot(short.velocity ~ blood.glucose, ylab = "fasting blood glucose
> [mmol/l]",
> xlab = "mean circumferential shortening velocity [%/s]",
> main = "Ventricular shortening velocity", pch = 20)
> abline(fit.LS, lwd = 2, col = myCol[1])
> abline(fit.M, lwd = 2, col = myCol[2])
> abline(fit.MM, lwd = 2, col = myCol[3])
> abline(fit.LTS, lwd = 2, col = myCol[4])
> lines(c(1, c(blood.glucose,25)), ALc1[1] + ALc1[2]*c(1,c(blood.glucose,25)),
> col = myCol[5], lwd = 2)
> lines(c(1, c(blood.glucose,25)), AL1[1] + AL1[2]*c(1,c(blood.glucose,25)),
> col = myCol[6], lwd = 2)
> legend("topleft", legend = c("LS", "M", "MM", "LTS", "ALc", "AL"),
> fill = myCol, ncol = 2)
> detach(thuesen)
>
> ###############################################################################
> ## Example 2
> ###############################################################################
> data(phones)
> attach(phones)
>
> ## LS estimator
> fit2.LS <- lm(calls ~ year)
> ## M estimator
> fit2.M <- rlm(calls ~ year, maxit = 50)
> ## MM estimator
> fit2.MM <- lmrob(calls ~ year)
> ## LTS estimator
> fit2.LTS <- ltsReg(calls ~ year)
>
> ## AL estimators
> ## regressor distribution: design measure
> K <- DiscreteMVDistribution(cbind(1,year))
> ## ALc estimator: only y contaminated
> system.time(IC1 <- rgsOptIC.ALc(r = 0.5, K = K, theta = fit2.M$coeff, scale
> = fit2.M$s))
> ## takes about 9 sec. on my system
> ALc2 <- oneStepEstimator(cbind(1,year,calls), IC1, c(fit2.M$coeff,
> fit2.M$s))
> ## AL estimator: x and y contaminated
> system.time(IC2 <- rgsOptIC.AL(r = 0.5, K = K, theta = fit2.MM$coeff, scale
> = fit2.MM$s))
> ## takes about 80 sec. on my system
> AL2 <- oneStepEstimator(cbind(1,year,calls), IC2, c(fit2.MM$coeff,
> fit2.MM$s))
>
> ## Plot
> plot(calls ~ year, ylab = "phone calls [Mio.]", xlab = "year",
> main = "Belgium Phone Calls 1950-1973", pch = 20)
> abline(fit2.LS, lwd = 2, col = myCol[1])
> abline(fit2.M, lwd = 2, col = myCol[2])
> abline(fit2.MM, lwd = 2, col = myCol[3])
> abline(fit2.LTS, lwd = 2, col = myCol[4])
> lines(c(1, c(year,75)), ALc2[1] + ALc2[2]*c(1,c(year,75)), col = myCol[5],
> lwd = 2)
> lines(c(1, c(year,75)), AL2[1] + AL2[2]*c(1,c(year,75)), col = myCol[6], lwd
> = 2)
> legend("topleft", legend = c("LS", "M", "MM", "LTS", "ALc", "AL"),
> fill = myCol, ncol = 2)
>
> #########################################################
>
> Like in case of RobLox there will be changes to RobRex in the near future
> which will clearly speed up these computations ...
>
> By the way, I made some further tests ...
> Using ROptEst version 0.5 (from CRAN archive) and the recent versions of
> distr (2.0.3), distrEx (2.0.3) and RandVar (0.6.6), the recent versions
> (0.6.1, available from RForge) of ROptRegTS and RobRex install and check
> without problems on my system (R version 2.8.0 Patched (2008-11-26 r47022),
> i686-pc-linux-gnu). I'm just about to submit the new versions of our
> distr-family as well as the new versions of RandVar, RobAStBase, ROptEst and
> RobLox to CRAN. Because of the incompatibility of ROptRegTS and RobRex with
> the new implementation one has to download the current versions of ROptRegTS
> and RobRex from R-Forge.
>
> Best,
> Matthias
>
> Johannes Graumann wrote:
>>
>> Yes, fresh R session (see below). Whether it's working? I was looking for
>> a function straight forward to use - along the lines of "roblox" - but
>> failed. Can you give a hint on how to actually set up something mimicking
>> library(robustbase)
>> (regression <- lmrob(
>> y~x
>> ))
>>
>> I have had very good results using "roblox" in the past and was interested
>> in trying the regression based on the same methodology, but being a user and
>> no specialist the package has me intimidated.
>>
>> Joh
>>
>>
>>
>>>
>>> sessionInfo()
>>>
>>
>> R version 2.8.0 (2008-10-20) x86_64-pc-linux-gnu
>> locale:
>> LC_CTYPE=en_US.UTF-8; <CUT>
>>
>> attached base packages:
>> [1] stats4 tcltk tools stats graphics grDevices utils
>> [8] datasets methods base
>> other attached packages:
>> [1] ROptRegTS_0.6.0 distrMod_2.0.2 MASS_7.2-44
>> [4] ROptEst_0.5.0 RandVar_0.6.5 distrEx_2.0.2
>> [7] evd_2.2-4 distr_2.0.2 SweaveListingUtils_0.1
>> [10] sfsmisc_1.0-6 startupmsg_0.5.2 robustbase_0.4-3
>> [13] RColorBrewer_1.0-2 geneplotter_1.20.0 annotate_1.20.1
>> [16] xtable_1.5-4 lattice_0.17-17 MaxQuantUtils_1.17
>> [19] tkrplot_0.0-18 TeachingDemos_2.3 plotrix_2.5
>> [22] gsubfn_0.3-7 proto_0.3-8 AnnotationDbi_1.3.12
>> [25] RSQLite_0.7-1 DBI_0.2-4 Biobase_2.1.7
>> [28] rkward_0.4.9
>> loaded via a namespace (and not attached):
>> [1] grid_2.8.0 KernSmooth_2.22-22 RobAStBase_0.1.2
>> On Wednesday 26 November 2008 19:32:47 Matthias Kohl wrote:
>>
>>>
>>> Dear Joh,
>>>
>>> was this a fresh R session?
>>>
>>> I just removed version 0.6 of package ROptEst and let the rest of the
>>> packages unchanged. Then, installed version 0.5 of ROptEst and after
>>> that installed version 0.5 of package ROptRegTS (from CRAN). Despite
>>> some warnings and notes during installation, all seems to work fine ...
>>>
>>> Does ROptRegTS work on your system?
>>>
>>> Best,
>>> Matthias
>>>
>>> Johannes Graumann wrote:
>>>
>>>>
>>>> ROptEst 0.5! Not ROptRegTS! Sorry.
>>>>
>>>> Joh
>>>>
>>>> Johannes Graumann wrote:
>>>>
>>>>>
>>>>> Hello Matthias,
>>>>>
>>>>> ROptRegTS 0.5 installs fine on top of the most recent rest of the
>>>>> bundle
>>>>> from R-Forge. When loading it I get the warnings below. Whether the
>>>>> functionality is impeded will be figured out later today.
>>>>>
>>>>> Thanks for your help!
>>>>>
>>>>> Joh
>>>>>
>>>>>
>>>>>>
>>>>>> warnings()
>>>>>>
>>>>>
>>>>> Warning messages:
>>>>> 1: In fun(...) ... : no valid postscript previewer found; consider
>>>>> setting options("eps_view"= "....") yourself
>>>>> 2: Multiple methods tables found for 'optIC'
>>>>> 3: Multiple methods tables found for 'neighbor<-'
>>>>> 4: Multiple methods tables found for 'CallL2Fam<-'
>>>>> 5: Multiple methods tables found for 'generateIC'
>>>>> 6: Multiple methods tables found for 'checkIC'
>>>>> 7: Multiple methods tables found for 'clip'
>>>>> 8: Multiple methods tables found for 'clip<-'
>>>>> 9: Multiple methods tables found for 'cent'
>>>>> 10: Multiple methods tables found for 'cent<-'
>>>>> 11: Multiple methods tables found for 'stand'
>>>>> 12: Multiple methods tables found for 'stand<-'
>>>>> 13: Multiple methods tables found for 'lowerCase'
>>>>> 14: Multiple methods tables found for 'lowerCase<-'
>>>>> 15: Multiple methods tables found for 'neighborRadius'
>>>>> 16: Multiple methods tables found for 'neighborRadius<-'
>>>>> 17: Multiple methods tables found for 'clipLo'
>>>>> 18: Multiple methods tables found for 'clipLo<-'
>>>>> 19: Multiple methods tables found for 'clipUp'
>>>>> 20: Multiple methods tables found for 'clipUp<-'
>>>>> 21: Multiple methods tables found for 'L2deriv'
>>>>> 22: Multiple methods tables found for 'FisherInfo'
>>>>> 23: Multiple methods tables found for 'checkL2deriv'
>>>>>
>>>>> Matthias Kohl wrote:
>>>>>
>>>>>>
>>>>>> Dear Joh,
>>>>>>
>>>>>> we are currently restructuring and extending our packages.
>>>>>>
>>>>>> Some parts of the package ROptEst (version 0.5) are now included in
>>>>>> the
>>>>>> new packages distrMod and RobAStBase. Unfortunatelly, the structure of
>>>>>> the packages ROptRegTS and RobRex has not been updated yet. Hence, if
>>>>>> you want to use package ROptRegTS you have to use the old version
>>>>>> (i.e., version 0.5) of ROptEst from the archive
>>>>>> (http://cran.at.r-project.org/src/contrib/Archive/ROptEst).
>>>>>>
>>>>>> These old versions of the packages ROptEst and ROptRegTS should work
>>>>>> with R 2.8.0 and the recent versions of the packages distr, distrEx
>>>>>> and
>>>>>> RandVar. I will check this later today and will then give you more
>>>>>> details on this.
>>>>>>
>>>>>> Best
>>>>>> Matthias
>>>>>>
>>>>>> Johannes Graumann wrote:
>>>>>>
>>>>>>>
>>>>>>> Hello,
>>>>>>>
>>>>>>> Trying to say "library(ROptRegTS)", I get the error below. Can anyone
>>>>>>> help me pinpoint what's wrong?
>>>>>>>
>>>>>>> Thanks, Joh
>>>>>>>
>>>>>>> Error in loadNamespace(package, c(which.lib.loc, lib.loc),
>>>>>>> keep.source
>>>>>>> = keep.source) :
>>>>>>> in 'ROptRegTS' methods for export not found: neighbor, bound, Curve,
>>>>>>> CallL2Fam
>>>>>>> In addition: Warning message:
>>>>>>> In fun(...) : no valid postscript previewer found; consider setting
>>>>>>> options("eps_view"= "....") yourself
>>>>>>> Error: package/namespace load failed for 'ROptRegTS'
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-SIG-Robust@r-project.org mailing list
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-robust
>>>>>>>
>>>>
>>>> _______________________________________________
>>>> R-SIG-Robust@r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-robust
>>>>
>>
>>
>
> --
> Dr. Matthias Kohl
> www.stamats.de
>
>
_______________________________________________
R-SIG-Robust@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-robust
