Simulation studies using Mixedmodels

The most noticeable advantages of fitting mixed-effects models using the MixedModels package for Julia relative to, say, lme4 in R are speed and reliability.

For simple models fit to data sets of modest size, speed is not that much of an issue. If the model can be fit in less time than it takes to organize and plot the data to determine an appropriate model form then it is not too important if one implementation is faster than another.

However, in a simulation study that may involve tens or hundreds of thousands of model fits speed of fitting even simple models can be important. The purpose of this notebook is to show how to perform simulation studies in Julia quickly and efficiently. The emphasis will be on simulation of linear mixed-effects models but many of the techniques and software practices apply more generally.

Overall strategy

It is almost trivial to state that the purpose of a simulation study is to generate a set of simulated values that we can use to assess whether a claim is reasonable. By the way, it is important to realize that a simulation study can only refute a claim by showing that the observed results are inconsistent with that claim. You cannot establish a result through a simulation study. You can only show that your results in a particular case are non inconsistent with the claim.

As with any other experiment it helps to determine up front what the scope of the study should be and what information from the study should be retained. The design space for the study can be very large, even for a simple case such as assessing the characteristics (Type I error and power) of a particular test, as described below. Also the amount of information generated can be very large and the expense of retaining a huge amount of information must be balanced with information loss by summarizing. There is a great temptation to perform too much data reduction, such as recording only the proportion of p-values falling below a threshold in a test. Retaining the simulated p-values themselves provides a richer source of information.

A basic rule in Julia is that, because of the way just-in-time (JIT) compilation works, any substantial computation should take place within a function. Thus the early phases of the simulation study are exercises in designing functions. In the later phases these functions should be wrapped in a package.

A simple example

In Matuschek, Kliegl, Vasishth, Baayen and Bates (2017) we used a simulation study to explore "Balancing Type I Error and Power in Linear Mixed Models". The experimental design consisted of a single binary experimental factor, C, which varied within subjects (S) and within items (I). There are m items and n subjects.

Begin by loading the packages to be used. The MixedModels package is obviously needed to define and fit the models as is the DataFrames package for the data representation. The Distributions package provides the complementary cumulative distribution function (ccdf) to evaluate p-values from test statistics assuming a given reference distribution, and the Gadfly package is used for plotting.


In [1]:
using Distributions, DataFrames, Gadfly, Iterators, MixedModels


WARNING: Method definition ==(Base.Nullable{S}, Base.Nullable{T}) in module Base at nullable.jl:238 overwritten in module NullableArrays at /home/bates/.julia/v0.6/NullableArrays/src/operators.jl:128.

Define a mkdata function to generate a data frame according to the number of subjects, the number of items, and the coding of the levels of the covariates. (My coauthors in that paper overruled me and used a [-0.5, 0.5] encoding. Here I use the [-1, 1] encoding preferred by statisticians.)


In [2]:
function mkdata(m=10, n=30)
    chars = vcat('A':'Z', 'a':'z')
    tuples = collect(product(['N','Y'], ('a':'z')[1:m], chars[1:n]))
    vecs = push!(Any[pool(getindex.(tuples, i)) for i in 1:3], 
        zeros(length(tuples)))
    DataFrame(vecs, [:C, :I, :S, :Y])
end


Out[2]:
mkdata (generic function with 3 methods)

In [3]:
const dat = mkdata()


Out[3]:
CISY
1NaA0.0
2YaA0.0
3NbA0.0
4YbA0.0
5NcA0.0
6YcA0.0
7NdA0.0
8YdA0.0
9NeA0.0
10YeA0.0
11NfA0.0
12YfA0.0
13NgA0.0
14YgA0.0
15NhA0.0
16YhA0.0
17NiA0.0
18YiA0.0
19NjA0.0
20YjA0.0
21NaB0.0
22YaB0.0
23NbB0.0
24YbB0.0
25NcB0.0
26YcB0.0
27NdB0.0
28YdB0.0
29NeB0.0
30YeB0.0
&vellip&vellip&vellip&vellip&vellip

In [4]:
typeof.(dat.columns)  # check the column types


Out[4]:
4-element Array{DataType,1}:
 DataArrays.PooledDataArray{Char,UInt8,1}
 DataArrays.PooledDataArray{Char,UInt8,1}
 DataArrays.PooledDataArray{Char,UInt8,1}
 Array{Float64,1}                        

The response, Y, is initialized to zeros so that it is available to construct LinearMixedModel objects. The const modifier in the assignment indicates that the type of dat will not change, although the actual data values may change (but in this case they won't). Using const, if appropriate, when assigning global variables is a good practice.

The lmm function creates, but does not fit, a LinearMixedModel. It is useful to have a model that has not yet been fit because it is the structure that can then be used to simulate! a response and refit! the model.

The formula for the linear mixed-effects model can be constructed with the @formula macro or using an explicit call to the Formula constructor. For example


In [5]:
srand(1234321)
refit!(simulate!(lmm(@formula(Y ~ 1 + (1|S) +(1|I)), dat), 
        β=[2000.], σ=300., θ=fill(1/3, 2)))


Out[5]:
Linear mixed model fit by maximum likelihood
 Formula: Y ~ 1 + (1 | S) + (1 | I)
     logLik       -2 logLik         AIC            BIC      
 -4.3350137×10³  8.6700273×10³  8.6780273×10³   8.695615×10³

Variance components:
              Column    Variance   Std.Dev.  
 S        (Intercept)   14072.258 118.626547
 I        (Intercept)    4291.301  65.508023
 Residual              101303.824 318.282617
 Number of obs: 600; levels of grouping factors: 30, 10

  Fixed-effects parameters:
             Estimate Std.Error z value P(>|z|)
(Intercept)   1955.61   32.6657 59.8675  <1e-99

or


In [6]:
srand(1234321)
const m0 = refit!(simulate!(lmm(Formula(:Y,:(1+(1|S)+(1|I))), dat),
        β=[2000.], σ=300., θ=fill(1/3, 2)))


Out[6]:
Linear mixed model fit by maximum likelihood
 Formula: Y ~ 1 + (1 | S) + (1 | I)
     logLik       -2 logLik         AIC            BIC      
 -4.3350137×10³  8.6700273×10³  8.6780273×10³   8.695615×10³

Variance components:
              Column    Variance   Std.Dev.  
 S        (Intercept)   14072.258 118.626547
 I        (Intercept)    4291.301  65.508023
 Residual              101303.824 318.282617
 Number of obs: 600; levels of grouping factors: 30, 10

  Fixed-effects parameters:
             Estimate Std.Error z value P(>|z|)
(Intercept)   1955.61   32.6657 59.8675  <1e-99

We obtain the simulated response and fit the alternative model as


In [7]:
const y0 = model_response(m0)
const m1 = refit!(lmm(Formula(:Y, :(1+C+(1|S)+(1|I))), dat), y0)


Out[7]:
Linear mixed model fit by maximum likelihood
 Formula: Y ~ 1 + C + (1 | S) + (1 | I)
     logLik       -2 logLik         AIC            BIC      
 -4.3350099×10³  8.6700198×10³  8.6800198×10³  8.7020044×10³

Variance components:
              Column    Variance   Std.Dev. 
 S        (Intercept)   14072.327 118.62684
 I        (Intercept)    4291.324  65.50820
 Residual              101302.465 318.28048
 Number of obs: 600; levels of grouping factors: 30, 10

  Fixed-effects parameters:
             Estimate Std.Error    z value P(>|z|)
(Intercept)   1956.74   35.1552    55.6601  <1e-99
C: Y         -2.25397   25.9875 -0.0867329  0.9309

The likelihood ratio test statistic for H₀: β₂ = 0 versus Hₐ: β₂ ≠ 0 is the difference in the objective (negative twice the log-likelihood) values for the fitted models.


In [8]:
lrtval = objective(m0) - objective(m1)


Out[8]:
0.0075225519667583285

The p-value, assuming a Chisq(1) reference distribution, is


In [9]:
ccdf(Chisq(1), lrtval)


Out[9]:
0.9308840260438832

Notice that the conclusions from the likelihoood ratio test are essentially the same as those from the Z test on β₂ in the coefficient summary table

  Fixed-effects parameters:
             Estimate Std.Error  z value P(>|z|)
(Intercept)   1955.61   34.7294    56.31  <1e-99
C             5.81124    12.988 0.447432  0.6546

for the alternative model m1. This is because the square of the z value, which is the ratio of the estimate of β₂ to its standard error, is essentially the same as lrtval


In [10]:
abs2(0.447432)


Out[10]:
0.200195394624

and a Chisq(1) distribution is the square of a standard normal distribution.

Generally a hypothesis test based on fitting both the null model and the alternative model to the data (i.e. the likelihood ratio test) is more reliable than one based on fitting only the alternative model and inferring what the change to the null model fit will be (the Z-statistic). In this case they give essentially the same results.

Plotting the data

Although it is tempting to continue with writing simulation functions, this is a good time to examine the responses that have been simulated. What do we expect them to look like?

To make it more convenient to plot the data we copy the contents of y0 vector to the :Y column of the data frame then look for a systematic shift


In [11]:
copy!(dat[:Y], y0)
plot(dat, x = :Y, Geom.density, Guide.xlabel("Simulated responses"))


Out[11]:
Simulated responses -5×10³ -4×10³ -3×10³ -2×10³ -1×10³ 0 1×10³ 2×10³ 3×10³ 4×10³ 5×10³ 6×10³ 7×10³ 8×10³ 9×10³ -4.0×10³ -3.8×10³ -3.6×10³ -3.4×10³ -3.2×10³ -3.0×10³ -2.8×10³ -2.6×10³ -2.4×10³ -2.2×10³ -2.0×10³ -1.8×10³ -1.6×10³ -1.4×10³ -1.2×10³ -1.0×10³ -8.0×10² -6.0×10² -4.0×10² -2.0×10² 0 2.0×10² 4.0×10² 6.0×10² 8.0×10² 1.0×10³ 1.2×10³ 1.4×10³ 1.6×10³ 1.8×10³ 2.0×10³ 2.2×10³ 2.4×10³ 2.6×10³ 2.8×10³ 3.0×10³ 3.2×10³ 3.4×10³ 3.6×10³ 3.8×10³ 4.0×10³ 4.2×10³ 4.4×10³ 4.6×10³ 4.8×10³ 5.0×10³ 5.2×10³ 5.4×10³ 5.6×10³ 5.8×10³ 6.0×10³ 6.2×10³ 6.4×10³ 6.6×10³ 6.8×10³ 7.0×10³ 7.2×10³ 7.4×10³ 7.6×10³ 7.8×10³ 8.0×10³ -5×10³ 0 5×10³ 1×10⁴ -4.0×10³ -3.5×10³ -3.0×10³ -2.5×10³ -2.0×10³ -1.5×10³ -1.0×10³ -5.0×10² 0 5.0×10² 1.0×10³ 1.5×10³ 2.0×10³ 2.5×10³ 3.0×10³ 3.5×10³ 4.0×10³ 4.5×10³ 5.0×10³ 5.5×10³ 6.0×10³ 6.5×10³ 7.0×10³ 7.5×10³ 8.0×10³ -0.0020 -0.0015 -0.0010 -0.0005 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 -0.00150 -0.00145 -0.00140 -0.00135 -0.00130 -0.00125 -0.00120 -0.00115 -0.00110 -0.00105 -0.00100 -0.00095 -0.00090 -0.00085 -0.00080 -0.00075 -0.00070 -0.00065 -0.00060 -0.00055 -0.00050 -0.00045 -0.00040 -0.00035 -0.00030 -0.00025 -0.00020 -0.00015 -0.00010 -0.00005 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 0.00045 0.00050 0.00055 0.00060 0.00065 0.00070 0.00075 0.00080 0.00085 0.00090 0.00095 0.00100 0.00105 0.00110 0.00115 0.00120 0.00125 0.00130 0.00135 0.00140 0.00145 0.00150 0.00155 0.00160 0.00165 0.00170 0.00175 0.00180 0.00185 0.00190 0.00195 0.00200 0.00205 0.00210 0.00215 0.00220 0.00225 0.00230 0.00235 0.00240 0.00245 0.00250 0.00255 0.00260 0.00265 0.00270 0.00275 0.00280 0.00285 0.00290 0.00295 0.00300 -0.002 0.000 0.002 0.004 -0.0015 -0.0014 -0.0013 -0.0012 -0.0011 -0.0010 -0.0009 -0.0008 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0012 0.0013 0.0014 0.0015 0.0016 0.0017 0.0018 0.0019 0.0020 0.0021 0.0022 0.0023 0.0024 0.0025 0.0026 0.0027 0.0028 0.0029 0.0030

In [14]:
plot(dat, x = :Y, Geom.density, Guide.xlabel("Simulated responses"),
    color = :C)


Out[14]:
Simulated responses -5×10³ -4×10³ -3×10³ -2×10³ -1×10³ 0 1×10³ 2×10³ 3×10³ 4×10³ 5×10³ 6×10³ 7×10³ 8×10³ 9×10³ -4.0×10³ -3.8×10³ -3.6×10³ -3.4×10³ -3.2×10³ -3.0×10³ -2.8×10³ -2.6×10³ -2.4×10³ -2.2×10³ -2.0×10³ -1.8×10³ -1.6×10³ -1.4×10³ -1.2×10³ -1.0×10³ -8.0×10² -6.0×10² -4.0×10² -2.0×10² 0 2.0×10² 4.0×10² 6.0×10² 8.0×10² 1.0×10³ 1.2×10³ 1.4×10³ 1.6×10³ 1.8×10³ 2.0×10³ 2.2×10³ 2.4×10³ 2.6×10³ 2.8×10³ 3.0×10³ 3.2×10³ 3.4×10³ 3.6×10³ 3.8×10³ 4.0×10³ 4.2×10³ 4.4×10³ 4.6×10³ 4.8×10³ 5.0×10³ 5.2×10³ 5.4×10³ 5.6×10³ 5.8×10³ 6.0×10³ 6.2×10³ 6.4×10³ 6.6×10³ 6.8×10³ 7.0×10³ 7.2×10³ 7.4×10³ 7.6×10³ 7.8×10³ 8.0×10³ -5×10³ 0 5×10³ 1×10⁴ -4.0×10³ -3.5×10³ -3.0×10³ -2.5×10³ -2.0×10³ -1.5×10³ -1.0×10³ -5.0×10² 0 5.0×10² 1.0×10³ 1.5×10³ 2.0×10³ 2.5×10³ 3.0×10³ 3.5×10³ 4.0×10³ 4.5×10³ 5.0×10³ 5.5×10³ 6.0×10³ 6.5×10³ 7.0×10³ 7.5×10³ 8.0×10³ 0.5 0.0 -0.5 -1.0 1.0 Color -0.0020 -0.0015 -0.0010 -0.0005 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 -0.00150 -0.00145 -0.00140 -0.00135 -0.00130 -0.00125 -0.00120 -0.00115 -0.00110 -0.00105 -0.00100 -0.00095 -0.00090 -0.00085 -0.00080 -0.00075 -0.00070 -0.00065 -0.00060 -0.00055 -0.00050 -0.00045 -0.00040 -0.00035 -0.00030 -0.00025 -0.00020 -0.00015 -0.00010 -0.00005 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 0.00045 0.00050 0.00055 0.00060 0.00065 0.00070 0.00075 0.00080 0.00085 0.00090 0.00095 0.00100 0.00105 0.00110 0.00115 0.00120 0.00125 0.00130 0.00135 0.00140 0.00145 0.00150 0.00155 0.00160 0.00165 0.00170 0.00175 0.00180 0.00185 0.00190 0.00195 0.00200 0.00205 0.00210 0.00215 0.00220 0.00225 0.00230 0.00235 0.00240 0.00245 0.00250 0.00255 0.00260 0.00265 0.00270 0.00275 0.00280 0.00285 0.00290 0.00295 0.00300 -0.002 0.000 0.002 0.004 -0.0015 -0.0014 -0.0013 -0.0012 -0.0011 -0.0010 -0.0009 -0.0008 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0012 0.0013 0.0014 0.0015 0.0016 0.0017 0.0018 0.0019 0.0020 0.0021 0.0022 0.0023 0.0024 0.0025 0.0026 0.0027 0.0028 0.0029 0.0030

There is no evidence of a systematic shift according to the level of the condition, C, which is what we would expect because we simulated the data from the null model.

Next consider the effect of the subjects and items


In [15]:
plot(dat, x = :S, y = :Y, Geom.boxplot, Guide.xlabel("Subject"),
    Guide.ylabel("Simulated response"))


Out[15]:
Subject A B C D E F G H I J K L M N O P Q R S T U V W X Y Z a b c d -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 -1500 -1400 -1300 -1200 -1100 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 -2000 0 2000 4000 6000 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 Simulated response

In [16]:
plot(dat, x = :I, y = :Y, Geom.boxplot, Guide.xlabel("Item"),
    Guide.ylabel("Simulated response"))


Out[16]:
Item a b c d e f g h i j -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 -1500 -1400 -1300 -1200 -1100 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 -2000 0 2000 4000 6000 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 Simulated response

In [17]:
plot(dat, x = :C, y = :Y, xgroup=:I, Geom.subplot_grid(Geom.boxplot),
    Guide.ylabel("Simulated response"))


Out[17]:
C by I j i h g f e d c b a -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0