[adegenet-commits] r73 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 20 14:39:14 CET 2008
Author: jombart
Date: 2008-03-20 14:39:14 +0100 (Thu, 20 Mar 2008)
New Revision: 73
Modified:
pkg/R/chooseCN.R
pkg/R/spca.R
pkg/TODO
pkg/man/chooseCN.Rd
pkg/man/spca.Rd
pkg/man/spcaIllus.Rd
Log:
New option in chooseCN: inverse distances.
Modified: pkg/R/chooseCN.R
===================================================================
--- pkg/R/chooseCN.R 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/R/chooseCN.R 2008-03-20 13:39:14 UTC (rev 73)
@@ -1,8 +1,8 @@
#####################
-# Function .chooseCN
+# Function chooseCN
#####################
-chooseCN <- function(xy,ask=TRUE, type=1, result.type="nb", d1=NULL, d2=NULL, k=NULL,
- plot.nb=TRUE, edit.nb=FALSE){
+chooseCN <- function(xy,ask=TRUE, type=NULL, result.type="nb", d1=NULL, d2=NULL, k=NULL,
+ a=NULL, dmin=NULL, plot.nb=TRUE, edit.nb=FALSE){
if(is.data.frame(xy)) xy <- as.matrix(xy)
if(ncol(xy) != 2) stop("xy does not have two columns.")
@@ -20,12 +20,20 @@
x <- xy[,1]
y <- xy[,2]
temp <- table(x,y)
- if(any(temp>1)){
+ if(any(temp>1) & (!is.null(type) && type!=7)){ # coords need not be unique if type==7 (inverse distances)
xy <- jitter(xy)
warning("Random noise was added to xy as duplicated coordinates existed.")
}
- ##
+
+ ## handle type argument
+ if(!is.null(type)){
+ type <- as.integer(type)
+ if(type < 1 |type > 7) stop("type must be between 1 and 7")
+ ask <- FALSE
+ }
+ if(is.null(type) & !ask) { type <- 1 }
+
### begin large while ###
chooseAgain <- TRUE
while(chooseAgain){
@@ -45,17 +53,18 @@
cat("\t Minimum spanning tree (type 4)\n")
cat("\t Neighbourhood by distance (type 5)\n")
cat("\t K nearest neighbours (type 6)\n")
+ cat("\t Inverse distances (type 7)\n")
cat("Answer: ")
type <- as.integer(readLines(n = 1))
- temp <- type < 1 |type > 6
+ temp <- type < 1 |type > 7
if(temp) cat("\nWrong answer\n")
} # end while
}
##
## graph types
- # type 1: Delaunay
+ ## type 1: Delaunay
if(type==1){
if(!require(tripack, quiet=TRUE)) stop("tripack library is required.")
cn <- tri2nb(xy)
@@ -67,20 +76,20 @@
cn <- graph2nb(cn, sym=TRUE)
}
- # type 3: Relative neighbours
+ ## type 3: Relative neighbours
if(type==3){
cn <- relativeneigh(xy)
cn <- graph2nb(cn, sym=TRUE)
}
- # type 4: Minimum spanning tree
+ ## type 4: Minimum spanning tree
if(type==4){
if(!require(ade4, quiet=TRUE)) stop("ade4 library is required.")
cn <- mstree(dist(xy))
cn <- neig2nb(cn)
}
- # type 5: Neighbourhood by distance
+ ## type 5: Neighbourhood by distance
if(type==5){
if(is.null(d1) |is.null(d2)){
tempmat <- as.matrix(dist(xy))
@@ -101,7 +110,7 @@
cn <- dnearneigh(x=xy, d1=d1, d2=d2)
}
- # type 6: K nearests
+ ## type 6: K nearests
if(type==6){
if(is.null(k)) {
cat("\n Enter the number of neighbours: ")
@@ -110,6 +119,29 @@
cn <- knearneigh(x=xy, k=k)
cn <- knn2nb(cn, sym=TRUE)
}
+
+ ## type 7: inverse distances
+ if(type==7){
+ if(is.null(a)) {
+ cat("\n Enter the exponent: ")
+ a <- as.numeric(readLines(n = 1))
+ }
+ cn <- as.matrix(dist(xy))
+ if(is.null(dmin)) {
+ cat("\n Enter the minimum distance \n(range = 0 -", max(cn),"): ")
+ dmin <- as.numeric(readLines(n = 1))
+ }
+ thres <- mean(cn)/1e8
+ if(dmin > thres) dmin <- thres
+ cn[cn < dmin] <- dmin
+ cn <- 1/(cn^a)
+ diag(cn) <- 0
+ cn <- prop.table(cn,1)
+ plot.nb <- FALSE
+ edit.nb <- FALSE
+ result.type <- "listw"
+ } # end type 7
+
## end graph types
if(ask & plot.nb) {
@@ -122,14 +154,21 @@
plot(cn,xy)
chooseAgain <- FALSE
}
- else {chooseAgain <- FALSE}
+ else {chooseAgain <- FALSE}
}
### end large while
if(edit.nb) {cn <- edit(cn,xy)}
- if(result.type == "listw") {cn <- nb2listw(cn, style="W", zero.policy=TRUE)}
+ if(result.type == "listw") {
+ if(type!=7) {
+ cn <- nb2listw(cn, style="W", zero.policy=TRUE)
+ } else {
+ cn <- mat2listw(cn)
+ cn$style <- "W"
+ }
+ }
res <- cn
Modified: pkg/R/spca.R
===================================================================
--- pkg/R/spca.R 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/R/spca.R 2008-03-20 13:39:14 UTC (rev 73)
@@ -16,8 +16,9 @@
################
# Function spca
################
-spca <- function(obj, xy=NULL, cn=NULL, scale=FALSE, scannf=TRUE, nfposi=1, nfnega=1, type=1,
- ask=TRUE, plot.nb=TRUE, edit.nb=FALSE ,truenames=TRUE, d1=NULL, d2=NULL, k=NULL) {
+spca <- function(obj, xy=NULL, cn=NULL, scale=FALSE, scannf=TRUE, nfposi=1, nfnega=1, type=NULL,
+ ask=TRUE, plot.nb=TRUE, edit.nb=FALSE ,truenames=TRUE, d1=NULL, d2=NULL, k=NULL,
+ a=NULL, dmin=NULL) {
if(!inherits(obj,c("genind","genpop"))) stop("obj must be a genind or genpop object.")
invisible(validObject(obj))
@@ -34,14 +35,15 @@
# connection network
if(is.null(cn)) {
if(is.null(xy)) stop("'xy' and 'cn' are both missing")
- resCN <- chooseCN(xy=xy, ask=ask, type=type, plot.nb=plot.nb, edit.nb=edit.nb, result.type="nb", d1=d1, d2=d2, k=k)
+ resCN <- chooseCN(xy=xy, ask=ask, type=type, plot.nb=plot.nb, edit.nb=edit.nb,
+ result.type="listw", d1=d1, d2=d2, k=k, a=a, dmin=dmin)
} else {
- if(!inherits(cn,"nb")) stop("cn is not of class 'nb' (package spdep)")
- resCN <- cn
+ if(inherits(cn,"nb") & !inherits(cn,"listw")) { cn <- nb2listw(cn, style="W", zero.policy=TRUE) }
+ if(!inherits(cn,"listw")) stop("cn does not have a recognized class ('nb' or 'listw', package spdep)")
+ resCN <- cn
}
- cn.nb <- resCN
- xy <- resCN at xy
- cn.lw <- nb2listw(cn.nb, style="W", zero.policy=TRUE)
+
+ xy <- attr(resCN,"xy")
# prepare data
f1 <- function(vec){
@@ -63,7 +65,7 @@
# perform analyses
pcaX <- dudi.pca(X, center=TRUE, scale=scale, scannf=FALSE)
- spcaX <- multispati(dudi=pcaX, listw=cn.lw, scannf=scannf, nfposi=nfposi, nfnega=nfnega)
+ spcaX <- multispati(dudi=pcaX, listw=resCN, scannf=scannf, nfposi=nfposi, nfnega=nfnega)
nfposi <- spcaX$nfposi
nfnega <- spcaX$nfnega
@@ -72,7 +74,7 @@
rownames(spcaX$xy) <- rownames(spcaX$li)
colnames(spcaX$xy) <- c("x","y")
- spcaX$cn <- cn.nb
+ spcaX$lw <- resCN
spcaX$call <- appel
@@ -137,7 +139,7 @@
print(sumry)
cat("\n$xy: matrix of spatial coordinates")
- cat("\n$cn: connection network (class nb)")
+ cat("\n$lw: a list of spatial weights (class 'listw')")
cat("\n\nother elements: ")
if (length(names(x)) > 10)
@@ -204,12 +206,12 @@
dudi <- dudi.pca(X, center=TRUE, scale=FALSE, scannf=FALSE, nf=nfposi+nfnega)
## end of pca
- listw <- nb2listw(object$cn, style="W", zero.policy=TRUE)
+ lw <- object$lw
# I0, Imin, Imax
n <- nrow(X)
I0 <- -1/(n-1)
- L <- listw2mat(listw)
+ L <- listw2mat(lw)
eigL <- eigen(0.5*(L+t(L)))$values
Imin <- min(eigL)
Imax <- max(eigL)
@@ -230,7 +232,7 @@
eig <- dudi$eig[1:nf]
cum <- cumsum(dudi$eig)[1:nf]
ratio <- cum/sum(dudi$eig)
- w <- apply(dudi$l1,2,lag.listw,x=listw)
+ w <- apply(dudi$l1,2,lag.listw,x=lw)
moran <- apply(w*as.matrix(dudi$l1)*dudi$lw,2,sum)
res <- data.frame(var=eig,cum=cum,ratio=ratio, moran=moran)
row.names(res) <- paste("Axis",1:nf)
@@ -248,7 +250,7 @@
nfposimax <- sum(eig > 0)
nfnegamax <- sum(eig < 0)
- ms <- multispati(dudi=dudi, listw=listw, scannf=FALSE,
+ ms <- multispati(dudi=dudi, listw=lw, scannf=FALSE,
nfposi=nfposimax, nfnega=nfnegamax)
ndim <- dudi$rank
@@ -290,7 +292,14 @@
z <- x$ls[,axis]
nfposi <- x$nfposi
nfnega <- x$nfnega
- neig <- nb2neig(x$cn)
+ ## handle neig parameter - hide cn if nore than 100 links
+ nLinks <- sum(card(x$lw$neighbours))
+ if(nLinks < 100) {
+ neig <- nb2neig(x$lw$neighbours)
+ } else {
+ neig <- NULL
+ }
+
sub <- paste("Score",axis)
csub <- 2
@@ -305,7 +314,7 @@
box()
# 3
- if(n<30) {neig <- nb2neig(x$cn)} else {neig <- NULL}
+ if(n<30) {neig <- nb2neig(x$lw$neighbours)} else {neig <- NULL}
s.value(xy,z, include.ori=FALSE, addaxes=FALSE, clegend=0, csize=.6,
neig=neig, sub=sub, csub=csub, possub="bottomleft")
Modified: pkg/TODO
===================================================================
--- pkg/TODO 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/TODO 2008-03-20 13:39:14 UTC (rev 73)
@@ -36,7 +36,7 @@
=====================
* Implement "sep" argument in df2genind
* Implement a method to merge different markers for the same individuals
-*
+* Implement spatial weights derived from inverse spatial distances in chooseCN -- done (TJ)
# TESTING:
==========
Modified: pkg/man/chooseCN.Rd
===================================================================
--- pkg/man/chooseCN.Rd 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/man/chooseCN.Rd 2008-03-20 13:39:14 UTC (rev 73)
@@ -9,36 +9,48 @@
connection network either with classe \code{nb} or \code{listw}.
}
\usage{
-chooseCN(xy, ask = TRUE, type = 1, result.type = "nb", d1 = NULL,
- d2 = NULL, k = NULL, plot.nb = TRUE, edit.nb = FALSE)
+chooseCN(xy, ask = TRUE, type = NULL, result.type = "nb", d1 = NULL,
+ d2 = NULL, k = NULL, a=NULL, dmin=NULL, plot.nb = TRUE, edit.nb = FALSE)
}
\arguments{
\item{xy}{an matrix or data.frame with two columns for x and y coordinates.}
\item{ask}{a logical stating whether graph should be chosen
- interactively (TRUE,default) or not (FALSE).}
- \item{type}{an integer giving the type of graph (see details). Used if
- \code{ask=FALSE}}
+ interactively (TRUE,default) or not (FALSE). Set to FALSE if \code{type} is
+ provided.}
+ \item{type}{an integer giving the type of graph (see details).}
\item{result.type}{a character giving the class of the returned
object. Either "nb" (default) or "listw", both from \code{spdep}
- package.}
+ package. See details.}
\item{d1}{the minimum distance between any two neighbours. Used if
\code{type=5.}}
\item{d2}{the maximum distance between any two neighbours. Used if
\code{type=5}.}
\item{k}{the number of neighbours per point. Used if
\code{type=6}.}
+ \item{a}{the exponent of the inverse distance matrix. Used if
+ \code{type=7}.}
+ \item{dmin}{the minimum distance between any two distinct points. Used
+ to avoid infinite spatial proximities (defined as the inversed
+ spatial distances). Used if \code{type=7}.}
\item{plot.nb}{a logical stating whether the resulting graph should be
plotted (TRUE, default) or not (FALSE).}
\item{edit.nb}{a logical stating whether the resulting graph should be
edited manually for corrections (TRUE) or not (FALSE, default).}
}
-\details{There are 6 kinds of graphs proposed: \cr
+\details{There are 7 kinds of graphs proposed: \cr
Delaunay triangulation (type 1)\cr
Gabriel graph (type 2)\cr
Relative neighbours (type 3)\cr
Minimum spanning tree (type 4)\cr
Neighbourhood by distance (type 5)\cr
K nearests neighbours (type 6)\cr
+ Inverse distances (type 7)\cr
+
+ The last option (type=7) is not a true neighbouring graph: all sites are
+ neighbours, but the spatial weights are directly proportional to the
+ inversed spatial distances.\cr
+ Also not that in this case, the output of the function is always a
+ \code{listw} object, even if \code{nb} was requested.
}
\value{Returns a connection network having the class \code{nb} or
\code{listw}. The xy coordinates are passed as attribute to the
Modified: pkg/man/spca.Rd
===================================================================
--- pkg/man/spca.Rd 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/man/spca.Rd 2008-03-20 13:39:14 UTC (rev 73)
@@ -26,8 +26,9 @@
autocorrelation\cr
}
\usage{
-spca(obj, xy=NULL, cn=NULL, scale=FALSE, scannf=TRUE, nfposi=1, nfnega=1, type=1, ask=TRUE,
-plot.nb=TRUE, edit.nb=FALSE ,truenames=TRUE, d1=NULL, d2=NULL, k=NULL)
+spca(obj, xy=NULL, cn=NULL, scale=FALSE, scannf=TRUE, nfposi=1, nfnega=1, type=NULL, ask=TRUE,
+plot.nb=TRUE, edit.nb=FALSE ,truenames=TRUE, d1=NULL, d2=NULL, k=NULL,
+ a=NULL, dmin=NULL)
\method{print}{spca}(x, \dots)
@@ -53,7 +54,7 @@
\item{nfnega}{an integer giving the number of negative eigenvalues retained
('local structures').}
\item{type}{an integer giving the type of graph (see details in
- \code{chooseCN} help page). Used if \code{ask=FALSE}}
+ \code{chooseCN} help page). If provided, \code{ask} is set to FALSE.}
\item{ask}{a logical stating whether graph should be chosen
interactively (TRUE,default) or not (FALSE).}
\item{plot.nb}{a logical stating whether the resulting graph should be
@@ -68,6 +69,11 @@
\code{type=5}.}
\item{k}{the number of neighbours per point. Used if
\code{type=6}.}
+ \item{a}{the exponent of the inverse distance matrix. Used if
+ \code{type=7}.}
+ \item{dmin}{the minimum distance between any two distinct points. Used
+ to avoid infinite spatial proximities (defined as the inversed
+ spatial distances). Used if \code{type=7}.}
\item{x}{a \code{spca} object.}
\item{object}{a \code{spca} object.}
\item{printres}{a logical stating whether results should be printed on
@@ -101,7 +107,7 @@
sPCA axes.}
\item{call}{the matched call.}
\item{xy}{a matrix of spatial coordinates.}
- \item{cn}{a connection network of class \code{nb}.}
+ \item{lw}{a list of spatial weights of class \code{listw}.}
Other functions have different outputs:\cr
- \code{summary.spca} returns a list with 3 components: \code{Istat}
@@ -118,10 +124,6 @@
Revealing cryptic spatial patterns in genetic variability by a new
multivariate method. Submitted to \emph{Heredity}.
-Smouse, P. E. and Peakall, R. (1999) Spatial autocorrelation analysis of
-individual multiallele and multilocus genetic structure. \emph{Heredity},
-\bold{82}, 561--573.
-
Wartenberg, D. E. (1985) Multivariate spatial correlation: a method for
exploratory geographical analysis. \emph{Geographical Analysis},
\bold{17}, 263--283.
Modified: pkg/man/spcaIllus.Rd
===================================================================
--- pkg/man/spcaIllus.Rd 2008-03-18 17:43:45 UTC (rev 72)
+++ pkg/man/spcaIllus.Rd 2008-03-20 13:39:14 UTC (rev 73)
@@ -67,7 +67,7 @@
# an auxiliary function for graphics
plotaux <- function(x,analysis,axis=1,lab=NULL,...){
neig <- NULL
-if(inherits(analysis,"spca")) neig <- nb2neig(analysis$cn)
+if(inherits(analysis,"spca")) neig <- nb2neig(analysis$lw$neighbours)
xrange <- range(x$other$xy[,1])
xlim <- xrange + c(-diff(xrange)*.1 , diff(xrange)*.45)
yrange <- range(x$other$xy[,2])
More information about the adegenet-commits
mailing list