小编Pas*_*cal的帖子

我正在绘制ROC并测量部分AUC作为生态位模型质量的度量.当我在R工作时,我正在使用ROCR和pROC包.我会选择一个使用,但就目前而言,我只是想看看它们是如何表现的,以及是否能更好地满足我的需求.

令我困惑的一件事是,在绘制ROC时,轴如下:

ROCR

x axis: 'true positive rate' 0 -> 1
y axis: 'false positive rate', 0 -> 1

PROC

x axis: 'sensitivity' 0 -> 1
y axis: 'specificity' 1 -> 0. 

但如果我使用两种方法绘制ROC,它们看起来完全相同.所以我只想确认一下:

true positive rate = sensitivity

false positive rate = 1 - specificity.

这是一个可重复的例子:

obs<-rep(0:1, each=50)
pred<-c(runif(50,min=0,max=0.8),runif(50,min=0.3,max=0.6))
plot(roc(obs,pred))

ROCRpred<-prediction(pred,obs)
plot(performance(ROCRpred,'tpr','fpr'))

我试图在R中绘制一个表,列名相对于表的角度.我想添加行来分隔这些列名称,与文本的角度相同.但是,看起来text()函数中指定的角度与绘图的纵横比无关,而我在segments()函数中使用的角度取决于绘图的纵横比.

这是我的意思的一个例子:

nRows <- 5
nColumns <- 3
theta <- 30

rowLabels <- paste('row', 1:5, sep='')
colLabels <- paste('col', 1:3, sep='')

plot.new()
par(mar=c(1,8,5,1), xpd=NA)
plot.window(xlim = c(0, nColumns), ylim = c(0, nRows), asp = 1)
text(labels = rowLabels, x=0, y=seq(from=0.5, to=nRows, by=1), pos=2)
text(labels = colLabels, x = seq(from = 0.4, to = nColumns, by = 1), y = nRows + 0.1, pos = 4, srt = theta, cex = 1.1)
segments(x0 = c(0:nColumns), x1 = …

我正在R中运行一些地理处理任务,我正在尝试创建一些用于剪切环境信息栅格的多边形.我正在缓冲一些复杂的多边形,这留下了我想摆脱的小子几何.在ArcGIS中,我认为这将涉及将我的多边形从多部分转换为单部分(或沿着这些部分的东西)然后溶解,但我不知道如何在R中执行此操作.

这是一个说明问题的示例:

require(maptools)
require(rgeos)

data(wrld_simpl)
wrld_simpl[which(wrld_simpl@data$NAME=='Greece'),]->greece
proj4string(greece)<-CRS('+proj=lonlat +datum=WGS84')
gBuffer(greece,width=0.5)->buf
plot(buf)

我真正想要的是多边形的外边界,里面没有别的东西.有任何想法吗？

我想计算从一组点到这些点的质心的多元距离.马哈拉诺比斯距离似乎适合这一点.但是,我收到一个错误(见下文).

任何人都可以告诉我为什么我会收到此错误,如果有办法解决它？

如果下载坐标数据和关联的环境数据,则可以运行以下代码.

require(maptools)
occ <- readShapeSpatial('occurrences.shp')
load('envDat.Rdata')

#standardize the data to scale the variables
dat <- as.matrix(scale(dat))
centroid <- dat[1547,]  #let's assume this is the centroid in this case

#Calculate multivariate distance from all points to centroid
mahalanobis(dat,center=centroid,cov=cov(dat))

Error in solve.default(cov, ...) : 
  system is computationally singular: reciprocal condition number = 9.50116e-19

给定一组xy坐标,我如何选择n个点,使得这些n个点彼此距离最远？

对于大数据集可能不会做得太好的低效方法如下(确定最远的1000个中的20个点):

xy <- cbind(rnorm(1000),rnorm(1000))

n <- 20
bestavg <- 0
bestSet <- NA
for (i in 1:1000){
    subset <- xy[sample(1:nrow(xy),n),]
    avg <- mean(dist(subset))
    if (avg > bestavg) {
        bestavg <- avg
        bestSet <- subset
    }
}

给定一组点，我尝试选择在这组点上分布最均匀的 n 个点的子集。换句话说，我试图精简数据集，同时仍然在空间中均匀采样。

到目前为止，我有以下内容，但这种方法可能不适用于较大的数据集。也许有一种更智能的方法来首先选择点的子集...以下代码随机选择点的子集，并寻求最小化该子集内的点与该子集外的点之间的距离。

建议表示赞赏！

evenSubset <- function(xy, n) {

    bestdist <- NA
    bestSet <- NA
    alldist <- as.matrix(dist(xy))
    diag(alldist) <- NA
    alldist[upper.tri(alldist)] <- NA
    for (i in 1:1000){
        subset <- sample(1:nrow(xy),n)
        subdists <- alldist[subset,-subset]
        distsum <- sum(subdists,na.rm=T)
        if (distsum < bestdist | is.na(bestdist)) {
            bestdist <- distsum
            bestSet <- subset
        }
    }
    return(xy[bestSet,])
}

xy2 <- evenSubset(xy=cbind(rnorm(1000),rnorm(1000)), n=20)
plot(xy)
points(xy2,col='blue',cex=1.5,pch=20)

我一直在查看教程和文档,但还没有弄清楚如何为所有列中的一个现有行分配一个值向量data.table.

我从一个data.table已经有正确数量的列和行的空开始:

dt <- data.table(matrix(nrow=10, ncol=5))

现在我计算一行之外的一些值data.table并将它们放在一个向量中vec,例如:

vec <- rnorm(5)

我怎样才能将vec的值分配给例如第一行data.table同时实现良好的性能(因为我还想逐步填充其他行)？

我正在为一组物种生成生态位模型,我想使用AUC作为生态位生态质量的指标.开发Maxent的Steven Phillips在他的Maxent手册中提供了用于计算R中AUC的代码.但是,我正在阅读报告部分AUC比率作为更健壮和概念上合理的度量的论文.我想我理解如何使用ROCR R软件包计算部分AUC,但是如何计算AUC比率？

这是Phillips的教程脚本:

presence<-read.csv("bradypus_variegatus_samplePredictions.csv")
background<-read.csv("bradypus_variegatus_backgroundPredictions.csv")
pp<-presence$Logistic.prediction
testpp<-pp[presence$Test.or.train=="test"]
trainpp<-pp[presence$Test.or.train=="train"]
bb<-background$logistic

combined<-c(testpp,bb)
label<-c(rep(1,length(testpp)),rep(0,length(bb)))
pred<-prediction(combined,label)
perf<-performance(pred,"tpr","fpr")
plot(perf,colorize=TRUE)
performance(pred,"auc")@y.values[[1]] #RETURNS AUC

AUC<-function(p,ind){
    pres<-p[ind]
    combined<-c(pres,bb)
    label<-c(rep(1,length(pres)),rep(0,length(bb)))
    predic<-prediction(combined,label)
    return(performance(predic,'auc')@y.values[[1]])
}

b1<-boot(testpp,AUC,100) #RETURNS AUC WITH STANDARD ERROR
b1

任何建议或建议将不胜感激!谢谢.

我有一张平均值和四分位数范围的表格.我想创建一个点图,其中点将显示此平均值,并且条形图将延伸通过点,以显示四分位数范围.换句话说,点将位于条的中点,其长度将等于我的四分位数范围数据.我在R工作

例如,

labels<-c('a','b','c','d')
averages<-c(10,40,20,30)
ranges<-c(5,8,4,10)
dotchart(averages,labels=labels)

然后将范围作为条形添加到该图中.

有任何想法吗？

谢谢!

我广泛使用R包rgdal和rgeos(目前使用R v3.2.2).最近,在我的ubuntu机器(ubuntu v15.10)上,当我加载rgdal包时,我看到以下内容:

> library(rgdal)
Loading required package: sp
rgdal: version: 1.0-7, (SVN revision 559)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.1, released 2014/09/24
 Path to GDAL shared files: /usr/local/share/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
 Path to PROJ.4 shared files: (autodetected)
WARNING: no proj_defs.dat in PROJ.4 shared files
 Linking to sp version: 1.2-1 

我在最近发布的r-sig-geo帖子中看到这是一个已知问题,并且已发布rgdal软件包的更新以解决此问题.但是,尽管重新安装了rgdal软件包,还是使用以下命令重新安装gdal和proj4:

sudo apt-get update && sudo apt-get install libgdal-dev libproj-dev

我没看见有分别.当我使用这些包的功能时,此问题会生成数百条警告消息.

相比之下,在我的Mac上,如果我加载rgdal包,我明白了

> …

我在相同区域Behrmann投影中有一个光栅,我想把它投射到Mollweide投影和情节.

然而,当我使用以下代码执行此操作时,绘图似乎不正确,因为地图延伸到两侧,并且有各种陆地的轮廓,我不期望它们.此外,地图延伸到绘图之外窗口.

任何人都可以帮助我得到这个很好的情节？

谢谢!

使用的数据文件可以从此链接下载.

这是我到目前为止的代码:

require(rgdal)
require(maptools)
require(raster)

data(wrld_simpl)

mollCRS <- CRS('+proj=moll')
behrmannCRS <- CRS('+proj=cea +lat_ts=30')

sst <- raster("~/Dropbox/Public/sst.tif", crs=behrmannCRS)

sst_moll <- projectRaster(sst, crs=mollCRS)
wrld <- spTransform(wrld_simpl, mollCRS)

plot(sst_moll)
plot(wrld, add=TRUE)

我编写了一个函数来运行系统发育的广义最小二乘法,一切看起来应该可以正常工作,但由于某种原因,在脚本(W)中定义的特定变量不断变为未定义.我已经盯着这段代码好几个小时了,无法弄清楚问题出在哪里.

有任何想法吗？

myou <- function(alpha, datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    W<-diag(vcv.phylo(tree)) # Weights
    fm <- gls(Trait1 ~ Trait2, data=dat, correlation = corMartins(alpha, tree, fixed = TRUE),weights = ~ W,method = "REML")
    return(as.numeric(fm$logLik))
}

corMartins2<-function(datax, datay, tree){
    data.frame(datax[tree$tip.label,],datay[tree$tip.label,],row.names=tree$tip.label)->dat
    colnames(dat)<-c("Trait1","Trait2")
    result <- optimize(f = myou, interval = c(0, 4), datax=datax,datay=datay, tree = tree, maximum = TRUE)
    W<-diag(vcv.phylo(tree)) # Weights
    fm <- gls(Trait1 ~ Trait2, data = dat, correlation = corMartins(result$maximum, tree, fixed =T),weights = ~ W,method = "REML")
    list(fm, result$maximum)}




#test


require(nlme) …

我正在尝试在 Rcpp 中创建一个函数，它将接受成对数字矩阵以及向量列表作为输入，每个元素都是行/列名称的子集。我希望这个函数识别与这些名称匹配的矩阵子集，并返回值的平均值。

下面我生成了一些类似于我拥有的数据类型的虚拟数据，然后尝试使用 Rcpp 函数。

library(Rcpp)

dat <- c(spA = 4, spB = 10, spC = 8, spD = 1, spE = 5, spF = 9)
pdist <- as.matrix(dist(dat))
pdist[upper.tri(pdist, diag = TRUE)] <- NA

这里我有一个由 pdist 中行/列名称的各个子集的字符向量组成的列表

subsetList <- replicate(10, sample(names(dat), 4), simplify=FALSE)

对于每一组名称，我想确定成对矩阵的子集并取值的平均值

这是我到目前为止所拥有的，它不起作用，但我认为它说明了我想要达到的目标。

cppFunction('
    List meanDistByCell(List input, NumericMatrix pairmat) {
    int n = input.size();
    List out(n);

    List dimnames = pairmat.attr( "dimnames" );
    CharacterVector colnames = dimnames[1];

    for (int i = 0; i < n; i++) { …