Light and Momentary

The post you've been waiting for

Have you been fretting, oh my friends, about the difficulties of making a nice-looking figure in R? Fret no more, because I am here to help.

I think I felt so crummy yesterday because I was coming down with a cold, or maybe a cold-plus. Today I am sniffly and achy and wobbly, but I am optimistic about getting my groove back once my immune system is less overloaded. I'm working on that R&R I told you about, and my revisions table has about 30 items in it. Earlier this week I spent a looooong time on just one of them: making a figure suitable for reproduction in black and white.

One of the things people rave about with R is its graphics. And indeed, you don't have to do very much googling to see beautiful R graphics. Just between you and me and the internet? They are a big pain. The closest I came to throwing over the whole PhD idea was when I had to fix a figure in R. It would have been pretty dumb to quit at that point, since I had successfully defended my dissertation, but I certainly thought about it.

So. You have a figure that awaits your attention? Here's what you do. First put your numbers where R can find them, like so:

> bored.stats <- c(3.2, 1.8)
> bored.nostats <- c(2.8, 4.7)
> sobored.stats <- c(4.8, 1.2)
> unbored.nostats <-c(1.2, 3.8)

For this figure we are going to imagine that your reaction to this post varies, depending on how bored you were when you started reading it and whether or not you have some stats background. Let's suppose that you are rating your before-reading boredom and your after-reading boredom on a 5-point Likert scale, where 1 is SO EXCITED and 5 is CRUSHED BY ENNUI.

If you are planning to save your finished figure, then you need to tell R about your plans at the outset. If you just want to mess around until you get it right, you can skip this step: 

> png(filename = "bored.png", height = 300, width = 225, bg = "white")

Now you can start plotting. Lay down the first set of points, without the stupid default labels, and then add your scale and your labels for each axis. Like so:

> plot(bored.stats, ylim = c(0,5), axes = F, ann = F)
> axis(1, at = 1:2, lab = c("Before reading", "After reading")
> axis (2, at = 0:5)
> lines(bored.stats, type = "o")

Hey, presto, you have a start. Now let's add those other lines, with different lines and point markers instead of different colors.

> lines(bored.nostats, type = "o", pch = 22, lty = 2) # plain English: square points, dashed line
> lines(sobored.stats, type = "o", pch = 23, lty = 3) # diamond points, dotted line
> lines(unbored.nostats, type = "o", pch = 24, lty = 4) # triangle points, dashed/dotted line

All of this I figured out with only a modicum of facepalming. The thing that tripped me up was making a legend. Just for you, my friends*, I am sharing the legend knowledge I ferreted out today.

*okay, okay, it's actually for me. Guilty as charged.

Once you tell R what to label your items, like this--

> legend.labels<-c("Bored, stats background", "Bored, no stats background", "Very bored, stats background", "Not bored, no stats background")

--then the following line of code will give you a legend:

> legend(1.15, 1, legend.labels, cex = 0.4, pch = 21:24, lty = 1:4) 

Translation: at the point (1.15, 1), use my legend labels to make a small legend (that's the cex = 0.4 part) with the specified sequence of point types and line types.

If you did that optional line of code above so R would save the file, then you need one more line of code at the end:

> dev.off()

If you're not sure where R put your file, just use getwd() to find the directory in which it lurks.

Tada! If it looks like the dog's breakfast, do it on a Mac instead of a Windows box. There's a way to deal with anti-aliasing issues in Windows, but I...am pleased enough that I figured out how to make a legend.

Wrangler

Remember those posts I used to write about R, my stats program? This is a post about R. I am going to see if I can make this interesting, but you'll have to let me know if I succeed.

I am doing a presentation in June that involves a bigger sample than most of my recent work -- about 700 kids. As the database gets bigger, it's easier to get lost in it or to miss goofy things that would be more obvious in a smaller sample.

Weird fact about me: I love playing in data. I decided to get a PhD because I thought I would love teaching. And I do -- I had a blast in my class this morning, and I think they did too. (I mean, I know it's not really about having a blast. If I can cover the material and keep them laughing, I think that's a win-win.) I didn't know how much I would enjoy the research side of things.

I'm looking at a novel measure of language development that I devised for my dissertation. Like most measures of language development, it's sensitive to age. I wanted to see what was left if I pulled out the effects of age in my sample. In other words:

> newdata$stres <- rstandard(lm(new.measure ~ age, data = newdata, na.action = na.exclude))

That last little bit cost me some head-scratching: it tells R to plug in an "NA" for anybody who has missing data, so I can match up the right kid with his very own standardized residual. (Don't you hate it when somebody gives you the wrong standardized residual?)

Next up, I wanted to get each family's data on one line -- both twins, one row. I probably shouldn't feel like a genius that I figured this out, but here's the truth: if I had a roll of aluminum foil in my office I would make a little medal that said "CONQUEROR." Here's the trick: if you separate the firstborn twins--

> newdata.odd <- subset(newdata, ID %% 2 > 0)

--from the second-born twins--

> newdata.even <-subset(newdata, ID %% 2 == 0)

--then you can merge them back together by family ID! TADA!!

> newdata.merge <-merge(newdata.odd, newdata.even, by = "familyID", all = T)

In the SPSS version of this syntax that I was combing through, you have to rename all the variables yourself. But R knows better than to give you multiple columns of variables with the same name. You'll have newdata.merge$stres.x for co-twin A, and newdata.merge$stres.y for co-twin B. And if you want to estimate genetic vs. environmental influences on your shiny new measure, you're most of the way there.

I bet you can hardly stand the excitement! Preliminary results suggest that shared environment is more important for this conversational language measure than for other related measures. Extra! Extra! Read all about in two years when I get it written up and wangled through peer review!

If you are still reading, I appreciate it! (And if you know more stats than I do and you spotted something goofy, let me know.) If you skimmed down to this paragraph in hopes that it might contain something of actual interest, no hard feelings. Regularly scheduled programming will resume soon.

Who knew a post about statistics software could contain this many exclamation points?

Ask the English Major

Dear Jamie, is it wise to back up my data at frequent intervals? How should I do that? --Prudence

Dear Prudence, it is very, very prudent to back up your data regularly. Want to know how I know this? Because I accidentally deleted more than 15% of the year 1 data when I got cocky with a merge on Wednesday. Whoops! Let's hear it for backups. To save all your stuff, just do this:

> save.image(file="directory_name/prudent_very_prudent.RData")

Then, dear Prudence, you can come out and play. --Jamie

Dear Jamie, my environment is getting cluttered. How can I clean it out efficiently? --Marla C.

Dear FlyLady, when did you become an R user? It leaves me in purple puddles sometimes too. To do a 27-fling boogie on your environment, use ls() to show you the numbers associated with your objects. Then you can do some mighty mighty rm-ing with those object numbers, like so:

> rm(list=ls()[c(1:35,38:129,131:350)])

That example would delete everything but the objects numbered 36, 37, and 130 or anything beyond 350. (HT)  Might want to back up first (see Prudence's question above). While you're on the line, can I ask you to stop sending me Rubba Scrubba ads? --Jamie

Dear Jamie, I am fascinated by the phenomenon of brilliant redheads. I have data on hair color and data on IQ. How can I find out which redheads in my sample have genius-level IQs? --Tiziano Vecelli

Dear TV, you loved them for their brains, huh? Here's how to get the tally, assuming the women have the same ID numbers in both datasets:

> genius.redheads<-intersect(subset(data$ID, data$hair.color=="red"), subset(other.data$ID, other.data$IQ>150))

The object genius.redheads will list the ID numbers for the redheads whose IQs are above the cutoff you specify. When you find some, send them my way so I don't delete my data again. Thanks and I hope it's not too moldy for your computer in there! --Jamie

Dear Jamie, please help! I am using someone else's dataset and he saved a continuous variable as a factor variable. Whatever shall I do? --Stanford Binet

Dear Stan, fret not. Try this:

> as.numeric(levels(IQ))[IQ]

I still havent figured out how to parse it but it works. (HT) --Jamie

Really Simple But Still Useful Things A Person Can Do In R Even If She Was An English Major

(This was an attempt at a Christmas picture from last weekend. How big is baby Stella? Soooo big. Not big enough to lie down with her brothers and stay put, but soooo big.)

Oh, Google, sometimes you are stupid. You are sending people to my post called "More R Tricks" even though the tricks therein are not remotely tricksy to someone googling "R tricks." I should amend it so the first line says "Keep on googling, because you are in the wrong place. Unless you are an English major."

I do know that these posts have been useful to at least a couple of people, so I am going to keep writing them from time to time. One of the things that made me scratch my head when I was first contemplating my big messy dataset was the need to merge in variables from a couple of even larger and messier datasets. It turns out to be easy-peasy. You can find helpful merging directions here, but this is what I do.

First, I make a new object with just the ID variable and the variables of interest. I am working with a subset of a giant sample (600 kids, 1500 variables), and I could wind up with a boatload of superfluous data. So:

> attach(big.sample) [this lets me work with variable names directly, instead of typing big.sample$ before each one]

> subsample <- data.frame(ID1, IQ, digit.span, zygosity)

Now I have the variables of interest in a smaller data frame called "subsample." I could also do it by column number, if I wanted, as discussed in the post linked above, with a command like this:

>subsample <- big.sample[,c(1, 600, 602, 14)] [this would give me all the rows in those four specified columns]

Next I'm going to merge those three new variables into my smaller database. I'm going to tell R to match by participant ID and only include the IDs that appear in both datasets, like this:

> year1 <- merge(year1, subsample, by = "ID1", all = F)

If you specify "all = T," as in the examples I linked to above, you will wind up with everything from both datasets. Since the command as written tells R to overwrite your old dataset with the new one, this could induce a fair amount of teeth-gnashing. You might prefer to create a transitional object while you are getting the hang of it:

> year1.for.nervous.nellies <- merge(year1, subsample, by = "ID1", all = F)

If your identifying variables have different names in their respective datasets, there is no need to fret. Just do it like this:

> year1 <- merge(year1, subsample, by.x = "ID", by.y = "ID1", all = F) [by.x is the variable name in the first dataset, by.y its name in the second]

Now you can check and see if it worked:

> names(year1)

The last three items should be IQ, digit.span, and zygosity. Now suppose you want some frequency counts for your snazzy new data. Do the boys and the girls shake out equally as identical and fraternal twins? All you need to do is type this--

> table(year1$sex, year1$zygosity)

--and R will tell you all about it.

More R Tricks

Do you hate these posts, O readers who subscribed to this blog for the tales of life with small children and not for the tips on deriving z-scores? I find these posts very useful personally and I do hope that unhappy googlers trying to learn R will also appreciate them, but they are probably a snore for the rest of you.

Yesterday I figured out part of something that had perplexed me. I would slap tables of variables together on a regular basis, and sometimes they would have column names and sometimes they would not. Part of the answer is that dataframes have column names; matrices don't. You can't do much (any?) statistical analysis on a matrix. But if you have three vectors of equal length named x, y, and z, you can say this--

> more.data<-data.frame(x,y,z)

--and get yourself a dataframe called "more.data," with columns headed x, y, and z.

Suppose you'd like to compare kids' performance on measures x, y, and z, but their scales are very different. You can derive z-scores as easy as pie:

> z.data<-scale(more.data)

Now you have a new dataframe called z.data, in which your measures have been converted to z-scores.

Suppose, though, that you don't want everything converted to z-scores -- suppose you have a column for participant ID#, and a column for group membership. Obviously, you don't want those to be z-scores. Participant #0.549879 is less than illuminating if you assigned him ID #62. You could convert the variables individually and then stick them into a new dataframe, but you can also do it more easily.

It took me a long time to figure out something about R that's actually pretty straightforward. In web pages about manipulating data, I kept seeing sets of brackets with mysterious commas in them. For a long time I just followed instructions blindly, but I finally saw what was going on. You can tell R what to do with your data in terms of rows and columns inside square brackets -- rows first, columns second. So this--

> part.of.my.data<-my.data[c(2:5),c(1,3,5,7,9)]

--tells R to give you a new array of data containing rows 2-5 and the odd-numbered columns of the old dataset. Or you could do this--

> part.of.my.data<-my.data[-c(16), ]

--which means "give me all the rows except 16 [that's the "-c(16)"], and give me every column [that's what you're saying with the ", ]". which is perhaps a little confusing]."

So if you want z-scores for your experimental measures x, y, and z (in columns 3-5) but not for the first two columns (participant ID and group membership), all you have to do is this:

> z.data<-cbind(more.data[ ,c(1,2)], scale(more.data[ ,c(3:5)]))

If you wanted to add a column with the mean of each participant's three z-scores, you could do it like this:

> z.data.with.means<-cbind(z.data, rowMeans(z.data[ ,c(3:5)]))

Extra points if you get all your brackets and parentheses lined up on the first try. You probably already know that if you say "rowmeans" instead of "rowMeans," R will blink blankly at you, waving its little "Klingon only" placard. You might not have considered that you must restrict the row means to the last three columns or else R will average in the participant ID numbers, leading you to wonder how your last participant had a mean z-score of 73. (That's like getting a 7800 on your SAT Verbal -- yes, the one with a maximum score of 800.)

Or maybe you have considered it and I am just dim.

Here's the requisite cute baby picture. A few days ago I said, "Darling baby mine, please stop trying to grab my socks in progress to see if you can tickle your cochlea with a size 1 knitting needle. Here, look at this nice fat book instead." Whereupon she flung the nice fat book to the ground and climbed upon it, the better to reach my knitting. She's a determined girl, that baby Stella.

In Which I Am Feeling Clever

Ha! I just discovered that even an English major can write an R function.

Here is a picture of all five of my lovely children (the oldest is wading, wearing the orange shirt) for those of you who are thinking, "Is she talking about her stats program AGAIN?" But perhaps you might also like to see the pretty picture I just made.

I am looking longitudinally at kids' language abilities, based on transcripts collected at three visits, each a year apart. In the first year, they were in first or second grade; by the last visit, they were in third or fourth grade. In addition to looking at group differences at each point, I am also looking at growth trajectories.

So-- I wanted to superimpose density curves. What did they look like in year 1 vs. year 2 vs. year 3? How were they changing? I figured out how to make these superimposed curves show up on my own computer screen, but I couldn't figure out how to save them in a usable format. I figured my committee would frown on a dissertation that said "Stop by my house and I will show you the pretty pictures I can't get into this document."

If you would like to make your very own superimposed density curves, it's actually not hard. Open a text editor and type the stuff in black.

progress<-function(x,y,z) [Hey, R, make me a function called "progress" where the expected input is three vectors.]
{
    plot(density(x),col=1)
    lines(density(y),col=2)
    lines(density(z),col=3)
}

Then save your function as "progress.R" and at the R command line type

> source("directory_name/progress.R")

Now we're all set to create a save-able image:

> png("image_name.png")
> progress(year1,year2,year3)
> dev.off()

And tada! This image shows the changes in the density of elaborated noun phrases across the three home visits. How much more often are kids using rich descriptive phrases as they mature? (More details here if you are pining away for want of knowledge about elaborated noun phrases.)

It's not perfect -- I'm going to change the title and the labels for the axes (which I know how to do) and ideally the scaling (which I have no idea how to do, but which would make it easier to compare curves across samples -- the numbers change each time due to attrition). But hey-- I wrote an R function and it did what I meant it to do! Hurray!

I Still Don't Speak Klingon

This morning I opened a new Google document that I called "Useful R Tips," but then it occurred to me that I might as well blog them. There must be other people out there in the world who are attempting to use R for stats and are equally baffled by the documentation. I thought I might have been exaggerating in this post, but it's really bad.

 Here is a cute baby picture for the rest of you, and a promise to be more interesting next time.

If you have a big monster spreadsheet and you want to look at descriptive stats for some of your variables, here's what you can do. Use cbind() to stick the variables of interest together into a table, like this:

> densities2<-cbind(adv2,mor2,lf2,ccnj2,ml2)

This gives me a mini-table, called densities2, with five columns containing all my data for just the five variables I want to know more about. If you want to stick together rows and not columns, use rbind().

You can also manipulate your variables while you are sticking them into that table. I am looking at things like conjunction density, which means all my values are less than 1. I told R to multiply them by 100 so I wasn't dealing with a million tiny numbers:

> densities2<-cbind(adv2,mor2,lf2,ccnj2,ml2)*100

You can use the pastecs package to make a handy table of descriptive statistics. There is information here in English, not Klingon, on how to install the package and use the stat.desc() function.

Now suppose you want to copy your table into a word processor, or just save it to read later. Use the write.table() function, like this:

> write.table(object.name.from.stat.desc.goes.here, "directory/pick_your_filename.csv", sep=",",row.names=TRUE)

Now when you look in the directory you specified, you will see pick_your_filename.csv, which you can open up in Excel or its Open Office clone to manipulate as you please.

I Hate Personal Growth

Two years ago I was taking a class in regression analysis. During the week that my second son turned eight and I needed to produce Halloween costumes for four children, I also had a big stats assignment due. I was doing this assignment with a freebie student version of SAS when the computer on which it lived died a sudden and yet harrowing death.

A Couple of Tips for Newbie Users of R

(Note to the usual crowd: sorry for the boring post. Halloween pix soon.)

In looking for helpful documentation on R, I was surprised by the unwelcoming tone of many of the resources ostensibly intended to help new R users. On one mailing list, for instance, I found a higher incidence of RTFM and STFW than I've ever seen before. Believe me, I wanted to RTFM. Unfortunately, the FM seems to be written in --I don't know, Klingon. Some language that uses an English sound system but is utterly unintelligible to me, some language spoken by a clever but apparently cruel race of beings who scorn the idea of intelligent life outside their demesnes.

To anyone who arrives here via Google, similarly battered by attempts to learn R, I suggest that you start with this friendly tutorial. If you want to copy a graphic from R to another program (like the Word document in which you are chronicling your maiden voyage), do it like this:

> png("file_name_here.png")
> plot(x1,x2) [or hist, or whatever else]
> dev.off()

I hope this saves you the many hours it took me to figure that one out. Oh, bummer, I was going to credit the web page where I finally found the magic sequence but I closed that tab and I have no hope of finding it in the moment. I am being instructed to hop off the computer so that Pete can watch Thomas the Tank Engine videos on YouTube.

Before I do: swimmermom asked me about my results at Freerice.com and I have gone back and forth about whether to reply because was it boasting to put that on the web? And maybe it is, but maybe an admission of weakness (see: the rest of this post) cancels out a boast. So I will say that I can get to level 50 pretty easily: the last time I was there, I missed 2 out of the first 47 words I was given. Then I stopped, because I wasn't sure about the word on the screen, and I have this compulsion about finishing at level 50, AND I have a statistics assignment that is (have I mentioned?) due tomorrow. So in sum: English vocabulary strong, Klingon vocabulary dismal.