16  Working with units and experimental errors

16.1 Working with units

It is easy to work with units in R thanks to the package units (see vignette).

Working with the units package can prove a very good idea to avoid conversion errors in your data treatment…

Here is the gist of it:

# Load the 'units' library
library(units)
t <- seq(0.1,1,length=3)
# attribute a unit, here 'seconds':
t <- set_units(t, "s") 
t

#> Units: [s]
#> [1] 0.10 0.55 1.00

And it works with the pipe too:

t <- seq(0.1,1,length=3) |> set_units("s")
t

#> Units: [s]
#> [1] 0.10 0.55 1.00

• When possible, automatic units conversion is performed. Also, it is possible to attribute units to columns of tibbles and data.frames:

d1    <- seq(1,2,length=3) |> set_units("m")
(tib1 <- tibble(t=t, d=d1, speed=d1/t))

#> # A tibble: 3 × 3
#>      t   d speed
#>    [s] [m] [m/s]
#> 1 0.1  1   10   
#> 2 0.55 1.5  2.73
#> 3 1    2    2

d2    <- seq(0.1, .3, length=3) |> set_units("cm")
(tib2 <- tibble(t=t, d=d2, speed=d2/t))

#> # A tibble: 3 × 3
#>      t    d  speed
#>    [s] [cm] [cm/s]
#> 1 0.1   0.1  1    
#> 2 0.55  0.2  0.364
#> 3 1     0.3  0.3

bind_rows(tib1,tib2)

#> # A tibble: 6 × 3
#>      t     d    speed
#>    [s]   [m]    [m/s]
#> 1 0.1  1     10      
#> 2 0.55 1.5    2.73   
#> 3 1    2      2      
#> 4 0.1  0.001  0.01   
#> 5 0.55 0.002  0.00364
#> 6 1    0.003  0.003

• You can moreover convert between units systems using set_units(vector, "unit") or units(vector) <- "unit":

T <- set_units(1, "fs")
set_units(1/T, "THz")

#> 1000 [THz]

1/T |> set_units(THz)

#> 0.001 [1/THz]

units(T) <- "ns"; T

#> 1e-06 [ns]

• Sometimes, you may define the unit outside of the set_units() call and want to retrieve it, or use a unit from another variable in a new variable. For this, use the mode="standard" option:

UNIT <- "m"
set_units(1:10, UNIT)

#> Error: In 'UNIT', 'UNIT' is not recognized by udunits.
#> 
#> See a table of valid unit symbols and names with valid_udunits().
#> Custom user-defined units can be added with install_unit().
#> 
#> See a table of valid unit prefixes with valid_udunits_prefixes().
#> Prefixes will automatically work with any user-defined unit.

set_units(1:10, UNIT, mode="standard")

#> Units: [m]
#>  [1]  1  2  3  4  5  6  7  8  9 10

set_units(1:10, units(T), mode="standard")

#> Units: [ns]
#>  [1]  1  2  3  4  5  6  7  8  9 10

• You can give units to table columns:

library(tidyverse)
starwars

#> # A tibble: 87 × 14
#>    name     height  mass hair_color skin_color eye_color birth_year sex   gender
#>    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
#>  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
#>  2 C-3PO       167    75 <NA>       gold       yellow         112   none  mascu…
#>  3 R2-D2        96    32 <NA>       white, bl… red             33   none  mascu…
#>  4 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
#>  5 Leia Or…    150    49 brown      light      brown           19   fema… femin…
#>  6 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
#>  7 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
#>  8 R5-D4        97    32 <NA>       white, red red             NA   none  mascu…
#>  9 Biggs D…    183    84 black      light      brown           24   male  mascu…
#> 10 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
#> # ℹ 77 more rows
#> # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
#> #   vehicles <list>, starships <list>

starwars |> 
    mutate(height = set_units(height,"cm"),
           mass   = set_units(mass,"kg"))

#> # A tibble: 87 × 14
#>    name     height mass hair_color  skin_color eye_color birth_year sex   gender
#>    <chr>      [cm] [kg] <chr>       <chr>      <chr>          <dbl> <chr> <chr> 
#>  1 Luke Sk…    172   77 blond       fair       blue            19   male  mascu…
#>  2 C-3PO       167   75 <NA>        gold       yellow         112   none  mascu…
#>  3 R2-D2        96   32 <NA>        white, bl… red             33   none  mascu…
#>  4 Darth V…    202  136 none        white      yellow          41.9 male  mascu…
#>  5 Leia Or…    150   49 brown       light      brown           19   fema… femin…
#>  6 Owen La…    178  120 brown, grey light      blue            52   male  mascu…
#>  7 Beru Wh…    165   75 brown       light      blue            47   fema… femin…
#>  8 R5-D4        97   32 <NA>        white, red red             NA   none  mascu…
#>  9 Biggs D…    183   84 black       light      brown           24   male  mascu…
#> 10 Obi-Wan…    182   77 auburn, wh… fair       blue-gray       57   male  mascu…
#> # ℹ 77 more rows
#> # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
#> #   vehicles <list>, starships <list>

• You can plot using units and easily convert between units while plotting:

p <- starwars |> 
    mutate(height = set_units(height,"cm"),
           mass   = set_units(mass,"kg")) |> 
    filter(sex != "hermaphroditic") |>
    ggplot(aes(x=height, y=mass, color=sex))+
        geom_point(size=2)+
        labs(x="Height", y="Mass")
p

p + ggforce::scale_x_unit(unit = "inches") +
    ggforce::scale_y_unit(unit = "pounds")

16.2 Working with experimental errors

When working in experimental science, you have to account for measurement errors and error propagation all along your data treatment. This is made really easy thanks to the quantities package that gathers the error and units packages. Most importantly, this allows you propagating the errors in the proper way. So, you input your experimental error once, and you don’t have to think about it anymore. Neat, isn’t it?

Here is the gist of it:

library(quantities)
options(errors.notation="plus-minus", errors.digits=4)
a <- set_errors(1, 0.1)
b <- 2 |> set_errors(0.2)
a+b

#> 3.0000 ± 0.2236

a*b

#> 2.0000 ± 0.2828

a^3

#> 1.0000 ± 0.3000

errors(a)

#> [1] 0.1

errors_min(a)

#> [1] 0.9

errors_max(a)

#> [1] 1.1

It thus becomes easy to plot the error bars from your experimental data. I recommend using the ggforce library to make ggplot2 work better with quantities:

library(ggforce)
options(errors.notation="parenthesis", errors.digits=1)
starwars |> 
    mutate(height=set_quantities(height,"cm",height*.05),
           mass=set_quantities(mass,"kg",mass*.05)
           )

#> # A tibble: 87 × 14
#>    name         height     mass hair_color skin_color eye_color birth_year sex  
#>    <chr>      (err) [… (err) [… <chr>      <chr>      <chr>          <dbl> <chr>
#>  1 Luke Skyw…   172(9)    77(4) blond      fair       blue            19   male 
#>  2 C-3PO        167(8)    75(4) <NA>       gold       yellow         112   none 
#>  3 R2-D2         96(5)    32(2) <NA>       white, bl… red             33   none 
#>  4 Darth Vad…  200(10)   136(7) none       white      yellow          41.9 male 
#>  5 Leia Orga…   150(8)    49(2) brown      light      brown           19   fema…
#>  6 Owen Lars    178(9)   120(6) brown, gr… light      blue            52   male 
#>  7 Beru Whit…   165(8)    75(4) brown      light      blue            47   fema…
#>  8 R5-D4         97(5)    32(2) <NA>       white, red red             NA   none 
#>  9 Biggs Dar…   183(9)    84(4) black      light      brown           24   male 
#> 10 Obi-Wan K…   182(9)    77(4) auburn, w… fair       blue-gray       57   male 
#> # ℹ 77 more rows
#> # ℹ 6 more variables: gender <chr>, homeworld <chr>, species <chr>,
#> #   films <list>, vehicles <list>, starships <list>

starwars |> 
    mutate(height=set_quantities(height,"cm",height*.05),
           mass=set_quantities(mass,"kg",mass*.05)
           ) |> 
    filter(sex!="hermaphroditic") |> 
    ggplot(aes(x=height, y=mass, color=sex))+
        geom_point(size=2)+
        labs(x="Height", y="Mass")+
        geom_errorbar(aes(ymin=errors_min(mass),
                          ymax=errors_max(mass)))+
        geom_errorbarh(aes(xmin=errors_min(height),
                           xmax=errors_max(height)))

16.3 Using physical constants

Now that we know how to work with units, it becomes super easy to use physical constants in our calculations. The constants package comes with a lot of physical constants defined with their units and errors:

library(constants)
library(units)
# Boltzmann constant
syms_with_units$k

#> 1.380649e-23 [J/K]

syms_with_units$k |> set_units(meV/K)

#> 0.0861733 [meV/K]

# Planck constant
syms_with_units$h

#> 6.62607e-34 [J/Hz]

syms_with_units$hbar

#> 1.054572e-34 [J*s]

# Avogadro number: find the symbol
lookup("Avogadro")

#>    symbol          quantity                             type        value
#> 36     na Avogadro constant Physico-chemical, Adopted values 6.022141e+23
#>    uncertainty  unit
#> 36           0 1/mol

syms_with_units$na

#> 6.022141e+23 [1/mol]