Control how the interaction graph is displayed in action
I am trying to construct an interaction term from a linear mixed effects model using an effects graph. See example below:
library(nlme)
fitA <- lme(PEE ~ Pupper*max_depth,
random=~1 + Pupper|ref, data=m4,
cor=corAR1(), method="ML")
Pupper
- continuous variable, and max_depth
- factorial variable (5 levels - 400 500 600 700 800).
When I plot the output of the model on a graph of effects with the term interaction, I can show how the relationship between PEE
and changes Pupper
depending on different levels of factors max_depth
:
library(effects)
plot(effect("Pupper*max_depth",fitA),
xlab=expression(paste("d"[-5]," P"[upper]," (m"^" -1",")")),
ylab=expression(paste("PEE rate (h"^" -1",")")),
factor.names=FALSE, layout=c(5,1), alternating=FALSE,
main="A", ticks.x=list(Pupper=list(at=seq(0.4,1.0,0.2))),
##layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
rotx=45, more=FALSE, grid=FALSE, lwd=1)
However, sometimes when I draw a model like this, the effects graph changes how the interaction is displayed (see below). In the effect plot below, the plot "decides" to display the relationship between PEE
and max_depth
and how it changes according to some arbitrary division of a continuous variable Pupper
(labeled Plower
in the code and 'k150' in the plot shown below):
fitB <- lme(PEE ~ Plower*max_depth,
random=~1 + Plower|ref, data=m4,
cor=corAR1(), method="ML")
plot(effect("Plower*max_depth",fitB),
xlab=expression(paste("d"[-5]," P"[lower]," (m"^" -1",")")),
ylab=expression(paste("PEE rate (h"^" -1",")")),
factor.names=FALSE, layout=c(5,1), alternating=FALSE,
main="A", ticks.x=list(Plower=list(at=seq(0.4,1.0,0.2))),
##layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
rotx=45, more=FALSE, grid=FALSE, lwd=1)
However, I am interested in how the factor max_depth
affects the relationship between PEE
and Plower
(as in the first picture above). I can't figure out why the effect function displays the same interaction term in two different ways. I would like to know how to control how the interaction term is presented in the effects plot, as this issue continues to cause an ugly head.
Below is a subset of my dataset:
structure(list(ref = c("2012-3corrige", "2011-28", "2011-26",
"2011-21", "2011-21", "2013-7", "2012-1corrige", "2012-6corrige",
"2013-4", "2011-21", "2013-10", "2013-4", "2013-13", "2011-26",
"2013-11", "2012-3corrige", "2013-1", "2012-14corrige", "2013-1",
"2011-27", "2012-6corrige", "2011-18", "2011-26", "2010-18",
"2012-14corrige", "2011-21", "2013-6", "2013-11", "2011-27",
"2011-18", "2012-16corrige", "2013-5", "2013-13", "2011-21",
"2012-14corrige", "2013-5", "2013-18", "2012-16corrige", "2011-28",
"2010-18", "2011-21", "2013-2", "2012-2corrige", "2013-4", "2013-5",
"2013-11", "2011-21", "2013-6", "2011-28", "2013-6", "2010-18",
"2011-21", "2013-18", "2011-16", "2012-11corrige", "2011-28",
"2011-27", "2012-3corrige", "2012-2corrige", "2013-3", "2012-1corrige",
"2012-14corrige", "2012-14corrige", "2013-10", "2012-6corrige",
"2010-18", "2012-11corrige", "2013-7", "2013-2", "2012-16corrige",
"2013-1", "2013-18", "2012-16corrige", "2013-6", "2012-4corrige",
"2013-4", "2013-10", "2013-3", "2013-2", "2011-16", "2012-1corrige",
"2011-21", "2011-21", "2013-18", "2013-3", "2011-26", "2010-18",
"2013-13", "2012-6corrige", "2013-3", "2012-16corrige", "2012-15corrige",
"2011-28", "2012-6corrige", "2012-6corrige", "2012-11corrige",
"2013-1", "2013-11", "2012-11corrige", "2013-6"), Pupper = c(0.861958207287982,
0.824829924556841, 0.958739109455748, 0.935401831656677, 0.955566680038604,
0.948368978826279, 0.745071680369673, 0.827539122942233, 0.726448658429027,
0.943103302931338, 0.858445846226439, 0.784802718309937, 0.881010495586365,
0.911770168408684, 0.90971638692581, 1.02155421458351, 0.851778844536538,
1.1553118943962, 0.887452083213511, 0.8218157295485, 0.871777265131409,
0.829892474962871, 1.01579427707254, 0.715539162683171, 1.12624787680155,
0.713105471394893, 0.802478037082636, 0.773243110590944, 0.762028205952159,
0.785089166910358, 0.844285844170484, 0.887514023676371, 0.870367623478723,
0.820303824472643, 0.636099278958915, 0.953776661488422, 0.816485694234068,
0.861493535070196, 0.787945463425822, 0.918041865421543, 0.877275056321815,
0.624152855209897, 0.971197595182818, 0.769613695304339, 0.941443459091764,
0.929070549770906, 1.031203743205, 0.692597025693873, 0.846978945035432,
0.72446749179426, 0.541564092852052, 0.744921803502444, 0.917786983273715,
0.702051561892398, 0.975310403563878, 0.808819367281032, 0.858040403089116,
0.741495941398947, 0.698143566897239, 0.979366380200314, 0.992046903013047,
0.995331870590213, 0.804437082665078, 0.8307554779262, 0.878549524310762,
0.654725061889849, 0.93024953667308, 0.654611447094126, 0.689696271315618,
0.77302453480932, 0.916283427766758, 0.894114399839305, 0.840205756601608,
0.767235548359607, 0.831544386468135, 0.685089269122402, 0.860269828471148,
0.895228365651283, 0.785946885397904, 0.812567650516969, 0.797256286689962,
0.800979891549511, 0.684467773772683, 0.846228645225391, 0.801015938251751,
0.964375424821682, 0.783654311543043, 0.951249150678552, 0.847095453102345,
0.782862048298847, 0.897798965949478, 0.79591714811698, 0.954852044385237,
0.885914708711347, 0.789575506205708, 1.10814372714012, 0.875651148193922,
0.851523408695002, 0.963324355206144, 0.795071091161036), Plower = c(0.705132769215998,
0.667302197075824, 0.629978835623335, 0.632452896796802, 0.641619045851976,
0.634150350206216, 0.521875889886134, 0.69048678481199, 0.620155894379255,
0.72673011955379, 0.644805071164551, 0.691418831100224, 0.561990510002912,
0.702502669034076, 0.5885329988032, 1.06019049650942, 0.610499795249761,
0.863589408611907, 0.671649710290516, 0.7008237216939, 0.613070958372683,
0.52121652570373, 0.743727100487806, 0.619214556245787, 1.0217109832694,
0.653199816289013, 0.653255947797901, 0.629436452185155, 0.621227279933305,
0.581484689776476, 0.605084016204913, 0.670828674932066, 0.694246594037978,
0.732994239783339, 0.728155423409921, 0.657673931367209, 0.681945582710071,
0.656113353702447, 0.55299186250794, 0.589741939797023, 0.760512984767519,
0.550684422635445, 0.888934443277143, 0.615143614667881, 0.736486026117717,
0.616589579139919, 0.640405340389975, 0.618517043688639, 0.612475849864031,
0.681245183469212, 0.642477842246546, 0.683125578173995, 0.636702442275825,
0.568741300299764, 0.681781639762194, 0.58956858049858, 0.697614984548545,
0.773372843818268, 0.599073358520381, 0.653548263966276, 0.846172099647715,
0.946825538132655, 0.635504629303462, 0.61980005655224, 0.594418483337567,
0.610786378368084, 0.809715477703094, 0.596886365511337, 0.601414998150196,
0.680138336678131, 0.672368946338244, 0.693205917779446, 0.736742863266092,
0.636678882954351, 0.611664395999418, 0.630585706572337, 0.6554630468205,
0.617362130357864, 0.615793526002561, 0.688748462389895, 0.733587834625896,
0.715468455706547, 0.695921451506322, 0.649384323802169, 0.654685675216232,
0.675344356606317, 0.617759392212426, 0.620895052860519, 0.652138022200822,
0.638494322605926, 0.798451637031189, 0.884450865414997, 0.895823643358446,
0.661857655055493, 0.743487278528243, 0.88302132854573, 0.660494764046872,
0.638155299450374, 0.515272975866573, 0.636047692132176), PEE = c(1.49625537031302,
0.579708304983786, 1.09755665230733, 0.79999598579792, 0.366971323405136,
0.534519852186464, 0.892172302701565, 0.764300080784289, 0.162161584516302,
1.05547854644252, 1.75994722974226, 0.502090519778478, 0.813556191910923,
0.72071830101183, 0.124737712452804, 0.24096278221797, 2.11763754191128,
0.0970872140009704, 0.214668546888839, 0.997227687637828, 0.449413221941473,
0.445533213208998, 0.719422276286327, 0.311417756794472, 0.0799998795735146,
0.836011454943211, 0.217381544231536, 0.131863894834852, 1.1881086717854,
0.562478645146688, 2.13755423670725, 0.260855398500945, 0.769228719926564,
0.792077729637109, 0.46631964160662, 0.727219913477435, 0.234599414042217,
1.24135448496734, 1.73912823566756, 0.437157161658986, 1.18491673172712,
1.57894265236866, 0.325374033367569, 0.133629870488068, 0.260855348600893,
0.279711624960117, 1.04650548272943, 1.13790951622142, 0.512819441159373,
2.51301278595252, 0.948078086013639, 0.183485515251287, 0.521708195407091,
0.834371581292816, 0.907354231373586, 0.263735732207326, 0.94736384877553,
0.865382874911045, 0.162378076290205, 1.80685106084338, 1.07131194190618,
1.20567188480079, 1.01009693910426, 0.352933736835024, 0.315767760495837,
0.901500577761704, 1.08481956174672, 0.553504151294972, 1.81542854475215,
2.23136249824668, 0.14018646847029, 0.58250584995577, 1.74754206600435,
0.404021461283339, 1.0507621403718, 3.1578487818322, 1.23592560921063,
0.428569941841892, 2.59927399028359, 0.462953155929056, 1.60334686111772,
0.610996390428844, 0.93749693604517, 0.374210416022193, 0.596024133599949,
1.07142175386991, 0.233917466116807, 0.773637607411201, 0.733915433670645,
0.693195231932592, 0.699678270730694, 0.75104196328333, 1.1707299559812,
0.376558572007052, 1.5725384212365, 0.17424722659426, 0.481925512179189,
0.127383975172354, 0.449990814000021, 0.828701950628209), max_depth = structure(c(5L,
1L, 5L, 5L, 3L, 2L, 3L, 2L, 2L, 2L, 2L, 5L, 3L, 3L, 3L, 4L, 5L,
1L, 5L, 1L, 5L, 2L, 4L, 3L, 3L, 3L, 3L, 2L, 4L, 2L, 2L, 5L, 2L,
2L, 4L, 5L, 4L, 2L, 1L, 5L, 2L, 1L, 4L, 4L, 4L, 3L, 4L, 1L, 2L,
2L, 5L, 5L, 4L, 1L, 3L, 1L, 3L, 2L, 4L, 1L, 5L, 1L, 5L, 4L, 4L,
1L, 2L, 4L, 1L, 1L, 4L, 4L, 2L, 2L, 2L, 4L, 2L, 5L, 1L, 2L, 1L,
3L, 4L, 3L, 1L, 5L, 1L, 4L, 4L, 3L, 3L, 3L, 1L, 3L, 1L, 1L, 4L,
1L, 4L, 2L), .Label = c("400", "500", "600", "700", "800"), class = "factor"),
fangle = structure(c(2L, 3L, 1L, 4L, 3L, 4L, 3L, 3L, 3L,
3L, 2L, 4L, 2L, 4L, 2L, 4L, 4L, 4L, 4L, 1L, 4L, 2L, 2L, 4L,
3L, 4L, 4L, 2L, 2L, 2L, 4L, 3L, 3L, 4L, 1L, 3L, 4L, 2L, 3L,
3L, 4L, 3L, 2L, 3L, 3L, 3L, 2L, 1L, 2L, 2L, 3L, 3L, 2L, 2L,
2L, 3L, 2L, 4L, 1L, 3L, 2L, 4L, 3L, 3L, 1L, 2L, 2L, 4L, 2L,
1L, 4L, 3L, 4L, 2L, 4L, 3L, 4L, 4L, 3L, 2L, 3L, 3L, 2L, 2L,
4L, 2L, 4L, 4L, 3L, 4L, 3L, 4L, 1L, 4L, 4L, 4L, 2L, 2L, 3L,
3L), .Label = c("0", "20", "40", "60"), class = "factor")), .Names = c("ref",
"Pupper", "Plower", "PEE", "max_depth", "fangle"), row.names = c(26297L,
18163L, 13367L, 10757L, 10813L, 43605L, 22984L, 27608L, 39808L,
11220L, 32882L, 39987L, 35960L, 13719L, 34174L, 25877L, 31747L,
19402L, 31394L, 14990L, 28537L, 9023L, 13684L, 1781L, 19411L,
9964L, 41834L, 33800L, 15277L, 8673L, 21864L, 40681L, 35425L,
11590L, 19901L, 40867L, 36845L, 21698L, 18302L, 470L, 11459L,
37414L, 24555L, 40026L, 40578L, 33627L, 9525L, 41816L, 17695L,
42057L, 294L, 9972L, 37137L, 8304L, 19086L, 15817L, 15351L, 26097L,
24896L, 39059L, 23703L, 20110L, 19937L, 32121L, 28556L, 13L,
19157L, 42865L, 37922L, 21887L, 31638L, 37008L, 21905L, 41848L,
26621L, 39864L, 32870L, 39107L, 37721L, 7969L, 23826L, 11903L,
12024L, 36500L, 38488L, 13287L, 462L, 36245L, 28096L, 38611L,
21500L, 20565L, 17140L, 27772L, 27773L, 18897L, 30992L, 34564L,
18553L, 41312L), class = "data.frame")
source to share
The data subset you supplied is not enough for models, but I think I can answer your question. Selection of the variable on the horizontal axis is controlled by the argument x.var
to plot.eff()
; from ?plot.eff
:
x.var
: index (number) or cited name of covariate or coefficient place on the horizontal axis of each panel of the effect graph. default is the predictor with the most levels or values.
In turn, the covariance values at which the effect is evaluated are controlled by the xlevels
before argument Effect()
, which in your case is called Effect()
; from ?Effect
:
xlevels
: This argument is used to set the number of levels for any focal predictor that is not a factor. If xlevels = NULL, the default is the number and level values for any numeric predictor as defined by grid.pretty. If xlevels = n is an integer, then each numeric predictor is represented by n equal bins. More typically, xlevels can be a named list of values, with a numeric predictor for each. For example, xlevels = list (x1 = c (2, 4, 7), x2 = 5) will use values 2, 4 and 7 for levels x1, 5 equally spaced levels for levels x2, and use the default for any other numeric predictors ... If partial residuals are calculated, then focus is a predictor,which should appear on the horizontal axis of the graph, the effect is evaluated at 100 equally spaced values across its entire range and, by default, other numeric predictors are evaluated in quantiles of the quantiles specified in the argument, unless their values are explicitly specified in xlevels.
Curiously, the second plot you are showing includes a variable k150
that is not displayed in the model.
I hope this helps,
John
source to share