# R/TangleInferenceGauss.R In GeoGenetix: Quantification of the effect of geographic versus environmental isolation on genetic differentiation

```TangleInferenceGauss <-
function (data, e, g, ...)
{
nc = ncol(data)
thetaOpt = function(fre, e, g, theta0 = c(0, 10, 10, 1),
type = c(1, 1, 1, 1)) {
cor = cor(fre)
D_E = as.matrix(dist(e)/max(dist(e)))
D_G = as.matrix(dist(g)/max(dist(g)))
n = ncol(cor)
I = diag(n)
type = as.logical(type)
infind = theta0 == Inf
theta0[infind] = 10^8
type[infind] = FALSE
ssc = function(theta) {
th = theta0
th[type] = theta
d = th[1]
be = th[2]
bg = th[3]
g = th[4]
P = D_E/be + D_G/bg
tcor = (1 - d) * exp(-P^g) + d * I
ss = sum((cor - tcor)^2)
return(ss)
}
dssc = function(theta) {
th = theta0
th[type] = theta
d = th[1]
be = th[2]
bg = th[3]
g = th[4]
P = D_E/be + D_G/bg
ds1 = (sum(2 * (cor - (1 - d) * exp(-P^g) - d * I) *
(exp(-P^g) - I)))
ds2 = -2 * (cor - (1 - d) * exp(-P^g) - d * I) *
(1 - d) * P^g * g * D_E * exp(-P^g)/(be^2 * P)
ds2 = (sum(ds2[!I]))
ds3 = -2 * (cor - (1 - d) * exp(-P^g) - d * I) *
(1 - d) * P^g * g * D_G * exp(-P^g)/(bg^2 * P)
ds3 = (sum(ds3[!I]))
ds4 = 2 * (cor - (1 - d) * exp(-P^g) - d * I) * (1 -
d) * P^g * log(P) * exp(-P^g)
ds4 = (sum(ds4[!I]))
dss = c(ds1, ds2, ds3, ds4)[type]
return(dss)
}
lb = c(0, 1e-07, 1e-07, 0)[type]
ub = c(1, Inf, Inf, ifelse(sum(infind) > 0, 2, 1))[type]
op = optim(theta0[type], ssc, dssc, lower = lb, upper = ub,
method = "L-BFGS-B", control = list(factr = 1e-10))
theta = theta0
theta[type] = op\$par
theta[infind] = Inf
names(theta) = c("delta", "beta_E", "beta_G", "gamma")
return(theta)
}
predError = function(data, e, g, theta, pred = 1) {
order = 1:nc
order[c(1, pred)] = c(pred, 1)
data = data[, order]
e = e[order]
g = g[order, ]
D_E = as.matrix(dist(e)/max(dist(e)))
D_G = as.matrix(dist(g)/max(dist(g)))
delta = theta[1]
beta_E = theta[2]
beta_G = theta[3]
gamma = theta[4]
S = (1 - delta) * exp(-(D_E/beta_E + D_G/beta_G)^gamma) +
delta * diag(nc)
x2 = t(as.matrix(data[, -1]))
S11 = S[1, 1]
S12 = S[1, -1]
S21 = S[-1, 1]
S22 = S[-1, -1]
mu1 = mu2 = mean(x2)
x1 = mu1 + S12 %*% solve(S22) %*% (x2 - mu2)
error = sum((x1 - data[, 1])^2)/nrow(data)
return(error)
}
MSE = rep(0, 3)
names(MSE) = c("e+g", "e", "g")
for (model in 1:3) {
g0 = ifelse(model == 2, Inf, 50)
e0 = ifelse(model == 3, Inf, 50)
SE = 0
for (k in 1:nc) {
theta = thetaOpt(data[, -k], e[-k], g[-k, ], theta0 = c(0,
e0, g0, 1))
error = predError(data, e, g, theta, pred = k)
SE = SE + error
}
MSE[model] = SE/nc
}
return(MSE)
}
```

## Try the GeoGenetix package in your browser

Any scripts or data that you put into this service are public.

GeoGenetix documentation built on May 1, 2019, 10:56 p.m.