R-bloggers

Update to PSID panel builder for R: psidR

leisuretronic — Sun, 19 May 2013 15:49:00 +0000

(This article was first published on plausibel, and kindly contributed to R-bloggers)

I just pushed the most recent version of the PSID panel data builder introduced a little while ago. Got some user feedback and made some improvements. The package is hosted on github.

News:

I added a reproducible example using artificial data which you can run by calling 'example(build.panel)'. This means you can try out the package before bothering to download anything and it provides a simple test of the main function.
I've included a suggestion to use the R survey package to analyse this dataset and made it explicit in the examples how to obtain the desired weights for each wave. Note that your results are invalid in the majority of cases if you ignore the survey design (i.e. the weights).
I got some useful comments from Anthony Damico (thanks!) and integrated the SAScii package. (check out his tutorials at http://www.asdfree.com/). This allows one to download the data directly from the PSID server into R, thereby removing any dependency on Stata or SAS to preprocess the raw data. (As is common with large datasets, the raw data come in ASCII format that needs to be fixed up into rows and columns.) The downside is that downloading directly takes a rather long time: downloading FAM1985ER, FAM1986ER and the index IND2009ER took 3 and a half hours.

Hopefully I can get another round of feedback (particularly from a windows user: I could not test that all the paths are written correctly on a unix system) before submitting to CRAN.

To leave a comment for the author, please follow the link and comment on his blog: plausibel.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Sharing my R notes

tylerrinker — Sun, 19 May 2013 00:56:48 +0000

(This article was first published on TRinker's R Blog » R, and kindly contributed to R-bloggers)

I started working with R 2 1/2 years ago. I remember opening R closing it and thinking it was the dumbest thing ever (command line to a non programmer is not inviting). Now it’s my constant friend. From the beginning I took notes to remind myself all of the things I learned and relearned. They’ve been invaluable to me in learning. They are not particularly well arranged nor do they credit sources properly. There are likely bad or outdated practices in there but I figured they may be helpful to others learning the language and so I’m sharing.

Note that :

1) they are poorly arranged
2) they may have mistakes
3) they don’t credit others work properly or at all

They were for me but now I think maybe others will find them useful so here they are:

click here

*Note that the file is larger ~7000KB and 274 pages worth.

To leave a comment for the author, please follow the link and comment on his blog: TRinker's R Blog » R.

R (Web Server) Solutions – Amplifying Artichokes

Pradeep Mavuluri — Sat, 18 May 2013 16:25:00 +0000

(This article was first published on Econometrics_Help, and kindly contributed to R-bloggers)

Every month I see one or more new R based web server solutions coming into the market, sight seeing some of them thought of sharing one of my old architecture map manifested to the client long back in early 2009 (good to see quick spreading of scalable and customizable open source statistical computing tool in the market).

To leave a comment for the author, please follow the link and comment on his blog: Econometrics_Help.

What is probabilistic truth?

Corey Chivers — Sat, 18 May 2013 15:34:40 +0000

(This article was first published on bayesianbiologist » Rstats, and kindly contributed to R-bloggers)

I am currently working on a validation metric for binary prediction models. That is, models which make predictions about outcomes that can take on either of two possible states (eg Dead/not dead, heads/tails, cat in picture/no cat in picture, etc.) The most commonly used metric for this class of models is AUC, which assesses the relative error rates (false positive, false negative) across the whole range of possible decision thresholds. The result is a curve that looks something like this:

Where the area under the curve (the curve itself is the Receiver Operator Curve (ROC)) is some value between 0 and 1. The higher this value, the better your model is said to perform. The problem with this metric, as many authors have pointed out, is that a model can perform very well in terms of AUC, but be completely miscalibrated in terms of the actual probabilities placed on each outcome.

A model which distinguishes perfectly between positive and negative cases (AUC=1) by placing a probability of 0.01 on positive cases and 0.001 on negative cases may be very far off in terms of the actual probability of a positive case. For instance, positive cases may actually occur with probability 0.6 and negative cases with 0.2. In most real situations, our models will predict a whole range of different probabilities with a unique prediction for each data point, but the general idea remains. If your goal is simply to distinguish between cases, you may not care whether the probabilities are not correct. However, if your model is purporting to quantify risk then you very much want to know if you are placing the probabilistically true predictions on cases that are yet to be observed.

Which begs the question: What is probabilistic truth?

This questions appears, at least at first, to be rather simple. A frequentist definition would say that the probability is correct, or true, if the predicted probability is equal to the long run outcomes. Think of a dice rolled over and over counting the number of times a one is rolled. We would compare this frequency to our predicted probability of rolling a one (1/6 for a fair six-sided die) and would say that our predicted probability was true if this frequency matched 1/6.

But what about situations where we can’t re-run an experiment over and over again? How then would we evaluate the probabilistic truth of our predictions?

I’ll be working through this problem in a series of posts in the coming weeks. Stay tuned!

To leave a comment for the author, please follow the link and comment on his blog: bayesianbiologist » Rstats.

Recent Changes to caret

Max Kuhn — Sat, 18 May 2013 12:19:35 +0000

(This article was first published on Blog - Applied Predictive Modeling, and kindly contributed to R-bloggers)

Here is a summary of some recent changes to caret.

Feature Updates:

train was updated to utilize recent changes in the gbm package that allow for boosting with three or more classes (via the multinomial distribution)
The Yeo-Johnson power transformation was added. This is very similar to the Box-Cox transformation, but it does not require the data to be greater than zero.

New models referenced by train:

Maximum uncertainty linear discriminant analysis (Mlda) and factor-based linear discriminant analysis (RFlda) from the HiDimDA package were added.
The kknn.train model in the kknn package was added. This is basically a more intelligent K-nearest neighbors model that can use distance weighting, non-Euclidean distances (via the o Minkowski distance) and a few other features.
The extraTrees function in the package of the same name was added. This generalizes the random forest model by adding randomness to the predictors and the split values that are evaluated at each split point.

Numerous bugs were also fixed in the last few releases.

The new version is 5.16-04. Feel free to email me at mxkuhn@gmail.com if you have any feature requests or questions.

To leave a comment for the author, please follow the link and comment on his blog: Blog - Applied Predictive Modeling.

Wiekvoet 2013-05-18 02:26:00

Wingfeet — Sat, 18 May 2013 08:26:00 +0000

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

I was reading Paul Hiemsta's blogpost on Much more efficient bubble sort in R using the Rcpp and inline packages, went back to his first post Bubble sort implemented in pure R and thought, surely we can do it better in pure R. So I cleaned inner loops, removed a recursion made some variations and finally made a big improvement as I vectorized it. I don't need to try Rcpp to know that is and will be faster, but surely closed the gap a bit.

Original code

This is how it was originally in the blog in pure R. If there are changes it is layout.

larger = function(pair) {
if(pair[1] > pair[2]) return(TRUE) else return(FALSE)
}
swap_if_larger = function(pair) {
if(larger(pair)) {
return(rev(pair))
} else {
return(pair)
}
}
swap_pass = function(vec) {
for(i in seq(1, length(vec)-1)) {
vec[i:(i+1)] = swap_if_larger(vec[i:(i+1)])
}
return(vec)
}
bubble_sort = function(vec) {
new_vec = swap_pass(vec)
if(isTRUE(all.equal(vec, new_vec))) {
return(new_vec)
} else {
return(bubble_sort(new_vec))
}
}

Improve inside of loops

This is essentially the same code, but cleaned the first two functions.

larger1 = function(pair) pair[1] > pair[2]

swap_if_larger1 = function(pair) {

if(larger1(pair)) rev(pair)

else pair

}

swap_pass1 = function(vec) {

for(i in seq(1, length(vec)-1)) {

vec[i:(i+1)] = swap_if_larger1(vec[i:(i+1)])

}

return(vec)

}

bubble_sort1 = function(vec) {

new_vec = swap_pass1(vec)

if(isTRUE(all.equal(vec, new_vec))) {

return(new_vec)

} else {

return(bubble_sort1(new_vec))

}

Improve outside loops

I cleaned the outside loop, then decided that the two inside functions were not needed any more, because it was reduced to one statement.

swap_pass2 = function(vec) {

for(i in 1:(length(vec)-1)) {

if (vec[i] > vec[i+1] ) vec[c language="(i,i+1)"][/c] <- vec[c language="(i+1,i)"][/c]

}

return(vec)

}

bubble_sort2 = function(vec) {

new_vec = swap_pass2(vec)

if(identical(vec, new_vec)) new_vec

else bubble_sort2(new_vec)

}

No recursion

This was actually not intended a speed update, but R was complaining about too deep recursion. The same swapp_pass2 as before is used.

bubble_sort3 = function(vec) {

new_vec = swap_pass2(vec)

while (!identical(vec, new_vec)) {

vec <- new_vec

new_vec = swap_pass2(vec)

}

new_vec

}

A different setup

The way I understood bubblesort is different from what is programmed so far. So far the sort is completed when no improvements are made, which is checked via a vector comparison. I always understood you don't check, the first bubble goes to the end, at which point the last element is the maximum. The second bubble then stops at end-1 etc. The final bubble is only the first two elements, after which the algorithm is finished.

swap_pass4 = function(vec,iend) {

for(i in 1:iend) {

if (vec[i] > vec[i+1] ) vec[c language="(i,i+1)"][/c] <- vec[c language="(i+1,i)"][/c]

}

return(vec)

}

bubble_sort4 = function(vec) {

for (iend in (length(vec)-1):1) vec <- swap_pass4(vec,iend)

vec

}

Tuning Paul's original algorithm

I was a bit disappointed with the improvements from using the true bubble sort. Apparently there is some gain by checking if the process is completed. Which is logical since the bubbles do some sorting of the intermediate elements. So, what about bubbles that go up and go down combined with a check on no improvements?

swap_pass3b = function(vec) {

for(i in length(vec):2) {

if (vec[i] < vec[i-1] ) vec[c language="(i,i-1)"][/c] <- vec[c language="(i-1,i)"][/c]

}

return(vec)

}

bubble_sort3b = function(vec) {

new_vec = swap_pass2(vec)

while (!identical(vec, new_vec)) {

new_vec = swap_pass2(vec)

vec <- swap_pass3b(new_vec)

}

vec

}

Vectorizing the bubbles

When you look at the bubblesort without checking, it seems it is only assumed that the first step brings the highest element at the last position. This can be achieved by just pulling the highest element and placing it at the end, without the intermediate swaps.

bubble_sort5 = function(vec) {

wm <- which.max(vec)

vec <- c(vec[-wm],vec[wm])

for (iend in ((length(vec)-1):2)) {

wm <- which.max(vec[1:iend])

vec <- c(vec[1:iend][-wm],vec[1:iend][wm],vec[(iend+1):length(vec)])

}

vec

}

Results

First a demonstration that they work;

test_vec = round(runif(10, 0, 100))

cbind(

bubble_sort(test_vec),

bubble_sort1(test_vec),

bubble_sort2(test_vec),

bubble_sort3(test_vec),

bubble_sort4(test_vec),

bubble_sort3b(test_vec),

bubble_sort5(test_vec)

)

[,1] [,2] [,3] [,4] [,5] [,6] [,7]

[1,] 28 28 28 28 28 28 28

[2,] 41 41 41 41 41 41 41

[3,] 42 42 42 42 42 42 42

[4,] 46 46 46 46 46 46 46

[5,] 52 52 52 52 52 52 52

[6,] 68 68 68 68 68 68 68

[7,] 89 89 89 89 89 89 89

[8,] 93 93 93 93 93 93 93

[9,] 98 98 98 98 98 98 98

[10,] 100 100 100 100 100 100 100

Speed

Now the timing. Cleaning the outer loops saves a third of the time. The bubble as I remembered it and the up-and-down bubble improve a bit. We are now at 13% of the original time. The final version, which uses the vectorization, blows them all out of the water.

test_vec = runif(1000, 0, 100)

system.time(bubble_sort(test_vec))

# user system elapsed

# 25.48 0.02 25.69

system.time(bubble_sort1(test_vec))

# user system elapsed

# 17.88 0.00 17.91

system.time(bubble_sort2(test_vec))

# user system elapsed

# 4.75 0.00 4.75

system.time(bubble_sort3(test_vec))

# user system elapsed

# 4.75 0.00 4.74

system.time(bubble_sort4(test_vec))

# user system elapsed

# 3.45 0.00 3.44

system.time(bubble_sort3b(test_vec))

# user system elapsed

# 3.57 0.00 3.59

system.time(bubble_sort5(test_vec))

# user system elapsed

# 0.06 0.00 0.06

How about sort()?

Sort is much faster and scales much better, even more so after subtracting the time for the sample() call.. But I don't need to do the exercise to know that a not optimal algorithm in an interpreted language cannot compete with modern algorithms in a compiled language.

library(microbenchmark)

vector_size = 100

print(microbenchmark(bubble_sort5(sample(vector_size)),

sort(sample(vector_size)),

sample(vector_size)) )

Unit: microseconds

expr min lq median uq

bubble_sort5(sample(vector_size)) 1730.041 1776.957 1806.2795 1866.758

sort(sample(vector_size)) 104.096 114.359 126.4540 180.334

sample(vector_size) 5.865 8.797 9.8965 11.729

max neval

9978.521 100

233.849 100

21.992 100

To leave a comment for the author, please follow the link and comment on his blog: Wiekvoet.

Interfacing XTide and R

Luke Miller — Sat, 18 May 2013 03:02:18 +0000

(This article was first published on lukemiller.org » R-project, and kindly contributed to R-bloggers)

XTide is an open-source program that predicts tide heights and current speeds for hundreds of tide and current stations around the United States. It can be used to produce tide predictions in the past and future for a site at your chosen interval (down to the minute), as well as producing sunrise and sunset times, [...]

To leave a comment for the author, please follow the link and comment on his blog: lukemiller.org » R-project.

Unit conversion in R

Karsten W. — Fri, 17 May 2013 22:22:00 +0000

(This article was first published on factbased, and kindly contributed to R-bloggers)

Last weekend I submitted an update of my R package datamart to CRAN. It has been more than a half year since the last update, however there are only minor advances. The package is still in its early stages, and very experimental.

One new feature is the function uconv. Think iconv, but instead of converting character vectors between different encodings, this function converts numerical vectors between different units of measurements. Now if you want to know how many centimeters one horse length is, you can write in R:

> #install.packages("datamart")
> library(datamart)
> uconv(1, "horse length", "cm")

and you will get the answer 240. I had the idea for this function when I had to convert between various energy units, including natural units of energy fuels like cubic metres of natural gas. The uconv function supports this, using common constants for the conversion.

> uconv(1, "Mtoe", "PJ")
[1] 41.88
> uconv(1, "m³ NG", "kWh")
[1] 10.55556

These conversions may be ambigious. For instance, the last one combines a volume and an energy dimension. An optional parameter allows the specification of the context, or unitset:

> uconv(1, "Mtoe", "PJ", uset="Energy")

The currently available unit sets and units therein can be inspected with

> uconvlist()

The first argument can be a numerical vector:

> set.seed(13)
> uconv(37+2*rnorm(5), "°C", "°F", uset="Temperature")
[1] 100.59558  97.59102 104.99059  99.27435 102.71309

To leave a comment for the author, please follow the link and comment on his blog: factbased.

Chutes & ladders: How long is this going to take?

Karl Broman — Fri, 17 May 2013 21:05:42 +0000

(This article was first published on The stupidest thing... » R, and kindly contributed to R-bloggers)

I was playing Chutes & Ladders with my four-year-old daughter yesterday, and I thought, “How long is this going to take?”

I saw an interesting mathematical analysis of the game a few years ago, but it seems to be offline, though you can read it via the wayback machine.

But that didn’t answer my specific question, namely, “How long is this going to take?”

So I wrote a bit of R code to simulate the game.

Here’s the distribution of the number of spins to complete the game, by number of players:

With two players, the average number of spins is 52, with a 90th percentile of 88.

If you add a third player, the average increases to 65, and the 90th percentile increases to 103. You’re playing fewer rounds, but each round is three times as long. If you add a fourth player, the average is 76 and the 90th percentile is 117.

So, in trying to minimize the agony, it seems best to not encourage my eight-year-old son to join us in the game. If he plays with us, there’s a 63% chance that it will take longer.

And that’s particularly true because then the chance of my daughter winning drops from about 1/2 to about 1/3.

That raises another question: if I let her go first, what advantage does that give her? Not much. The chance that the person who goes first winning is 50.9%, 34.4%, and 25.9%, respectively, when there are 2, 3, and 4 players. So not a noticeable amount. Thus I cheat (on her behalf). Really, thought, I’m cheating in order to shorten the game as much as to ensure that she wins.

Note: There’s a close connection between this problem and my work on the multiple-strain recombinant inbred lines. (See this and that.) I’m tempted to play around with it some more.

To leave a comment for the author, please follow the link and comment on his blog: The stupidest thing... » R.

Which Torontonians Want a Casino? Survey Analysis Part 2

inkhorn82 — Fri, 17 May 2013 19:52:55 +0000

(This article was first published on Data and Analysis with R, at Work, and kindly contributed to R-bloggers)

In my last post I said that I would try to investigate the question of who actually does want a casino, and whether place of residence is a factor in where they want the casino to be built. So, here goes something:

The first line of attack in this blog post is to distinguish between people based on their responses to the third question on the survey, the one asking people to rate the importance of a long list of issues. When I looked at this list originally, I knew that I would want to reduce the dimensionality using PCA.

library(psych)
issues.pca = principal(casino[,8:23], 3, rotate="varimax",scores=TRUE)

The PCA resulted in the 3 components listed in the table below. The first component had variables loading on to it that seemed to relate to the casino being a big attraction with lots of features, so I named it “Go big or Go Home”. On the second component there seemed to be variables loading on to it that related to technical details, while the third component seemed to have variables loading on to it that dealt with social or environmental issues.

	Go Big or Go Home	Concerned with Technical Details	Concerned with Social/Environmental Issues or not	Issue/Concern
Q3_A	0.181	0.751		Design of the facility
Q3_B	0.366	0.738		Employment Opportunities
Q3_C	0.44	0.659		Entertainment and cultural activities
Q3_D	0.695	0.361		Expanded convention facilities
Q3_E		0.701	0.346	Integration with surrounding areas
Q3_F	0.808	0.266		New hotel accommodations
Q3_G	-0.117		0.885	Problem gambling & health concerns
Q3_H			0.904	Public safety and social concerns
Q3_I		0.254	0.716	Public space
Q3_J	0.864	0.218		Restaurants
Q3_K	0.877	0.157		Retail
Q3_L	0.423	0.676	-0.1	Revenue for the city
Q3_M	0.218	0.703	0.227	Support for local businesses
Q3_N	0.647	0.487	-0.221	Tourist attraction
Q3_O		0.118	0.731	Traffic concerns
Q3_P	0.497	0.536	0.124	Training and career development

Once I was satisfied that I had a decent understanding of what the PCA was telling me, I loaded the component scores into the original dataframe.

casino[,110:112] = issues.pca$scores
names(casino)[110:112] = c("GoBigorGoHome","TechnicalDetails","Soc.Env.Issues")

In order to investigate the question of who wants a casino and where, I decided to use question 6 as a dependent variable (the one asking where they would want it built, if one were to be built) and the PCA components as independent variables. This is a good question to use, because the answer options, if you remember, are “Toronto”, “Adjacent Municipality” and “Neither”. My approach was to model each response individually using logistic regression.

casino$Q6[casino$Q6 == ""] = NA
casino$Q6 = factor(casino$Q6, levels=c("Adjacent Municipality","City of Toronto","Neither"))

adj.mun = glm(casino$Q6 == "Adjacent Municipality" ~ GoBigorGoHome + TechnicalDetails + Soc.Env.Issues, data=casino, family=binomial(logit))
toronto = glm(casino$Q6 == "City of Toronto" ~ GoBigorGoHome + TechnicalDetails + Soc.Env.Issues, data=casino, family=binomial(logit))
neither = glm(casino$Q6 == "Neither" ~ GoBigorGoHome + TechnicalDetails + Soc.Env.Issues, data=casino, family=binomial(logit))

Following are the summaries of each GLM:
Toronto:

Adjacent municipality:

Neither location:

And here is a quick summary of the above GLM information:

Judging from these results, it looks like those who want a casino in Toronto don’t focus on the big social/environmental issues surrounding the casino, but do focus on the flashy and non-flashy details and benefits alike. Those who want a casino outside of Toronto do care about the social/environmental issues, don’t care as much about the flashy details, but do have a focus on some of the non-flashy details. Finally, those not wanting a casino in either location care about the social/environmental issues, but don’t care about any of the details.

Here’s where the issue of location comes into play. When I look at the summary for the GLM that predicts who wants a casino in an adjacent municipality, I get the feeling that it’s picking up people living in the down-town core who just don’t think the area can handle a casino. In other words, I think there might be a “not in my backyard!” effect.

The first inkling that this might be the case comes from an article from the Martin Prosperity Institute (MPI), who analyzed the same data set, and managed to get a very nice looking heat map-map of the responses to the first question on the survey, asking people how they feel about having a new casino in Toronto. From this map, it does look like people in Downtown Toronto are feeling pretty negative about a new casino, whereas those in the far east and west of Toronto are feeling better about it.

My next evidence comes from the cities uncovered by geocoding the responses in the data set. I decided to create a very simple indicator variable, distinguishing those for whom the “City” is Toronto, and those for whom the city is anything else. I like this better than the MPI analysis, because it looks at peoples’ attitudes towards a casino both inside and outside of Toronto (rather than towards the concept of a new Casino in Toronto). If there really is a “not in my backyard!” effect, I would expect to see evidence that those in Toronto are more disposed towards a casino in an adjacent municipality, and that those from outside of Toronto are more disposed towards a casino inside Toronto! Here we go:

As you can see here, those from the outside of Toronto are more likely to suggest building a casino in Toronto compared with those from the inside, and less likely to suggest building a casino in an adjacent municipality (with the reverse being true about those from the inside of Toronto).

That being said, when you do the comparison within city of residence (instead of across it like I just did), those from the inside of Toronto seem equally likely to suggest that the casino be built in our outside of the city, whereas those outside are much more likely to suggest building the casino inside Toronto than outside. So, depending on how you view this graph, you might only say there’s evidence for a “not in my backyard!” effect for those living outside of Toronto.

As a final note, I’ll remind you that although these analyses point to which Torontonians do want a new casino, the fact from this survey remains that about 71% of respondents are unsupportive of a casino in Toronto, and 53% don’t want a casino built in either Toronto or and adjacent municipality. I really have to wonder if they’re still going to go ahead with it!

To leave a comment for the author, please follow the link and comment on his blog: Data and Analysis with R, at Work.