MSS - 6 months ago 36

R Question

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).

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:

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.

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")
```