The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
library(fastLaplace)
if(requireNamespace("mgcv")){
sigma2 = 1
phi = 0.2
beta.true = c(1,1)
n = 400
n.pred = 100
coords.all<- matrix(runif((n+n.pred)*2),ncol=2,nrow=n+n.pred) # simulate data locations
X.all <- matrix(runif((n+n.pred)*2),ncol=2,nrow=(n+n.pred))
suppressMessages(require(fields))
dist.all <- fields::rdist(coords.all,coords.all) # compute distance matrix
matern <- function(phi,mat.dist){
K = (1+sqrt(5)/phi*mat.dist+ 5/(3*phi^2)*mat.dist^2)*exp(-sqrt(5)/phi*mat.dist)
return(K)
} # matern 2.5
V.all <- sigma2*matern(phi,dist.all) # compute covariance matrix
set.seed(1)
r.e.all <- mgcv::rmvn(1,rep(0,nrow(coords.all)),V.all) # simulate random effects
pi.all <- X.all%*%beta.true + r.e.all # linear model
p.all <- exp(pi.all)/(1+exp(pi.all)) # compute the probability of z = 1 for binary process
Y.all <- sapply(p.all, function(x) sample(0:1, 1, prob = c(1-x, x))) # simulate binary observations
} else{
stop("Package \"mgcv\" needed to generate a simulated data set")
}
#> Loading required namespace: mgcv
Y <- as.matrix(Y.all[1:n],nrow = n)
X <- X.all[1:n,]
coords <- coords.all[1:n,]
data <- data.frame(cbind(Y,X))
colnames(data) = c("Y","X1","X2")
mod.glm <- glm(Y~-1+X1+X2,family="binomial",data=data)
mod.glm.esp <- predict(mod.glm,data, type="response")
mod.glm.s2 <- var(Y - mod.glm.esp)
mod.glm.phi <- 0.1*max(dist(coords))
startinit <- c(mod.glm$coef,log(mod.glm.s2),log(mod.glm.phi))
names(startinit) <- c("X1","X2","logsigma2","logphi")
result.bin <- fsglmm(Y~-1+X1+X2, kappa=2.5, inits = startinit, data = data,coords = coords, family = "binomial", ntrial = 1, offset = NA,method.optim = "CG", method.integrate = "NR", control = list(maxit=1000,ndeps=rep(1e-2,4),reltol=0.01),rank = 50)
result.bin$summary
#> Maximum likelihood estimation
#>
#> Call:
#> bbmle::mle2(minuslogl = nlikSGLMM, start = inits, method = method.optim,
#> data = list(Y = Y, X = X, mat.dist = mat.dist, ntrial = ntrial,
#> family = family, method = method.integrate, kappa = kappa,
#> offset = offset, U1 = U1, rank = rank), vecpar = TRUE,
#> skip.hessian = TRUE, control = control)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> X1 1.89994 0.44509 NA NA
#> X2 0.48805 0.41668 NA NA
#> logsigma2 -0.29190 0.53395 NA NA
#> logphi -1.63619 0.37567 NA NA
#>
#> -2 log L: 401.2075
X.pred <- X.all[(n+1):(n+n.pred),]
coords.pred <- coords.all[(n+1):(n+n.pred),]
pred.sglmm(result.bin,data=X.pred,coords=coords.pred)
#> [1] 0.7173744 0.9011669 0.8174600 0.7143686 0.6964783 0.9054776 0.7366131
#> [8] 0.7031458 0.9269817 0.9051901 0.8931835 0.4344861 0.8930961 0.7268835
#> [15] 0.6794283 0.4742220 0.8197480 0.8347397 0.9134523 0.3202977 0.8738509
#> [22] 0.6427999 0.7457062 0.5924379 0.9058852 0.8466672 0.6078437 0.7995875
#> [29] 0.5693686 0.6750436 0.9437651 0.7544525 0.7625658 0.9218002 0.6513120
#> [36] 0.8064205 0.9295038 0.9182620 0.6105966 0.8203830 0.5307209 0.8577690
#> [43] 0.7125925 0.7156856 0.8385185 0.9215306 0.9043965 0.6113454 0.6571189
#> [50] 0.5723820 0.8917714 0.7052825 0.6184694 0.7408175 0.8327119 0.7398930
#> [57] 0.9234421 0.6905265 0.7738696 0.9242138 0.9136343 0.8689911 0.7954075
#> [64] 0.4641804 0.8459859 0.4840722 0.8223483 0.5877907 0.8536701 0.8729155
#> [71] 0.8922757 0.7047227 0.7456766 0.9367009 0.8174338 0.6571998 0.9012756
#> [78] 0.7099313 0.8327465 0.6874149 0.9082839 0.9439709 0.5455138 0.9463924
#> [85] 0.9189615 0.5550945 0.8069886 0.9413874 0.5549568 0.7051886 0.9503762
#> [92] 0.9232483 0.8562200 0.6091460 0.7334006 0.9171320 0.6331991 0.5884551
#> [99] 0.8798014 0.8443174
if(requireNamespace("mgcv")){
sigma2 = 1
phi = 0.2
prec = 2
beta.true = c(1,1)
n = 400
n.pred = 100
coords.all<- matrix(runif((n+n.pred)*2),ncol=2,nrow=n+n.pred) # simulate data locations
X.all <- matrix(runif((n+n.pred)*2),ncol=2,nrow=(n+n.pred))
suppressMessages(require(fields))
dist.all <- fields::rdist(coords.all,coords.all) # compute distance matrix
matern <- function(phi,mat.dist){
K = (1+sqrt(5)/phi*mat.dist+ 5/(3*phi^2)*mat.dist^2)*exp(-sqrt(5)/phi*mat.dist)
return(K)
} # matern 2.5
V.all <- sigma2*matern(phi,dist.all) # compute covariance matrix
set.seed(1)
r.e.all <- mgcv::rmvn(1,rep(0,nrow(coords.all)),V.all) # simulate random effects
mu.all <- exp(X.all%*%beta.true+r.e.all)
Y.all <- rnbinom(length(mu.all), mu=mu.all,size=prec)
} else {
stop("Package \"mgcv\" needed to generate a simulated data set")
}
if(requireNamespace("MASS")){
Y <- as.matrix(Y.all[1:n],nrow = n)
X <- X.all[1:n,]
coords <- coords.all[1:n,]
data <- data.frame(cbind(Y,X))
colnames(data) = c("Y","X1","X2")
mod.glm <- MASS::glm.nb(Y~-1+X1+X2,data=data)
mod.glm.esp <- predict(mod.glm, data)
mod.glm.s2 <- var( log(Y+1) - mod.glm.esp)
mod.glm.phi <- 0.1*max(dist(coords))
mod.glm.prec <- mod.glm$theta
startinit <- c(mod.glm$coef,log(mod.glm.s2),log(mod.glm.phi),log(mod.glm.prec))
names(startinit) <- c("X1","X2","logsigma2","logphi","logprec")
} else {
stop("Package \"MASS\" needed to set the initial parameters")
}
result.nb <- fsglmm(Y~-1+X1+X2, kappa=2.5, inits = startinit, data = data,coords = coords, family = "negative.binomial", offset = NA,method.optim = "CG", method.integrate = "NR", control = list(maxit=1000,ndeps=rep(1e-2,5),reltol=0.01),rank = 50)
result.nb$summary
#> Maximum likelihood estimation
#>
#> Call:
#> bbmle::mle2(minuslogl = nlikSGLMM, start = inits, method = method.optim,
#> data = list(Y = Y, X = X, mat.dist = mat.dist, ntrial = ntrial,
#> family = family, method = method.integrate, kappa = kappa,
#> offset = offset, U1 = U1, rank = rank), vecpar = TRUE,
#> skip.hessian = TRUE, control = control)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> X1 0.71409 0.18226 NA NA
#> X2 0.79465 0.17990 NA NA
#> logsigma2 0.16537 0.41144 NA NA
#> logphi -1.48301 0.20974 NA NA
#> logprec 0.61702 0.11935 NA NA
#>
#> -2 log L: 1973.482
X.pred <- X.all[(n+1):(n+n.pred),]
coords.pred <- coords.all[(n+1):(n+n.pred),]
pred.sglmm(result.nb,data=X.pred,coords=coords.pred)
#> [1] 1.0796784 2.5511352 6.2148177 6.4954555 1.6994834 1.1658579
#> [7] 3.4941327 2.7992210 1.9913884 13.2394374 11.4312459 9.5249584
#> [13] 5.6843326 3.2474477 1.3924241 11.1788155 6.3858692 15.6510471
#> [19] 30.5209369 32.5874261 33.8961389 2.3121162 1.0648111 1.4463879
#> [25] 1.6027887 15.2776189 1.6313418 4.2806683 7.8745649 0.4147299
#> [31] 3.3787403 26.5047551 2.2723055 3.3217126 2.9620488 3.7096393
#> [37] 2.1696730 1.0133735 11.9323421 6.8693548 2.8557717 1.1899519
#> [43] 3.1238762 6.8792812 5.6720593 7.1592524 2.3571801 1.1341541
#> [49] 23.1115252 0.7456393 19.5059490 1.2288857 9.0283915 4.2649047
#> [55] 1.4327406 1.2139899 1.9578175 6.5373381 0.8973056 8.8158868
#> [61] 2.1944846 4.5888146 2.2190761 0.7039620 2.2549576 1.4654028
#> [67] 3.0896944 0.8932928 3.5688911 0.9714466 3.9925341 3.0443346
#> [73] 12.2336647 4.8261169 0.9185217 4.7886871 0.4901198 3.0649941
#> [79] 14.9865994 2.4464208 9.9471592 1.5416763 1.2350628 1.6442479
#> [85] 0.5339071 2.8120669 9.0368075 1.5598224 3.3950542 4.6703657
#> [91] 1.1676898 3.2353459 2.0608186 1.5531116 11.2147363 18.4602356
#> [97] 0.8459434 3.4027514 3.8329304 12.2404823
if(requireNamespace("ngspatial")&
requireNamespace("mgcv")){
n = 30
A = ngspatial::adjacency.matrix(n)
Q = diag(rowSums(A),n^2) - A
x = rep(0:(n - 1) / (n - 1), times = n)
y = rep(0:(n - 1) / (n - 1), each = n)
X = cbind(x, y) # Use the vertex locations as spatial covariates.
beta = c(1, 1) # These are the regression coefficients.
P.perp = diag(1,n^2) - X%*%solve(t(X)%*%X)%*%t(X)
eig = eigen(P.perp %*% A %*% P.perp)
eigenvalues = eig$values
q = 400
M = eig$vectors[,c(1:q)]
Q.s = t(M) %*% Q %*% M
tau = 6
Sigma = solve(tau*Q.s)
set.seed(1)
delta.s = mgcv::rmvn(1, rep(0,q), Sigma)
lambda = exp( X%*%beta + M%*%delta.s )
Z = c()
for(j in 1:n^2){Z[j] = rpois(1,lambda[j])}
Y = as.matrix(Z,ncol=1)
data = data.frame("Y"=Y,"X"=X)
colnames(data) = c("Y","X1","X2")
} else {
stop("Packages \"ngspatial\" and \"mgcv\" are needed to generate a simulated data set")
}
#> Loading required namespace: ngspatial
linmod <- glm(Y~-1+X1+X2,data=data,family="poisson") # Find starting values
linmod$coefficients
#> X1 X2
#> 1.044368 1.041890
starting <- c(linmod$coefficients,"logtau"=log(1/var(linmod$residuals)) )
result.pois.disc <- fsglmm.discrete(Y~-1+X1+X2, inits = starting, data=data,family="poisson",ntrial=1, method.optim="BFGS", method.integrate="NR", rank=50, A=A)
result.pois.disc$summary
#> Maximum likelihood estimation
#>
#> Call:
#> bbmle::mle2(minuslogl = nlikSGLMM.discrete, start = inits, method = method.optim,
#> data = list(Y = Y, X = X, family = family, method = method.integrate,
#> ntrial, offset = offset, M = M, MQM = MQM, rank = rank),
#> vecpar = TRUE, skip.hessian = FALSE, control = list(maxit = 1000))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> X1 0.991924 0.050628 NA NA
#> X2 1.031817 0.050312 NA NA
#> logtau 1.724322 0.346589 NA NA
#>
#> -2 log L: 3480.376
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.