More on linear regression... take care of initial point for robust regression

 It is better to supply a initial point for robust regression !


x <- c(0.428571428571429, 0.2, 0.3, 0, 0, 0.2, 0, 0, 0.1, 0, 0.1, 

  0.1, 0, 0, 0, 0.3, 0.2, 0.3, 0.2, 0.2, 0, 0, 0, 0.2, 0.222222222222222, 

  0.1, 0, 0.4, 0.3, 0.5, 0, 0.4, 0.5, 0.8, 0.3, 0.1, 0.2, 0, 0.1, 

  0, 0.1, 0.1, 0.4, 0, 0, 0, 0, 0, 0, 0.333333333333333, 0.444444444444444, 

  0.2, 0.222222222222222, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

  0, 0, 0, 0, 0.4, 0, 0, 0, 0, 0.111111111111111, 0, 0, 0, 0, 0, 

  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

  0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.2, 0, 0, 0, 0, 0, 0, 0)

y <- c(0.369999388816617, 0.152260672623962, 0.392636473518975, 0.107742543910461, 

       0.105802942749025, 0.147888875829182, 0.0180958542177892, 0.00376865991773073, 

       0.00852273837059186, 0.00857162539273257, 0.112051915806605, 

       0.0755798935035902, 0.0268356635188166, 0.00419682981147497, 

       0.00525227149049588, 0.120478282490271, 0.433740156299308, 0.504343565678989, 

       0.0206028265235165, 0.151280541814684, 0.0103854597483139, 0.0361491733598324, 

       0.0472461159232045, 0.617197541010412, 0.341487842962594, 0.301148642757091, 

       0.196289575076544, 0.581773949785228, 0.231922322336444, 0.0794048724091231, 

       0.00537870061961155, 0.098836566232596, 0.91476603065251, 0.616117348565266, 

       0.630629550337848, 0.0299893115473008, 0.144611328935455, 0.00316949810200154, 

       0.00364313171075452, 0.00841171463044064, 0.029740302359379, 

       0.0035511361061388, 0.00584692618954925, 0.00100829236005141, 

       0.000763929818703641, 0.000489764007946297, 0.0420444966611161, 

       0.0377279236240397, 0.0755041088036021, 0.0702674671437556, 0.490342044264684, 

       0.0126472157965609, 0.0198354059018282, 0.00302651586313457, 

       0.00280412047780873, 0.00206798894186777, 0.000689349234298032, 

       0.0118392137557611, 0.00294598087896767, 0.00360192399011589, 

       0.000177536212304921, 0.000372443888258384, 0.0360944855490864, 

       0.124415609631037, 0.00935459188019089, 0.00378842359914763, 

       0.00685095396765276, 0.0103847747243049, 0.00134801747284256, 

       0.000591917764291339, 0.0212501726144608, 0.0115559930482827, 

       0.0339710160832701, 0.0173148151553399, 0.00433684377042257, 

       0.00331646545220465, 0.000129380544070732, 0.00168666576557622, 

       0.00202293274563766, 0.000127617791462482, 0.000126412903174323, 

       0.0050245051731696, 0.0029260069948355, 0.000546063655789477, 

       9.47601490327725e-05, 0.00526659632136166, 0.000108655976769835, 

       0.000443878503125248, 0.000170564355810093, 0.00140919430563582, 

       0.000899982324830479, 0.000637780267140836, 0.00112316932596192, 

       0.000405537765870411, 0.00360093692558572, 0.003761521337125, 

       0.0125607769366116, 0.00172106229606284, 0.000912577233480949, 

       0.00523099962479612, 0.000536960099118963, 0.0165559987979115, 

       0.000306510798257877, 0.00229079171256135, 0.000184489822114836, 

       0.000924333005338229, 0.00404542140403508, 0.00140001440983028, 

       0.00799258757598758, 0.000705978887791876, 0.00560785783318743, 

       0.00275163168909473, 0.0140781508814297, 0.00138060195094072, 

       0.00499557472754309, 0.0160496188963424, 0.00986865781418728, 

       0.000623581108784438, 0.00547952864148045, 0.000973862265134229, 

       0.00678300356505618, 0.00456741578551256, 0.0563840117464262, 

       0.12167652045372, 0.0133864818580382, 0.0130352885016593, 0.00052780369438233, 

       0.00024670139960277, 0.000635307193894255, 0.000319077056548938, 

       0.000887158541859849, 0.00556904105860274, 0.000197728927761402, 

       0.773494870178593, 0.0232722370275186, 0.000265508467810429, 

       0.000411718635448104, 0.00218648841943213, 0.000671109445425718, 

       0.00553495122778798, 0.00921147407566352, 0.000298669678588454)


plot(x, y)


l <- lm(y ~ x)

abline(l, col= "green", lty=1, lwd=2) 

text(x = 0.8, y=0.8, labels = "Ordinary least-square", col="green")



library("robustbase")

# Default version... very wrong

lr <- lmrob(y ~ x, setting="KS2014")

abline(lr, col= "red", lty=2, lwd=2) 

text(x = 0.8, y=0.1, labels = "'Robust'", col="red")


# Using the lm default: it solves the problem

lr <- lmrob(y ~ x, setting="KS2014", init = list(coefficients=l$coefficients, scale = 1))

abline(lr, col= "red", lty=2, lwd=1) 

text(x = 0.8, y=0.5, labels = "Robust", col="red")


Commentaires

Enregistrer un commentaire

Posts les plus consultés de ce blog

Standard error from Hessian Matrix... what can be done when problem occurs

An Hessian matrix is a square matrix of partial second order derivative. From the square-root of the inverse of the diagonal, it is possible to estimate the standard error of parameters. Let take an example: val=rnorm(100, mean=20, sd=5) # Return -ln L of values in val in Gaussian distribution with mean and sd in par fitnorm<-function(par, val) {   -sum(dnorm(val, par["mean"], par["sd"], log = TRUE)) } # Initial values for search p<-c(mean=20, sd=5) # fit the model result <- optim(par=p, fn=fitnorm, val=val, method="BFGS", hessian=TRUE) # Inverse the hessian matrix to get SE for each parameters mathessian <- result$hessian inversemathessian <- solve(mathessian) res <- sqrt(diag(inversemathessian)) # results data.frame(Mean=result$par, SE=res)   library(MASS) (r<-fitdistr(val, "normal")) It works well. However, the inverse of the matrix cannot be calculated if the second order derivative for some paramete
Lire la suite

stepAIC from package MASS with AICc

Due to a bug in dropterm() and addterm() in package MASS, it is impossible to use AICc with stepAIC() . Here is a solution. Enter these commands and now you can use: stepAIC(g0, criteria="AICc") or stepAIC(g0, criteria="BIC") For example: datax <- data.frame(y=rnorm(100), x1=rnorm(100),                     x2=rnorm(100), x3=rnorm(100), x4=rnorm(100)) g0 <- glm(y ~ x1+x2+x3+x4, data=datax) gx_AICc <- stepAIC(g0, criteria="AICc") gx_AIC <- stepAIC(g0, criteria="AIC") gx_AIC <- stepAIC(g0, criteria="BIC")
Lire la suite

Install treemix in ubuntu 20.04

Informations are available here  https://bitbucket.org/nygcresearch/treemix/downloads/ In terminal, enter for  Ubuntu 20.04: wget https://bitbucket.org/nygcresearch/treemix/downloads/treemix-1.13.tar.gz tar -xvf treemix-1.13.tar.gz sudo apt-get install libgsl-dev sudo apt install libboost-dev sudo apt install libboost-all-dev cd treemix-1.13/ ./configure make sudo make install And it works: $ treemix TreeMix v. 1.13 $Revision: 231 $ Options: -h display this help -i [file name] input file -o [stem] output stem (will be [stem].treeout.gz, [stem].cov.gz, [stem].modelcov.gz) -k [int] number of SNPs per block for estimation of covariance matrix (1) -global Do a round of global rearrangements after adding all populations -tf [file name] Read the tree topology from a file, rather than estimating it -m [int] number of migration edges to add (0) -root [string] comma-delimited list of populations to set on one side of the root (for migration) -g [vertices file name] [edges file name] read the g
Lire la suite