1
$\begingroup$

Im working on a dataset of individuals with the covid-19 infection. The dataset consist of 2010 individuals and four columns:

  • deceased: if the person died of corona (1:yes, 0:no)
  • sex: male/female
  • age: age of the person (ranging from 2 to 99 years old)
  • country: which country the person is from (France, Japan, Korea, Indonesia)

The dataset can be loaded with the following code:

id = "1CA1RPRYqU9oTIaHfSroitnWrI6WpUeBw"
d.corona = read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download",id),header = T)

Im trying to answer the following questions:

  • Is age a greater risk factor for males than for females?
  • Is age a greater risk factor for the French population than for the Korean population?
  • It seems like the French population is at much higher risk of dying from Covid-19 than the other countries (just from inspecting data). Does this seem like a valid result, or could it be influenced by the way the data were collected?

I was thinking it would be a good idea to fit some appropriate models to answer the questions. I already fit a logistic regression model with glm but I don't see how I could use that to answer these specific questions. For the very last question I believe the mean age in France is significantly higher than in the other countries, which would explain why it seems like the probability of dying from the virus is higher in France based on the data. But Im sure there could be a better explanation.

$\endgroup$
2
  • 2
    $\begingroup$ There are a number of subtle issues at play here, especially in terms of data biases. I would caution you from drawing conclusions based on a simple model without some serious consideration of how those influence patterns in your data. $\endgroup$ – mkt - Reinstate Monica Apr 24 '20 at 12:37
  • $\begingroup$ @mkt-ReinstateMonica Im not doing professional statistical analysis, this exercise is just for educational purposes. Many of the conclusions made based on this data won't agree with reality due to the limitations of the dataset. $\endgroup$ – Pame Apr 24 '20 at 13:25
2
+50
$\begingroup$

Your question did not provide an extensive example and it is not clear to me. However, for what I understand, you would like to model a simple logistic regression.
It is not clear which kind of plot do you like to produce.

I think you need to be more specific regarding your goal. However, you've a simple code to work on:

id = "1CA1RPRYqU9oTIaHfSroitnWrI6WpUeBw"
d.corona = read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download",id),header = T)

results <- data.frame()

# cycle for country
for (i in levels(d.corona$country)) 
{
  index <- i
  df <- subset(d.corona, country == index)

fit <- glm(deceased ~ sex + age, family = binomial(link = "logit"), data = df)
or <- round(exp(cbind("OR" = coef(fit), confint.default(fit, level = 0.95))), 2)
or <- data.frame(or[2:3,])
or$country <- index
results <- rbind(results, or)
rm(index, df, or, fit)
}
rm(i)

results contains:

enter image description here

As you can see, you'll have Odds Ratios. Please be more specific regarding which probabilities you're looking for (predicted probabilities?). A Good reading regarding this topic here.

Please note Japan estimation is not computed due to no events in female:

df <- subset(d.corona, country == "japan")
table(df$deceased, df$sex)

   female male
  0    120  171
  1      0    3

As you can see, age is generally the most associated variable with the outcome deceased, while sex did show same features. Notice that Indonesia has different trend in GLM estimates


Edits for generating plots

Based on your request here a very basic R example for generating plots you need:

library(visreg)
library(ggplot2)
library(gridExtra)

id = "1CA1RPRYqU9oTIaHfSroitnWrI6WpUeBw"
d.corona = read.csv(sprintf("https://docs.google.com/uc?id=%s&export=download",id),header = T)

# Setting needed spaces
results <- data.frame()
p <- list()
x <-1

# cycle for country
for (i in levels(d.corona$country)) 
{
  index <- i
  df <- subset(d.corona, country == index)
  fit <- glm(deceased ~ sex + age, family = binomial(link = "logit"), data = df)

  # Generating plot for each country
  p[[x]] <- visreg(fit, "age", by="sex", overlay=TRUE, ylab="Log odds", gg=TRUE) + ggtitle(index)

  rm(index, df, fit)
  x <- x+1
}
rm(i, x)
do.call(grid.arrange,p)

Producing these:

enter image description here

I would suggest to read visreg documentation for improving the plots.

Moreover, I would suggest to think on the possibility to dichotomize (i.e < >65yrs) or categorize (i.e. <45; 45-65; >65yrs) your age variable.


EDT III: Fatality prediction example

For the description of the following code reffering here.

# Subsetting France
df <- subset(d.corona, country == "France")
mod <- glm(deceased ~ sex + age, family = binomial(link = "logit"), data = df)

# Producing data for prediction
preddata <- data.frame(sex=rep(c("male", "female"), each = length(seq(min(df$age), max(df$age)))), 
                       age = seq(min(df$age), max(df$age)))

# Making prediction
preds <- predict(mod, preddata, type = "link", se.fit = TRUE)

# Confidence interval on the linear predictor is:
critval <- 1.96 ## approx 95% CI
upr <- preds$fit + (critval * preds$se.fit)
lwr <- preds$fit - (critval * preds$se.fit)
fit <- preds$fit

# Now for fit, upr and lwr we need to apply the inverse of the link function to them
fit2 <- mod$family$linkinv(fit)
upr2 <- mod$family$linkinv(upr)
lwr2 <- mod$family$linkinv(lwr)

library(ggplot2)
preddata$lwr <- lwr2 
preddata$upr <- upr2 
ggplot(data=df, mapping=aes(x=age,y=deceased, color = sex)) + geom_point() +         
  stat_smooth(method="glm", method.args=list(family=binomial), se = FALSE)+
  geom_line(data=preddata, mapping=aes(x = age, y=upr), linetype = 3) + 
  geom_line(data=preddata, mapping=aes(x = age, y=lwr), linetype = 3) +
  labs(x = "Age (Yrs)", y = "Fatality", title = "France")+
  theme_light()

enter image description here

$\endgroup$
6
  • $\begingroup$ If you ignore the plotting part, I think it should be fairly straight forward to find out whether age is a greater risk factor for males than for females, and whether age is a greater risk factor for the French than for the Korean. $\endgroup$ – Pame Apr 22 '20 at 12:00
  • $\begingroup$ Simply use subset() on glm. fit <- glm(deceased ~ age, family = binomial(link = "logit"), data = subset(df, sex == "male")) as suggested here $\endgroup$ – Borexino Apr 22 '20 at 13:20
  • $\begingroup$ However, I would suggest to define your main model including all the covariates related to your outcome w/o splitting your dataset. I.e. male had more than 3 times the risk of the outcome compare to female in Korea in the example I provided. $\endgroup$ – Borexino Apr 22 '20 at 13:35
  • $\begingroup$ Okay, so using the subset method you suggested and fitting two separate models for male and female, I get that the age coefficient for females is slightly higher than the age coefficient for males. So Im assuming this means age is a greater risk factor for females than for males? Since an increase in age leads to a greater odds increase of dying for females compared to males. $\endgroup$ – Pame Apr 23 '20 at 9:10
  • $\begingroup$ Also, how would you do the plot if you're plotting predicted probability of dying as a function of age separately for the two sexes and each country? $\endgroup$ – Pame Apr 23 '20 at 10:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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