text:
# Introduction to R for Excel Users

As the saying goes, when all you have is a hammer, everything looks like a nail. Excel was designed to do simple financial analyses and to craft financial statements. Though its capabilities have been expanded over the years, it was never designed to perform the sort of data analysis that industry scientists, engineers and Six Sigma belts need to perform on a daily basis.

Most data analyses performed in Excel look more like simple financial spreadsheets rather than actual data analysis, and this quality of work translates into bad—or at least sub-optimal—business decisions. There are alternatives to Excel, and the free, open-source data analysis platform R is one of them.

Unfortunately, R has a steep learning curve. I’m offering, for free, a short primer on R [PDF] where I’ve sought to make that learning curve a little less painful for engineers and scientists who normally work in Excel.

### Background

A couple of years ago, I was developing a short course to teach R to scientists and engineers in industry who normally used Excel. The goal was to help them transition to a more capable tool. My course design notes morphed into a handout, and when plans for the course fell through, that handout grew into a self-study guide, which I later adapted into this seventy-page, stand-alone introduction for Excel users.

### Organization

The primer walks the reader through the basics of R, starting with a brief overview of capabilities, then diving into installation, basic operations, graphical analysis and basic statistics. I believe that a picture is worth a thousand words, so it’s light on text and heavy on examples and visuals.

The end of the book rounds out with a look at some of the most useful add-ons, the briefest of introductions to writing your own, custom functions in R, and a cross-reference of common Excel functions with their equivalents in R.

The text is broken up into chapters and fully indexed so that it can be used either as a walk-through tutorial or as a quick reference.

# Can We do Better than R-squared?

If you're anything like me, you've used Excel to plot data, then used the built-in “add fitted line” feature to overlay a fitted line to show the trend, and displayed the “goodness of fit,” the r-squared (R2) value, on the chart by checking the provided box in the chart dialog.

The R2 calculated in Excel is often used as a measure of how well a model explains a response variable, so that “R2 = 0.8” is interpreted as “80% of the variation in the 'y' variable is explained by my model.” I think that the ease with which the R2 value can be calculated and added to a plot is one of the reasons for its popularity.

There's a hidden trap, though. R2 will increase as you add terms to a model, even if those terms offer no real explanatory power. By using the R2 that Excel so helpfully provides, we can fool ourselves into believing that a model is better than it is.

Below I'll demonstrate this and show an alternative that can be implemented easily in R.

### Some data to work with

First, let's create a simple, random data set, with factors a, b, c and response variable y.

head(my.df)

##       y a       b      c
## 1 2.189 1 -1.2935 -0.126
## 2 3.912 2 -0.4662  1.623
## 3 4.886 3  0.1338  2.865
## 4 5.121 4  1.2945  4.692
## 5 4.917 5  0.1178  5.102
## 6 4.745 6  0.4045  5.936

Here is what this data looks like:

### Calculating R-squared

What Excel does when it displays the R2 is create a linear least-squares model, which in R looks something like:

my.lm <- lm(y ~ a + b + c, data = my.df)

Excel also does this when we call RSQ() in a worksheet. In fact, we can do this explicitly in Excel using the Regression analysis option in the Analysis Pack add-on, but I don't know many people who use this, and Excel isn't known for its reliability in producing good output from the Analysis Pack.

In R, we can obtain R2 via the summary() function on a linear model.

summary(my.lm)

##
## Call:
## lm(formula = y ~ a + b + c, data = my.df)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -1.2790 -0.6006  0.0473  0.5177  1.5299
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    2.080      0.763    2.72    0.034 *
## a             -0.337      0.776   -0.43    0.679
## b             -0.489      0.707   -0.69    0.515
## c              1.038      0.817    1.27    0.250
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.1 on 6 degrees of freedom
## Multiple R-squared:  0.833,  Adjusted R-squared:  0.75
## F-statistic:   10 on 3 and 6 DF,  p-value: 0.00948

Since summary() produces a list object as output, we can grab just the R2 value.

summary(my.lm)$r.squared  ## [1] 0.8333  Normally, we would (somewhat loosely) interpret this as telling us that about 83% of the variation in the response y is explained by the model. Notice that there is also an "adjusted r-squared” value given by summary(). This tells us that only 75% of the variation is explained by the model. Which is right? ### The problem with R-squared Models that have many terms will always give higher R2 values, just because more terms will slightly improve the model fit to the given data. The unadjusted R2 is wrong. The calculation for adjusted R2 is intended to partially compensate for that “overfit,” so it's better. It's nice that R shows us both values, and a pity that Excel won't show the adjusted value. The only way to get an adjusted R2 in Excel is to run the Regression analysis; otherwise, we have to calculate adjusted R2 manually. Both R2 and adjusted R2 are measures of how well the model explains the given data. However, in industry we usually want to know something a little different. We don't build regression models to explain only the data we have; we build them to think about future results. We want R2 to tell us how well the model predicts the future. That is, we want a predictive R2. Minitab has added the ability to calculate predictive R2 in Minitab 17, and has a nice blog post explaining this statistic. ### Calcuting predictive R-squared Neither R nor Excel provide a means of calculating the predictive R2 within the default functions. While some free R add-on packages provide this ability (DAAG, at least), we can easily do it ourselves. We'll need a linear model, created with lm(), for the residuals so we can calculate the “PRESS” statistic, and then we need the sum of squares of the terms so we can calculate a predictive R2. Since the predictive R2 depends entirely on the PRESS statistic, we could skip the added work of calculating predictive R2 and just use PRESS, as some authors advocate. The lower the PRESS, the better the model is at fitting future data from the same process, so we can use PRESS to compare different models. Personally, I'm used to thinking in terms of R2, and I like having the ability to compare to the old R2 statistic that I'm familiar with. To calculate PRESS, first we calculate the predictive residuals, then take the sum of squares (thanks to (Walker’s helpful blog post) for this). This is pretty easy if we already have a linear model. It would take a little more work in Excel. pr <- residuals(my.lm)/(1 - lm.influence(my.lm)$hat)
PRESS <- sum(pr^2)
PRESS

## [1] 19.9

The predictive R2 is then (from a helpful comment by Ibanescu on LikedIn) the PRESS divided by the total sum of squares, subtracted from one. The total sum of squares can be calculated directly as the sum of the squared residuals, or obtained by summing over Sum Sq from an anova() on our linear model. I prefer using the anova function, as any statistical subtleties are more likely to be properly accounted for there than in my simple code.

# anova to calculate residual sum of squares
my.anova <- anova(my.lm)
tss <- sum(my.anova$"Sum Sq") # predictive R^2 pred.r.squared <- 1 - PRESS/(tss) pred.r.squared  ## [1] 0.5401  You'll notice that this is smaller than the residual R2, which is itself smaller than the basic R2. This is the point of the exercise. We don't want to fool ourselves into thinking we have a better model than we actually do. One way to think of this is that 29% (83% – 54%) of the model is explained by too many factors and random correlations, which we would have attributed to our model if we were just using Excel's built-in function. When the model is good and has few terms, the differences are small. For example, working through the examples in Mitsa's two posts, we see that for her model 3, R2 = 0.96 and the predictive R2 = 0.94, so calculating the predictive R2 wasn't really worth the extra effort for that model. Unfortunately, we can't know, in advance, which models are “good.” For Mitsa's model 1 we have R2 = 0.95 and predictive R2 = 0.32. Even the adjusted R2 looks pretty good for model 1, at 0.94, but we see from the predictive R2 that our model is not very useful. This is the sort of thing we need to know to make correct decisions. ### Automating In R, we can easily wrap these in functions that we can source() and call directly, reducing the typing. Just create a linear model with lm() (or an equivalent) and pass that to either function. Note that pred_r_squared() calls PRESS(), so both functions have to be sourced. pred_r_squared <- function(linear.model) { lm.anova <- anova(linear.model) tss <- sum(lm.anova$"Sum Sq")
# predictive R^2
pred.r.squared <- 1 - PRESS(linear.model)/(tss)
return(pred.r.squared)
}

PRESS <- function(linear.model) {
pr <- residuals(linear.model)/(1 - lm.influence(linear.model)$hat) PRESS <- sum(pr^2) return(PRESS) }  Then we just call the function to get the result: pred.r.squared <- pred_r_squared(my.lm) pred.r.squared  ## [1] 0.5401  I've posted these as Gists on GitHub, with extra comments, so you can copy and paste from here or go branch or copy them there. ### References and further reading # Normality and Testing for Normality Many of our statistical tests make assumptions about the distribution of the underlying population. Many of the most common—ImR (XmR) and XbarR control charts, ANOVA, t-tests—assume normal distributions in the underlying population (or normal distributions in the residuals, in the case of ANOVA), and we’re often told that we must carefully check the assumptions. At the same time, there’s a lot of conflicting advice about how to test for normality. There are the statistical tests for normality, such as Shapiro-Wilk or Anderson-Darling. There’s the “fat pencil” test, where we just eye-ball the distribution and use our best judgement. We could even use control charts, as they’re designed to detect deviations from the expected distribution. We are discouraged from using the “fat pencil” because it will result in a lot of variation from person to person. We’re often told not to rely too heavily on the statistical tests because they are not sensitive with small sample sizes and too sensitive to the tails. In industrial settings, our data is often messy, and the tails are likely to be the least reliable portion of our data. I’d like to explore what the above objections really look like. I’ll use R to generate some fake data based on the normal distribution and the t distribution, and compare the frequency of p-values obtained from the Shapiro-Wilk test for normality. ### A Function to test normality many times First, we need to load our libraries library(ggplot2) library(reshape2)  To make this easy to run, I’ll create a function to perform a large number of normality tests (Shapiro-Wilk) for sample sizes n = 5, 10 and 1000, all drawn from the same data: #' @name assign_vector #' @param data A vector of data to perform the t-test on. #' @param n An integer indicating the number of t-tests to perform. Default is 1000 #' @return A data frame in "tall" format assign_vector <- function(data, n = 1000) { # replicate the call to shapiro.test n times to build up a vector of p-values p.5 <- replicate(n=n, expr=shapiro.test(sample(my.data, 5, replace=TRUE))$p.value)
p.10 <- replicate(n=n, expr=shapiro.test(sample(my.data, 10, replace=TRUE))$p.value) p.1000 <- replicate(n=n, expr=shapiro.test(sample(my.data, 1000, replace=TRUE))$p.value)
#' Combine the data into a data frame,
#' one column for each number of samples tested.
p.df <- cbind(p.5, p.10, p.1000)
p.df <- as.data.frame(p.df)
colnames(p.df) <- c("5 samples","10 samples","1000 samples")
#' Put the data in "tall" format, one column for number of samples
#' and one column for the p-value.
p.df.m <- melt(p.df)
#' Make sure the levels are sorted correctly.
p.df.m <- transform(p.df.m, variable = factor(variable, levels = c("5 samples","10 samples","1000 samples")))
return(p.df.m)
}

### Clean, random data

I want to simulate real-word conditions, where we have an underlying population from which we sample a limited number of times. To start, I’ll generate 100000 values from a normal distribution. To keep runtimes low I’ll have assign_vector() sample from that distribution when performing the test for normality.

n.rand <- 100000
n.test <- 10000
my.data <- rnorm(n.rand)
p.df.m <- assign_vector(my.data, n = n.test)

We would expect that normally distributed random data will have an equal probability of any given p-value. i.e. 5% of the time we’ll see p-value ≤ 0.05, 5% of the time we’ll see p-value > 0.05 and ≤ 0.10, and so on through > 0.95 and ≤ 1.00. Let’s graph that and see what we get for each sample size:

ggplot(p.df.m, aes(x = value)) +
geom_histogram(binwidth = 1/10) +
facet_grid(facets=variable ~ ., scales="free_y") +
xlim(0,1) +
ylab("Count of p-values") +
xlab("p-values") +
theme(text = element_text(size = 16))

Histogram of p-values for the normal distribution, for sample sizes 5, 10 and 1000.

This is, indeed, what we expected.

Now let’s compare the normal distribution to a t distribution. The t distribution would pass the “fat pencil” test—it looks normal to the eye:

ggplot(NULL, aes(x=x, colour = distribution)) +
stat_function(fun=dnorm, data = data.frame(x = c(-6,6), distribution = factor(1)), size = 1) +
stat_function(fun=dt, args = list( df = 20), data = data.frame(x = c(-6,6), distribution = factor(2)), linetype = "dashed", size = 1) +
scale_colour_manual(values = c("blue","red"), labels = c("Normal","T-Distribution")) +
theme(text = element_text(size = 12),
legend.position = c(0.85, 0.75)) +
xlim(-4, 4) +
xlab(NULL) +
ylab(NULL)

Starting with random data generated from the t-distribution:

my.data <- rt(n.rand, df = 20)

Histogram of p-values for the t distribution, for sample sizes 5, 10 and 1000.

The tests for normality are not very sensitive for small sample sizes, and are much more sensitive for large sample sizes. Even with a sample size of 1000, the data from a t distribution only fails the test for normality about 50% of the time (add up the frequencies for p-value > 0.05 to see this).

### Testing the tails

Since the t distribution is narrower in the middle range and has longer tails than the normal distribution, the normality test might be failing because the entire distribution doesn’t look quite normal; we haven’t learned anything specifically about the tails.

To test the tails, we can construct a data set that uses the t distribution for the middle 99% of the data, and the normal distribution for the tails.

my.data <- rt(n.rand, df = 20)
my.data.2 <- rnorm(n.rand)
# Trim off the tails
my.data <- my.data[which(my.data < 3 & my.data > -3)]
# Add in tails from the other distribution
my.data <- c(my.data, my.data.2[which(my.data.2 < -3 | my.data.2 > 3)])

Histogram of p-values for sample sizes 5, 10 and 1000, from a data set constructed from the t distribution in the range -3 to +3 sigmas, with tails from the normal distribution below -3 and above +3.

Despite 99% of the data being from the t distribution, this is almost identical to our test with data from just the normal distribution. It looks like the tails may be having a larger impact on the normality test than rest of the data

Now let’s flip this around: data that is 99% normally-distributed, but using the t distribution in the extreme tails.

my.data <- rnorm(n.rand)
my.data.2 <- rt(n.rand, df = 20)
# Trim off the tails
my.data <- my.data[which(my.data < 3 & my.data > -3)]
# Add in tails from the other distribution
my.data <- c(my.data, my.data.2[which(my.data.2 < -3 | my.data.2 > 3)])

Histogram of p-values for sample sizes 5, 10 and 1000, from a data set constructed from the normal distribution in the range -3 to +3 sigmas, with tails from the t-distribution below -3 and above +3.

Here, 99% of the data is from the normal distribution, yet the normality test looks almost the same as the normality test for just the t-distribution. If you check the y-axis scales carefully, you’ll see that the chance of getting p-value ≤ 0.05 is a bit lower here than for the t distribution.

To make the point further, suppose we have highly skewed data:

my.data <- rlnorm(n.rand, 0, 0.4)

This looks like:

For small sample sizes, even this is likely to pass a test for normality:

### What have we learned?

• With small sample sizes, everything looks normal.
• The normality tests are, indeed, very sensitive to what goes on in the extreme tails.

In other words, if we have enough data to fail a normality test, we always will because our real-world data won’t be clean enough. If we don’t have enough data to reliably fail a normality test, then there’s no point in performing the test, and we have to rely on the fat pencil test or our own understanding of the underlying processes.

Don’t get too hung up on whether your data is normally distributed or not. When evaluating and summarizing data, rely mainly on your brain and use the statistics only to catch really big errors in judgement. When attempting to make predictions about future performance, e.g. calculating Cpk or simulating a process, recognize the opportunities for errors in judgment and explicitly state you assumptions.

# Scientific Consensus

I have seen a number of comments lately that “consensus” has no place in science, and that claims that there is a “scientific consensus” are just thinly veiled political double-speak. I have to take issue with such criticisms. In fact, consensus is one of the cornerstones of science.

This is not quite the same consensus used in politics or everyday life. Consensus means “general agreement” and is usually achieved through some form of discussion and negotiation. Consensus is therefore agreement over opinions, and is often agreement over a course of action despite differing opinions.

In science, consensus is derived from data and independent replication of experiments. It is the consensus that an idea—specifically, a testable hypothesis—is correct. It is the expression of scientists that a hypothesis is (a) scientifically testable and falsifiable, (b) that it has not been falsified and (c) that it explains the universe better than competing hypotheses. It is a consensus derived from the replication of observations or tests by other researchers.

Part of the problem, here, is that science has a habit of taking everyday words and developing very specific meanings around them. This happens because scientists need to communicate clearly and exactly at times, and language is messy and full of fuzzy concepts. The same thing happens in a lot of occupations. For example, accountants also develop specific meanings for everyday words.

Scientists no longer argue over the validity of Newton’s hypothesis on the gravitational force because there is broad consensus that the hypothesis is correct (as far as it goes). Objects attract each other according to their mass and the distance between them, and there have been plenty of independent experiments confirming the specific relationship, F = -G M1m2 / r2. The consensus is so strong that it’s referred to as Newton’s law of gravity. Likewise, scientists no longer argue over the geocentric model of the universe because there is broad consensus—derived from data collected over centuries by many independent researchers—that the Earth is not at the center of the solar system, let alone the universe.

Conversely, there is very little consensus when it comes to the accelerating expansion of the universe. Cosmologists agree that the universe is expanding faster than can be explained by our current understanding of the universe, but there are many, conflicting hypotheses about the causes. There is consensus over the fact of the accelerating expansion but the data does not yet support consensus on the underlying physical processes or mechanism causing it, and so there is no consensus about the physical process.

In one sense, scientific consensus is stronger, or more robust, than we are used to thinking about in politics and everyday life, precisely because it is based on observation and careful analysis by independent groups. It’s not just consensus based on what we think might be true, or what we want to be true, but based on what careful observation tells us must be true. Scientific consensus is not subject to whim.

From another perspective, though, scientific consensus is much weaker than we are used to. In science, there is no downside to abandoning or overturning a consensus when the data points in a different direction. In fact, there is significant benefit to being the person who can overturn a previous consensus; we remember Galileo, Newton, Darwin, Einstein and others precisely because their work, collecting and analyzing data, was so pivotal in altering the scientific consensus. In everyday life, if you back out of a consensus agreement, it’s likely that others party to the agreement will feel betrayed. There can be a significant social cost to pay for backing out of a consensus, even when you are convinced that you are right. Scientists may sometimes feel this same social pressure, but the scientific method provides clear guidance for adopting or abandoning consensus, and it doesn’t focus on the people involved but rather on external, objective observations of how the universe works. Scientific consensus is, perhaps, more readily changed than conventional consensus.

So consensus does exist is science and it plays an important role in science. We should be careful to distinguish, though, between consensus based on independent replication of results and consensus based on preconceptions and social negotiation.

# And You Thought Physics Was *YAWN*

Part of my day job involves monitoring the renewable energy market, and particularly keeping abreast of storage technologies. It’s like combining my hobby with my job.

A new wind turbine design was recently announced, the SeaTwirl. It’s an off-shore turbine design using vertical blades. The key technological advance here is that it includes a method to store energy, so that it can continue to produce electricity when the wind stops blowing. This ability to deliver a constant output is important, because, as you may have heard, wind energy is intermittent; you only get electricity when the wind blows, and only to the degree that it’s blowing. Demand, unfortunately, doesn’t follow wind’s intermittancy—nobody stops to check that the wind is blowing before they turn on their lights—and the utilities, transmission system operators and distribution system operators all have to supply electricity to meet demand.

The SeaTwirl stores energy by integrating an unusual wind turbine design with a pumped hydro system. It has the turbine blades on a large circular ring, which rotates parallel to the water’s surface kind of like a hula-hoop. This ring is hollow (more or less), and is filled up with water when the wind blows. When the wind stops blowing, the momentum of the water keeps the tube spinning, generating electricity, and the water is allowed to drain back out, spinning a hydro turbine to generate electricity.

The nice thing about this is that the storage can always be “recharged” and it’s “free.” Or at least it seems to be. If you’ve ever held a bucket while spinning around, you know that spinning with an empty bucket takes a lot less effort than spinning around with a full bucket. In part, this is because of a property known as the moment of inertia. The heavier or larger a spinning object gets, the more it resists changes to its rate of rotation (or rpm).

If the SeaTwirl is filling up this horizontal “hula-hoop” with water, then the weight of the tube is increasing and so is the moment of inertia. As the moment of inertia increases, the energy needed to reach a given rpm increases. Wind turbines normally generate electricity in proportion to the wind speed, because the rpm of the blades is proportional to the wind speed. Increase the moment of inertia and you decrease the rpm, which means you generate less electricity for a given wind speed.

Now comes the physics. For something shaped like a hula-hoop, the moment of inertia, I, is calculated from the mass, M, and the radius of the hoop, R, according to:

$I = M R^{2}$

The energy, E, in a spinning object is equal to the moment of inertia times the speed of rotation, ω, according to

$E=\frac{1}{2} I \omega ^{2}$

If we know the energy (because we know the wind speed), then we can calculate the speed of rotation, ω, by rearranging that equation to get ω on the left-hand side:

$\omega = \sqrt{\frac{2E}{I}}$

We can then replace I with mass and radius from the first equation to get

$\omega = \sqrt{\frac{2E}{MR^{2}}}$

So we can see that, if we don’t change the energy E (or don’t change the wind speed), and don’t change the radius R of the spinning hoop, then increasing the mass M results in a slower rate of rotation.

From SeaTwirl’s website and press releases, we can estimate how big the SeaTwirl is, which will let us estimate how much slower a full SeaTwirl will spin than an empty one, and therefore how much less electricity must be generated. We can calculate this by taking the ratio of ω full to ω empty, so that the parts that we don’t have to know E and R.

$\frac{\omega_{full}}{\omega_{empty}} = \frac{\sqrt{\frac{2E}{M_{full}R^{2}}}}{\sqrt{\frac{2E}{M_{empty}R^{2}}}} = \sqrt{\frac{\frac{2E}{R^{2}}}{\frac{2E}{R^{2}}}}\sqrt{\frac{M_{empty}}{M_{full}}}=\sqrt{\frac{M_{empty}}{M_{full}}}$

For the SeaTwirl, we now have to find out what R is, and estimate M for both filled and empty

The whole turbine assembly is made of composite materials, which probably have a density, $\rho_{c}$, of around 2500 kilograms per cubic meter (similar to fiberglass). Water has a density, $\rho_{w}$, of near 1000 kilograms per cubit meter (depending on temperature). The diameter of the turbine will be near 180 meters, so the radius, R, of our “hula-hoop” is half that, or 90 meters. From the pictures, it looks like the thickness of that hula-hoop is a few percent of the total diameter of the turbine, so we can figure an outside diameter of the “hula-hoop” of about 2 meters, for a radius, r, of 1 meter. Figure that at least ten percent of this is composite, and the rest is the hollow, water-filled portion.

To estimate the weight of the water in the “hula-hoop,” we can approximate the water as being a cylinder of radius $r_{w} = 0.9r$ and length equal to the circumference of the “hula-hoop,” $l = 2\pi R$. The volume of such a cylinder is equal to the cross-sectional area of the water column, $A_{w}=\pi r_{w}^{2}$ times the length of the column, l. The total mass of the water, $m_{w}$ is the density times this volume.

$m_{w} = \rho_{w}\pi r_{w}^{2}2\pi R = \rho_{w}3\pi (0.9r)^{2} R$

Plugging in our estimates for the above values gives us

$m_{w} = 1000 \cdot 3\pi (0.9)^{2} 90 = 687000$ kg

That’s a lot of water.

Now for the empty “hula-hoop.” We can treat it in the same way: a cylinder of material of radius r, length $l = 2\pi R$. However, we don’t want to calculate for a solid cylinder of composite; we have to subtract out the hollow part with radius $r_{w}$. So the mass of the composite is

$m_{c} = \rho_{c} 2 \pi R ( \pi r^{2} - \pi r^{2}_{w} )$

$m_{c} = 2500 \cdot 3 \pi 90 (1^{2} - 0.9^{2} ) = 403000$ kg

So the water more than doubles the weight of the hoop.

From the picture, you can see that there’s another hoop at the top, and the two hoops are connected by the turbines, which combined are probably worth at least another hoop in weight, so we can further assume that the mass of this bottom hoop, empty, is roughly one-third of the total mass of the movable parts of turbine.

The mass of the turbine, empty, is therefore about 1200000 kg, or 1200 tons. Filled with water, this goes up to about 1900000 kg, or 1900 tons. Empty, that’s maybe twice the largest off-shore turbine currently in existence, but this thing is easily twice as big as any current turbine, so our estimate appears to be in the right neighborhood.

Now we go back to our equation for the ratio of the rotational velocities, ω, and plug in these weights:

$\frac{\omega_{full}}{\omega_{empty}} = \sqrt{\frac{M_{empty}}{M_{full}}} =\sqrt{\frac{1200}{1900}} = 0.8$

So we get about 80% as much electricity from a water-filled turbine as from an empty one, when the wind blows. This is a direct efficiency loss due to the storage of energy in the spinning-water-hoop.

In addition, there’s the efficiency losses in loading water into the hoop, or “charging” the hoop, and the efficiency losses of “discharging” the hoop, running the water back out through a hydro turbine. Pumped hydro is usually about 72% efficient, or less, in each direction, so the total round-trip efficiency of storage + discharge is about 50% efficient. There are a lot of other storage technologies that do at least this good, if not better.

These two figures, the 80% efficiency loss of just operating the turbine and the 50% storage-discharge efficiency, can be used to directly compare SeaTwist with other wind turbine + storage technology solutions. Any storage technology that has at least a 50% round-trip efficiency and increases the total system cost by less than 20% over the system’s operating lifetime will outperform SeaTwist in terms of return on investment.

# R Function Reference

Updated below

The R Function Reference is a mind map that I created as a guide for novice and intermediate users of the R statistics language. When you first open it, I suggest that you collapse all the nodes by clicking on the “Expand/Collapse all nodes” button in the bottom left of the screen to make the map easier to navigate. You can also adjust the zoom level with the slider next to that button.

The top-level nodes of the R Function Reference

The mind map is arranged in eight sections, or main branches, arranged by task. What do you want to do? Each branch covers a general set tasks, such as learning to use R, running R, working with data, statistical analysis or plotting data. The end of each string of nodes is generally a function and example. The Reference provides code fragments, rather than details of the function or complete reproducible code blocks. Once you’ve followed the Reference and have an idea of how to accomplish something, you can look up the details in R’s help system (e.g. “?read.csv” to learn more about using the read.csv() function), or search Google or the online R-Help mailing list archives for answers using the function name.

There are a lot of useful nodes and examples, especially in the “Graphs” section, but the mind map is not complete; some trails end before you get to a useful function reference. I am sorry for that, but it’s a work in progress, and will be slowly updated over time.

### Update 1

In comments, several users reported problems opening the mind map. With a little investigation, it appears that the size of the mind map is the problem. To try to fix the problems , I have split the mind map out into several small mind maps, all linked together.

The new main mind map is the R Function Reference, Main. The larger branches on this main map no longer expand to their own content, but contain a link to a “child” mind map. The link looks like a sheet of paper with an arrow pointing to the right, click on it and little cartoon speech bubble will pop up with a link that you have to click on to go to the child mind map. Likewise, the central nodes on the child mind maps contain a link back to the main mind map.

Due to load times and the required extra clicks, this may slightly reduce usability for users who didn’t have a problem with the all-in-one version, but will hopefully make the mind map accessible to a broader audience.

I have to offer praise to the developers of Mind42. Though I couldn’t directly split branches off into their own mind maps or duplicate the mind map, it was very easy to export the mind map as a native Mind42 file and then import it multiple times, editing the copies without any loss of data or links. The ability to link directly between mind maps within Mind42 was also a key enabling feature. Considering that this is a free web app, its capabilities are most impressive. They were also quick to respond when I posted a call for help on the Mind42 forum.

Please let me know how the new, “improved” version works.

The old mind map, containing everything, is still available, but I will not update it.

# Graphing Highly Skewed Data

Recently Chandoo.org posted a question about how to graph data when you have a lot of small values and a few larger values. It’s not the first time that I’ve come across this question, and I’ve seen a lot of answers, many of them really bad. While all solutions involve trade-offs for understanding and interpreting graphs, some solutions are better than others.

Data graphs tell stories by revealing patterns in complex data. Good data graphs let the data tell the story by revealing the patterns, rather than trying to impose patterns on the data.

As William Cleveland discusses in The Elements of Graphing Data and his 1993 paper A Model for Studying Display Methods of Statistical Graphics, there are two basic visual operations that people employ when looking at and interpreting graphs: pattern perception and table look-up. Pattern perception is where we see the geometric patterns in a graph: groupings; relative differences (larger/smaller); or trends (straight/curved or increasing/decreasing). Table look-up is where we explore details of values and names on a graph. These two operations are distinct and complimentary, and it is through these two operations that the data’s story is told.

month sales
1 Feb 09 200
2 Mar 09 300
3 Apr 09 200
4 May 09 300
5 Jun 09 200
6 Jul 09 300
7 Aug 09 350
8 Sep 09 400
9 Oct 09 450
10 Nov 09 1200
11 Dec 09 100000
12 Jan 10 85000
13 Feb 10 450

So suppose that we have some data like that at right, where we are interested in the patterns of smaller, individual values, but there are also a few extremely large values, or outliers. We describe such data as being skewed. How do we plot this data? First, for such a small data set, a simple table is the best approach. People can see the numbers and interpret them, there aren’t too many numbers to make sense of and the table is very compact. For more complicated data sets, though, a graph is needed. There’s a few basic options:

• Graph as-is;
• Graph with a second axis;
• Graph the logarithm of the data;
• Use a scale break.
• Plot the data multiple times.

### Graph As-Is

A bar chart with all data, including outliers, plotted on the same scale.

This is the simplest solution, and if you’re only interested in knowing about the outliers (Dec ’09 and Jan ’10) then it will do. However, it completely hides whatever is happening in the rest of the months. Pattern recognition tells us that two months near the end of the series have the big numbers. Table-lookup tells us the approximate values and that these months are around December ’09 and February ’10, but the way the labels string together and overlap the tick marks, it’s not clear exactly what the labels are, let alone which label applies to which bar (which months are those, precisely? Is that “09 Dec” and “09 Feb?” Do the numbers even go with the text, or are they separate labels?).

For all but the simplest of messages, this rendition defeats both pattern recognition and table look-up. We definitely need a better solution.

### Use a Secondary Axis

Excel gives us an easy solution: break the data into two columns (“small” numbers in one and “large” numbers in the other) and plot them on separate axes. Now we can see all the data, including the patterns in all the months.

Bar chart, with outliers plotted using a secondary axis.

Unfortunately, pattern recognition tells us that the big-sales months are about the same as all the other months. It’s only the table look-up that tells us how big of a difference there is between the two blue columns and the rest of the data. This is why I’ve added data labels to the two columns: to aid table look-up.

Even if we tweaked around with the axes to set the outliers off from the rest of the data, we’d still have the same basic problem: pattern recognition would tell us that there is a much smaller difference than there actually is. By using a secondary axis, we’ve set up a basic conflict between pattern recognition and table look-up. Worse, it’s easy to confuse the axes; which bars go with which axis? Reproduction in black and white or grayscale would make it impossible to correctly connect bars to the correct axis. Some types of color blindness would similarly make it difficult to interpret the graph. Table look-up is easily defeated with secondary axes.

The secondary axis presents so many problems that I always advise against using it. Stephen Few, author of Show Me The Numbers and Information Dashboard Design, calls graphs with secondary axes “dual-scaled graphs.” In his 2008 article Dual-Scaled Axes in Graphs, he concludes that there is always a better way to display data than by using secondary axes. Excel makes it easy to create graphs like this, but it’s always a bad idea.

### Take the Logarithm

In scientific applications, skewed data is common, and the usual solution is to plot the logarithm of the values.

Bar chart plotting skewed with logarithmic axis.

With the logarithm, it is easy to plot, and see, all of the data. Trends in small values are not hidden. Pattern perception immediately tells us the overall story of the data. Table look-up is easier than with secondary axes, and immediately tells us the scale of the differences. Plotting the logarithm allows pattern perception and table look-up to compliment each other.

Below, I’ve created the same graph using a dot plot instead of a bar chart. Dot plots have many advantages over bar charts: most obviously, dot plots provide a better arrangement for category labels (e.g. the months); also, dot plots provide a clearer view of the data by plotting the data points rather than filling in the space between the axis and the data point. There are some nice introductions to dot plots, including William Cleveland’s works and a short introduction by Naomi Robbins. The message is clear: any data that you might present with a bar chart (or pie chart) will be better presented using dot plots.

Skewed data plotted on a dot plot using a logarithmic scale.

### Use a Scale Break

Another approach, which might be better for audiences unfamiliar with logarithmic scales, is to use a scale break, or broken axis. With some work, we can create a scale break in Excel or OpenOffice.org.

Bar chart with outliers plotted by introducing a subtle scale break on the y-axis.

There are plenty of tutorials for how to accomplish this in Excel. For this example, I created the graph in OpenOffice.org Spreadsheet, using the same graph with the secondary axis, above. I adjusted the two scales, turned off the labels for both y-axes and turned off the tick marks for the secondary y-axis. Then I copied the graph over to the OpenOffice.org Draw application and added y-axis labels and the break marks as drawing objects.

That pretty much highlights the first problem with this approach: it takes a lot of work. The second problem is that those break marks are just too subtle; people will miss them.

The bigger problem is with interpretation. As with the secondary axis, this subtle scale break sets up a basic conflict between the two basic operations of graph interpretation. Pattern recognition tells us that the numbers are comparable; it’s only table look-up that tells us what a large difference there is.

Cleveland’s recommendation, when the logarithm won’t work, is to use a full-panel scale break. In this way, pattern recognition tells that there are two distinct groups of data, and table look-up tells us what they are.

Dot plot with a full scale break to show outliers.

The potential disadvantage of this approach is that pattern perception might be fooled. While the scale break visually groups the “large” values from the “small” ones, the scale also changes, so that the broader panel on the left actually represents a much narrower range of values (about 1100 dollars range) than the narrower panel on the right (about 17000 dollars range). Our audience might have difficulties interpreting this correctly.

### Small Multiples

Edward Tufte has popularized the idea of small multiples, the emphasis of differences by repeating a graph or image with small changes from one frame to the next. In this case, we could show the full data set, losing fidelity in the smaller values, and then repeat the graph while progressively zooming in on a narrower and narrower slice with each repetition.

The full data, with outliers, is plotted on the left. On the right, a zoomed view showing detail in the smaller values.

This shares many similarities to Cleveland’s full scale break, but provides greater flexibility. With this data, there are two natural ranges: 0 – 100000 and 0 – 1200. If there were more data between 1200 and 85000, we might repeat the graph several times, zooming in more with each repetition to show lower levels of detail.

I think there are two potential pitfalls. As with the full scale break, the audience might fail to appreciate the effect of the changes to scale. Worse, the audience might be fooled into thinking that each graph represented a different set of data, rather than just a different slice of the same data. Some care  in preparing such graphs will be needed for successful communication.

### Summary

When presenting data that is, like the data above, arranged by category, use a dot plot instead of bar charts. When your data is heavily skewed, the best solution is to graph the logarithm of the data. However, if your audience will be unable to correctly interpret the logarithm, try a full scale break or small multiples.


text:
Example of algebraic field not honoring the Total Order relation?

I've started studying Calculus. I've stumbled upon the definition of an 'Ordered field'. One of the requirements are for the field to honor the Total Order relation. Meaning that for every a,b, a could either be:

• a > b
• a < b
• a = b

While this sounds trivial, I'm trying for intuition purposes to find an counter example of a field not honoring the Total Order relation. Can you give me an example?

Thanks

-

$\Bbb C$, the field of complex numbers. If $\Bbb C$ were an ordered field, ordered by some linear order $\le$, the set $P=\{z\in\Bbb C:z>0\}$ would have the following properties:

1. For every $z\in\Bbb C$, exactly one of the following is true: $z\in P$, $-z\in P$, $z=0$.
2. $P$ is closed under multiplication: if $w,z\in P$, then $wz\in P$.

Now $i\ne 0$, so by (1) either $i\in P$ or $-i\in P$. If $i\in P$, then $-1=i^2\in P$ by (2). If, on the other hand, $-i\in P$, then $-1=(-i)^2\in P$. Thus, we can be sure that $-1\in P$. But then $1=(-1)^2\in P$, so both $-1$ and $1$ are in $P$, violating (1). Thus, there is no linear order that makes $\Bbb C$ an ordered field.

-
ah thanks! wasn't thinking of imaginary numbers. –  vondip Nov 5 '12 at 7:53
@vondip: You’re welcome! –  Brian M. Scott Nov 5 '12 at 7:55

It is rather immediate that if a field $K$ admits a total ordering, then $-1$ is not a sum of squares.

Indeed, if we have a total ordering on $K$, then $1>0$, thus $-1<0$. Also for every $a\in K$ we have $a^2>0$. Thus $-1\neq \sum_{i=1}^n a_i^2$, since one side is positive and the other is negative.

From this observation one immediately concludes that fields of positive characteristic (e.g. finite fields) do not admit total orderings. (A field $K$ has positive characteristic if there exists $n>0$ such that $\underbrace{1+\cdots+1}_{n\ \textrm{times}}=0$ in $K$.)

It also gives a way to construct fields, say, sub-fields of $\mathbb{C}$ that admit no total ordering.

For example the cyclotomic fields: Let $n>2$ and let $\zeta = e^{2\pi i/n}$ be a primitive $n$-th root of unity. Then $\mathbb{Q}(\zeta)$, the minimal sub-field of $\mathbb{C}$ that contains $\zeta$, admits no total ordering. To see this we separate between two cases. In the first case $n$ is even, and then $-1=\zeta^{n/2}$ is a square.

In the second case $n$ is odd. Then the sets $\{\zeta^{k} \mid k=0, \ldots, n-1\}$ and $\{\zeta^{2k} \mid k=0, \ldots,n-1\}$ coincide. So using the formula for the sum of terms in geometric series we get that

$\displaystyle\zeta^{n-1} + \zeta^{n-2} + \cdots + \zeta +1 = \frac{1-\zeta^n}{1-\zeta} =0,$

so

$\displaystyle -1=\zeta^{n-1} + \zeta^{n-2} + \cdots + \zeta = \zeta^{2(n-1)} + \zeta^{2(n-2)} + \cdots + \zeta^2 = (\zeta^{n-1})^2 + (\zeta^{n-2})^2 + \cdots + \zeta^2$

is a sum of squares.

-

In fact, a field $k$ can be ordered iff $-1$ is not a sum of squares.

If $-1$ is a sum of squares in an ordered field, then $-1$ is positive whence $1$ is negative; but $1$ is a square and so positive. It is a contradiction.

Let us suppose that $-1$ is not a sum of squares. Then set of field $K \supset k$ in which $-1$ is not a sum of squares is inductive; by Zorn's lemma, there exists a maximal field $\tilde{k}$ containing $k$ in which $-1$ is not a sum of squares.

Let $\gamma \in \tilde{k}$ not be a square. Then $\tilde{k}\subsetneq \tilde{k}(\sqrt{\gamma})$ where $\tilde{k}(\sqrt{\gamma})$ is a rupture field of $X^2-\gamma$. By definition of $\tilde{k}$, there exist $a_i$ and $b_i$ such that $-1= \sum\limits_{i=1}^n (a_i+ \sqrt{\gamma}b_i)= \sum\limits_{i=1}^n (a_i^2+ \gamma b_i^2) + 2 \sqrt{\gamma}\sum\limits_{i=1}^n a_ib_i$. Because $\sqrt{\gamma} \notin \tilde{k}$, $\sum\limits_{i=1}^na_ib_i=0$ and because $-1$ is not a sum of squares in $\tilde{k}$, $\gamma$ is not a sum of squares in $\tilde{k}$.

So, in $\tilde{k}$, a sum of squares is a square. Because $-1$ is not a sum of squares in $\tilde{k}$, $\sum\limits_{i=1}^n b_i^2 \neq 0$. Then we deduce that $\displaystyle - \gamma = \frac{1^2+ \sum\limits_{i=1}^n a_i^2}{\sum\limits_{i=1}^n b_i^2}$ is a square. So, for all $\gamma \in \tilde{k}$, $\gamma$ or $-\gamma$ is a square. Finally, define on $k$ the following order: $y \leq x$ iff $x-y$ is a square in $\tilde{k}$ and prove that $(k,\leq)$ is an ordered field.

For example, $\mathbb{C}$ and $\mathbb{F}_{p^n}$ cannot be ordered because $-1=i^2$ in $\mathbb{C}$ and $-1= \underset{p^n-1}{\underbrace{1^2+...+1^2}}$ in $\mathbb{F}_{p^n}$.

Reference: Algebra, van der Waerden.

-

format: plain


text:
# self-study

1. ### What is the smallest value of r that gives a test of size not exceeding 0.05?

0down votefavorite I've been given a question that I don't understand. Question: Consider a test which rejects H0 if the number of observed successes is ≥r≥r , where r is some number. What is the smallest value of r that gives a test of size not exceeding 0.05? We have been given: o=50 ,for...
2. ### CDF of a PDF with a negative power

I need to solve the following problem: For a marine hull insurance policy, a claim occurs with probability 0.DD, where DD represents the two date-digits of your birthday (DD.MM.YYYY). If a claim occurs, it has a probability density given by fB(x)=c·x^(−1.5), 1<x<YYYY, where YYYY are the...
3. ### Distribution of X+Y of a bivariate normally distributed (X,Y)

I need to give the distribution of $X+Y$ and $X-Y$ knowing that $(X,Y)$ is bivariate normally distributed with marginal means 1, marginal variances 1 and correlation $p=2$.Is it right that the marginal mean are simply $E(X)=1$ and $E(Y)=1$ ?I don't see where I should start knowing that.
4. ### Computer game log normally distributed

I'm playing a computer game and I want to boost the high score. I assume that attempts are independent (say, you do not improve by playing), and in each game my achieved points are lognormally distributed with mean 500 and variance 1’000’000. How often do I have to play so that your high score...
5. ### Sudokua problem using related binomial distribution

It is given: If you succeed to solve a Sudoku with probability 0.4 and you try many times, how likely is it that you succeed twice before you fail 5 times? What is the expected number of failures before two successes? I don't understand which one between the binomial, geometric, negative...
6. ### Resulting expected value of a bounded exponential pdf

I need to determine the distribution function of an exponential random variable with mean 2, given that its outcome is larger than 2 and give the resulting expected value? I know that the mean is 1/lambda, so lambda should be 1/2. I get that way 0.5*e^(-0.5x) >= 2. How can I deal with the...
7. ### Distribution function of an exponential random variable

Let Y ⇠ Exp(2) be the size of an insurance loss. If the insurance company pays min(Y, 4), what is the variance of that payment? I know that lambda =2 If I put that in the pdf I get: 2*e^(-2x). I don't understand how to use the information given by min(Y, 4)
8. ### Definition of random variable with respect to event space

It is given that: We toss two dice with sample space Ω = {(i, j), 1 ≤ i, j ≤ 6} and the σ-algebra is generated by the events Ak = {(i,j) : max(i,j) = k} (k = 1, . . . , 6). Show that whereas X1(i, j) = max(i, j) is a random variable in the corresponding probability space, X2(i, j) = i + j is...
9. ### Prove: Birth-death chain has stationary distribution -> satisfy detail bal. equations

I'm working my way through a course in stochastic processes but came across this Lemma in my coursebook where they didn't show the working (Lemma 8.5). I'm not familiar with doing proofs and it isn't taught in any of my courses but I can read/understand them, how did they reach the conclusion of...
10. ### Calculating Overall Relative Bias

I am in some trouble to understand how is to calculate the overall relative bias. In this link, there are results of overall relative bias in "Parameter estimates" sub-section under "Results" section. There they mentioned that : Could you please explain me how did they calculate it? I had...
11. ### Project Issue.

I am a student of statistics discipline. I have to conduct a project which i can apply in society or nation. My preferable topics are agricultural/technological/environmental statistics. But i am not understanding how can i specialize my topic, make it applicable and can have an ample online...
12. ### Obtaining sample variance from grouped data for goodness of fit test

This is a practice question I came across when dong some goodness of fit test examples.A company sells cloths by mail order.The size of clothes is defined by hip size; thus the height of customers of a particular customer may vary considerably. Data set of heights sent in by customers of size 18...
13. ### Proof of constistency of Maximum Likelihood Estimators (MLE)

Hi all, I would appreciate some help comprehending a logical step in the proof below about the consistency of MLE. It comes directly from Introduction to Mathematical Statistics by Hogg and Craig and it is slightly different than the standard intuitive one that makes use of the Weak Law of...
14. ### Need Help Choosing a Textbook

Hi, I'm looking for a textbook/book to let me self-study regression analysis. Some background: I have a research psych background academically (I've done many means-tests, chi-squareds, pearson's r, but always with random subjects/experiments) so I'm not totally new to the material, but...


text:
## Thomas' Calculus 13th Edition

Rewriting, $\displaystyle \frac{z^{2}}{1^{2}}+\frac{y^{2}}{1^{2}}-\frac{x^{2}}{1^{2}}=1\qquad$ this is of the form of a hyperboloid of one sheet, (compare to Table 12.1 and interchange x and z). It opens along the x-axis. In planes parallel to the yz plane ($x=k)$, the cross sections are circles. At $x=0$, the radius is 1 and at $x=\pm 2$, the radius is $\sqrt{5}$ (use for sketching). The cross sections with planes parallel to xy and xz are hyperbolas.


text:
## Monday, April 21, 2008

### [AD 11] Solution of the PL-Kurtosis puzzle - The SD model

The question if powerlaw diffusion is compatible with positive kurtosis and, in particular, with an exponential step width distribution, could finally be settled by directly constructing a discrete time series with the requested properties:

We consider a simple stochastic process, which produces a time series of values x_k=x(t_k) for discrete time points t_k = k*dt:

$x_n=x_{n-1}+s_n\;d_n$

where the s_n are sign factors

$s_n&space;\in&space;\left\{&space;-1,+1&space;\right\}$

with equal distribution function

$P(s)=\frac{1}{2}\delta_{s,-1}+\frac{1}{2}\delta_{s,+1}$

and a powerlaw auto-correlation:

$C_{ss}(m)=\left<&space;s_n\;&space;s_{n+m}&space;\right>=\delta_{m,0}\;\sigma_s^2&space;+&space;(1-\delta_{m,0})\;m^{-p}$

The d_n are positive step widths

$d_n\in&space;\left[&space;0,\infty&space;&space;\righ]$

with an arbitrary distribution function of finite variance (Gaussian, Exponential, Poisson...)

$P(d>0)&space;\;\;\mbox{arbitrary,&space;with}\;\;&space;\sigma_d^2=\left<&space;(d\!-\!\overline{d})^2&space;\right>&space;<\infty$

but without any temporal correlations (white noise):

$C_{dd}(m)=\left<&space;(d_n\!-\!\overline{d})\;(d_{n+m}\!-\!\overline{d})&space;\right>=$
$=\left<&space;d_n\;d_{n+m}&space;\right>&space;-\overline{d}^2&space;=\sigma_d^2\;\delta_{m,0}$

which means that

$\left<&space;d_n\;d_{n+m}&space;\right>&space;=&space;\overline{d}^2\;\;\mbox{for}\;\;m\neq&space;0$

We assume there are no cross-correlations between the sign-factors and the step widths:

$C_{sd}(m)=\left<&space;s_n\;d_{n+m}&space;\right>=0$

The velocity of the nth step is

$v_n=\frac{x_n-x_{n+1}}{\Delta&space;t}=\frac{s_n\;d_n}{\Delta&space;t}$

For the auto-correlation function of the velocity we obtain

$C_{vv}(m)=\left<&space;v_n\;v_{n+m}&space;\right>=\frac{1}{\Delta&space;t^2}\left<&space;s_n&space;d_n\;s_{n+m}d_{n+m}&space;\right>$

which can be factorized because the s and d and uncorrelated:

$C_{vv}(m)=\frac{1}{\Delta&space;t^2}\left<&space;s_n\;s_{n+m}&space;\right>\left<&space;d_n\;d_{n+m}&space;\right>$

which yields for all m except 0:

$C_{vv}(m\!&space;\neq&space;\!0)=\frac{\overline{d}^2}{\Delta&space;t^2}\;\;m^{-p}$

This means that the process has powerlaw VAC (and PSD and MSD) for arbitrary step widths distributions.

By using separate processes x_n and y_n of the above kind for the different spatial directions, the method can be generalized to multiple dimensions.


text:
## Introduction

Daily, lunar and seasonal cycles of natural light have been key forms of environmental information for organisms since shortly after the first emergence of life1. Artificial light at night (ALAN) has rapidly increased since the 19th century and has been introduced outdoors in places, at times, spectra and intensities at which it does not naturally occur at night2,3. Global simulations of skyglow (one form of light pollution) on clear nights demonstrate the dramatic extent to which artificial night sky brightness (NSB) has altered nightscapes worldwide4, with for example 88% of the EU experiencing light pollution on clear nights. We have recently shown that ALAN is growing in area and intensity by more than 2% per year on a global scale, with individual countries having increases of more than 10% per year5.

ALAN thus has become an anthropogenic stressor - light pollution, that can affect humans6, animals7,8, plants9 and microorganisms10. Recent findings include disruption of pollination11 or the disturbance of reproduction of mammals12. Despite the many studies regarding the impact of ALAN on the environment, a comprehensive quantitative understanding of the amount of artificial and natural light in the nocturnal environment for all weather conditions and how it is perceived by organisms is still lacking. In particular, while reference values for sky luminance and illuminance exists for clear nocturnal skies4, these reference points have not been defined for cloudy conditions. However, knowledge of such reference points is needed for the design of experiments and for ecological modeling to fully evaluate the environmental impacts of ALAN7,8.

Atmospheric conditions can change light levels rapidly during the day13 and the night14. In a natural setting with no ALAN, clouds should darken the night sky in most cases. For large distances around today’s urban areas however, clouds usually amplify the effects of ALAN, creating a complete reversal of natural conditions14,15,16,17,18. NSB data on rural sites with low ALAN levels, however, is sparse and in general few studies with imaging devices and clouds exist.

One reason for these knowledge gaps in light pollution research is due to the fragmentation of expertise across very different disciplines and the lack of widely applicable commercial measurement tools until recently. For decades, NSB monitoring was mainly of interest for astronomers as they were first to study light pollution19,20,21,22. For logical reasons, astronomers mainly concentrated their research on clear skies, but also performed comprehensive observations at several important (often very remote) astronomical sites, for example at high altitudes23,24,25 or latitudes26. The advent of small photometers customized for NSB measurements (mainly the Sky Quality Meter: SQM) has sparked NSB research beyond professional astronomy. Nowadays scientists from different fields, amateur astronomers, and citizen scientists obtain NSB data from local to global scale3,4,16,27. However, these single channel sensors usually provide zenith NSB data in a single spectral band, which provides only limited information of the nocturnal night light environment and might require cross-calibration28,29. Commercial digital cameras with fisheye lenses are a promising tool for all-sky nighttime photometry to fill this gap, yielding spatially resolved NSB data of nearly the full hemisphere in three spectral channels18,30,31,32,33. A more comprehensive discussion of state-of-the-art NSB measurement techniques can be found in a recent review34.

Recently, we reported on lowered zenith luminance values for cloudy nights using an SQM at Lake Stechlin, Germany35 while others observed both, brightening and darkening by clouds in SQM datasets at Montsec Astronomical Park, Spain17 and in Austria36. All-sky imaging data of cloud amplification were investigated in Madrid, Spain37 and near Montsec Astronomical Park, Spain18. However, imaging data on the darkening by clouds in the context of ALAN has not been published, yet. To set a lower limit for nighttime luminance and illuminance such data is essential as it provides information of the NSB across the full hemisphere.

In this work, we use differential all-sky photometry38 with a commercial DSLR camera to investigate the changes that clouds have upon the NSB of both the entire sky, and segments of the sky. We examine two rural locations: one 70 km north of Berlin, Germany (LakeLab, Lake Stechlin, Brandenburg) and one 180 km east of Cape Town, South Africa (Night Sky Caravan Farm, Bonnievale). The calculated all-sky luminance maps and the derived color corrected temperature (CCT) from the three spectral channels of the DSLR camera enable us to study the NSB distribution in detail. The spatially resolved data allows to identify and study the minimally polluted parts of the sky for locations that are not completely free of artificial skyglow. Our results fill a research gap in the context of ecological light pollution and we hope this work triggers more such spatially resolved NSB measurements on a global scale to increase understanding of the impacts of ALAN.

## Results

### Overcast sky at LakeLab, Lake Stechlin, Germany

The first observations were taken at the LakeLab, situated about 70 km north of the city of Berlin. Data was taken during a clear night on August 11th 2016 at 01:29 local time and the following overcast night on August 12th 2016 at 01:27 local time (+2 hours GMT). The cloud base height for the overcast night was inferred to be hCB ≈ 850 m ± 200 m. More information on the site and the surroundings is summarized in the methods section.

#### Full color images, LakeLab

The upper row in Fig. 1 shows all-sky images in full color (RGB) obtained at LakeLab on Lake Stechlin for a) a clear sky and b) an overcast sky. The upper part of the image is pointing towards north, the lower towards the south, west is on the right-hand side and east on the left-hand side, respectively. In Fig. 1(a), the Milky Way is prominent and some skyglow is apparent near the horizon. In Fig. 1(b), the Milky Way is obstructed by clouds, and a bright glow is apparent in the west. A change in color from blue-green (clear sky) to yellow (clouds) is notable. Three regions of interest for further investigation are indicated in Fig. 1(a,b). The red region points to a sky segment with almost no apparent skyglow from nearby settlements, roughly to the north, centered at an azimuth angle of 350° spanning over a width of 10°. The blue region points roughly to the south, centered at an azimuth angle of 160° (Berlin direction) spanning over a width of 30°. The green region points roughly to the west, centered at an azimuth angle of 280° (Rheinsberg power plant direction) spanning over a width of 10°.

#### Luminance maps

From the all-sky images shown in Fig. 1(a,b), the luminance (similar to “brightness” with respect to human vision) was calculated using the software Sky Quality Camera (SQC, see methods) for each pixel. The resulting luminance maps are shown in Fig. 1(c) for the clear night and d) the overcast night.

For the clear night (Fig. 1c), the luminance ranges from near-natural values of about 0.25 mcd/m2 around the zenith (excluding the Milky Way) to 1.0–2.0 mcd/m2 in the direction of nearby and distant settlements (natural sky is ca. 0.20–0.25 mcd/m24). These include the city of Berlin (azimuth ca. 160°) in the south-southeast (bottom left) and the towns of Neustrelitz (azimuth ca. 5°) and Neubrandenburg (azimuth ca. 10°), both to the north (top of image). There is a small skyglow belt near the horizon with an average luminance reaching values on the order of 0.7 mcd/m2. See Fig. 5 in the methods section showing maps of the region.

For the overcast night (Fig. 1d), the luminance of 0.20 mcd/m2 near zenith is slightly lower than for the clear night, and the skyglow belt near the horizon has vanished. The horizon is dark in most directions, with the notable exception of bright areas in the direction of nearby light sources. Most striking of these is the former nuclear power plant of Rheinsberg at ca. 2.8 km distance to the west (azimuth ca. 280°, right hand side). Others are Fürstenberg (azimuth ca. 60°) and Rheinsberg (azimuth ca. 240°). All three local sites are not very apparent in the clear sky data.

The difference in luminance between the two nights was acquired by subtracting the clear night data from the cloudy night data. The result is shown in Fig. 1(e), with green regions being darker with clouds and red regions brighter. The subtracted image shows that the decrease in NSB at the zenith is mainly due to the attenuation of the Milky Way’s light by clouds. Furthermore, clouds obviously darken the sky near the horizon with the exceptions of the power plant (right side, west) and the two skyglow regions at the very edge of the image (upper right, northeast and lower left, southwest). Strikingly, on clear nights the skyglow from distant sources (e.g. from Berlin, Neustrelitz, Neubrandenburg) is dominant, while the overcast skyglow is dominated by local light sources (Power plant, Rheinsberg and Fürstenberg).

#### CCT maps

The lower row in Fig. 1(f–h) shows CCT maps calculated from the RGB images at LakeLab. Figure 1f) shows the clear night and g) the overcast night data, while Fig. 1h) shows the difference where the data from the clear night was subtracted from the cloudy night. During the clear night (Fig. 1f), the CCT near zenith (Milky Way up) is 4700 K and the CCT values decrease towards the horizon. For the overcast night (Fig. 1g), the CCT values around zenith (10° circle) are on the order of 3400 K. The skyglow regions to the south east, north east and south west reach low CCT values on the order of 2000 K to 2500 K. The skyglow from the power plant (west) has about 3000 K CCT, hardly distinguishable from the background. In the CCT difference map, shown in Fig. 1(h), red and yellow colors indicate a decrease in CCT due to clouds, while green and blue regions indicate an increase in CCT by clouds. Overall, only a decrease in CCT by clouds is observed.

#### Angular luminance distribution

From the spatially resolved data shown in Fig. 1(c,d), we can derive the angular distribution of luminance from zenith to the horizon. This can be done for the whole image, and for segments of certain azimuth angles (see Fig. 1(a,b) and text above. The angular luminance distribution is shown in Fig. 2 for a) the average over all azimuth angles, b) the sky segment facing the north, c) the sky segment facing the south and d) the sky segment facing the west. Each graph shows clear (blue diamonds) and overcast (black circles) data. The region in the north shown in Fig. 2(b) was chosen as it is least affected by skyglow, the one in the south (Berlin direction), shown in Fig. 2(c), was chosen as it is most affected by distant skyglow and the one in the west (power plant direction), shown in Fig. 2(d), was chosen as it is most affected by nearby skyglow.

All plots show a darkening by clouds near zenith. For the region in the north and the south as well as the average over the whole image, the luminance of the cloudy night is lower than that of the clear night for all zenith angles. For the region in the west, the luminance of the cloudy night is lower only for zenith angles below 40°.

The general trend for the clear sky data is that the near-zenith values are elevated due to the presence of the Milky Way and that apart from that, the luminance rises towards the horizon. This rise depends on the amount of artificial skyglow present at the individual direction. No such general trend is obvious in the cloudy sky data. However, for the north and south, the luminance of the cloudy sky decreases towards the horizon and then peaks at different zenith angles. Minima occur at about 71° for the north reaching 0.11 mcd/m2 and at 60° for the south reaching 0.15 mcd/m2. The peaks occur at about 86° for the north reaching 0.24 mcd/m2 and at 81° for the south reaching 0.30 mcd/m2, respectively.

The decrease in luminance towards the horizon on the cloudy night is almost overcompensated for the direction towards the west. However, luminance values remain almost constant between 0.20 and 0.21 mcd/m2 up to 30°. A distinct peak with a maximum luminance of 1.1 mcd/m2 is apparent at a zenith angle of 72°. We attribute this peak to the lighting system of the former nuclear power plant of Rheinsberg, which is not visible in the clear sky data.

#### Cloud attenuation factors

We calculated cloud attenuation factors for each investigated sky segment from the angular luminance distribution shown in Fig. 2(a–d) by dividing the clear sky data by the cloudy sky data. The result is plotted in the lower row of Fig. 2(e). Black squares represent the northern segment, diamonds the southern segment, crosses the average and circles the western segment. Near zenith, a darkening by almost a factor of 2 is observed. Above 30°, the western segment shows a brightening by clouds, while the other show darkening. In the north (black squares) a maximum attenuation factor of 5.2 is observed at 78°, in the south a maximum attenuation factor of 4.2 is observed at 82° and in the averaged data, a maximum attenuation factor of 3.8 is observed at 84°. In the west the attenuation factor reaches a minimum of 0.4 (i.e. a brightening factor of 2.5) at 72°.

#### Illuminance

Horizontal and scalar illuminance (the total luminous flux incident on a surface, see methods) is calculated by the SQC software by integrating over the whole image. It is further possible to estimate (a hypothetical) illuminance assuming a certain angular luminance distribution for zenith or azimuth angles. For example, it is common to assume a constant luminance for all azimuth and zenith angles to estimate the illuminance from zenith measurements taken with single channel devices like an SQM (see methods). We can further assume a non homogeneous luminance distribution varying with zenith angle, but with a uniform value for all azimuth angles. In that case, the average luminance distribution shown in Fig. 2(a) would produce equal illuminance as that obtained via integration over all data points. Furthermore, we can infer a hypothetical illuminance under minimal skyglow conditions by assuming that the angular luminance distribution of the northern sky segment (least affected by skyglow) would be valid for all azimuth angles. In Table 1, the horizontal and scalar illuminance from the full image is compared to the inferred illuminance values from the northern sky segment and for illuminance inferred from an SQM handheld photometer (see methods).

### Cloudy sky in remote area near Cape Town, South Africa (Night Sky Caravan Park, Bonnievale)

The second set of observations were obtained at Night Sky Caravan Farm near Bonnievale, South Africa for a nearly cloud-free sky on March 15th 2018 at 22:44 local time and a cloudy sky on March 14th 2018 at 21:50 local time (+2 hours GMT). The cloud base height for the cloudy night was inferred to be hCB ≈ 1.5 km ± 0.2 km. The site lies in the Cape Wine district about 180 km east of Cape Town (see methods for details).

#### Full color images, Bonnievale

The upper row in Fig. 3 shows all-sky RGB images obtained at Night Sky Caravan Park near Bonnievale, South Africa for a) an almost clear sky and b) a cloudy sky. The upper part of the image is pointing towards north, the lower towards the south, west is on the right-hand side and east on the left-hand side, respectively. In Fig. 3(a), the well structured Milky Way is prominent and some skyglow is apparent near the horizon, particularly in the northeast. Some small patchy clouds can be seen in the eastern direction. In Fig. 3(b), the Milky Way is obstructed by clouds and a bright glow is apparent in the northeast. The southwestern direction however, is very dark. As for the photos from Germany, a change in color from blueish (clear sky) to yellowish (clouds) is also notable.

Two regions of interest for further investigation are indicated in Fig. 3(a,b). The red region points to a region with almost no apparent skyglow from nearby settlements, roughly to the south-west, centered at an azimuth angle of 216° spanning over a width of 10°. The blue region points roughly to the region with the highest skyglow near the horizon, centered at an azimuth angle of 36° (Bonnievale direction) spanning over a width of 10°.

#### Luminance maps

The middle row in Fig. 3(c–e) shows luminance maps c) for the clear night and d) the overcast night, whereas e) shows the difference from subtracted data. For the clear night (Fig. 3c), the luminance ranges from near-natural values at zenith of about 0.20 mcd/m2 around the zenith (Milky Way is not at zenith) to 2.0 mcd/m2 in the north-eastern direction. This skyglow originates from the small town of Bonnievale (12 km distance). Skyglow from other nearby towns like Robertson (25 km north), Swellendam (40 km east) or Worcester (65 km north-west) are not apparent in the luminance map as the skyglow is confined due to the nearby mountain ranges or obstructed by trees. See Fig. 6 in the methods section for a map of the region.

For the overcast night (Fig. 3d), the luminance around the zenith (10° circle) is slightly lower than for the clear night, reaching values of about 0.19 mcd/m2. The skyglow of Bonnievale is clearly amplified by clouds, reaching luminance values up to 10 mcd/m2. The sky in the opposite direction, however, is darker with the presence of clouds, reaching minima of about 0.05 mcd/m2.

The difference in luminance between the two nights is shown in Fig. 3e) with green regions being darker with clouds and red regions brighter with clouds. The skyglow in the north-east is amplified (upper left) and the sky light in the south-west attenuated by the clouds, while the luminance directly at zenith is similar under both conditions.

#### CCT maps

The CCT maps obtained from the RGB data for Bonnievale are shown in Fig. 3(f) for the clear night and g) the overcast night, while Fig. 3(h) shows the CCT difference from subtracted data. During the clear night, the CCT values near zenith (10° circle) are on the order of 4000 K, with the Milky Way having slightly higher values on the order of 4500 K compared to the sky background. The CCT values decrease towards the horizon. The skyglow of Bonnievale near the horizon has values of 2200 K. For the overcast night, the CCT values around zenith (10° circle) are on the order of 3400 K. Over most of the sky, the CCT values remain relatively constant. The skylow regions to the southeast, northeast and southwest reach low CCT values on the order of 2000 K to 2500 K. The difference in CCT between the two nights is shown in Fig. 3(h), where the data from the clear night was subtracted from the cloudy night. Red and yellow colors indicate a decrease in CCT due to clouds, while green and blue regions indicate an increase in CCT by clouds. Overall, only a decrease in CCT by clouds is observed.

#### Angular luminance distribution

The angular luminance distribution is shown in the two upper graphs in Fig. 4 for a), the sky segment facing the southwest (red area in Fig. 3) b) the sky segment facing the northeast (blue area in Fig. 3). Each graph shows clear night (blue diamonds) and overcast night (black circles) data.

For the region in the southwest (low skyglow) shown in Fig. 4(a), the luminance of the cloudy night is lower than that of the clear night for all zenith angles. At zenith, the clear sky has near-natural luminance values on the order of 0.2 mcd/m2. The luminance rises with zenith angle, peaking at around 30° due to the Milky Way and then rising to about 0.5 mcd/m2 near the horizon. The luminance of the cloudy sky is also close to 0.2 mcd/m2 near the zenith and gradually decreases towards the horizon. The luminance minimum of 0.05 mcd/m2 is reached at about 75°. For zenith angles above 85° the luminance drops to values near zero, as this part of the sky is shielded by nearby objects.

For the region in the northeast (Bonnievale direction) shown in Fig. 4(b), the luminance of the cloudy night is higher than that of the clear night for almost all zenith angles. The clear sky has near-natural luminance values near zenith and the luminance rises towards the horizon to values of 2.5 mcd/m2 at 86°, before it drops sharply. The luminance of the cloudy sky is also close to 0.2 mcd/m2 at zenith, remaining below 0.9 mcd/m2 up to angles of 67°. Above that angle, it rapidly rises before it peaks at 9.1 mcd/m2 at about 83°.

#### Cloud attenuation factors

We calculated the factor for brightening (amplification) or darkening (attenuation) due to the presence of clouds for the two investigated sky segments from the angular luminance distribution shown in Fig. 4(a,b). The result, obtained by dividing the clear sky data by the cloudy sky data, is plotted in Fig. 4(c). Grey circles show the northeastern segment and black squares represent the sky segment facing the southwest. The attenuation by clouds observed in the southwestern segment (grey circles) is near unity at zenith rising slowly to values of about 3 at 65°, with a local maximum of 4.8 at 37° due to the Milky Way. Above 65° it rises to a maximum of 11.3 at 86°. In the other direction (northeast, black squares) the attenuation is also near unity at zenith and then drops slowly to values of about 0.5 at 70°. Above 70° it drops to a minimum of 0.15 at 83°, representing a brightening by a factor of about 6.7.

#### Illuminance

As described above, horizontal and scalar illuminance is calculated by the SQC software and it is further possible to estimate (a hypothetical) illuminance assuming a certain angular luminance distribution. In Table 2, the horizontal and scalar illuminance from the full image is provided and compared to the inferred illuminance values from the two selected sky segments, and for zenith luminance measured an SQM handheld photometer (see methods).

## Summary and Discussion

### Reference points for natural nocturnal light environments

ALAN influences a large portion of the natural nocturnal landscapes worldwide4, and its use is growing in area and intensity5. The NSB for clear skies is relatively well understood and a global model for NSB has been published recently4. A natural clear starlit sky has a zenith luminance around 0.2–0.3 mcd/m2 (depending on the position of the Milky Way and the amount of airglow), and that the clear sky horizontal illuminance is around 0.6–0.9 mlx34, setting lower limits for moonless clear nights.

While some studies regarding the brightening of the night sky by clouds in urban areas exist, very few data on the attenuation (darkening) by clouds has been published (but see17,35,36 for data from single channel devices). A commonly cited value of 0.1 mlx for illuminance for cloudy nights7,39,40 is based on extensive (daytime and nighttime) illuminance and luminance measurements performed by D.R.E. Brown at multiple locations several decades ago41. Brown’s seminal report has produced excellent daylight and moonlight lookup charts that have been proven relatively robust. However, the accuracy of the Luckiesh-Taylor brightness meter42 used for the nighttime measurements at very low light levels remains unclear. Furthermore, the attenuation by clouds at night were not measured but estimated. Knowledge of the lower limit of nocturnal darkness for cloudy nights, however, is an important reference point to estimate and model impacts of ALAN on species level or even for whole ecosystems. At the moment, the bright values for cloudy nights in urban areas are compared to the reference values of a typical clear night (see e.g.16 for a brightening as large as a factor of 2200). However, cloudy night values obtained with ALAN should be better compared with the cloudy night values without ALAN as the amplification is much higher. The reference point idea can be put forward to other weather phenomena such as fog and extended to nights with moonlight, as darkening and brightening of moonlit skies will be much more complex and will depend on cloud height, type and optical thickness.

Reference points for illuminance under natural conditions are also useful for planning the experimental design of laboratory experiments, to mimic natural (daytime and nocturnal) conditions of specific habitats indoor. Right now, many laboratory experiments use simply “darkness” as a dark control level43,44 (a value of “0 lx” is often reported, although we note that this presumably indicates the illuminance was simply below sensor threshold, which in the case of inexpensive lightmeters could mean as much as 0.5 lx). Sometimes experiments use (often constant) moonlight levels on the order of 0.1–1 lx45 which does not represent a natural nocturnal light environment.

### Cloud attenuation

At the rural site at LakeLab in Germany, we observed an attenuation by clouds at zenith by a factor of 1.5 and a reduction of horizontal illuminance by a factor of 1.6, when taking into account the full data-set. Off-zenith cloud attenuation factors of more than 5 were observed for the sky segment least affected by ALAN and cloud amplification factors of more than 4 were obtained for the sky segment most affected by ALAN. The minimum off-zenith sky luminance was determined to be 0.11 mcd/m2. From the sky segment least affected by skyglow, a hypothetical horizontal illuminance of 1.04 mlx for clear sky and 0.51 mlx for overcast sky was inferred from the measured angular luminance distribution.

At the South African site, we observed an attenuation by clouds at zenith by a factor of 1.25, and a reduction of horizontal illuminance by a factor of 1.25. Off-zenith attenuation factors of more than 11 were obtained for the sky segment least affected by ALAN, while a brightening by a factor of more than 6 was observed in the opposite direction for the sky segment most affected by ALAN. The minimum off-zenith sky luminance was determined to be 0.05 mcd/m2. By inspecting the sky segment least affected by skyglow, a hypothetical horizontal illuminance of 0.93 mlx for clear sky and 0.34 mlx for overcast sky was inferred.

### Color shift

The multi-spectral information from the individual red, green and blue color channels of the camera indicates a spectral change with changing weather conditions. The CCT data shows an overall reduction (indicating a spectral shift towards red) in color temperature for overcast conditions at both sites, in agreement with previously reported results using SQMs and color filters in an urban context15,46. We are aware of the deficiencies of the CCT method to characterize artificial light sources that deviate from a perfect black body radiator47,48. At this point, however, CCT is used widely in the lighting industry and the public awareness for CCT is growing. Furthermore, if the raw camera data is stored and archived properly, the three color channels can be investigated with a different method at a later stage. Please be aware that a standard DSLR is still limited for ecological studies as it lacks ultraviolet or near infrared sensitivity. A modification of digital cameras to extend the measurement range is not straightforward, but possible49.

The ability to observe changes in the nighttime spectrum over time is important not only for astronomical observations, but also for the evaluation of possible ecological impacts of ALAN. Due to the transition towards solid-state lighting with LEDs, it is difficult to judge whether a site has become brighter or not50 when devices are used that are not matching photometric visual band such as the SQM (see the comprehensive discussion on this topic by Sanchez de Miguel et al.51 or recent work by Bouroussis52) or even nighttime satellites, as the DMSP and VIIRS sensors lack sensitivity at short visible wavelengths5,53.

### Distance shift

The presence of clouds changes the ALAN sources that affect an individual site. While for clear skies distant ALAN sources can have a major impact on the NSB, these distant sources can be attenuated due to atmospheric conditions and local light sources can dominate with their contribution to the NSB. This effect will also depend on the cloud height17. In our study, this effect is exemplified in the LakeLab data obtained near Berlin (best seen in the subtracted image Fig. 1(e). There, the long distance skyglow (e.g. from Berlin, Neustrelitz, Neubrandenburg) dominates the clear sky while local light sources (e.g. Rheinsberg and Fürstenberg) dominate the overcast sky. Most strikingly of the local light sources is the former nuclear power plant of Rheinsberg at ca. 2.5 km distance to the west. Due to the spatial information, we were able to identify the former power plant as the main source of (local) skyglow for the overcast conditions, while the distant skyglow from e.g. Berlin was attenuated by clouds. During the clear night, the light emission from the power plant and other local sites was not dominant in the NSB observations while the Berlin skyglow was prominent. However, the power plant light emission had only a small impact on the zenith luminance for both clear and cloudy night. This shift in source would not have been apparent via observations of zenith luminance with a single channel sensor devices such as an SQM or illuminance with a luxmeter. This means that a site that seems free from ALAN under conditions of clear skies (and therefore potentially fulfills requirements for a dark sky park32), is not necessarily free from ecological light pollution.

We have used a commercial DSLR camera with a fisheye lens for our study, as the easy-to-use method30,31,32 has proven to be applicable for field work18,33,38. The data was processed with a commercial software, SkyQualityCamera, SQC (see methods) that has an intuitive graphical user interface and a comprehensive set of features. SQC allows the analysis of such data without any programming skills or knowledge of photometry. An alternative to this software that requires some basic programming knowledge is the free software DiCaLum54, which we used in an earlier publication18.

The all-sky multi-spectral imaging data is superior to single channel devices as the luminance maps not only provide a zenith NSB value, but also allow a comprehensive investigation of the angular luminance distribution and the color temperature. The spatially resolved data enables us to identify the (nearly) unpolluted parts of the sky and extract cloud attenuation factors for these parts, as well as to perform a direct comparison with polluted skies, using simultaneously acquired data. Although it is feasible to attach single channel devices to translation stages to scan the different segments of the sky55 this is time consuming and requires very stable atmospheric conditions, which is problematic for cloud investigations.

Overall, the method applied here has potential to be widely used by amateurs as well as professionals from astronomy, life sciences or interdisciplinary fields to complement or replace SQM measurements.

### Current obstacles and future directions

At each site, two consecutive nights, one clear and one with clouds were investigated to ensure that artificial illumination and natural sky background are comparable. A comparison between different months or even seasons is not advisable, as lighting practice and natural background can change considerably. We note that it is rare to have the opportunity of being in an area that is almost unpolluted by ALAN, and then have the weather conditions change from clear to overcast from one night to the next, with both occurring during the monthly window of moonless nights. Despite this, we were able to acquire two data sets at two different rural locations. The amount of data available for further analysis could be increased either by setting up permanently measuring systems, or increasing the number of people acquiring data.

We therefore urge amateur and professional astronomers or interested citizen scientists as well as ALAN researchers from other fields to obtain such cloudy night data, and send it to light pollution researchers, or to analyze and publish it by themselves. A first step could be to approach us or the authors of the DiCaLum code54. Perhaps astrophotographers who travel long distances to very remote sites tend not to take images if the sky is overcast. However, such photographers could help building a data set that improves our understanding of light pollution. Furthermore, more of such studies on a systematic basis are required to cover different ecosystems (e.g. marine systems or polar regions), and different climate zones or other weather conditions like fog, snow and ice in the atmosphere.

## Conclusions

Understanding of the Earth’s illuminance on cloudy nights under natural conditions is a major knowledge gap. This reference point is essential for developing new experiments and performing ecological modeling in ALAN research. We have investigated the spectral color shift and the attenuation (or amplification) of the sky brightness by nocturnal clouds. Our main results can be listed as follows:

• Luminance levels as low as 0.05 mcd/m2 (off-zenith) were observed for cloudy nights. Our results suggest an upper limit on natural horizontal illuminance under overcast conditions of 0.34 mlx. We expect more remote areas to have even lower illuminance.

• Clouds strongly shift the distance scale of skyglow; sites that appear to be nearly unpolluted on clear nights might be quite strongly affected by (local) ALAN on cloudy nights.

• A color shift to the red spectrum was observed for cloudy nights, which matches earlier observations in urban contexts15,46.

These results were obtained using a simple and robust method: differential all-sky photometry with a commercial digital single lens reflex camera and a fisheye lens. The results represent a first step towards discovering the missing reference points for natural nocturnal darkness in the context of ecological light pollution. Because both our sites were affected to some extent by artificial light, further investigations of the impact of clouds are necessary to fully understand the potential consequences of ecological light pollution on the environment. We urge that researchers and interested photographers not only measure zenith NSB on clear nights, but also use imaging sensors under different atmospheric conditions to quantify the nocturnal light environment.

## Methods

### Study sites

#### LakeLab, Lake Stechlin, Germany

The first set of observations were made at the LakeLab (53° 8′ 36.5″N, 13° 1′42.7″E), a floating ecological research platform installed in the southwestern bay of Lake Stechlin, a clearwater lake with a surface area of 4.25 km2 and a maximum depth of 69.5 m. The lake is situated about 75 km north of the city of Berlin in a forested area. The floating construction consists of one large enclosure with 30 m diameter and 24 smaller enclosures that are 9 m in diameter. The entire structure is anchored to the lake bottom to prevent drift. For more details about the LakeLab, see our previous report of the zenith brightness over several months35. The new World Atlas of artificial NSB4 predicts a clear sky luminance at zenith of 0.22 mcd/m2 at this site56.

Lake Stechlin is part of a nature reserve belonging to one of the areas of Germany least affected by skyglow (according to Falchi et al.4). The nearest village is Neuglobsow (53° 8′48 ″N, 13°2′49″E) with less than 300 inhabitants and the largest nearby major settlement is the city of Berlin (distance ca. 75 km, azimuth ca. 162°, about 3.5 Mio. inhabitants) in the south-south-east. Other settlements contributing to skyglow at LakeLab are the towns of Neustrelitz (ca. 25 km distance, azimuth ca. 5°, ca. 20,000 inhabitants) and Neubrandenburg (ca. 50 km distance, azimuth ca. 10°, ca. 65,000 inhabitants), both to the north, Fürstenberg (ca. 9 km distance, azimuth ca. 60°, ca. 6000 inhabitants) in the northeast and Rheinsberg (ca. 10 km distance, azimuth ca. 243°, ca. 8000 inhabitants) in the southwest. The former nuclear power plant of Rheinsberg is located at ca. 2.8 km distance to the west (azimuth ca. 280°). For locations and distances see maps in Fig. 5.

#### Bonnievale, South Africa

The second set of observations were made at the Night Sky Caravan Park (34°0′44.8″S, 19°59′41.4″E), an astro-tourism site near Bonnievale, South Africa. The site lies in the Cape Winelands district, part of the western Cape province about 180 km east of Cape Town, see Fig. 6 for a map of the region. Nearby settlements are the towns of Bonnievale (12 km distance, azimuth ca. 12°, ca. 9000 inhabitants) to the northeast, Robertson (25 km, azimuth ca. 335°, ca. 28.000 inhabitants) to the north, Swellendam (40 km, azimuth ca. 90°, 17.000 inhabitants) to the east and Worcester (65 km, azimuth ca. 310°, ca. 97.000 inhabitants) to the northwest.

### DSLR camera and image processing software

Images were obtained with a commercial DSLR camera, a Canon EOS 6D, which has a full-frame 20.2 Megapixel (5496 × 3670) CMOS sensor. The camera was mounted on a tripod and operated with a circular fisheye lens (Sigma EX DG with 8 mm focal length) always at full aperture of 3.5. To acquire all-sky images, the camera was aligned with the center of the lens oriented towards the zenith. At Lake Stechlin, Germany, ISO 3200 and shutter speed of 30 s was used. At Bonnievale in South Africa, the camera was set to ISO 1600 and a shutter speed of 30 s.

For image processing we used the commercial “Sky Quality Camera” software (Version 1.8, Euromix, Ljubljana, Slovenia). The camera is calibrated by the software manufacturer. The manufacturer’s calibration includes photometric calibration using the green channel of the camera, as well as correction of optical aberrations like vignetting. The software processes the luminance Lv,sky of the sky for each pixel of the camera (luminance is commonly referred to as “brightness” referenced to human vision). The software calculates the cosine corrected illuminance Ev,cos in the imaging plane:

$${E}_{v,cos}={\int }_{0}^{\frac{\pi }{2}}\,{\int }_{0}^{2\pi }\,{L}_{v,sky}(\theta ,\varphi )sin\theta cos\theta d\varphi d\theta ,$$
(1)

and the scalar illuminance for the imaging hemisphere Ev,scal,hem without cosine correction57:

$${E}_{v,scal,hem}={\int }_{0}^{\frac{\pi }{2}}\,{\int }_{0}^{2\pi }\,{L}_{v,sky}(\theta ,\varphi )sin\theta d\varphi d\theta \mathrm{.}$$
(2)

In the equations, Lv, sky is the sky luminance, θ is the zenith angle and ϕ is the azimuth angle. For all-sky images, i.e. when imaging in the horizontal plane, Ev,cos is usually termed horizontal illuminance.

From the three spectral channels, color correlated temperature (CCT) is calculated by the software. CCT describes the temperature of a thermal source that most closely matches the color spectrum of the observed light. Candle light has a low CCT (around 2000 K), the sun has a CCT of about 5800 K, and even higher temperatures appear bluish-white. In urban areas, a spectral shift from blue to red (i.e. reduction in color temperature) has been reported for overcast compared to clear skies15,46. The software provides luminance and CCT maps, allows the subtraction of one image from another, to analyze sectors of the image, and to do cross sections.

### Sky quality meter (SQM)

For comparison to the imaging results, we used a sky quality meter with a lens (SQM-LU). It measures the sky radiance LSQM in units magSQM/arcsec2 in a cone with a full width half maximum of 20° in a single spectral channel only approximately resembling human photopic vision (see review comparing devices34, or51 for a detailed discussion of the spectral sensitivity of the SQM). Due to the spectral mismatch, the sky radiance LSQM measured with an SQM can only approximately be converted to a sky luminance at zenith Lv,sky,zen using:

$${L}_{v,sky,zen}\,[{\mathrm{cd}/{\rm{m}}}^{2}]\approx 10.8\times {10}^{4}\times {10}^{-0.4\ast {L}_{{\rm{SQM}}}}\mathrm{.}$$
(3)

Please be aware that this commonly used equation overestimates the precision achievable with the SQM. The interested reader is referred to some recent work by Bará on the conversion58. Using SQM data obtained at zenith, we inferred illuminance values assuming a homogeneous brightness distribution over the hemisphere according to Ev,hor ≈ πLv,sky,zen59, which is a rough estimation commonly done in the literature to link SQM measurements with illuminance measurements often done in field work. While the angular NSB distribution will vary from site to site and from season to season, this crude estimate provides a reasonable approximation. A discussion on the relation between zenith sky brightness and horizontal illuminance is outside the scope of this work, but can be found in59.

### Estimating cloud base height from the all-sky imaging data

The information of cloud base height is important for studies of ALAN on overcast nights. Ribas et al. used laser ceilometer data in their study17, but such data is not always easily accessible. With the all-sky imaging technique used here, the cloud base height information can be inferred from the imaging data itself. If the distance to a light source is known it is possible to estimate the cloud base height hCB from the projection of the light on the cloud base by triangulation using the distance d and the elevation angle ε, with tan(ε) = hCB/d. This method was used in the past to estimate cloud heights in optical ceilometers using a narrow beam of light60 and recently demonstrated with night-time imaging of light pollution61. For the German site at LakeLab, the distance to the power plant is ca. 2.8 km with ε ≈ 17° ± 1° resulting in a cloud base height of hCB ≈ 850 m ± 100 m. This value can be confirmed by checking data for Fürstenberg. The distance here is ca. 8.9 km with ε ≈ 6° ± 1° resulting in a cloud base height of hCB ≈ 950 m ± 150 m. For the South African site, the distance to the town of Bonnievale is ca. 12 km with ε ≈ 7° ± 1°, which results in a cloud base height of hCB ≈ 1.5 km ±  0.2 km. Please be aware that this method does not work at large distances due to the curvature of the Earth and lens aberrations for small angles.

## Data Accessibility

The datasets generated during and/or analysed during the current study are available from the corresponding author on request.