Reputation: 307
I would like to reproduce some calculations from a book (logit regression). The book gives a contingency table and the results.
Here is the Table:
.
example <- matrix(c(21,22,6,51), nrow = 2, byrow = TRUE)
#Labels:
rownames(example) <- c("Present","Absent")
colnames(example) <- c(">= 55", "<55")
It gives me this:
>= 55 <55
Present 21 22
Absent 6 51
But to use the glm()-function the data has to be in the following way:
(two colums, one with "Age", and one with "Present", filled with 0/1)
age <- c(rep(c(0),27), rep(c(1),73))
present <- c(rep(c(0),21), rep(c(1),6), rep(c(0),22), rep(c(1),51))
data <- data.frame(present, age)
> data
present age
1 0 0
2 0 0
3 0 0
. . .
. . .
. . .
100 1 1
Is there a simple way to get this structure from the table/matrix?
Upvotes: 4
Views: 1809
Reputation: 2743
So, glm
is not quite that inflexible. In part ?glm
reads
For ‘binomial’ and ‘quasibinomial’ 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.
I'll assume you want to test the effect of age on Present/Absent
.
The key is for to specify the response like (in psueudo-code) c(success, failure)
.
So you need data like data.frame(Age= ..., Present = ..., Absent)
. The easiest way to do this from your example
is to transpose, then coerce to data.frame
, and add a column:
example_t <- as.data.frame(t(example))
example_df <- data.frame(example_t, Age=factor(row.names(example_t)))
which gives you
Present Absent Age
>= 55 21 6 >= 55
<55 22 51 <55
Then, you can run the glm:
glm(cbind(Present, Absent) ~ Age, example_df, family = 'binomial')
to get
Call: glm(formula = cbind(Present, Absent) ~ Age, family = "binomial",
data = example_for_glm)
Coefficients:
(Intercept) Age<55
1.253 -2.094
Degrees of Freedom: 1 Total (i.e. Null); 0 Residual
Null Deviance: 18.7
Residual Deviance: -1.332e-15 AIC: 11.99
Addendum
You could also get here via the answer by @therimalaya. But it's just the first step
as.data.frame(as.table(example))
(only gets you part way there)
Var1 Var2 Freq
1 Present >= 55 21
2 Absent >= 55 6
3 Present <55 22
4 Absent <55 51
but to actually have a column of successes and failures, you need to do something more. You could use tidyr
to get there
as.data.frame(as.table(example)) %>% tidyr::spread(Var1, Freq)
is similar to my example_df
above
Var2 Present Absent
1 >= 55 21 6
2 <55 22 51
Upvotes: 1
Reputation: 4592
reshape2::melt(example)
This will give you,
Var1 Var2 value
1 Present >= 55 21
2 Absent >= 55 6
3 Present <55 22
4 Absent <55 51
which you can easily use for glm
Upvotes: 2
Reputation: 11419
The code below might look long but only the group_by()
and do()
instruction deal with expanding the data. All the rest is about changing the data in long format and encoding character variables as 0 and 1. I tried to start from the exact matrix you gave in your question.
Load data manipulation packages
library(tidyr)
library(dplyr)
Create a matrix as in your example, but avoid ">" signs in column names
example <- matrix(c(21,22,6,51), nrow = 2, byrow = TRUE)
rownames(example) <- c("Present","Absent")
colnames(example) <- c("above55", "below55")
Convert the matrix to a data frame
example <- data.frame(example) %>%
add_rownames("chd")
Or simply create a data frame directly
data.frame(chd = c("Present", "Absent"),
above55 = c(21,6),
below55 = c(22,51))
data2 <- example %>%
gather(age, nrow, -chd) %>%
# Encode chd and age as 0 or 1
mutate(chd = ifelse(chd=="Present",1,0),
age = ifelse(age=="above55",1,0)) %>%
group_by(chd, age) %>%
# Expand each variable by nrow
do(data.frame(chd = rep(.$chd,.$nrow),
age = rep(.$age,.$nrow)))
head(data2)
# Source: local data frame [6 x 2]
# Groups: chd, age [1]
#
# chd age
# (dbl) (dbl)
# 1 0 0
# 2 0 0
# 3 0 0
# 4 0 0
# 5 0 0
# 6 0 0
tail(data2)
# Source: local data frame [6 x 2]
# Groups: chd, age [1]
#
# chd age
# (dbl) (dbl)
# 1 1 1
# 2 1 1
# 3 1 1
# 4 1 1
# 5 1 1
# 6 1 1
table(data2)
# age
# chd 0 1
# 0 51 6
# 1 22 21
Same as your example except for the age encoding issue mentioned in my comment above.
Upvotes: 1
Reputation: 2362
I would go for:
library(data.table)
tab <- data.table(AGED = c(1, 1, 0, 0),
CHD = c(1, 0, 1, 0),
Count = c(21, 6, 22, 51))
tabExp <- tab[rep(1:.N, Count), .(AGED, CHD)]
Edit: Quick explanation, as it took me some time to figure it out:
In data.table
objects .N
stores the number of rows of a group (if grouped with by
) or just the number of rows of the whole data.table
, so in this example:
tab[rep(1:.N, Count)]
and
tab[rep(1:4, Count)]
and finally
tab[rep(1:4, c(21, 6, 22, 51)]
are equivalent.
Same with base R:
tab2 <- data.frame(AGED = c(1, 1, 0, 0),
CHD = c(1, 0, 1, 0),
Count = c(21, 6, 22, 51))
tabExp2 <- tab2[rep(1:nrow(tab2), tab2$Count), c("AGED", "CHD")]
Upvotes: 1
Reputation: 24198
You could perhaps use the countsToCases
function as defined here.
countsToCases(as.data.frame(as.table(example)))
# Var1 Var2
#1 Present >= 55
#1.1 Present >= 55
#1.2 Present >= 55
#1.3 Present >= 55
#1.4 Present >= 55
#1.5 Present >= 55
# ...
You can always recode the variables to numeric afterwards, if you prefer.
Upvotes: 2