The term bivariate analysis refers to the analysis of two variables. You can remember this because the prefix “bi” means “two.”
The purpose of bivariate analysis is to understand the relationship between two variables
There are three common ways to perform bivariate analysis:
1. Scatterplots
2. Correlation Coefficients
3. Simple Linear Regression
The following example shows how to perform each of these types of bivariate analysis using the following dataset that contains information about two variables: (1) Hours spent studying and (2) Exam score received by 20 different students:
#create data frame df <- data.frame(hours=c(1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 8), score=c(75, 66, 68, 74, 78, 72, 85, 82, 90, 82, 80, 88, 85, 90, 92, 94, 94, 88, 91, 96)) #view first six rows of data frame head(df) hours score 1 1 75 2 1 66 3 1 68 4 2 74 5 2 78 6 2 72
1. Scatterplots
We can use the following syntax to create a scatterplot of hours studied vs. exam score in R:
#create scatterplot of hours studied vs. exam score plot(df$hours, df$score, pch=16, col='steelblue', main='Hours Studied vs. Exam Score', xlab='Hours Studied', ylab='Exam Score')
The x-axis shows the hours studied and the y-axis shows the exam score received.
From the plot we can see that there is a positive relationship between the two variables: As hours studied increases, exam score tends to increase as well.
2. Correlation Coefficients
A Pearson Correlation Coefficient is a way to quantify the linear relationship between two variables.
We can use the cor() function in R to calculate the Pearson Correlation Coefficient between two variables:
#calculate correlation between hours studied and exam score received
cor(df$hours, df$score)
[1] 0.891306
The correlation coefficient turns out to be 0.891.
This value is close to 1, which indicates a strong positive correlation between hours studied and exam score received.
3. Simple Linear Regression
Simple linear regression is a statistical method we can use to find the equation of the line that best “fits” a dataset, which we can then use to understand the exact relationship between two variables.
We can use the lm() function in R to fit a simple linear regression model for hours studied and exam score received:
#fit simple linear regression model fit <- lm(score ~ hours, data=df) #view summary of model summary(fit) Call: lm(formula = score ~ hours, data = df) Residuals: Min 1Q Median 3Q Max -6.920 -3.927 1.309 1.903 9.385 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 69.0734 1.9651 35.15 < 2e-16 *** hours 3.8471 0.4613 8.34 1.35e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.171 on 18 degrees of freedom Multiple R-squared: 0.7944, Adjusted R-squared: 0.783 F-statistic: 69.56 on 1 and 18 DF, p-value: 1.347e-07
The fitted regression equation turns out to be:
Exam Score = 69.0734 + 3.8471*(hours studied)
This tells us that each additional hour studied is associated with an average increase of 3.8471 in exam score.
We can also use the fitted regression equation to predict the score that a student will receive based on their total hours studied.
For example, a student who studies for 3 hours is predicted to receive a score of 81.6147:
- Exam Score = 69.0734 + 3.8471*(hours studied)
- Exam Score = 69.0734 + 3.8471*(3)
- Exam Score = 81.6147
Additional Resources
The following tutorials provide additional information about bivariate analysis:
An Introduction to Bivariate Analysis
5 Examples of Bivariate Data in Real Life
An Introduction to Simple Linear Regression
I am doing my analysis where i have two binary outcome variables . What R code can i use to run bivariate copula joint modelling using gjrm package or any other package?
Hi Joachim…Hi,
I’d be happy to help you perform bivariate copula joint modeling with two binary outcome variables using R. The `GJRM` package (Generalized Joint Regression Modeling) is well-suited for this task. Below, I’ll guide you through the process and provide example code.
—
### **Overview**
Bivariate copula models allow you to model the joint distribution of two correlated outcome variables by specifying their marginal distributions and a copula function that captures the dependency structure between them. This is particularly useful when dealing with binary outcomes where traditional multivariate distributions may not be appropriate.
### **Using the `GJRM` Package**
The `GJRM` package in R provides a flexible framework for fitting a wide range of joint regression models, including those using copulas for dependence modeling.
#### **Installation**
First, install and load the `GJRM` package if you haven’t already:
“`R
install.packages(“GJRM”)
library(GJRM)
“`
#### **Preparing Your Data**
Assuming you have a dataset `mydata` with two binary outcome variables `Y1` and `Y2`, and predictor variables `X1`, `X2`, etc.
“`R
# Example data
mydata <- data.frame( Y1 = rbinom(200, 1, 0.5), Y2 = rbinom(200, 1, 0.5), X1 = rnorm(200), X2 = rnorm(200) ) ``` #### **Specifying Marginal Models** Define the formulas for each of the binary outcome variables: ```R # Marginal models for Y1 and Y2 f1 <- Y1 ~ X1 + X2 f2 <- Y2 ~ X1 + X2 ``` #### **Fitting the Bivariate Copula Model** Use the `gjrm()` function to fit the joint model: ```R # List of formulas formula_list <- list(f1, f2) # Fit the model fit <- gjrm( formula = formula_list, data = mydata, margins = c("probit", "probit"), Model = "Bivariate", BivD = "Frank" # Choose the copula family ) ``` **Parameters Explained:** - `formula`: A list containing the formulas for the marginal models. - `data`: The dataset containing all variables. - `margins`: Specifies the distribution and link function for the marginals; `"probit"` is commonly used for binary outcomes. - `Model`: Set to `"Bivariate"` to indicate a bivariate model. - `BivD`: Specifies the copula family (e.g., `"Frank"`, `"Clayton"`, `"Gumbel"`, `"Joe"`, `"AMH"`, `"FGM"`, `"Galambos"`, `"Hüsler–Reiss"`, `"Tawn"`, `"Plackett"`). #### **Examining the Results** ```R # Summarize the fitted model summary(fit) ``` The summary will provide estimates for: - The coefficients of the marginal models. - The copula parameter, which indicates the dependency between `Y1` and `Y2`. --- ### **Interpreting the Output** - **Marginal Effects**: Look at the coefficients for `X1`, `X2` in both marginals to understand their effects on `Y1` and `Y2`. - **Copula Parameter**: A significant copula parameter suggests a dependency between the two outcomes. --- ### **Alternative Copula Families** You can try different copula families to see which fits your data best. ```R # Using the Clayton copula fit_clayton <- gjrm( formula = formula_list, data = mydata, margins = c("probit", "probit"), Model = "Bivariate", BivD = "Clayton" ) # Compare models using AIC AIC(fit, fit_clayton) ``` --- ### **Predicting and Plotting** #### **Predictions** ```R # Predict probabilities for the marginals predictions <- predict(fit, newdata = mydata, type = "response") ``` #### **Plotting the Results** ```R # Load plotting package library(ggplot2) # Extract fitted values fitted_values <- data.frame( Y1 = mydata$Y1, Y2 = mydata$Y2, Fitted_Y1 = fit$fitted.values[[1]], Fitted_Y2 = fit$fitted.values[[2]] ) # Plot observed vs. fitted for Y1 ggplot(fitted_values, aes(x = Fitted_Y1, y = Y1)) + geom_jitter(height = 0.1, width = 0) + geom_smooth(method = "loess") + labs(title = "Observed vs. Fitted for Y1") # Plot observed vs. fitted for Y2 ggplot(fitted_values, aes(x = Fitted_Y2, y = Y2)) + geom_jitter(height = 0.1, width = 0) + geom_smooth(method = "loess") + labs(title = "Observed vs. Fitted for Y2") ``` --- ### **Using Other Packages** While `GJRM` is highly recommended for its flexibility and ease of use, there are other packages you might consider: #### **1. `VGAM` Package** The `VGAM` package can handle some bivariate models, but it's less flexible for copula models with binary outcomes. #### **2. `copula` Package** The `copula` package is more general-purpose and requires manual specification of the likelihood, which can be complex for binary data. --- ### **Complete Example** Putting it all together: ```R # Load necessary packages install.packages("GJRM") library(GJRM) library(ggplot2) # Example data mydata <- data.frame( Y1 = rbinom(200, 1, 0.5), Y2 = rbinom(200, 1, 0.5), X1 = rnorm(200), X2 = rnorm(200) ) # Define marginal models f1 <- Y1 ~ X1 + X2 f2 <- Y2 ~ X1 + X2 formula_list <- list(f1, f2) # Fit the bivariate copula model fit <- gjrm( formula = formula_list, data = mydata, margins = c("probit", "probit"), Model = "Bivariate", BivD = "Frank" ) # Summarize the results summary(fit) # Predict probabilities predictions <- predict(fit, newdata = mydata, type = "response") # Plot observed vs. fitted values fitted_values <- data.frame( Y1 = mydata$Y1, Y2 = mydata$Y2, Fitted_Y1 = fit$fitted.values[[1]], Fitted_Y2 = fit$fitted.values[[2]] ) ggplot(fitted_values, aes(x = Fitted_Y1, y = Y1)) + geom_jitter(height = 0.1, width = 0) + geom_smooth(method = "loess") + labs(title = "Observed vs. Fitted for Y1") ggplot(fitted_values, aes(x = Fitted_Y2, y = Y2)) + geom_jitter(height = 0.1, width = 0) + geom_smooth(method = "loess") + labs(title = "Observed vs. Fitted for Y2") ``` --- ### **Additional Tips** - **Missing Data**: Ensure your dataset has no missing values in the variables used. - **Copula Selection**: Try different copula families to find the best fit. Use information criteria like AIC for comparison. - **Diagnostics**: Check residuals and perform diagnostics to assess model fit. - **Variable Selection**: Include relevant predictors in the marginal models to improve accuracy. --- ### **References** - **GJRM Package Documentation**: [CRAN - Package GJRM](https://cran.r-project.org/package=GJRM) - **Marra, G., & Radice, R. (2017)**. Semi‐parametric copula-based regression models for longitudinal data. *Scandinavian Journal of Statistics*, 44(3), 647-664. ---