1 Introduction
The purpose of this section is to give users an easy introduction,
starting from simple linear regression. Based on simulations we will
explain the main functions, the input and the output. First we load the
LMMsolver and ggplot2 packages:
1.1 Linear Regression
We will start with a simple example where the true function is linear in variable \(x\):
Using this function we simulate data and add normal distributed noise:
set.seed(2016)
n <- 25
x <- seq(0, 1, length = n)
sigma2e <- 0.04
y <- f1(x) + rnorm(n, sd = sqrt(sigma2e))
dat1 <- data.frame(x = x, y = y)We can fit the data using the LMMsolve function:
We can make predictions using the predict()
function:
newdat <- data.frame(x = seq(0, 1, length = 300))
pred1 <- predict(obj1, newdata = newdat, se.fit = TRUE)
# adding the true values for comparison
pred1$y_true <- f1(pred1$x)Note that for this linear model we could have used the standard
lm() function, which will give the same result.
The following plot gives the simulated data with the predictions, and pointwise standard-error bands. The true value is plotted as dashed red line.
ggplot(data = dat1, aes(x = x, y = y)) +
geom_point(col = "black", size = 1.5) +
geom_line(data = pred1, aes(y=y_true), color = "red",
linewidth = 1, linetype = "dashed") +
geom_line(data = pred1, aes(y = ypred), color = "blue", linewidth = 1) +
geom_ribbon(data = pred1, aes(x=x,ymin = ypred-2*se, ymax = ypred+2*se),
alpha = 0.2, inherit.aes = FALSE) +
theme_bw()
1.2 Fitting a non-linear function
In this section we will use the following non-linear function for the simulations:
The simulated data is generated by the following code
set.seed(12)
n <- 150
x <- seq(0, 1, length = n)
sigma2e <- 0.04
y <- f2(x) + rnorm(n, sd = sqrt(sigma2e))
dat2 <- data.frame(x, y)We can use the spline argument to fit the non-linear
trend:
where spl1D(x, nseg = 50) defines a mixed model
P-splines with 50 segments.
The model fit can be summarized in terms of effective dimensions:
summary(obj2)
#> Table with effective dimensions and penalties:
#>
#> Term Effective Model Nominal Ratio Penalty
#> (Intercept) 1.00 1 1 1.00 0.0
#> lin(x) 1.00 1 1 1.00 0.0
#> s(x) 11.28 53 51 0.22 0.0
#> residual 136.72 150 148 0.92 30.3
#>
#> Total Effective Dimension: 150The intercept and the slope lin(x) define the linear (or
fixed) part of the model, the non-linear (or random) part is defined by
s(x), with effective dimension 11.28.
Making predictions on the interval \([0,1]\) and plotting can be done in the same way as for the linear regression example:
newdat <- data.frame(x = seq(0, 1, length = 300))
pred2 <- predict(obj2, newdata = newdat, se.fit = TRUE)
pred2$y_true <- f2(pred2$x)
ggplot(data = dat2, aes(x = x, y = y)) +
geom_point(col = "black", size = 1.5) +
geom_line(data = pred2, aes(y = y_true), color = "red",
linewidth = 1, linetype ="dashed") +
geom_line(data = pred2, aes(y = ypred), color = "blue", linewidth = 1) +
geom_ribbon(data= pred2, aes(x=x, ymin = ypred-2*se, ymax = ypred+2*se),
alpha=0.2, inherit.aes = FALSE) +
theme_bw() 
1.3 Non-Gaussian data
The LMMsolver package can also be used for non-gaussian
data, using the family argument, with default
family = gaussian. As an example we use count data using
the Poisson distribution, defined by \[
Pr(X=k) = \frac{\lambda^k e^{-\lambda}}{k!} \;,
\] with parameter \(\lambda >
0\) and \(k\) is the number of
occurrences. More general, the value of the parameter \(lambda\) can depend on another variable
time \(x\), for example time. Here we
will assume that \(x\) is defined on
the interval \([0,1]\) and defined
by:
\[ \lambda(x) = 4 + 3x + 4 \sin(7 x) \] Using this function we simulate the following data
set.seed(1234)
n <- 150
x <- seq(0, 1, length=n)
fun_lambda <- function(x) { 4 + 3*x + 4*sin(7*x) }
x <- seq(0, 1, length = n)
y <- rpois(n = n, lambda = fun_lambda(x))
dat3 <- data.frame(x = x, y = y)Now we fit the data with the argument
family = poisson():
obj3 <- LMMsolve(fixed = y ~ 1,
spline = ~spl1D(x, nseg = 50),
family = poisson(),
data = dat3)
summary(obj3)
#> Table with effective dimensions and penalties:
#>
#> Term Effective Model Nominal Ratio Penalty
#> (Intercept) 1.00 1 1 1.00 0
#> lin(x) 1.00 1 1 1.00 0
#> s(x) 6.54 53 51 0.13 0
#> residual 141.46 150 148 0.96 1
#>
#> Total Effective Dimension: 150Making predictions and plotting the data is similar to the Gaussian data we showed before:
newdat <- data.frame(x = seq(0, 1, length = 300))
pred3 <- predict(obj3, newdata = newdat, se.fit = TRUE)
pred3$y_true <- fun_lambda(pred3$x)
ggplot(data = dat3, aes(x = x, y = y)) +
geom_point(col = "black", size = 1.5) +
geom_line(data = pred3, aes(y = y_true), color = "red",
linewidth = 1, linetype ="dashed") +
geom_line(data = pred3, aes(y = ypred), color = "blue", linewidth = 1) +
geom_ribbon(data= pred3, aes(x=x, ymin = ypred-2*se, ymax = ypred+2*se),
alpha=0.2, inherit.aes = FALSE) +
theme_bw() 
1.4 Smoothing combining two experiments
In this section we will give a bit more complicated example, to show
some further options of LMMsolver. Suppose there are two
experiments, A and B, with the same true unobserved non-linear function
f2(x) as defined before.
The simulated data is given by the following code:
set.seed(1234)
nA <- 50
nB <- 100
mu_A <- 0.10
mu_B <- -0.10
sigma2e_A <- 0.04
sigma2e_B <- 0.10
x1 <- runif(n = nA)
x2 <- runif(n = nB)
y1 <- f2(x1) + rnorm(nA, sd = sqrt(sigma2e_A)) + mu_A
y2 <- f2(x2) + rnorm(nB, sd = sqrt(sigma2e_B)) + mu_B
Experiment <- as.factor(c(rep("A", nA), rep("B", nB)))
dat4 <- data.frame(x = c(x1, x2), y = c(y1,y2), Experiment = Experiment)Before analyzing the data in further detail a boxplot gives some insight:
ggplot(dat4, aes(x = Experiment, y = y, fill = Experiment)) +
geom_boxplot() +
geom_point(position = position_jitterdodge(), alpha = 0.3) 
Comparing the two experiments we can see that:
- There is a clear difference in mean/median between the two
experiments. This can be corrected for by adding the argument
random = ~Experiment. - The variance in experiment A is smaller than in B. This implies that
is important to allow for heterogeneous variances which can be modelled
by defining
residual = ~Experiment.
The model in LMMsolve() is given by:
obj4 <- LMMsolve(fixed= y ~ 1,
spline = ~spl1D(x, nseg = 50),
random = ~Experiment,
residual = ~Experiment,
data = dat4)The table of effective dimensions is given by:
summary(obj4)
#> Table with effective dimensions and penalties:
#>
#> Term Effective Model Nominal Ratio Penalty
#> (Intercept) 1.00 1 1 1.00 0.00
#> lin(x) 1.00 1 1 1.00 0.00
#> Experiment 0.93 2 1 0.93 77.98
#> s(x) 7.89 53 51 0.15 0.00
#> Experiment_A!R 43.65 50 50 0.87 32.15
#> Experiment_B!R 95.52 100 100 0.96 9.01
#>
#> Total Effective Dimension: 150And making the predictions:
newdat <- data.frame(x=seq(0, 1, length = 300))
pred4 <- predict(obj4, newdata = newdat, se.fit = TRUE)
pred4$y_true <- f2(pred4$x)
ggplot(data = dat4, aes(x = x, y = y, colour = Experiment)) +
geom_point(size = 1.5) +
geom_line(data = pred4, aes(y = y_true), color="red",
linewidth = 1, linetype = "dashed") +
geom_line(data = pred4, aes(y = ypred), color = "blue", linewidth = 1) +
geom_ribbon(data = pred4, aes(x = x,ymin = ypred-2*se, ymax = ypred+2*se),
alpha = 0.2, inherit.aes = FALSE) +
theme_bw()
The estimated random effects for Experiment can be obtained using the
coef() function:
The sum of the effects is equal to zero, as expected for a standard random term.
