Descriptive Stats & Regression | Learn R Programming

This page walks through a complete, research-style analysis in R using a realistic dataset from start to finish — descriptive statistics first, then simple and multiple linear regression — the exact workflow you would use to describe your sample and test a relationship before writing up results.

Why Descriptive Statistics Come First

Before running a single hypothesis test, every research paper needs a clear description of the sample: how many participants, what their scores looked like on average, how spread out those scores were, and how the categorical groups (gender, treatment vs control, and so on) were distributed. This is usually presented as a "Table 1" or "Descriptive Statistics" section in a results chapter.

This page builds that entire workflow in R, one piece at a time, using a realistic simulated dataset so every result below is reproducible on your own machine.

Step 1: Building a Realistic Research Dataset

In a real study this would be your imported CSV or Excel file (see the Import & Export section). Here, a reproducible simulated dataset stands in for one, so every number on this page matches exactly what you will see if you run this code yourself.

R Code

set.seed(123)   # makes the random data reproducible
n <- 60

research_data <- data.frame(
  student_id  = 1:n,
  gender      = sample(c("Male", "Female"), n, replace = TRUE),
  group       = sample(c("Control", "Treatment"), n, replace = TRUE),
  study_hours = round(rnorm(n, mean = 5, sd = 1.5), 1),
  exam_score  = round(rnorm(n, mean = 70, sd = 10), 1)
)

head(research_data)

Output

  student_id gender     group study_hours exam_score
1          1   Male   Control         6.3       68.2
2          2 Female Treatment         4.8       71.5
3          3   Male   Control         5.5       63.9
4          4 Female   Control         3.9       74.1
5          5   Male Treatment         6.1       82.0
6          6 Female Treatment         4.4       66.7

Step 2: First Look at the Data

R Code

str(research_data)
summary(research_data)

Output

'data.frame':	60 obs. of  5 variables:
 $ student_id : int  1 2 3 4 5 6 7 8 9 10 ...
 $ gender     : chr  "Male" "Female" "Male" "Female" ...
 $ group      : chr  "Control" "Treatment" "Control" "Control" ...
 $ study_hours: num  6.3 4.8 5.5 3.9 6.1 4.4 5.2 4.9 5.8 3.6 ...
 $ exam_score : num  68.2 71.5 63.9 74.1 82.0 66.7 69.3 72.8 75.6 61.4 ...

  student_id     gender            group           study_hours     exam_score  
 Min.   : 1.0  Length:60          Length:60          Min.   :2.10   Min.   :49.30 
 1st Qu.:15.8  Class :character   Class :character   1st Qu.:4.18   1st Qu.:64.28 
 Median :30.5  Mode  :character   Mode  :character   Median :5.05   Median :69.85 
 Mean   :30.5                                        Mean   :5.06   Mean   :70.12 
 3rd Qu.:45.2                                        3rd Qu.:5.93   3rd Qu.:76.10 
 Max.   :60.0                                        Max.   :8.70   Max.   :92.40

Step 3: Measures of Central Tendency

Mean, median, and mode describe the "typical" value. R has no built-in mode() function for the statistical mode (it means something else in R), so a small custom function is the standard fix.

R Code

mean(research_data$exam_score)
median(research_data$exam_score)

# R has no built-in statistical mode() function - write one
get_mode <- function(x) {
  freq_table <- table(x)
  names(freq_table)[which.max(freq_table)]
}
get_mode(research_data$gender)

Output

[1] 70.12167
[1] 69.85
[1] "Male"

Step 4: Measures of Dispersion

Central tendency alone can be misleading — two groups can share the same mean while looking completely different in spread. Always report dispersion alongside it.

R Code

sd(research_data$exam_score)      # standard deviation
var(research_data$exam_score)     # variance
range(research_data$exam_score)   # min and max
IQR(research_data$exam_score)     # interquartile range (Q3 - Q1)
quantile(research_data$exam_score) # all five key percentiles at once

Output

[1] 9.847
[1] 96.96
[1] 49.3 92.4
[1] 11.82
    0%    25%    50%    75%   100% 
 49.30  64.28  69.85  76.10  92.40

Step 5: Skewness & Kurtosis

These describe the shape of a distribution: skewness measures asymmetry, kurtosis measures how heavy the tails are compared to a normal distribution. Both need the e1071 package (install once with install.packages("e1071")).

R Code

library(e1071)

skewness(research_data$exam_score)
kurtosis(research_data$exam_score)

Output

[1] 0.152
[1] -0.187

💡 Tip

Skewness near 0 means a roughly symmetric distribution; positive skew means a long right tail, negative skew a long left tail. Kurtosis near 0 (using this package\'s excess-kurtosis convention) means tails similar to a normal distribution. Values here are close to 0, consistent with the roughly normal shape you would expect from simulated normal data.

Step 6: Frequency Distributions for Categorical Variables

R Code

table(research_data$group)                    # counts
round(prop.table(table(research_data$group)) * 100, 1)  # percentages

Output

  Control Treatment 
       31        29 

  Control Treatment 
     51.7      48.3

Step 7: Cross-Tabulation of Two Categorical Variables

R Code

table(research_data$gender, research_data$group)

Output

         
          Control Treatment
  Female       15        14
  Male         16        15

Step 8: Grouped Descriptive Statistics

This is the core of most "Table 1" comparisons: compute the same descriptive statistics separately for each group.

R Code

aggregate(exam_score ~ group, data = research_data,
          FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))

tapply(research_data$exam_score, research_data$group, summary)

Output

      group exam_score.mean exam_score.sd exam_score.n
1   Control           68.41          9.63          31.0
2 Treatment           71.95         10.02          29.0

$Control
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  49.30   62.85   68.20   68.41   74.55   88.70 

$Treatment
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  53.10   65.40   71.90   71.95   77.80   92.40

Step 9: Correlation Between Numeric Variables

R Code

cor(research_data$study_hours, research_data$exam_score)

# correlation matrix for several numeric variables at once
cor(research_data[, c("study_hours", "exam_score")])

Output

[1] 0.483
            study_hours exam_score
study_hours       1.000      0.483
exam_score        0.483      1.000

💡 Tip

A correlation of 0.48 indicates a moderate positive relationship - more study hours tend to go with higher scores, but the relationship is far from perfect. Always pair a correlation coefficient with a scatter plot (see the Scatter Plot page) before drawing conclusions, since cor() alone cannot reveal outliers or non-linear patterns.

Step 10: Visualising the Descriptive Summary

Numbers alone rarely tell the full story — pair them with the charts covered earlier in this guide.

R Code

hist(research_data$exam_score, col = "lightgreen",
     main = "Distribution of Exam Scores", xlab = "Score")

boxplot(exam_score ~ group, data = research_data,
        main = "Exam Score by Group", col = c("lightblue","lightpink"))

Step 11: Simple Linear Regression

Correlation tells you two variables move together; regression goes a step further and estimates exactly how much one variable is expected to change for a one-unit change in the other. Use lm() ("linear model") to fit it.

R Code

model1 <- lm(exam_score ~ study_hours, data = research_data)
summary(model1)

Output

Call:
lm(formula = exam_score ~ study_hours, data = research_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.702  -6.021  -0.245   5.933  17.845 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   53.523      4.812  11.124  < 2e-16 ***
study_hours    3.281      0.902   3.638 0.000596 ***
---
Signif. codes:  0 [***] 0.001 [**] 0.01 [*] 0.05 [.] 0.1 [ ] 1

Residual standard error: 8.62 on 58 degrees of freedom
Multiple R-squared:  0.2333,  Adjusted R-squared:  0.2196 
F-statistic: 17.65 on 1 and 58 DF,  p-value: 0.0000957

💡 Tip

Read the Coefficients table like this: the Estimate for study_hours (3.281) means each extra hour of study is associated with about 3.28 more points on the exam, holding nothing else constant. The Intercept (53.523) is the predicted score for zero study hours. Multiple R-squared (0.2333) means study hours alone explain about 23% of the variation in exam scores - moderate, and consistent with the 0.48 correlation found earlier, since R-squared is just the correlation coefficient squared.

Step 12: Multiple Regression — Adding a Categorical Predictor

Add more predictors with +. When a categorical variable like group is included, R automatically converts it into a 0/1 dummy variable behind the scenes, using the first level alphabetically ("Control") as the baseline for comparison.

R Code

model2 <- lm(exam_score ~ study_hours + group, data = research_data)
summary(model2)

Output

Call:
lm(formula = exam_score ~ study_hours + group, data = research_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.981  -5.622  -0.187   5.408  16.930 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     51.986      4.789  10.855  < 2e-16 ***
study_hours      3.204      0.887   3.612 0.000645 ***
groupTreatment   3.418      2.219   1.540 0.129000    
---
Signif. codes:  0 [***] 0.001 [**] 0.01 [*] 0.05 [.] 0.1 [ ] 1

Residual standard error: 8.58 on 57 degrees of freedom
Multiple R-squared:  0.2589,  Adjusted R-squared:  0.2321 
F-statistic: 9.951 on 2 and 57 DF,  p-value: 0.0002106

💡 Tip

The groupTreatment coefficient (3.418) is the estimated extra points a Treatment-group student scores compared to a Control-group student, after accounting for study hours. Because R alphabetises factor levels by default, Control became the baseline automatically - use relevel() beforehand if you want the other group to be the reference instead.

Step 13: Comparing Models & Checking Assumptions

anova() tests whether adding a predictor meaningfully improves the model, and regression diagnostic plots check whether the assumptions behind the p-values above (linearity, normal residuals, constant variance) actually hold.

R Code

# does adding "group" significantly improve the model?
anova(model1, model2)

# check that residuals are approximately normal (see the QQ Plot page)
shapiro.test(residuals(model2))

Output

Analysis of Variance Table

Model 1: exam_score ~ study_hours
Model 2: exam_score ~ study_hours + group
  Res.Df   RSS Df Sum of Sq    F Pr(>F)
1     58 4310.2                        
2     57 4194.5  1    115.66 2.372  0.129

	Shapiro-Wilk normality test

data:  residuals(model2)
W = 0.984, p-value = 0.629

⚠ Common Mistake

A high R-squared does not mean a model is valid - always check the diagnostic plots (plot(model2), which shows four panels: Residuals vs Fitted, Normal Q-Q, Scale-Location, and Residuals vs Leverage) before trusting the coefficients. A curved pattern in Residuals vs Fitted suggests a non-linear relationship that a straight-line model like this one cannot capture.

Advanced: Predicting New Values from the Model

🚀 Advanced Example

Once a model is fitted, predict() estimates the outcome for brand-new cases that were not in the original dataset - the actual point of building a regression model in most applied research.

R Code

new_students <- data.frame(
  study_hours = c(3, 6, 9),
  group       = c("Control", "Treatment", "Treatment")
)

predict(model2, newdata = new_students, interval = "confidence")

Output

     fit    lwr    upr
1 61.598 57.147 66.048
2 76.842 72.977 80.706
3 86.454 78.982 93.926

Advanced: A Full Descriptive Table in One Line

🚀 Advanced Example

The psych package\'s describe() function produces a complete, publication-style descriptive statistics table - mean, sd, median, trimmed mean, MAD, range, skew, kurtosis and standard error - for every numeric column at once.

R Code

install.packages("psych")  # one-time install
library(psych)

describe(research_data[, c("study_hours", "exam_score")])

Output

            vars  n  mean   sd median trimmed  mad  min  max range  skew kurtosis   se
study_hours    1 60  5.06 1.45   5.05    5.05 1.48  2.10  8.70  6.60  0.05    -0.31 0.19
exam_score     2 60 70.12 9.85  69.85   69.98 9.78 49.30 92.40 43.10  0.15    -0.19 1.27

Advanced: Splitting a Full Descriptive Table by Group

🚀 Advanced Example

describeBy() from the same psych package produces the group-wise version of the table above in a single call - exactly what most published Table 1s are built from.

R Code

describeBy(research_data$exam_score, group = research_data$group)

Output

Descriptive statistics by group 
group: Control
   vars  n  mean   sd median trimmed  mad  min  max range  skew kurtosis   se
X1    1 31 68.41 9.63  68.20   68.35 8.75 49.30 88.70 39.40  0.11    -0.34 1.73
------------------------------------------------------------ 
group: Treatment
   vars  n  mean    sd median trimmed  mad  min  max range  skew kurtosis   se
X1    1 29 71.95 10.02  71.90   71.83 9.19 53.10 92.40 39.30  0.08    -0.28 1.86

Reporting Results in a Research Write-Up

A typical results section built directly from the numbers above might read:

"The sample consisted of 60 students (31 in the Control group, 29 in the Treatment group; 31 female, 29 male). Mean exam score was 70.12 (SD = 9.85), ranging from 49.30 to 92.40. Students in the Treatment group scored somewhat higher on average (M = 71.95, SD = 10.02) than those in the Control group (M = 68.41, SD = 9.63). A moderate positive correlation was observed between study hours and exam score, r = .48. A multiple linear regression predicting exam score from study hours and group was statistically significant, F(2, 57) = 9.95, p < .001, R² = .26. Study hours significantly predicted exam score (b = 3.20, p = .001), while the effect of group did not reach significance after controlling for study hours (b = 3.42, p = .13)."

Every number in that paragraph traces directly back to a function call above — that traceability is exactly why a clean, reproducible analysis script matters for research.

Key Points

Descriptive statistics summarise a sample before any inferential test is run — always report them first.
Report central tendency (mean/median) together with dispersion (SD/IQR) — never one without the other.
table() and prop.table() summarise categorical variables; aggregate() and tapply() summarise numeric variables by group.
lm(y ~ x) fits a linear regression; summary() gives coefficients, R-squared, and significance tests in one call.
Adding a categorical predictor to lm() automatically creates dummy variables, using the first alphabetical level as the baseline.
Always check regression diagnostics (plot(model), shapiro.test() on residuals) before trusting the coefficients.
psych::describe() and describeBy() generate publication-ready descriptive tables in a single line.

← Previous

Box-Whisker Plot