1

Background: I am running a little A/B test, with 2x2 factors (foreground's black and background's white, off-color vs normal color), and Analytics reports the number of hits for each of the 4 conditions and at what rate they 'converted' (a binary variable, which I define as spending at least 40 seconds on page). It's easy enough to do a little editing and get in a nice R dataframe:

rates <- read.csv(stdin(),header=TRUE)
Black,White,N,Rate
TRUE,FALSE,512,0.2344
FALSE,TRUE,529,0.2098
TRUE,TRUE,495,0.1919
FALSE,FALSE,510,0.1882

Naturally, I'd like to look at a logistic regression on something like Rate ~ Black * White but R's glm wants a dataframe of 2046 rows each reporting a TRUE or FALSE conversion value & the values of Black and White. This... is a little more tricky. I googled around and checked SO but while I found some clunky code on how to convert a table of contingency counts to a dataframe, I didn't find anything about percentages/rates.

After a lot of trouble, I came up with a loop over the 4 conditions in which I repeat a dataframe rate * n times with the relevant condition values and the result True and then do the same thing but for (1 - rate) * n and the result False, and then stitch together all 8 dataframes into one giant dataframe:

ground <- NULL
for (i in 1:nrow(rates)) {
        x <- rates[i,]
        y <- do.call("rbind", replicate((x$N * x$Rate),     data.frame(Black=c(x$Black),White=c(x$White),Conversion=c(TRUE)),  simplify = FALSE))
        z <- do.call("rbind", replicate((x$N * (1-x$Rate)), data.frame(Black=c(x$Black),White=c(x$White),Conversion=c(FALSE)), simplify = FALSE))
        ground <- rbind(ground,y,z)
}

The resulting dataframe ground looks right:

sum(rates$N)
[1] 2046
nrow(ground)
[1] 2042
# the missing 4 are probably from the rounding-off of the reported conversion rate
summary(ground); head(ground, n=20)
   Black           White         Conversion     
 Mode :logical   Mode :logical   Mode :logical  
 FALSE:1037      FALSE:1020      FALSE:1623     
 TRUE :1005      TRUE :1022      TRUE :419      
 NA's :0         NA's :0         NA's :0        
   Black White Conversion
1   TRUE FALSE       TRUE
2   TRUE FALSE       TRUE
3   TRUE FALSE       TRUE
4   TRUE FALSE       TRUE
5   TRUE FALSE       TRUE
6   TRUE FALSE       TRUE
7   TRUE FALSE       TRUE
8   TRUE FALSE       TRUE
9   TRUE FALSE       TRUE
10  TRUE FALSE       TRUE
11  TRUE FALSE       TRUE
12  TRUE FALSE       TRUE
13  TRUE FALSE       TRUE
14  TRUE FALSE       TRUE
15  TRUE FALSE       TRUE
16  TRUE FALSE       TRUE
17  TRUE FALSE       TRUE
18  TRUE FALSE       TRUE
19  TRUE FALSE       TRUE
20  TRUE FALSE       TRUE

And likewise, the logistic regression spits out a sane-looking answer:

g <- glm(Conversion ~ Black*White, family=binomial, data=ground); summary(g)
...
Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.732  -0.683  -0.650  -0.643   1.832  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)           -1.472      0.114  -12.94   <2e-16
BlackTRUE              0.291      0.154    1.88    0.060
WhiteTRUE              0.137      0.156    0.88    0.381
BlackTRUE:WhiteTRUE   -0.404      0.220   -1.84    0.066

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2072.7  on 2041  degrees of freedom
Residual deviance: 2068.2  on 2038  degrees of freedom
AIC: 2076

Number of Fisher Scoring iterations: 4

So my question is: is there any more elegant way of turning my Analytics's rate data into glm input than that awful loop?

2
  • It's not true that glm requires one line per case. There is a method to use aggregated data. See ?glm.
    – IRTFM
    Commented Sep 13, 2013 at 19:12
  • I see, so it can take a matrix of counts? I would still need to convert to counts though, which I don't know how to do elegantly either without using some of the tricks suggested here.
    – gwern
    Commented Sep 14, 2013 at 19:47

3 Answers 3

1
rates$counts <- rates$N*rates$Rate
rates$counts <- round(rates$counts,0)
 rates
#----------
  Black White   N   Rate counts
1  TRUE FALSE 512 0.2344    120
2 FALSE  TRUE 529 0.2098    111
3  TRUE  TRUE 495 0.1919     95
4 FALSE FALSE 510 0.1882     96

> rates$failures <-rates$N -rates$counts    s
> glm(cbind(counts,failures)~Black*White, data=rates, family="binomial")

Call:  glm(formula = cbind(counts, failures) ~ Black * White, family = "binomial", 
    data = rates)

Coefficients:
        (Intercept)            BlackTRUE            WhiteTRUE  
            -1.4615               0.2777               0.1356  
BlackTRUE:WhiteTRUE  
            -0.3894  

Degrees of Freedom: 3 Total (i.e. Null);  0 Residual
Null Deviance:      4.104 
Residual Deviance: -7.461e-14   AIC: 33.05 
4
  • I see. Use the matrix input, and then it's easy to extract the count of failure & success from N and Rate - just think in columns/vectors. It looks so obvious in retrospect once you know how to do the input as a matrix of counts...
    – gwern
    Commented Sep 14, 2013 at 20:04
  • It is "so obvious" if the first GLM capable program you used was GLIM v4.
    – IRTFM
    Commented Sep 14, 2013 at 23:38
  • I've never used GLIM, so I wouldn't know.
    – gwern
    Commented Sep 17, 2013 at 3:17
  • Yeah, it was an old geezer comment.
    – IRTFM
    Commented Sep 17, 2013 at 5:53
1

One thing is how to convert your data. Another is why. From ?glm: "[f]or binomial [...] famil[y] the response can [...] be specified as a factor (when the first level denotes failure and all others success) or as a two-column matrix with the columns giving the numbers of successes and failures.". The first way corresponds to your "R's glm wants a dataframe of 2046 rows each reporting a TRUE or FALSE conversion". The second way basically corresponds to the your original data set where the "successes" easily can be calculated from Rate and N. A third way would be to use the proportion of successes per treatment combination as response variable, in which case the number of trials must be supplied as the weights argument.

set.seed(1)
 # one row per observation
 df1 <- data.frame(x = sample(c("yes", "no"), 40, replace = TRUE),
                 y = sample(c("yes", "no"), 40, replace = TRUE),
                 z = rbinom(n = 40, size = 1, prob = 0.5))
df1

library(plyr)
# aggregated data with one row per treatment combination
df2 <- ddply(.data = df1, .variables = .(x, y), summarize,
             n = length(z),
             rate = sum(z)/n,
             success = n*rate,
             failure = n - success)  
df2

# three different ways to specify the models,
# which all give the same parameter estimates for x, y and x*y
mod1 <- glm(z ~ x * y, data = df1, family = binomial) 
mod2 <- glm(cbind(success, failure) ~ x * y, data = df2, family = binomial)
mod3 <- glm(rate ~ x * y, data = df2, weights = n, family = binomial)

summary(mod1)
summary(mod2)
summary(mod3) 
5
  • The weights is a cool-looking approach, even simpler than the cbind solution, but I'm a little bit troubled that these apparently-equivalent approaches aren't turning in the same result. It's the same data in the same regression, isn't it? The parameters may be the same, but other values differ. For example, the complexity measure: look at mod1$aic; mod2$aic; mod3$aic and you'll see the long format is turning in a very different AIC.
    – gwern
    Commented Sep 14, 2013 at 19:59
  • @gwern, I'm sorry, but I don't understand what you mean by "it's the same data in the same regression". AIC values of the three models I present are not directly comparable, because they have different response variables and are fit using two different data sets.
    – Henrik
    Commented Sep 14, 2013 at 22:00
  • Alright, then that's a problem: I don't see why 3 ostensibly equivalent approaches are using 'two different data sets'. If that's true, doesn't that seem like a big deal? I wouldn't use an algorithm that randomly made up data for no reason, why would I use data conversion commands that made up data...
    – gwern
    Commented Sep 16, 2013 at 0:22
  • @gwern. Please read ?glm again. Especially the section I referred to in my answer: "For binomial [...] families the response can also be specified as a factor (when the first level denotes failure and all others success) or as a two-column matrix with the columns giving the numbers of successes and failures.". "For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes". Please also re-read @DWIN's first comment to your question and his answer. The data used is certainly not just "made up".
    – Henrik
    Commented Sep 16, 2013 at 0:40
  • @gwern, please also see here. My mod2 is specified in the same way as the model in DWIN's answer, both based on aggregated data. In addition, I created a data set with one observation per row (df1), on which the aggregated data is based, with the intention that you should see the link between the two data sets and the results they produced.
    – Henrik
    Commented Sep 16, 2013 at 0:50
0

Not quite clear what you're converting, but if all you need is n rows for each value in column N, then EDIT -- I was very sloppy. First thing- convert all factors in your original file to numeric or character as appropriate. then,

# just put in placeholder values
newdf<-data.frame(Black="n",White="n",Rate=0,stringsAsFactors=FALSE) 
newdf[1:rates[1,3],]<-rates[1,c(1,2,4)]
    newdf[4:rates[2,3],] <- rates[2,c(1,2,4)]

and so on for each row in your original rates dataframe.

1
  • I don't understand what you're trying to do, but it doesn't seem to work when I try it. Your definition of newdf breaks for me (Error in data.frame(Black, White, Rate) : object 'Black' not found) and when I replace it with data.frame(Black="",White="",Rate="") and then run newdf[1:rates[1,3],]... it yields a bunch of warnings (Warning messages: 1: In '[<-.factor'('*tmp*', iseq, value = c(TRUE, TRUE, TRUE, TRUE, : invalid factor level, NA generated) and then 512 rows of NAs.
    – gwern
    Commented Sep 14, 2013 at 19:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.