[view raw Rmd]

Stevens’s measurement scales

S.S. Stevens’s classical 1946 paper On the Theory of Scales of Measurement tells us there are four measurement scales:

  • nominal,
  • ordinal,
  • interval, and
  • ratio.

R is pretty good at representing the first two by using factor and ordered,

(f = factor(c("d", "a", "b", "c", "a", "b")))

## [1] d a b c a b
## Levels: a b c d

(o = ordered(c("d", "a", "b", "c", "a", "b")))

## [1] d a b c a b
## Levels: a < b < c < d

which give warnings about meaningless operations, like

(e = f * 2)

## Warning in Ops.factor(f, 2): '*' not meaningful for factors

## [1] NA NA NA NA NA NA

and R combines interval and ratio into numeric variables. Having different representations between these different measurement scales has, in my opinion, always been a major advantage of R. It prevents you from doing things that are statistically not meaningful.

Why physical units?

In physics class, we learned that every physical quantity has a measurement unit. If a represents speed, with unit m/s, we can’t add it meaningfully to b which has unit seconds, but we can add it to c measured in km/h after proper unit conversion. Dimensional analysis tracks units of measurements of variables when computations are performed. It is used to determine the unit of measure of the result, but also catches computations that aren’t physically meaningful. Can this be automated?

Physical unit databases, and conversion software

The Unified Code for Units of Measure, or UCUM, is based on the ISO 80000: 2009 Quantities and Units standards series that specify the use of System International (SI). UCUM comes with a BNF grammar and a machine-readable (XML) document with all the units, or all those that are useful – the amount of derivable units is infinite.

Being rather formal, and close to ISO, it is no surprise that UCUM has been recommended for encoding units of measures by many open geospatial consortium standards for spatial data.

A more pragmatic and practical approach is taken by udunits, developed by the geo/atmospheric scientists of UCAR/unidata. Udunits not only consists of an XML file with all the units, their names and symbols, but also of a software library that can validate units, check whether they are convertible (like km/h and m/s) and carry out this conversion. James Hiebert wrote an R package, udunits2, which interfaces to this software library, but does little more than exposing its functions as R functions.

Using physical units in R: the units package

I have always wondered why R has no support for dimensions built in, or at least have a package that does this. Date and POSIXt objects have implicit units (1 day, 1 second), but only time difference difftime objects have explicit, and modifiable units:

t = Sys.time() + 0:3 * 3600
(deltat = diff(t))

## Time differences in hours
## [1] 1 1 1

units(deltat) = "mins"
deltat

## Time differences in mins
## [1] 60 60 60

When I discovered the udunits2 R package, I couldn’t resist writing the units R package, which works similarly to difftime, but for all physical units supported by udunits2. Thus, after

library(units)
(a = as.units(1:5, "m/s"))

## Units: m/s
## [1] 1 2 3 4 5

we can do simple arithmetic:

2 * a

## Units: m/s
## [1]  2  4  6  8 10

a + a

## Units: m/s
## [1]  2  4  6  8 10

a * a

## Units: (m/s)*(m/s)
## [1]  1  4  9 16 25

but also automatic unit conversion

b = as.units(1:5, "km/h")
a + b

## Units: m/s
## [1] 1.277778 2.555556 3.833333 5.111111 6.388889

b + a

## Units: km/h
## [1]  4.6  9.2 13.8 18.4 23.0

a / b

## Units: (m/s)/(km/h)
## [1] 1 1 1 1 1

a * b

## Units: (m/s)*(km/h)
## [1]  1  4  9 16 25

as you can see, units are propagated and converted to that of the first argument when needed, but not simplified. Wrong units trigger an error:

s = as.units(1:5, "s")
e = try(x <- a + s)
attr(e, "condition")[[1]]

## [1] "cannot convert s into m/s"

We can also do comparison and apply basic functions, subset, or concatenate

signif(a^2.5, 3)

## Units: (m/s)^2.5
## [1]  1.00  5.66 15.60 32.00 55.90

a[2:4]

## Units: m/s
## [1] 2 3 4

c(a,b)

## Units: m/s
##  [1] 1.0000000 2.0000000 3.0000000 4.0000000 5.0000000 0.2777778 0.5555556
##  [8] 0.8333333 1.1111111 1.3888889

c(b,a)

## Units: km/h
##  [1]  1.0  2.0  3.0  4.0  5.0  3.6  7.2 10.8 14.4 18.0

Conversion to and from difftime, use in data.frames or matrix objects is illustrated in the package vignette.

The further future

When dealing with measurement unit rigorously, the output of linear regression of two variables, zinc with units ppm and dist with units m would ideally contain them:

> library(sp)
> data(meuse)
> summary(lm(zinc ~ dist, meuse))

Call:
lm(formula = zinc ~ dist, data = meuse)

Residuals (ppm):
    Min      1Q  Median      3Q     Max 
-475.20 -189.94  -52.94  120.15 1088.80 

Coefficients:
            Estimate Units Std. Error t value Pr(>|t|)    
(Intercept)   756.70   ppm      35.66   21.22   <2e-16 ***
dist        -1195.67 ppm/m     114.84  -10.41   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 281.7 ppm on 153 degrees of freedom
Multiple R-squared:  0.4147,    Adjusted R-squared:  0.4109 
F-statistic: 108.4 on 1 and 153 DF,  p-value: < 2.2e-16

I’m convinced this would help understand what residuals, regression coefficient estimates, and standard errors mean.

Getting output like this automatically may not happen any time soon: when solving the normal equations, each entry of the cross product matrix \(X’X\) would need to store its own physical unit, and matrix product and solve routines would need to propagate them.

The near future

It is of course good to know whether R variables are stored as factor or character, as integer or double, but it doesn’t prevent you from adding apples and oranges. Verifying compatibility of physical units does. Dimensional analysis helps here, and helps understanding and verifying meaningfulness of results.

I would be more than happy to hear of any use cases using the units package, be it for educational or operational projects, and also for suggestions how (or pull requests) to improve this package. My wish list right now:

  • add units by default to axis labels of plots
  • support log-units handling of udunits
  • integrate with spatial and temporal reference systems in R
  • link this to our work on meaningful spatial statistics and provenance of data generation