Thwarted Again

The ol’ blog doesn’t allow comments and discourages search engine indexing now so basically putting things up here is somewhat pointless but I’m going to keep doing it.

I was talking about R and statistics and stuff and I ran into an interesting problem. Let’s say I have a csv file with data for a dependent, response variable W and three independent variables X, Y, and Z. Additionally, I have a second csv file with some additional values for the input variables. My goal is to produce a model for W in terms of X, Y and Z from the first data file and then see what the model predicts for the values of X, Y and Z in the second data file. If you had measured values for W in the second file then you could compare your predictions with the actual and get an idea of how good your model is. This is called out-of-sample model validation. Anyway, I load the data files like so:

> data <- read.csv("datafile.csv") > data2 <- read.csv("datafile2.csv")

However, if you don't have a handy data set you can just get R to fake one for you. Note that I've chosen different sizes here.

> data <- data.frame(W=rnorm(10), X=rnorm(10), Y=rnorm(10), Z=rnorm(10)) > data2 <- data.frame(W=rnorm(5), X=rnorm(5), Y=rnorm(5), Z=rnorm(5))

It turns out that there are two ways to make your linear model. The right way and the wrong way. The wrong way goes like this:

> lmfit <- lm(data$W ~ data$X + data$Y + data$Z) > print(lmfit$fitted.values)
1 2 3 4 5 6 7 8 9 10
-0.8562631 1.1639021 1.1048939 0.4355315 0.7493099 0.1484507 -0.8458088 -0.3278646 0.3195509 1.3771322


When you try to predict the values of W using data2 it will complain at you in a kind of unhelpful way.

> predW <- predict(lmfit, data2) Warning message:
'newdata' had 5 rows but variables found have 10 rows
> print(predW)
1 2 3 4 5 6 7 8 9 10
-0.8562631 1.1639021 1.1048939 0.4355315 0.7493099 0.1484507 -0.8458088 -0.3278646 0.3195509 1.3771322


R is being kind of unhelpful here. I found a fairly long discussion of this that ultimately helped me solve my problem but it goes around the long way to get there. http://faustusnotes.wordpress.com/2012/02/16/problems-with-out-of-sample-prediction-using-r/

The short version is that when we built our model we said we were specifically looking at "data$W", not just "W". So when we ask it to predict based on inputs from data2 it says, "I can't find data$W in here." and then it just stuffs the original predictions for 'data' into predW. That's not very helpful and possibly misleading. However, this can all be solved by giving our model variable names less scope and providing data via the data argument.

> lmfit_new <- lm(W ~ X + Y + Z, data=data) > print(lmfit_new$fitted.values)
1 2 3 4 5 6 7 8 9 10
-0.8562631 1.1639021 1.1048939 0.4355315 0.7493099 0.1484507 -0.8458088 -0.3278646 0.3195509 1.3771322
> predW_new <- predict(lmfit_new, data2) > print(predW_new)
1 2 3 4 5
-1.3766995 -1.8670918 1.1370810 2.8103150 0.6789566


All better. And then we can get the Root Mean Square of the difference between our prediction and our (completely made up) data2 W:

> sqrt(sum((data2$W - predW_new)^2))
[1] 3.655083


For a while I’ve used JMP for data analysis professionally. Recently we’ve been talking about experiment design and statistical power at work and during some of these discussions I decided it might be good if I knew how to do some R. The learning curve is OK but there’s a lot of kind of magical behavior that I don’t really get yet.

So far I’m working on getting a handle on linear models. JMP has a dialog that is Analyze->Fit Model. It then gives you a nice little browsing interface and lets you build your model with the crosses and nesting and all that. Then you tell it what kind of model you want and what kind of report you’re interested in seeing and let it go. I have a sample dataset that has nothing to do with work that has a response variable “delta”, a couple of continuous variables that effect delta and then a nominal variable “coding” which is essentially, control / experimental. I’ve fit several different models to this data but knowing what I know about how I built the dataset, it makes sense to nest the continuous variables within “coding”. The report it produces has a lot of nicely formatted data and you can run additional tests and reports from within that.

R is a much more programmy kind of tool. It’s possible to get the same model in one line of code.
> lmfit <- lm(data$delta ~ data$coding + data$dt*data$coding + data$dt2*data$coding, read.csv("derp.csv"))
But then you have a result object that you interrogate for information. I'm still struggling with the data structures R presents. The way you index into them to get at specific data really hasn't clicked yet and I find myself fumbling around in the terminal a lot. This is what it looks like when you ask it for the linear model and a summary thereof.
> lmfit
Call:
lm(formula = data$delta ~ data$coding + data$dt * data$coding +
data$dt2 * data$coding, data = read.csv("derp.csv"))

Coefficients:
(Intercept) data$codingb data$dt
2.15264 0.01773 0.75516
data$dt2 data$codingb:data$dt data$codingb:data$dt2
0.51304 1.67854 0.96067

> summary(lmfit)
Call:
lm(formula = data$delta ~ data$coding + data$dt * data$coding +
data$dt2 * data$coding, data = read.csv("derp.csv"))

Residuals:
Min 1Q Median 3Q Max
-2.25889 -0.22361 0.03567 0.39381 2.18111

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.15264 0.37886 5.682 2.23e-06 ***
data$codingb 0.01773 0.55084 0.032 0.97452
data$dt 0.75516 0.21063 3.585 0.00104 **
data$dt2 0.51304 0.20672 2.482 0.01818 *
data$codingb:data$dt 1.67854 0.29369 5.715 2.02e-06 ***
data$codingb:data$dt2 0.96067 0.29090 3.302 0.00226 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.018 on 34 degrees of freedom
Multiple R-squared: 0.9205, Adjusted R-squared: 0.9088
F-statistic: 78.75 on 5 and 34 DF, p-value: < 2.2e-16

It's not as pretty but all of the information is there. Also you can save everything as a script and rerun it if you get more data or a different dataset with the same columns. JMP also supports scripting but I honestly have never used it.

I do find R's form of the fit equation to be easier to understand but either way you have to make sense of the idea that "if experimental use these coefficients, if control use these others":
R:

(Intercept) 2.15264
data$codingb 0.01773
data$dt 0.75516
data$dt2 0.51304
data$codingb:data$dt 1.67854
data$codingb:data$dt2 0.96067

vs
JMP:

5.476 + Match( :coding, "a", (:dt - 1.3) * 0.755, "b", (:dt - 1.3) * 2.434)
+ Match( :coding, "a", (:dt2 - 1.25) * 0.513, "b", (:dt2 - 1.25) * 1.474)
+ Match( :coding, "a", -1.700, "b", 1.700, . )


reddit take us such nice places.