MSS MSS - 3 months ago 17
R Question

Change axis and mirror graph to complete non-linear surface in R

I hope you can help me to solve this issue, I've been trying different things but nothing work so far:

I have a 3D graph that is using a squared term (x2) on the x axis (values go from 0 to 100). The original x has positive and negative values (values go from -10 to 10). In x2 and therefore in the X axis of my 3D graph the values are all positive. Thinking that x2=100 is the value obtained from x=-10^2 and x=10^, x2=25 comes from x= -5^2 and x=5^2 and so on. I have only "half" of the graph and I would like to:
1) Have the graph with the original scale going from -10 to 10 on the X axis.
2) Complete the other half of the graph to have the non-linear relationship (i.e. to complete the surface that corresponds from -10 to 0, which I assume should be a mirror of the one I have right now).

enter image description here

Using different colours you can see better the nonlinear relationship, but I didn't include them here to simplify the code.

Since it is not possible to get the negative values back to plot x because the square root will be always positive, I duplicated the data in Excel. I added negative values (values now go from -100 to 100), I made an R list again. This is not a solution because it still has the same scale of x2, but anyways it doesn’t work.

This is how I plot the graph:

Data: https://www.dropbox.com/s/fv943jf35eqtkd8/NSSH.csv?dl=0

link function code:

logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta) days * plogis(eta)^(days-1) *
.Call("logit_mu_eta", eta, PACKAGE = "stats")
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta, name = link),
class = "link-glm")
}


The 3D graph:

library(akima)
x <- NSSH$reLDM
x2<- x^2
y <- NSSH$yr
y2 <-y^2
n <-NSSH$AgeDay1
z <- NSSH$survive
m <- glm(z~x2+y+y2+x2:y+n,family=binomial(link=logexp(NSSH$exposure)))

# interaction
i <- 25
xtemp <- seq(min(x),max(x),length.out=i)
xrange <- rep(xtemp,times=i)
x2temp <- seq(min(0),max(100),length.out=i)
x2range <- rep(x2temp,times=i)
ytemp <- seq(min(y),max(y),length.out=i)
yrange <- rep(ytemp,each=i)
y2temp <- seq(min(y2),max(y2),length.out=i)
y2range <- rep(y2temp,each=i)
ntemp <- rep(mean(n),times=i)
nrange <- rep(ntemp,times=i)

newdata <- data.frame(x2=x2range,y=yrange,y2=y2range,n=nrange)
zhat <- predict(m,newdata=newdata)
NS <- zhat^27
xyz <- interp(x2range,yrange,NS)


quartz()
persp(xyz,
theta = 35, phi = 50,col="blue", border="grey40", ticktype = "detailed", zlim=c(0,1)) -> res2


Is there a way I can copy the "half" graph I have as a “mirror” and put it next to the part I already have and use the original scale from x?

Thanks a lot for your help!

UPDATE:

The 3D graph is perfect!

But when I use the "half graph" to make a contour plot it looks like this:
enter image description here

And now with the new graph it looks like this, I wonder why the origin around 0 next to the value 0.7 (area in the red circle) doesn't look the same as the first contour plot. Do you have any idea? is it possible to fix it? Thanks again.

enter image description here

this is the code of the contour plot:

image(xyz2,col = "white")
contour(xyz2,add=T)

Answer

I think you don't have to worry about the small things with the exception that make X and Y increase and dim(Z) are c(length(X), length(Y)) .

xyz2 <- interp(sqrt(x2range), yrange, NS)   # change scale before interpolate
xyz2$x <- c(rev(xyz2$x)*-1, xyz2$x)         # reverse and combine
xyz2$x[41] <- 1.0E-8                    # because [40] = [41] = 0 (40 is interp's nx value)
xyz2$z <- rbind(apply(xyz2$z, 2, rev), xyz2$z) # reverse and combine

persp(xyz2,xlab="Relative laying date",ylab="Year",zlab="Nest success",
      theta = 35, phi = 50,col="blue", border="grey40", ticktype = "detailed")

enter image description here