FAQs

R-related FAQ

What version of R is needed?

It is good practice to keep R up to date, so it is recommended that you have the most recent version of R installed. In any case, you should at least have the R version 4.0 or newer installed. You can check the version by running “version” or “sessionInfo()” in R.

What packages are needed?

The packages needed in each unit are loaded with library() in the course book chapter and in the practical. The main ones are:

tidyverse (contains ggplot2, readr, dplyr, tidyr, forcats, and others)
ggfortify (diagnostic plots with autoplot())
car (variance inflation factors with vif())
multcomp (Tukey multiple comparisons)
GGally (ggpairs())
lme4 and lmerTest (mixed models)
vegan (ordination)

On the UZH RStudio Server most of these are already installed. On your own computer, install each one once with install.packages().

I am having issues installing a package

Problems with installing packages are somewhat common. Often, these problems be avoided by using

install.packages( "NAME", type = "binary")

instead of simply

install.packages( "NAME", )

The reason is that quite often, packages are available a source packages, which must be compiled locally, which is often not possible on standard R machines. Therefore, by specifying type = “binary”, you tell R to install the already compiled version, even though it might be an older version.

My output of summary(), str() etc. looks different from older materials

Older materials (books, web pages, older videos) sometimes show output that looks a bit different from what you get. Usually the reason is that, since R version 4.0.0, text columns are read in as character (<chr>) variables rather than factors, and functions such as summary(), glimpse() and str() show character variables differently from factors. Package updates can also change how output looks. In BIO144 we use read_csv() from the readr package, which reads text columns as character variables. That is fine: linear models work with character variables as well as factors (see “Should I convert character variables to factors?” below).

What does rm(list=ls()) do? Should I use it?

rm(list=ls()) deletes all the objects in R’s environment (ls() lists their names and rm() removes them). You will see it at the top of many older scripts, including some older BIO144 materials.

In BIO144 we now recommend something better: restart R (Session > Restart R in RStudio) and run your script from the top. Restarting gives a completely fresh R session: it also unloads packages and resets options, which rm(list=ls()) does not. If your script then runs from top to bottom without errors, you know it does not depend on anything you did outside the script. Working in an RStudio project (see the Unit 1 practical) helps too.

How do I set the number of decimal places shown?

Use the round() function to tell R how many decimal places to report and what rounding to do. You otherwise cannot know for sure how a number displayed in the console has been rounded or not.

What does the predict() function do?

In general, the predict() function takes as argument a model object, for instance one fitted with lm(). It then predicts/estimates/calculates the response variable based on the values of the covariates (explanatory variables) and the estimated coefficients according to the fitted model. You can do this by hand (good for understanding).

If you do not specify the data to use for the prediction, it will automatically use the data that was used in fitting the model. You can also specify the dataset used, which is useful, for instance when you have multiple covariates and you want to visualize the relation of one specific explanatory variable with the response variable. In that case you create a new dataset in which you keep every covariate at constant values (e.g. their mean) with the exception for the one in you are interested in. Then you can use this dataset in the predict() function and then, if you want, you can plot the relation of that variable with the response.

Note: predict() does not fit a model, it evaluates a model with some given data

What does the expand.grid() function do?

The expand.grid() function is a handy function to create new data, as it returns every possible combinations of the input provided. For instance expand.grid(a = 1:2, b = 4:5) returns

We often use this function to create new data because we want to visualize the effect of one particular variable on the response variable according to the fitted model. When we have a model with more than 1 explanatory variable, these variables can vary together and if they do the effect of one particular variable on the response variable is more difficult to see. Hence to best show the importance of a particular variable, we keep all other variables at constant values. These constant values should ideally be sensible values, i.e. keeping the mass at 0 does not make sense, but it makes more sense to keep at the mean value.

Which diagnostic plots does autoplot() show?

By default, autoplot() (from the ggfortify package) shows four diagnostic plots for a linear model: residuals vs fitted values, the normal QQ-plot, the scale–location plot, and residuals vs leverage. You can choose which plots to show with the which argument, as with plot() for linear models:

autoplot(my_model, which = 1:6)   # all six plots
autoplot(my_model, which = 4)     # Cook's distance for each observation

When a model has only categorical explanatory variables, all observations in a group have the same leverage, so the residuals vs leverage plot is replaced by a plot of residuals against the factor levels. That’s fine: the other plots are the more important ones.

Note that autoplot() does not work for mixed models (lmer()): use plot(my_model) for the residuals vs fitted values, and qqnorm(resid(my_model)) for a QQ-plot.

gather() vs pivot_longer() and spread() vs pivot_wider()

Some older materials (and many older web pages) use the tidyr functions gather() and spread(). These have been replaced by pivot_longer() and pivot_wider(), which are easier to use. In BIO144 we use pivot_longer() and pivot_wider() (see the course book, Chapter 2).

unique() vs levels()

The function levels() only works with factors. read_csv() reads text variables as characters, not factors, so levels() gives NULL for them. Use unique() instead of levels() to see the different values of a character variable.

Should I convert character variables to factors?

read_csv() reads text variables (e.g. country, continent) as characters (chr), not factors. Linear models work with both, so you do not have to convert them. The main reason to convert a variable to a factor is to set the order of its levels (e.g. for a graph, or to choose the reference level; see the next question). If you do want to convert one variable:

dd <- dd |> mutate(continent = factor(continent))

And to convert all character variables at once (not usually needed):

dd <- dd |> mutate(across(where(is.character), as.factor))

How do I change the reference level of a categorical variable?

In summary() of a linear model, the coefficients of a categorical variable are differences from the reference level, which by default is the first level in alphabetical order. To use a different reference level, for example “control”, use fct_relevel() from the forcats package (loaded with the tidyverse), then fit the model again:

dd <- dd |>
  mutate(treatment = fct_relevel(treatment, "control"))

The model is the same model (the fitted values and the anova() table do not change), but the coefficients and their t-tests in summary() now compare each level with “control”. See the course book, Chapters 5 and 7.

Plotting one regression line when the points are coloured by group

Question: I map colour = continent in ggplot(), and geom_smooth() then draws a separate regression line for each continent. How do I get one line for all the data, with the points still coloured by continent?

Answer: Aesthetic mappings given in ggplot() are inherited by all the geoms, so geom_smooth() also uses colour = continent. Put the mapping only in the geom that needs it:

ggplot(health_data, aes(x = log10_expenditure, y = log10_mortality)) +
  geom_point(aes(colour = continent)) +
  geom_smooth(method = "lm", se = FALSE)

(Alternatively, keep the mapping in ggplot() and remove it from the smooth with geom_smooth(method = "lm", aes(colour = NULL)).)

ggpairs() shows too many graphs. Can I show only some variables?

Yes: select the variables you want first, then make the graph:

bodyfat_data |>
  select(bodyfat, age, abdomen, height, weight) |>
  ggpairs()

How do I find the extreme values of each variable?

To see the row with the largest value of one variable, use slice_max() (and slice_min() for the smallest):

dd |> slice_max(abdomen, n = 1)

Note that the largest value is not necessarily an extreme value: it may be perfectly in line with the other values. One common definition of an extreme value (used for the points shown separately in box plots) is a value more than 1.5 times the interquartile range below the first quartile or above the third quartile. Graphs (histograms, box plots, scatterplots) are usually the best way to spot extreme values.

Statistics-related FAQ

What is the difference between fitted and predicted values?

Fitted values: The response values you get when you use the observed values of the explanatory variables, i.e. the ones you used to fit the model: the original data

Predicted values: In general, the response values you get when you use new observations of the explanatory variables (i.e. new data). Note: if the new data is the same as the original data, then predicted values == fitted values.

This being said, it happens that the terms fitted values and predicted values are used interchangeably, even though there is a (small) difference.

Interpretation of p-values

The p-value is best interpreted as:

The p-value is the probability of having a test statistics value at least as extreme as the value actually observed (e.g. the t-value of the t-test or F-value of the F-test) given that the Null hypothesis is correct (often H0: parameter=0). If the p-value is below 0.05 we reject H0 and go with H1 (often H1: parameter =/= 0). The p-value is not the probability of H0 being correct.

Interpretation of confidence intervals

Let’s assume we have a simple linear regression with body-fat as the response and BMI as the explanatory variable. We find that the slope of BMI is 1.75 with a corresponding 95% confidence interval [1.55,1.97]. This is interpreted as follows:

“All true values for beta1 (slope for BMI) between 1.55 and 1.97 are compatible with the observed data”.
To rephrase this slightly, this means: the true slope for BMI could have any value between 1.55 and 1.97 and we would not be surprised to see the data that we have (and thus, if the true slope was outside the interval, it would be surprising to see the data that we have: the estimated slope based on the data does not match the true slope.

Note: I want to stress that it is not correct to say that there is a 95% probability that the true slope is between 1.55 and 1.96 (see below).

So, with the interpretation given above, if the confidence interval of the slope does not cover 0, then a slope of 0 is not compatible with the data. In fact, in the case of the 95% confidence interval this means that the p-value is below 0.05, which means that we reject the Null hypothesis H0 that the slope for BMI is 0. Note: this is neither “good” nor “bad”, it is just what it is: we are just interested in the truth.

There is an alternative interpretation of confidence intervals which is arguably more precise (but you can use the above for this course). To explain it, here are a couple of things that I personally think are good to keep in mind:

The parameter of interest (e.g. the slope of BMI) is estimated based on the available data, meaning that the estimate is a function of the data. If for instance the experiment is repeated the data will look a little bit different and the estimate will be a little bit different. However: the true value of the parameter of interest (the slope) is the same in both experiments!
The same thing goes for confidence intervals: they are estimated based on the data, so if the data changes a bit, the upper and the lower bound of the confidence interval will be a bit different as well.
An example of wrong interpretation: continuing the example above, let’s say we know the true slope of BMI and it is 2. What is the probability that the confidence interval includes the true slope? 0 (!), because 2 is not in [1.55,1.97]. And what would the probability be if the true slope is 1.6? 1 (!), because 1.6 is included in [1.55,1.97]. This is why it is incorrect to say that there is a 95% probability that the estimated confidence interval contains the true slope: the estimated confidence interval either contains or it does not contain the true slope.

Now, because of probability theory reasons, a confidence interval is best interpreted like this:

“If we repeat the entire experiment many times (and thus collect slightly different data each time) and for each experiment we calculate the 95% confidence interval, then 95% of the calculated confidence intervals will contain the true parameter value”.
This definition is a bit harder to get, but it basically tells us that it is a measure of uncertainty of the estimated slope.

Model diagnostics residuals versus fitted values

The central assumption of linear regression is that the errors are independent (from each other) and identically distributed (i.i.d.) following a normal distribution with mean = 0 and variance = sigma^2. We cannot measure the errors directly, but we have estimates of them: the residuals. This means that we can check these assumptions based on the residuals, because given that they are estimates of the errors and that the model is valid the residuals have to be i.i.d. following a normal distribution as well. In the figure residuals vs. fitted values we can check for some of these assumptions:

Mean = 0: the mean has to be zero across all fitted values. So the points in the figure have to be spread approximately symmetrically above and below 0, so that the mean is approximately 0 across all fitted values. If for instance we see that for larger fitted values the mean of the residuals increases or decreases clearly away from 0, then this assumption might be violated.
Variance == sigma^2: constant variance. Similarly to above, the variance has to be independent from the fitted values, meaning that the spread of point below and above 0 needs to be approximately the same across all fitted values. If this is not the case this results in a pattern (e.g. larger fitted values result in a higher variance in the residuals) and this assumption might be violated.

Difference between binomial and count data

For count data we use Poisson regression and for binomial data we use binomial regression (commonly called logistic regression)

Count data come from questions like “how many times have you voted in your life?” and the answer can be a positive integer (including 0), while binary data come from questions like “did you vote last year?” to which the answer can only be yes (1) or no (0).
We get binomial data if we have n (=number of persons) independent Bernoulli trials (i.e. binary trials), so we would basically count how many people have voted last year out of the n we asked. In other words, if we have 7 yes and 3 no to the question “did you vote last year?”, they come from 10 different people (sample size n=10), while if to the question “how many times have you voted in your life?” we get an answer like “I’ve voted 10 times in my life” this comes from a single person (sample size n=1), so keep this in mind.

Why is the dispersion reported by a quasi-Poisson model different from residual deviance / residual df?

In BIO144 we check for overdispersion with the ratio of the residual deviance to the residual degrees of freedom. When you then fit a quasi-Poisson (or quasibinomial) model, summary() reports a “dispersion parameter” that is usually a bit different from that ratio. This is because R estimates the dispersion from the Pearson residuals, not the deviance residuals:

sum(residuals(my_quasi_model, type = "pearson")^2) / df.residual(my_quasi_model)

Both are reasonable ways to estimate the dispersion; they just use slightly different kinds of residuals. For example, in the Unit 8 grouse tick data the deviance ratio is about 8.6 and the quasi-Poisson dispersion is about 12.5. Either way, the data are clearly overdispersed.

How large does the dispersion have to be before I should worry?

Values above 1 mean overdispersion, but a little overdispersion is not a problem. The course book (Chapters 8 and 9) uses the rough guide that a ratio of residual deviance to residual degrees of freedom above about 1.5–2 suggests overdispersion that should be dealt with (for example with a quasi-Poisson or quasibinomial model). Like all rules of thumb, this is only a rough guide: how much overdispersion matters also depends on the sample size. Note that this check does not work for binary (0/1) data.

ANOVA or multiple regression?

ANOVA and (multiple) regression are both linear models, fitted in R with lm(). ANOVA is the special case in which all explanatory variables are categorical. The difference is mostly in which output we look at first:

For a categorical explanatory variable with k levels, we first look at the anova() table: one F-test of whether the group means differ at all. The summary() table instead gives k − 1 separate t-tests, each comparing one level with the reference level, which is usually not the question we want to answer first.
For a continuous explanatory variable, the summary() table directly gives the estimate of the slope and its test.

When a model has both kinds of explanatory variable, use both tables.

Degrees of freedom in a one-way ANOVA

With n observations in g groups: the total sum of squares has n − 1 degrees of freedom; the groups (the model) use g − 1; the residuals have n − g. The F-test uses g − 1 and n − g degrees of freedom. For example, 20 observations in 4 groups give F with 3 and 16 degrees of freedom. (Counting the intercept as well, the model estimates g parameters in total, which is why the residual df is n − g.) See “Degrees of freedom” in the course book, Chapters 3 and 5.

Does a t-test work when the two groups have different numbers of observations?

Yes. A t-test (or a linear model) does not need equal group sizes, as long as its assumptions are met and each group has enough observations. Unequal group sizes slightly reduce the precision of the estimate of the difference, compared with the same total split equally. With unequal variances as well, use var.equal = FALSE (Welch’s t-test), which estimates the variance separately in each group.

It is often more important to ask why the group sizes differ. In the class reaction time data, for example, if people of one sex at birth were more willing to take part than the other, the sample may not represent the class, whatever the test says.

What does “a distribution with a long tail of high values” mean?

It means that large values are more frequent than they would be in a symmetric distribution (such as the normal distribution): the distribution is right-skewed. The mean is pulled towards the high values more than the median is, so in such a distribution the mean is larger than the median.

Organization-related FAQ

How to ask a question so we can easily help?

Is practical attendance mandatory?

Can I change which practical afternoon I attend?

Are the sections named Extras part of the exam?

Device with enough battery for exam

About the Discussion Forum