3

I am trying to plot an interaction term from a linear mixed effects model using the effects plot. See example below: 

library(nlme)
fitA <- lme(PEE ~ Pupper*max_depth,
  random=~1 + Pupper|ref, data=m4,
  cor=corAR1(), method="ML")

Pupper is a continuous variable and max_depth is a factorial variable (5 levels - 400,500,600,700,800).

When I plot model output in an effects plot with an interaction term I am able to show how the relationship between PEE and Pupper changes according to the different factors levels of 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)

plot1

However, sometimes when I plot a similar model the effects plot will change how I want the interaction to be displayed (see below). In the effects plot below, the plot 'decides' to display the relationship between PEE and max_depth and how this changes according to some arbitrary division of the continuous variable Pupper (labeled Plower in 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)

plot2

However, I am interested in how the factor max_depth influences the relationship between PEE and Plower (similar to the first figure above). I cannot work out why the effect function displays the same interaction term in two different ways. I would love to know how to control the way an interaction term is represented in an effect plot as this problem keeps popping its ugly head up.

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")
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
mal1981
  • 31
  • 4
  • 3
    Please provide a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input data so we can run the code and see the results ourselves. Even if you don't have enough rep to inline images, you can upload them to imgur.com and include the link. A higher rep user can include them for you. – MrFlick Jul 20 '15 at 21:04
  • Thanks MrFlick. I have given a sample of input data to run the code. Hopefully this helps. Cheers – mal1981 Jul 20 '15 at 22:05

1 Answers1

1

The data subset you've supplied are insufficient to fit the models but I think that I can answer your question. The selection of the variable on the horizontal axis is controlled by the x.var argument to plot.eff(); from ?plot.eff:

x.var: the index (number) or quoted name of the covariate or factor to place on the horizontal axis of each panel of the effect plot. The default is the predictor with the largest number of levels or values.

In turn, the covariate values at which an effect is evaluated is controlled by the xlevels argument to Effect(), which in your case is called by 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, then the number and values of levels for any numeric predictor is determined by grid.pretty. If xlevels=n is an integer, then each numeric predictor is represented by n equally spaced levels. More generally, xlevels can be a named list of values at which to set each numeric predictor. For example, xlevels=list(x1=c(2, 4, 7), x2=5) would use the values 2, 4 and 7 for the levels of x1, 5 equally spaced levels for the levels of x2, and use the default for any other numeric predictors. If partial residuals are computed, then the focal predictor that is to appear on the horizontal axis of an effect plot is evaluated at 100 equally spaced values along its full range, and, by default, other numeric predictors are evaluated at the quantiles specified in the quantiles argument, unless their values are given explicitly in xlevels.

Curiously, the second plot that you show appears to include a variable, k150, that doesn't appear in the model.

I hope this helps,

John

John Fox
  • 284
  • 1
  • 3