Save R objects, and other stuff

Yesterday, Christopher asked me how to store an R object, in order to save some time, when working on the project.

First, download the csv file for searches related to some keyword, via http://www.google.com/trends/, for instance “sunglasses“. Recall that csv files store tabular data (numbers and text) in plain-text form, with comma-separated values (where csv term comes from). Even if it is not a single and well-defined format, it is a text file, not an excel file!

Recherche sur le Web : intérêt pour sunglasses
Dans tous les pays; De 2004 à ce jour

Intérêt dans le temps
Semaine,sunglasses
2004-01-04 - 2004-01-10,48
2004-01-11 - 2004-01-17,47
2004-01-18 - 2004-01-24,51
2004-01-25 - 2004-01-31,50
2004-02-01 - 2004-02-07,52

(etc) The file can be downloaded from the blog,

> report=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/glasses.csv",
+ skip=4,header=TRUE,sep=",",nrows=464)

Then, we have run a function on this data frame, to transform it. It can be found from a source file

> source(“http://freakonometrics.blog.free.fr/public/code/H2M.R”)

In this source file, there is function that transforms a weekly series into a monthly one. The output is either a time series, or a numeric vector,

> sunglasses=H2M(report,lang="FR",type="ts")

Here, we asked for a time series,

> sunglasses
          Jan      Feb      Mar      Apr
2004 49.00000 54.27586 66.38710 80.10000
2005 48.45161 58.25000 69.93548 80.06667
2006 49.70968 57.21429 67.41935 82.10000
2007 47.32258 55.92857 70.87097 84.36667
2008 47.19355 54.20690 64.03226 79.36667
2009 45.16129 50.75000 63.58065 76.90000
2010 32.67742 44.35714 58.19355 70.00000
2011 44.38710 49.75000 59.16129 71.60000
2012 43.64516 48.75862 64.06452 70.13333
          May      Jun      Jul      Aug
2004 83.77419 89.10000 84.67742 73.51613
2005 83.06452 91.36667 89.16129 76.32258
2006 86.00000 92.90000 93.00000 72.29032
2007 86.83871 88.63333 84.61290 72.93548
2008 80.70968 80.30000 78.29032 64.58065
2009 77.93548 70.40000 62.22581 51.58065
2010 71.06452 73.66667 76.90323 61.77419
2011 74.00000 79.66667 79.12903 66.29032
2012 79.74194 82.90000 79.96774 69.80645
          Sep      Oct      Nov      Dec
2004 56.20000 46.25806 44.63333 53.96774
2005 56.53333 47.54839 47.60000 54.38710
2006 51.23333 46.70968 45.23333 54.22581
2007 56.33333 46.38710 44.40000 51.12903
2008 51.50000 44.61290 40.93333 47.74194
2009 37.90000 30.38710 28.43333 31.67742
2010 50.16667 46.54839 42.36667 45.90323
2011 52.23333 45.32258 42.60000 47.35484
2012 54.03333 46.09677 43.45833

that we can plot using

> plot(sunglasses)

Now we would like to store this time series. This is done easily using

> save(sunglasses,file="sunglasses.RData")

Next time we open R, we just have to use

> load("/Users/UQAM/sunglasses.RData")

to load the time series in R memory. So saving objects is not difficult.

Last be not least, for the part on seasonal models, we will be using some functions from an old package. Unfortunately, on the CRAN website, we see that

but nicely, files can still be found on some archive page. On Linux, one can easily install the package using (in R)

> install.packages(
+ "/Users/UQAM/uroot_1.4-1.tar.gz",
+ type="source")

With a Mac, it is slightly more complicated (see e.g. Jon’s blog): one has to open a Terminal and to type

R CMD INSTALL /Users/UQAM/uroot_1.4-1.tar.gz

On Windows, it is possible to install a package from a zipped file: one has to download the file from archive page, and then to spot it from R.

The package is now installed, we just have to load it to play with it, and use functions it contains to tests for cycles and seasonal behavior,

> library(uroot)
> CH.test
function (wts, frec = NULL, f0 = 1, DetTr = FALSE, ltrunc = NULL)
{
s <- frequency(wts)
t0 <- start(wts)
N <- length(wts)
if (class(frec) == "NULL")
frec <- rep(1, s/2)
if (class(ltrunc) == "NULL")
ltrunc <- round(s * (N/100)^0.25)
R1 <- SeasDummy(wts, "trg")
VFEalg <- SeasDummy(wts, "alg")

(etc)

The D.I.Y. strategy, and why academics should blog !

Last week, I went to the Econometrics seminar of Montréal, at UdM, where Alfred Galichon was giving a great talk on marriage market. Alfred is a former colleague (from France), a co-author, an amazing researcher, and above all, a friend of mine. And he has always be supportive about my blogging activities. So while we were having lunch, after the seminar, Alfred mentioned my blogging activity to the other researchers. I should say researchers in Econometrics (yes, with a capital E, since it is a Science, as mention in an old paper by David Hendry by the end of the 70’s). Usually, when I am involved in this kind of meeting, I start with some apologies, explaining that I do like theoretical econometrics (if not, I would not come to the seminar), but I do like my freakonometrics activity. I do like to use econometrics (or statistical techniques) to figure out (at least to try) why some things works the way they do. I try to find data, and then try to briefly analyze them to answer some simple questions. Or sometime, I just run simulations to answer more theoretical questions (or at least to give clues).

But above all, I like the fact that blogging gives me the opportunity to interact with people I would never meet without this activity. For instance, last May, I was discussing (on Twitter) with @coulmont, @joelgombin and @imparibus about elections in France. Then @coulmont asked me “yes, everyone knows that there should be some ecological fallacies behind my interpretation, but I am not so sure since I have data with a small granularity. As an econometrician, what do you think ?” Usually, I hate having a label, like “… I ask you since you’re a mathematician“, or “as an economist, what do you think of…“. Usually, when people ask me economic questions, I just claim being a mathematician, and vice-versa. But here, I even put on the front of my blog the word “econometrics” (more or less). So here, I could not escape… And the truth is, that while I was a student, I never heard anything about this “ecological fallacy“. Neither did I as a researcher (even if I have been reading hundreds of econometric articles, theoretical and applied). Neither did I as a professor (even if I have been teaching econometrics for almost ten years, and I have read dozens textbooks to write notes and handouts). How comes ? How come researchers in sociology and in political sciences know things in econometrics that I have never heard about ?

The main reason – from my understanding – is the following: if everyone talks about “interdisciplinarity” no one (perhaps a few) is really willing to pay the price of working on different (not to say many) areas. I tried, and trust me, I found it difficult. It is difficult to publish a paper in a climate journal when you’re not specialist in climate (and you just want to give your opinion as a statistician). It is difficult to assume that you might waste weeks (not to say months) reading articles in geophysics if you want to know more about earthquakes risks, going to seminars, etc. Research is clearly a club (“club as defined in Buchanan (1965)) story.

This week, I planned to go to some journal club in biology and physics, at McGill (kindly, a colleague there invited me, but we got a time misunderstanding)… this has nothing to do with my teaching, nor with my research activities. But I might learn something ! Yes, I do claim that I am paid just to have fun, to read stuff that I do find interesting, trying to understand the details of a proof, trying to understand how data were obtained. In most cases, it might (and should) be a complete waste of time, since I will not publish anything (anything serious, published in some peer reviewed journal) on that topic… but should I really care ? As I explained earlier (in French), I do also claim that I have a moral obligation to return everything I have seen, heard, read. And since I am not a big fan of lectures (and that I do not think I have skills for that) I cannot give my opinion, neither on economics facts (as @adelaigue or @obouba might do on their blogs) or on science results (as @tomroud does). But I think I can help on modeling and computational issues. My point being: never trust what you read (even on my blog) but please, try to do it yourself! You read that “90% of French executive think about expatriation” (as mentioned here)? Then try to find some data that should confront that statement. And see if you come up with the same conclusion… And since it might be a bit technical sometimes, here are some lines of code, to do it on your own… Academics have a legitimacy when they give their opinions on technical issues. At least they can provide with a list of reference everyone should read to get an overview of the topic. But academics can also help people read graphs, or data. To give them “numeracy” (or a culture in numbers) necessary to understand complex issues.

To conclude, I should mention that I understood what this “ecological fallacy” was from  Thomsen (1987) and many more documents could be found on Soren Thomsen’s page http://www.mit.ps.au.dk/srt/. But I got most of the information I was looking for from a great statistician, who happens to be also an amazing blogger: Andrew Gelman (see http://andrewgelman.com/). I will probably write a post someday about this, since I found the question extremely interesting, and important.

Actuariat 1, ACT2121, septième cours

Toujours dans le cadre de la préparation à l’examen P de la SOA, une série d’exercices. Comme il reste trois semaines (en plus de l’examen final), on va essayer de finir de revoir l’ensemble des notions. Les 50 exercices sont en ligne ici. Je mettrais en ligne (très bientôt) l’énoncé et la correction de l’examen intra de lundi (avec – comme la dernière fois – les statistiques de réponse par question)

Introduction aux processus SARIMA

Quelques transparents en plus, qui devraient correspondre aux deux prochains cours de séries temporelles, sur les processus autorégressifs (AR) et moyennes mobiles (MA), les ARMA, les ARIMA (intégrés) et les SARIMA (saisonniers). J’ai mis des notes sur les tests de racine unités, je rajouterais quelques transparents la semaine prochaine sur les tests de saisonnalité, et quelques exemples pratiques de prévision. Les transparents sont en ligne ici,

Actuariat en Afrique subsaharienne francophone

Depuis que le blog (ou les précédents) existe, je sais (par les messages que je reçois par courriel) qu’il est beaucoup lu en Afrique (francophone). J’ai reçu avant hier un livre publié par Aymric Kamega et Frédéric Planchet, tiré de la thèse de doctorat d’Aymric (en ligne sur http://halshs.archives-ouvertes.fr/) défendue en décembre dernier. J’étais alors rapporteur extérieur, et j’avais souligné l’intérêt des travaux d’Aymric dans le cadre de la volonté de la CIMA (autorité de contrôle régionale des marchés d’assurance pour l’Afrique subsaharienne francophone) de fournir aux assureurs de la région des outils adaptés. En particulier des tables de mortalité d’expérience, propres à la région. En effet, (comme le rappelle Aymric), , suite aux états généraux (de l’assurance vie)  en 2007, le principe de la construction de tables de mortalité d’expérience a été adopté en remplacement des tables de mortalité de la population générale française entre 1960 et 1964 (tables dites PM 60-64 et PF 60-64) jusqu’alors imposées. Aymric avait travaillé sur ce sujet, et a soutenu une thèse de doctorat qui avait fait l’unanimité. J’avais pris beaucoup de plaisir à lire la thèse, j’en pense que j’en aurais à lire le livre (disponible sur le site de l’éditeur http://seddita.com/).

Introduction aux séries temporelles

Nous allons commencer demain la modélisation des séries temporelles. Les transparents de la séance sont en ligne ici. Je rappelle que les notes de cours (complètes) sont également en ligne, . Je continuerai à poster régulièrement des billets contenant des commandes R.

Au programme cette semaine, l’utilisation des méthodes de régression pour extraire tendance et cycle, le lissage exponentielle (simple, double et saisonnier), et une présentation des notions importantes dans le cours (stationnarité, autocorrélations, bruit blanc, etc).

In less than 50 days: 2013, year of statistics

A couple of days ago, Jeffrey sent me an e-mail, encouraging me to write about the international year of statistics, hoping that I would participate to the event. The goals of Statistics2013 are (as far as I understood) to increase public awareness of the power and impact of statistics on all aspects of society, nuture Statistics as a profession and promote creativity and development in Statistics.

So I will try to write more posts (yes, I surely can try) on statistical methodology, with sexy applications… So see you in a few days on this blog, to celebrate statistics ! I can even invite guests to write here… all contributors are welcome !

 

Le mois des (boites à) moustaches

On célèbre cette année les 35 ans d’un livre qui a révolutionné la statistique, tant sur les aspects computationnels que graphiques: Exploratory Data Analysis par John Tukey. La légende prétend (on pourra relire a brief history of S par Richard Becker) que ce livre a inspiré ce qui allait devenir R,
Commercial software did not fit well into our research environment. It often used a ‘‘shotgun’’ approach — print out everything that might be relevant to the problem at hand because it could be several hours before the statistician could get another set of output. This was reasonable when computing was done in a batch mode. However, we wanted to be able to interact with our data, using Exploratory Data Analysis (Tukey, 1971) techniques. In addition, commercial statistical software usually didn’t compute what we wanted and was not set up to be modified.

En particulier, c’est dans ce livre qu’est introduit un outil graphique intéressant pour visualiser des données, le boxplot, appelé en français boite à moustaches (que tous les statisticiens se doivent d’utiliser en ce mois de Movembre). Dans la version originale, les boites sont verticales

mais mis horizontalement, on voit une boite, avec des moustaches. Cela dit, les moustaches ne sont pas toujours ce que l’on croit. Par exemple, on pense souvent (et moi le premier, il y a encore peu) que la boite à moustaches permettait de visualiser des valeurs extrêmes, voire des quantiles s’il y a trop de données (le maximum peut – en effet – être très très grand)

C’est malheureusement plus compliqué que ca… Afin de mieux comprendre ce qu’on trace, regardons le jeu de données suivant,

> set.seed(1)
> X=rlnorm(99)


(on est ici en échelle logarithmique). Pour commencer, on peut faire une boite à moustache avec une visualisation du minimum (à gauche) et du maximum (à droite)

> boxplot(X,horizontal=TRUE,log="x",range="0")


Et effectivement, on retrouve facilement tous les paramètres: ceux de la boite (en rouge)

> abline(v=median(X),lty=2,col="red")
> abline(v=quantile(X,.25),lty=2,col="red")
> abline(v=quantile(X,.75),lty=2,col="red")

et ceux de la moustache (en bleu)

> abline(v=min(X),lty=2,col="blue")
> abline(v=max(X),lty=2,col="blue")

Mais ce n’est pas la version standard de la boite à moustaches, qui est – sous R – la suivante

> boxplot(X,horizontal=TRUE,log="x")


Si la borne inférieure de la moustache est toujours le minimum (en tous cas ici), pour la partie supérieure… c’est plus compliqué. Malheureusement ce n’est pas un quantile standard (ceux à 90% et 95% sont représentés ci-dessus en mauve). Pour comprendre ce que c’est, il faut regarder le code: ici, on rajoute au quartile supérieur (la partie de droite de la boite) un pourcentage de la distance inter-quartile (la longueur de la boite), et ce pourcentage est le suivant

> M=boxplot(X)$stats
> (M[5]-M[4])/(M[4]-M[2])
[1] 1.350628
> 2*qnorm(.75)
[1] 1.34898

Plusieurs sites ou ouvrages suggèrent d’utiliser 1.5 fois la longueur interquartile, mais sous R, on utilise . Il faut le savoir (je l’ai appris un peu par hasard en faisant une formation sur la régression quantile au R Montréal Group au printemps passé). Et effectivement, on retrouve les bornes de la moustache,

> abline(v=quantile(X,.75)+2*qnorm(3/4)*IQR(X),lty=2,col="blue")
> abline(v=min(X),lty=2,col="blue")

et je laisse les quantiles pour montrer que (malheureusement) ce n’est pas ce qui est utilisé ici,

> abline(v=quantile(X,.9),lty=2,col="purple")
> abline(v=quantile(X,.95),lty=2,col="purple")

Cela dit, comme je l’évoquais lors de la formation, beaucoup de monde laisse planer une réelle ambiguité sur ce qui est représenté par ces boites à moustaches… Et le fait que ce ne soit pas un quantile me gene un peu… Par exemple, ici on avait la boite à moustaches suivante

Si je prends les 25 plus petites observations, et que je les divise par 4, on obtient

La partie de gauche bouge (c’est normal, c’est la partie qu’on a modifiée), et la médiane reste identique, ainsi que la partie droite de la boite (ce qui est normal, ce sont des quantiles au delà de 25%). Mais de manière un peu gênante, la borne supérieure de la moustache s’est ralongée…  C’est normal, car on utilise une distance inter-quartile. Mais on n’a pas touché aux grandes observations, c’est donc gênant de la voir bouger ainsi… non ?

sur les transformations dans un modèle linéaire

Je voulais prendre 5 minutes pour reprendre une question posée par courriel, qui me permettra de poursuivre sur des choses évoquées mercredi dernier en cours. Je vais reformuler la question, mais en gros, cela disait: la méthode de Box-Cox (évoquée en début de semaine, ici) avait pour objet de choisir entre un modèle linéaire (que l’on étudié dans tous les sens depuis le début du cours) et un modèle log-linéaire (évoqué ici, par exemple). L’idée était d’associer le cas linéaire à  (dans la transformée de Box-Cox) et à  pour le cas multiplicatif (ou log-linéaire). Mais que se passe-t-il si la valeur optimale rejette ces deux cas, et est proche (disons) de  ?
Considérons le cas suivant (exemple que je ne me lasse d’utiliser)

> reglm=lm(dist~speed,data=cars)
> library(MASS)
> boxcox(reglm)

La valeur optimale est effectivement proche de . En posant , on aurait envie de considérer un modèle de la forme

(comme le suggère la transformation de Box-Cox). On peut faire la régression,

> regsqrt=lm(sqrt(dist)~speed,data=cars)
+ summary(regsqrt)

Call:
lm(formula = sqrt(dist) ~ speed, data = cars)

Residuals:
Min      1Q  Median      3Q     Max
-2.0684 -0.6983 -0.1799  0.5909  3.1534

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.27705    0.48444   2.636   0.0113 *
speed        0.32241    0.02978  10.825 1.77e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.102 on 48 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.7034
F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14

et effectivement, le résultat est concluant. Mais on ne peut pas en rester là, comme pour pour le passage au logarithme, notre variable d’intéret reste .
La prédiction sur le modèle transformé était

et

Comme avec le logarithme, on pourrait prendre comme prédiction pour  mais

Bref, on a un estimateur biaisé. Et là encore, comme on est sur des variables positives, l’inégalité de Jensen doit meme nous garantir que l’une des quantités domine toujours l’autre (je laisse les amateurs de convexité l’écrire dans le bon sens). Comment faire ?
La première solution est de noter que

Donc de manière très simple, on a notre prédiction (en prenant le carré, comme intuité, mais en rajoutant une terme – positif – lié à la variance).
Une seconde consiste à noter que Z (conditionnellement à la variable explicative) était gaussienne. Donc en prenant le carré (moyennant quelques changement d’échelle), on devrait tomber sur une loi du chi-deux qui est une loi que l’on connait bien (sinon je peux renvoyer ici).
Essayons de creuser un peu ces idées (en particulier, pour dériver des intervalles de confiance). Si on visualise la prédiction de ce modèle sur la racine carrée de la distance de freinage, on obtient

> plot(speed,sqrt(dist))
> x=seq(0,30,by=.2)
> distsqrtp=predict(regsqrt,newdata=
+ data.frame(speed=x),interval="prediction")
> polygon(c(x,rev(x)),c(distsqrtp[,2],
+ rev(distsqrtp[,3])),col="yellow",border=NA)
> lines(x,distsqrtp[,1],lwd=2,col="red")

Pour passer d’une loi normale à une loi du chi-deux, il faut retrancher la moyenne (pour centrer la variable) et diviser par l’écart-type (pour avoir une variance unitaire),

> s=summary(regsqrt)$sigma
> mu=predict(regsqrt)
> distsqrtp01=(sqrt(dist)-mu)/s

On a ainsi une variable qui, conditionnellement à la variable explicative est supposé suivre une loi ,

> plot(speed,distsqrtp01,ylim=c(-3,3))
> abline(h=qnorm(c(.025,.975),0,1),lty=2,col="red")

On prend ensuite le carré pour obtenir notre loi , on on trace la bande de confiance (il faut que la valeur à la deux soit faible),

> plot(speed,distsqrtp01^2)
> abline(h=qchisq(.95,df=1),lty=2,col="red")

On a maitenant nos intervalles de confiance, en utilisant cette loi du chi-deux, en inversant notre transformation,

> mu=predict(regsqrt,newdata=data.frame(speed=x))
> distsup=(mu+s*sqrt(qchisq(.95,df=1)))^2
> distinf=(mu-s*sqrt(qchisq(.95,df=1)))^2

que l’on peut visualiser simplement (sur les données de base, avec en ordonnée la distance)

> plot(cars)
> lines(x,distsup,lty=2,col="red")
> lines(x,distinf,lty=2,col="red")

Et la valeur prédite ? On va utiliser notre relation liant la variance, l’espérance du carré, et le carré de l’espérance,

> distesp=mu^2+s^2
> lines(x,distesp,lwd=2,col="red")

Nice, n’est-ce pas ?
Sinon, comme je le disais en cours, si la transformation optimale sur  est de prendre , peut-etre pourrait-on envisager de prendre  (de manière un peu duale) comme variable explicative. Encore une fois, c’est une idée, car rien ne le garantit. En effet

En particulier, le second modèle est un modèle linéaire classique, avec de l’homoscédasticité: on aura une dispersion uniforme autour de notre parabole. En revanche, avec le premier modèle, on a une espèce d’hétéroscédasticité (comme tenu du double produit), avec une variance du terme d’erreur qui croit avec  (car les coefficients sont positifs). Bref, en terme d’intervalles de confiance, on devrait avoir des choses assez différentes.
Regardons la régression de la distance (cette fois) sur la vitesse, et le carré de la vitesse,

> reglm=lm(dist~speed+I(speed^2),data=cars)
> distp=predict(reglm,newdata=
+ data.frame(speed=x),interval="prediction")

Si on visualise la prédiction, avec un intervalle de confiance, on obtient

> plot(cars)
> polygon(c(x,rev(x)),c(distp[,2],
+ rev(distp[,3])),col="yellow",border=NA)
> lines(x,distp[,1],lwd=2,col="red")

Autrement dit, on suppose que notre modèle est homoscédastique. Les régions de confiance entre les deux approches sont clairement différentes… meme si les prédictions sont presque superposées (oui oui, il y a deux courbes, une rouge et une bleue).

On Box-Cox transform in regression models

A few days ago, a former student of mine, David, came back to me about Box-Cox tests in linear models. It made me look more carefully at the test, and I do not understand what is computed, to be honest. Let us start with something simple, like a linear simple regression, i.e.

https://latex.codecogs.com/gif.latex?Y_i=\beta_0+\beta_1%20X_i+\varepsilon_i

Let us introduced – as suggested in Box & Cox (1964) – the following family of (power) transformations

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}%20=%20\begin{cases}%20\dfrac{Y_i^\lambda-1}{\lambda}%20&\text{%20%20}%20(\lambda%20\neq%200)\\[8pt]%20\log{(Y_i)}%20%20&\text{%20}%20(\lambda%20=%200)%20\end{cases}

on the variable of interest. Then assume that

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}=\beta_0+\beta_1%20X_i+\varepsilon_i

As mentioned in Chapter 14 of Davidson & MacKinnon (1993) – in French – the log-likelihood of this model (assuming that observations are independent, with distribution https://latex.codecogs.com/gif.latex?\varepsilon_i\sim\mathcal{N}(0,\sigma^2)) can be written

https://latex.codecogs.com/gif.latex?\log%20\mathcal{L}=-\frac{n}{2}\log(2\pi)-n\log\sigma%20\\-\frac{1}{2\sigma^2}\sum_{i=1}^n\left[Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)\right]^2+(\lambda-1)\sum_{i=1}^n\log%20Y_i

We can then use profile-likelihood techniques (see here) to derive the optimal transformation.

This can be done in R extremely simply,

> library(MASS)
> boxcox(lm(dist~speed,data=cars),lambda=seq(0,1,by=.1))

we then get the following graph,

If we look at the code of the function, it is based on the QR decomposition of the https://latex.codecogs.com/gif.latex?\boldsymbol{X}matrix (since we assume that https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a full-rank matrix). More precisely, https://latex.codecogs.com/gif.latex?\boldsymbol{X}=QR where https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a https://latex.codecogs.com/gif.latex?n\times%202 matrix, https://latex.codecogs.com/gif.latex?Q is a https://latex.codecogs.com/gif.latex?n\times%202 orthonornal matrix, andhttps://latex.codecogs.com/gif.latex?R is a https://latex.codecogs.com/gif.latex?2\times2 upper triangle matrix. It might be convenient to use this matrix since, for instance, https://latex.codecogs.com/gif.latex?R\widehat{\boldsymbol{\beta}}=Q%27Y.  Thus, we do have an upper triangle system of equations.

> X=lm(dist~speed,data=cars)$qr

The code used to get the previous graph is (more or less) the following,

> g=function(x,lambda){
+ y=NA
+ if(lambda!=0){y=(x^lambda-1)/lambda}
+ if(lambda==0){y=log(x)}
+ return(y)} 
> n=nrow(cars)
> X=lm(dist~speed,data=cars)$qr
> Y=cars$dist
> logv=function(lambda){
+ -n/2*log(sum(qr.resid(X, g(Y,lambda)/
+ exp(mean(log(Y)))^(lambda-1))^2))}
> L=seq(0,1,by=.05)
> LV=Vectorize(logv)(L)
> points(L,LV,pch=19,cex=.85,col="red")

As we can see (with those red dots) we can reproduce the R graph. But it might not be consistent with other techniques (and functions described above). For instance, we can plot the profile likelihood function, https://latex.codecogs.com/gif.latex?\lambda\mapsto\log\mathcal{L}

> logv=function(lambda){
+ s=summary(lm(g(dist,lambda)~speed,
+ data=cars))$sigma
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ -n/2*log(2 * pi)-n*log(s)-.5/s^2*(sum(e^2))+
+ (lambda-1)*sum(log(Y))
+ }
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> (maxf=optimize(logv,0:1,maximum=TRUE))
$maximum
[1] 0.430591

$objective
[1] -197.6966

> abline(v=maxf$maximum,lty=2)

The good point is that the optimal value of https://latex.codecogs.com/gif.latex?\lambda is the same as the one we got before. The only problem is that the https://latex.codecogs.com/gif.latex?y-axis has a different scale. And using profile likelihood techniques to derive a confidence interval will give us different results (with a larger confidence interval than the one given by the standard function),

> ic=maxf$objective-qchisq(.95,1)
> #install.packages("rootSolve")
> library(rootSolve)
> f=function(x)(logv(x)-ic)
> (lower=uniroot(f, c(0,maxf$maximum))$root)
[1] 0.1383507
> (upper=uniroot(f, c(maxf$maximum,1))$root)
[1] 0.780573
> segments(lower,ic,upper,ic,lwd=2,col="red")

Actually, it possible to rewrite the log-likelihood as

https://latex.codecogs.com/gif.latex?\mathcal{L}=\star-\frac{n}{2}\log\left[\sum_{i=1}^n\left(\frac{Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)}{\dot{Y}^\lambda}\right)^2\right]

(let us just get rid of the constant), where

https://latex.codecogs.com/gif.latex?\dot{Y}=\exp\left[\frac{1}{n}\sum_{i=1}^n%20\log%20Y_i\right]

Here, it becomes

> logv=function(lambda){
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ elY=(exp(mean(log(Y))))
+ -n/2*log(sum((e/elY^lambda)^2))
+ }
>
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> optimize(logv,0:1,maximum=TRUE)
$maximum
[1] 0.430591

$objective
[1] -47.73436

with again the same optimal value for https://latex.codecogs.com/gif.latex?\lambda, and the same confidence interval, since the function is the same, up to some additive constant.

So we have been able to derive the optimal transformation according to Box-Cox transformation, but so far, the confidence interval is not the same (it might come from the fact that here we substituted an estimator to the unknown parameter https://latex.codecogs.com/gif.latex?\sigma.

Somewhere else, part 19

Interesting discussions, this week, on statistics, and confidence,

And alot of posts on Sandy (and consequences in the NYC area)

And a usual, a lot of very interesting posts here and there

avec toujours des billets en Français

Did I miss something ?