title: Add raw data points to jp.int (sjPlot)
body: <p>For my manuscript, I plotted a lme with an interaction of two continuous variables:</p>
  
  <p>Create data</p>
  
  <pre><code>mydata &lt;- data.frame( SID=sample(1:150,400,replace=TRUE),age=sample(50:70,400,replace=TRUE), sex=sample(c("Male","Female"),200, replace=TRUE),time= seq(0.7, 6.2, length.out=400), Vol =rnorm(400),HCD =rnorm(400))  
  mydata$time &lt;- as.numeric(mydata$time)
  </code></pre>
  
  <p>Run the model:</p>
  
  <pre><code>model &lt;- lme(HCD ~ age*time+sex*time+Vol*time, random=~time|SID, data=mydata)
  </code></pre>
  
  <p>Make plot:</p>
  
  <pre><code>sjp.int(model, swap.pred=T, show.ci=T, mdrt.values="meansd")
  </code></pre>
  
  <p>The reviewer now wants me to add the raw data points to this plot. How can I do this? I tried adding geom_point() referring to mydata, but that is not possible.
  Any ideas?</p>
  
  <p>Update:
  I thought that maybe I could extract the random slope of HCD and then residuals HCD for the covariates and also residuals Vol for the covariates and plot those two to make things easier (then I could plot the points in a 2D plot).
  So, I tried to extract the slopes and use these to fit a linear regression, but the results are different (in the reproducible example less significant, but in my data: the interaction became non-significant (and was significant in the lme)). Not sure what that means or whether this just shows that I should not try to plot it this way.</p>
  
  <p>get the slopes:</p>
  
  <pre><code>  model &lt;- lme(HCD ~ time, random=~time|SID, data=mydata)
    slopes &lt;- rbind(row.names(model$coefficients$random$SID),  model$coef$random$SID[,2])
    slopes2 &lt;- data.frame(matrix(unlist(slopes), nrow=144, byrow=T))
    names(slopes2)[1] &lt;- "SID"
    names(slopes2)[2] &lt;- "slopes"
  </code></pre>
  
  <p>(save the slopes2 and reopen, because somehow R sees it as a factor)</p>
  
  <p>Then create a cross-sectional dataframe and merge the slopes:</p>
  
  <pre><code>  mydata$time2 &lt;- round(mydata$time)
    new &lt;- reshape(mydata,idvar = "SID", timevar="time2", direction="wide")
    newdata &lt;- dplyr::left_join(new, slop, by="SID")
  </code></pre>
  
  <p>The lm:</p>
  
  <pre><code>modelw &lt;- lm(slop$slopes ~ age.1+sex.1+Vol.1, data=newdata)
  </code></pre>
  
  <p>Vol now has a p-value of 0.8 (previously this was 0.14)</p>