3C 279: Modelling Brightness and Variability Before vs. After the 2023 Dimming

A mixed-effects (GLMM/LMM) analysis of AAVSO magnitudes

Author

AAVSO — analysis for Rodney H

Published

July 1, 2026

# Utility functions ----------------------------------------------------------

#' Julian Date -> calendar Date (UTC).
#' JD 2440587.5 corresponds to 1970-01-01 00:00 UTC (the Unix epoch).
jd_to_date <- function(jd) as.Date(jd - 2440587.5, origin = "1970-01-01")

#' Standard error of the mean
se <- function(x) sqrt(var(x, na.rm = TRUE) / sum(!is.na(x)))

#' +/-3 SD summary (kept from the AAVSO template for a familiar readout)
sd3 <- function(x) {
  m <- round(mean(x, na.rm = TRUE), 3)
  s <- round(sd(x,   na.rm = TRUE), 3)
  tibble(
    label = c("-3sd", "-2sd", "-1sd", "Mean", "+1sd", "+2sd", "+3sd"),
    value = c(m - 3*s, m - 2*s, m - s, m, m + s, m + 2*s, m + 3*s)
  )
}

Abstract

In plain language. This report asks whether the quasar 3C 279 has really dimmed since 2023, and whether its brightness has become steadier or more erratic — using AAVSO amateur observations going back to about 1999.

The answer to both is yes, it changed — it is now fainter and quieter. Comparing the years before 2022 with 2023 onward (2022 is left out as a transition), and looking at the well-standardised Johnson V filter so we compare like with like, 3C 279 is about 1.7 magnitudes fainter than before (a magnitude is a factor of ~2.5 in brightness, so this is a large drop), and its night-to-night scatter has roughly halved. In other words, the object didn’t just dim — it also settled into a calmer, less variable state.

One important subtlety, spelled out in the analysis: if you lump all filters together, the “quieter” signal almost disappears. That is not because the source stayed variable — it is because different filters were used in different eras and they have different natural scatter. Once we hold the filter fixed (the V-only check), the drop in variability is clear. This is exactly the kind of confound the mixed-effects model is designed to handle.

Everything below shows the data cleaning, the statistical models, and the diagnostic checks behind these two conclusions.

1 Data

1.1 Load and clean

The raw file has columns JD, mag, I_Magnitude, band, val.

I_Magnitude is constant (all 0) and carries no information — dropped.
val is a free-text note field (“Z”, “V”, exposure strings) — not used in the model.
band is the photometric passband. Magnitudes are not comparable across bands, so band is a first-class effect in the model, and non-photometric / junk codes (N/A, TG, SR, CV) are removed.
Rows with a missing mag are dropped.

raw <- read_csv(DATA_CSV, show_col_types = FALSE)

# Photometric bands we trust for a magnitude comparison.
KEEP_BANDS <- c("V", "CR", "B", "R", "I", "Vis.")

dat <- raw |>
  transmute(
    jd    = as.numeric(JD),
    mag   = as.numeric(mag),
    band  = as.character(band)
  ) |>
  filter(!is.na(jd), !is.na(mag), band %in% KEEP_BANDS) |>
  mutate(
    date  = jd_to_date(jd),
    year  = year(date),
    band  = factor(band, levels = KEEP_BANDS)
  ) |>
  # Guard against obviously bad magnitudes (3C 279 lives ~11-18 mag).
  filter(mag > 8, mag < 20) |>
  arrange(date)

cat(sprintf("Rows: %d raw -> %d after cleaning\n", nrow(raw), nrow(dat)))

Rows: 7850 raw -> 7765 after cleaning

glimpse(dat)

Rows: 7,765
Columns: 5
$ jd   <dbl> 2451233, 2451249, 2451250, 2451252, 2451252, 2451253, 2451254, 24…
$ mag  <dbl> 14.80, 14.20, 15.27, 14.20, 15.24, 15.27, 14.20, 15.16, 14.20, 15…
$ band <fct> Vis., Vis., Vis., Vis., Vis., Vis., Vis., Vis., Vis., Vis., Vis.,…
$ date <date> 1999-02-23, 1999-03-11, 1999-03-12, 1999-03-14, 1999-03-14, 1999…
$ year <dbl> 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999,…

1.2 Define the epoch contrast

model_dat <- dat |>
  mutate(
    epoch = case_when(
      date <  PRE_END      ~ "pre_2022",
      date >= DIMMING_FROM ~ "dimming_2023on",
      TRUE                 ~ NA_character_        # 2022 transition -> excluded
    )
  ) |>
  filter(!is.na(epoch)) |>
  mutate(epoch = factor(epoch, levels = c("pre_2022", "dimming_2023on")))

# How much data, and which bands, land in each epoch?
model_dat |>
  count(epoch, band) |>
  pivot_wider(names_from = band, values_from = n, values_fill = 0) |>
  kable(caption = "Observation counts by epoch and band (2022 excluded as transition)")

Observation counts by epoch and band (2022 excluded as transition)
epoch	V	CR	B	R	I	Vis.
pre_2022	2782	2622	770	229	85	739
dimming_2023on	206	0	158	1	50	52

model_dat |>
  group_by(epoch) |>
  summarise(
    n         = n(),
    mean_mag  = mean(mag),
    sd_mag    = sd(mag),
    var_mag   = var(mag),
    .groups   = "drop"
  ) |>
  kable(digits = 3,
        caption = "Brightness (mean) and variability (SD/variance) by epoch")

Brightness (mean) and variability (SD/variance) by epoch
epoch	n	mean_mag	sd_mag	var_mag
pre_2022	7227	14.527	1.330	1.769
dimming_2023on	467	16.942	0.866	0.749

Reminder on magnitudes: larger mag = fainter. If the object has dimmed since 2023, we expect a higher mean magnitude in the dimming_2023on epoch.

2 The light curve

p_lc <- ggplot(dat, aes(x = date, y = mag, colour = band)) +
  geom_point(alpha = 0.35, size = 0.9) +
  # Shade the excluded 2022 transition band.
  annotate("rect", xmin = PRE_END, xmax = DIMMING_FROM,
           ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "grey20") +
  geom_vline(xintercept = c(PRE_END, DIMMING_FROM),
             linetype = "dashed", colour = "grey40") +
  scale_y_reverse() +                       # brighter (smaller mag) at the top
  labs(
    title    = "3C 279 — AAVSO light curve with epoch split",
    subtitle = "Shaded band = 2022 transition (excluded from the contrast)",
    x = NULL, y = "Magnitude (fainter →)", colour = "Band"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

p_lc

3 Distribution by epoch

ggplot(model_dat, aes(x = mag, fill = epoch)) +
  geom_density(alpha = 0.45) +
  scale_x_reverse() +
  scale_fill_manual(values = c(pre_2022 = "steelblue",
                               dimming_2023on = "firebrick")) +
  labs(title = "Magnitude distribution: pre-2022 vs. 2023-onward",
       x = "Magnitude (fainter →)", y = "Density", fill = "Epoch") +
  theme_minimal(base_size = 12)

4 Model 1 — Mean brightness (LMM)

We first model the mean magnitude. Magnitude is a continuous measurement, so this is a linear mixed model (Gaussian response), not a count GLM.

Fixed effect: epoch — the pre-vs-post contrast we care about.
Random effect: (1 | band) — each passband gets its own intercept, absorbing the systematic zero-point differences between filters.

\[\text{mag}_{ij} = \beta_0 + \beta_1\,\text{epoch}_i + b_j^{\text{band}} + \varepsilon_{ij}\]

m_mean <- lmer(mag ~ epoch + (1 | band), data = model_dat, REML = TRUE)
summary(m_mean)

Linear mixed model fit by REML ['lmerMod']
Formula: mag ~ epoch + (1 | band)
   Data: model_dat

REML criterion at convergence: 14681.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7096 -0.4497  0.1188  0.5756  4.8277 

Random effects:
 Groups   Name        Variance Std.Dev.
 band     (Intercept) 1.1621   1.0780  
 Residual             0.3922   0.6262  
Number of obs: 7694, groups:  band, 6

Fixed effects:
                    Estimate Std. Error t value
(Intercept)         14.99644    0.44030   34.06
epochdimming_2023on  1.40910    0.03116   45.23

Correlation of Fixed Effects:
            (Intr)
epchdm_2023 -0.008

# The epoch coefficient = change in mean magnitude, holding band constant.
tidy(m_mean, effects = "fixed", conf.int = TRUE) |>
  kable(digits = 4,
        caption = "Fixed effects. 'epochdimming_2023on' > 0 means fainter since 2023.")

Fixed effects. ‘epochdimming_2023on’ > 0 means fainter since 2023.
effect	term	estimate	std.error	statistic	conf.low	conf.high
fixed	(Intercept)	14.9964	0.4403	34.0595	14.1335	15.8594
fixed	epochdimming_2023on	1.4091	0.0312	45.2272	1.3480	1.4702

5 Model 2 — Variability (location–scale GLMM)

The central question is about variability, not just the average. A standard mixed model assumes one common residual variance. Here we let the residual dispersion depend on epoch via glmmTMB’s dispformula, so the model estimates a separate spread for each era and tests whether they differ.

m_scale <- glmmTMB(
  mag ~ epoch + (1 | band),
  dispformula = ~ epoch,          # <- variance is allowed to change by epoch
  data   = model_dat,
  family = gaussian()
)
summary(m_scale)

 Family: gaussian  ( identity )
Formula:          mag ~ epoch + (1 | band)
Dispersion:           ~epoch
Data: model_dat

      AIC       BIC    logLik -2*log(L)  df.resid 
  14684.0   14718.8   -7337.0   14674.0      7689 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 band     (Intercept) 0.9695   0.9846  
 Residual                 NA       NA  
Number of obs: 7694, groups:  band, 6

Conditional model:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         14.99618    0.40219   37.29   <2e-16 ***
epochdimming_2023on  1.40862    0.02978   47.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dispersion model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -0.464897   0.008348  -55.69   <2e-16 ***
epochdimming_2023on -0.055637   0.035598   -1.56    0.118    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Back-transform the dispersion model to a residual SD per epoch.
disp <- fixef(m_scale)$disp
sd_pre     <- sqrt(exp(disp[["(Intercept)"]]))
sd_dimming <- sqrt(exp(disp[["(Intercept)"]] + disp[["epochdimming_2023on"]]))

tibble(
  epoch          = c("pre_2022", "dimming_2023on"),
  residual_sd    = c(sd_pre, sd_dimming)
) |>
  kable(digits = 3,
        caption = "Model-estimated residual SD by epoch (the variability answer)")

Model-estimated residual SD by epoch (the variability answer)
epoch	residual_sd
pre_2022	0.793
dimming_2023on	0.771

cat(sprintf(
  "Variability change: residual SD goes from %.3f (pre-2022) to %.3f (2023 on) — a %.0f%% %s.\n",
  sd_pre, sd_dimming,
  100 * abs(sd_dimming - sd_pre) / sd_pre,
  ifelse(sd_dimming > sd_pre, "increase", "decrease")
))

Variability change: residual SD goes from 0.793 (pre-2022) to 0.771 (2023 on) — a 3% decrease.

5.1 Cross-check: simple variance test

An independent, model-free check on the same claim.

var.test(mag ~ epoch, data = model_dat) |>
  broom::tidy() |>
  kable(digits = 4, caption = "F test for equal variances between epochs")

F test for equal variances between epochs
estimate	num.df	den.df	statistic	p.value	conf.low	conf.high	method	alternative
2.3606	7226	466	2.3606	0	2.0593	2.6848	F test to compare two variances	two.sided

6 Diagnostics

diag_df <- tibble(
  fitted   = fitted(m_mean),
  residual = resid(m_mean)
)

p1 <- ggplot(diag_df, aes(fitted, residual)) +
  geom_point(alpha = 0.3, size = 1, colour = "steelblue") +
  geom_hline(yintercept = 0, colour = "firebrick") +
  geom_smooth(method = "loess", se = FALSE,
              colour = "grey30", linetype = "dashed") +
  labs(title = "Residuals vs Fitted", x = "Fitted", y = "Residual") +
  theme_minimal(base_size = 11)

p2 <- ggplot(diag_df, aes(sample = residual)) +
  stat_qq(alpha = 0.3, size = 1, colour = "steelblue") +
  stat_qq_line(colour = "firebrick") +
  labs(title = "Normal Q-Q", x = "Theoretical", y = "Sample") +
  theme_minimal(base_size = 11)

p1 + p2

7 +/-3 SD summary by epoch

model_dat |>
  group_by(epoch) |>
  group_modify(~ sd3(.x$mag)) |>
  pivot_wider(names_from = epoch, values_from = value) |>
  kable(digits = 3, caption = "Magnitude +/-3 SD summary by epoch")

Magnitude +/-3 SD summary by epoch
label	pre_2022	dimming_2023on
-3sd	10.537	14.344
-2sd	11.867	15.210
-1sd	13.197	16.076
Mean	14.527	16.942
+1sd	15.857	17.808
+2sd	17.187	18.674
+3sd	18.517	19.540

8 Sensitivity check — Johnson V band only

The pooled models above adjust for band with a random intercept, but they still pool the residual scatter across bands whose composition shifts over time (recent years are CR-heavy). To check whether the variability conclusion is an artifact of that shift, we repeat the analysis on the single most homogeneous, well-populated passband — Johnson V — where a like-for-like comparison needs no band effect at all.

v_dat <- model_dat |> filter(band == "V") |> droplevels()

v_dat |>
  group_by(epoch) |>
  summarise(n = n(), mean_mag = mean(mag), sd_mag = sd(mag), .groups = "drop") |>
  kable(digits = 3, caption = "V band only: brightness and variability by epoch")

V band only: brightness and variability by epoch
epoch	n	mean_mag	sd_mag
pre_2022	2782	15.280	0.808
dimming_2023on	206	17.024	0.376

# Full V-band history (all dates, so the 2022 transition is visible too).
v_all <- dat |>
  filter(band == "V") |>
  mutate(epoch = case_when(
    date <  PRE_END      ~ "pre_2022",
    date >= DIMMING_FROM ~ "dimming_2023on",
    TRUE                 ~ "2022 (transition)"
  ))

# Per-epoch mean magnitude, drawn as horizontal levels to show the dimming.
v_levels <- v_all |>
  filter(epoch != "2022 (transition)") |>
  group_by(epoch) |>
  summarise(
    xmin = min(date), xmax = max(date), mean_mag = mean(mag), .groups = "drop"
  )

ggplot(v_all, aes(date, mag)) +
  annotate("rect", xmin = PRE_END, xmax = DIMMING_FROM,
           ymin = -Inf, ymax = Inf, alpha = 0.08, fill = "grey20") +
  geom_point(aes(colour = epoch), alpha = 0.5, size = 1.1) +
  geom_segment(
    data = v_levels,
    aes(x = xmin, xend = xmax, y = mean_mag, yend = mean_mag),
    colour = "black", linewidth = 0.9, linetype = "solid"
  ) +
  geom_vline(xintercept = c(PRE_END, DIMMING_FROM),
             linetype = "dashed", colour = "grey40") +
  scale_y_reverse() +                       # brighter (smaller mag) at the top
  scale_colour_manual(values = c(
    "pre_2022"          = "steelblue",
    "2022 (transition)" = "grey55",
    "dimming_2023on"    = "firebrick"
  )) +
  labs(
    title    = "3C 279 — Johnson V light curve",
    subtitle = "Black lines = per-epoch mean magnitude; shaded = 2022 transition (excluded)",
    x = NULL, y = "V magnitude (fainter →)", colour = "Epoch"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

The single-band view makes the shift unambiguous: the 2023-onward cloud (red) sits both lower (fainter) and tighter (less scattered) than the pre-2022 cloud (blue), with the per-epoch mean levels showing the ~1.7 mag drop directly.

With a single band there is no random effect to include, so we fit a location–scale model with epoch in both the mean and the dispersion:

m_v <- glmmTMB(
  mag ~ epoch,
  dispformula = ~ epoch,
  data   = v_dat,
  family = gaussian()
)
summary(m_v)

 Family: gaussian  ( identity )
Formula:          mag ~ epoch
Dispersion:           ~epoch
Data: v_dat

      AIC       BIC    logLik -2*log(L)  df.resid 
   6899.3    6923.3   -3445.6    6891.3      2984 


Conditional model:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         15.27997    0.01533   997.0   <2e-16 ***
epochdimming_2023on  1.74380    0.03027    57.6   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dispersion model:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -0.21278    0.01341  -15.87   <2e-16 ***
epochdimming_2023on -0.76877    0.05106  -15.06   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Mean shift and per-epoch residual SD, back-transformed.
dv <- fixef(m_v)$disp
tibble(
  quantity = c("mean shift (mag, + = fainter)",
               "residual SD pre_2022",
               "residual SD dimming_2023on"),
  value = c(
    fixef(m_v)$cond[["epochdimming_2023on"]],
    sqrt(exp(dv[["(Intercept)"]])),
    sqrt(exp(dv[["(Intercept)"]] + dv[["epochdimming_2023on"]]))
  )
) |>
  kable(digits = 3, caption = "V-only: mean and variability change")

V-only: mean and variability change
quantity	value
mean shift (mag, + = fainter)	1.744
residual SD pre_2022	0.899
residual SD dimming_2023on	0.612

var.test(mag ~ epoch, data = v_dat) |>
  broom::tidy() |>
  kable(digits = 4,
        caption = "V-only F test for equal variances between epochs")

V-only F test for equal variances between epochs
estimate	num.df	den.df	statistic	p.value	conf.low	conf.high	method	alternative
4.6322	2781	205	4.6322	0	3.7552	5.6141	F test to compare two variances	two.sided

Interpretation. Restricted to a single band, the dimming (~1.7 mag fainter) persists and the variability drop that the pooled model washed out reappears clearly (residual SD roughly halves). This tells us the flat residual SD in the pooled model was driven by changing band composition, not by the source: measured like-for-like, 3C 279 has become both fainter and quieter since 2023.

9 Summary

Data: AAVSO magnitudes for 3C 279, cleaned to trusted photometric bands (V, CR, B, R, I, Vis.), spanning ~1999 to present.
Mean brightness: Model 1 estimates the average magnitude shift from the pre-2022 era to the 2023-onward era (see the epoch coefficient — positive = dimmer). The object is clearly fainter post-2023.
Variability: Model 2 (location–scale GLMM) estimates a separate spread for each epoch. On the pooled data the residual SD looks nearly unchanged, but that is a band-composition artifact — see the caveat below.
Sensitivity check (Section 8): restricting to Johnson V alone removes the band-mix confound and shows the variability genuinely dropped after 2023 (residual SD roughly halves), alongside the dimming. The like-for-like conclusion is: 3C 279 became both fainter and quieter.
Caveats: epochs are unbalanced across bands (e.g. CR dominates recent years); the (1 | band) term adjusts for filter zero-points but not for changing band composition over time — which is exactly why the single-band check matters. A useful next step is adding an observer effect if the raw AAVSO download includes an observer code. ```