# # This file is meant to be used in R http://r-project.org # Copyright (C) 2008-2010 William Briggs # matt@wmbriggs.com # Comments and corrections more than welcome # # This is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA # newdbinom = function(x, n.new, k.old, n.old){ # New observables for binomial distribution given old data and guess of how many new data points there will be n.new # beta-binomial a = k.old+1 b = n.old-k.old+1 (ans = exp(lchoose(n.new,x)+lbeta(x+a,n.new-x+b)-lbeta(a,b))) } newpbinom = function(x, n.new, k.old, n.old){ # New observables for binomial distribution given old data and guess of how many new data points there will be n.new # beta-binomial (ans = cumsum(newdbinom(0:x,n.new, k.old, n.old))) } newpnorm = function(x, x.old){ # New observables for normal distribution given old data x.old. # t fit = glm(x.old~1) newdata= data.frame(A=1) s1 = obs.glm(fit,newdata) (ans = pt((x-s1$central)/s1$scale,s1$df)) } glm.posterior = function(fit){ # approximate posterior probabilities for parameters in glm models # how to run # glm.posterior(fit) par(ask=TRUE) a = summary(fit) d = dim(a$coefficients) for (i in 1:d[1]){ b = a$coefficients[i,1:2] s = seq(b[1]-3*b[2],b[1]+3*b[2],length.out=1000) plot(s,dnorm(s,b[1],b[2]),xlab=row.names(a$coefficients)[i],ylab="Density",type='l') abline(v = 0, lty = 3) cat("Probability parameter ", row.names(a$coefficients)[i], "< 0 | x, E_N =",pnorm(0,b[1],b[2]),"\n") } par(ask=FALSE) } obs.glm = function(fit,newdata, nb = 10000){ x = model.matrix(fit) y = as.matrix(model.frame(fit)[,1]) n = length(y) tt <- terms(fit) Terms <- delete.response(tt) m <- model.frame(Terms, newdata, xlev = fit$xlevels) z <- model.matrix(Terms, m) if (fit$family$family == "gaussian"){ # flat prior constants b=0.0001;c=0.001;d=0.001; #theta1 = sum(x.old)/(b+length(x.old)) theta1<-solve(b+t(x) %*% x) %*% (t(x) %*% y) #fn1 = (b + length(x.new)+length(x.old)) fn1<-1-z %*% solve(t(z) %*% z+b+t(x) %*% x) %*% t(z) #beta1 = d+ beta1<-d+t(y-x %*% theta1) %*% y/2 #print( sqrt(as.numeric(fn1*(c+n/2)/as.numeric(beta1))) ) f <- list(central=as.numeric(z %*% theta1), scale=1/sqrt(as.numeric(fn1*(c+n/2)/as.numeric(beta1))),df=2*c+n,family=fit$family$family) #f <- list(central=as.numeric(z %*% theta1), scale=1/sqrt(as.numeric(fn1*(c+n/2)/beta1)),df=2*c+n,family=fit$family$family) } else { # most likely a binomial # crude approximation, waiting for a better man to compute a = summary(fit)$coefficients #s = matrix(0,nb,dim(a)[1]) #for (i in 1:dim(a)[1]){ # s[1:nb,i] = rnorm(nb,a[i,1],a[i,2]) #} v = vcov(fit) s = rmultnorm(nb,a[,1],v) s = s %*% t(z) s = 1/(1+exp(-s)) f <- list(s=s, newdata=newdata,family=fit$family$family) } } # how to run #newdata= data.frame(Campaign="A") #newdata= data.frame(Sex="M",n=1000) #s1 = obs.glm(fit,newdata) #s1 #newdata= data.frame(Campaign="B") #newdata= data.frame(Sex="F",n=1000) #s2 = obs.glm(fit,newdata) #s2 #obs.glm.prob(s1,s2) glm.pred = function(fit,p= c(0.1,.5,.9), nsd=3){ a = model.frame(fit) y = names(a[,1,drop=FALSE]) z = as.data.frame(a[,-1]) names(z) = names(a)[-1] # the most common values of 'x' (or z) stored in q q = data.frame(z[1,]) names(q) = names(a)[-1] for(i in 1:dim(z)[2]){ if(is.factor(z[,i])){ a=sort(table(z[,i]),decreasing=TRUE) q[1,i] = names(a)[1] } else { a=quantile(z[,i],probs=0.5,na.rm=TRUE) q[1,i] = a } } #print(y) #write.csv(q,row.names=FALSE) #cat('\n') cat('Variable',',','Level',',','Central ',y,',','Prob(',y,'>base level)','\n',sep='') for(i in 1:dim(z)[2]){ if(is.factor(z[,i])){ a=levels(z[,i]) m = length(a) newdata = q[rep(1,m),,drop=FALSE] newdata[,i] = a } else { m = length(p) a=as.numeric(quantile(z[,i],probs=p,na.rm=TRUE)) newdata = q[rep(1,m),,drop=FALSE] newdata[,i] = a } b=obs.glm(fit,newdata) central = b$central scale = diag(matrix(b$scale,m,m))^2 #df = a$df d=pnorm(0,central[1]-central,sqrt(scale[1] + c(0,scale[-1]))) w = data.frame(variable = names(z)[i], level = a, central = signif(central,nsd), prob.gt.bl = signif(d,nsd)) write.table(w,row.names=FALSE,col.names=FALSE,sep=',') } } skill = function(fit1,fit2,plot=FALSE){ # fit2 is assumed to be the worse, or ``null", model # these might be different because of different numbers of missing obs in the two fits n1 = length(fit1$y) n2 = length(fit2$y) ny = max(n1,n2) # find the missing values! # assumes no missing values for the y variable at this time b = as.data.frame(eval(attr(fit$terms,"variables"))) a = na.omit(b) b = as.integer(attr(a,"na.action")) crps = 0 score1 = rep(NA,ny) score2 = rep(NA,ny) k = 0 if (fit1$family$family == "gaussian"){ x = seq(.9*min(fit1$y),1.1*max(fit1$y),length.out = ny*10) for (i in 1:ny) { newdata = data.frame(fit1$model[i,]) s1 = obs.glm(fit1,newdata) newdata = data.frame(fit2$model[i,]) s2 = obs.glm(fit2,newdata) sd1 = pt((x-s1$central)/s1$scale,s1$df) sd2 = pt((x-s2$central)/s2$scale,s2$df) # rough approximation if (!is.na(fit1$y[i])) score1[i] = sum((sd1-(x>fit1$y[i]))^2)/length(x) if (!is.na(fit2$y[i])) score2[i] = sum((sd2-(x>fit2$y[i]))^2)/length(x) #if (!is.na(fit1$y[i])) score1 = score1 + sum((sd1-(x>fit1$y[i]))^2)/length(x) #if (!is.na(fit2$y[i])) score2 = score2 + sum((sd2-(x>fit2$y[i]))^2)/length(x) # crps #if (!is.na(fit1$y[i])) score1 = score1 + crps(fit1$y[i], c(s1$central,s1$scale)) #if (!is.na(fit2$y[i])) score2 = score2 + crps(fit2$y[i], c(s2$central,s2$scale)) } } else { # this takes the raw prob of obs and computes, lord help us, # the brier score; still waiting for better man for (i in 1:ny) { newdata = data.frame(fit1$model[i,], n=1) s1 = obs.glm(fit1,newdata) newdata = data.frame(fit2$model[i,], n=1) s2 = obs.glm(fit2,newdata) # rough approximation if (!(i %in% b)) { k = k + 1 score1[i] = sum((s1$s-fit1$y[i])^2)/length(s1$s) } score2[i] = sum((s2$s-fit2$y[i])^2)/length(s2$s) #if (!is.na(fit1$y[i])) score1[i] = sum((s1$s-fit1$y[i])^2)/length(s1$s) #if (!is.na(fit2$y[i])) score2[i] = sum((s2$s-fit2$y[i])^2)/length(s2$s) #if (!is.na(fit1$y[i])) score1 = score1 + sum((s1$s-fit1$y[i])^2)/length(s1$s) #if (!is.na(fit2$y[i])) score2 = score2 + sum((s2$s-fit2$y[i])^2)/length(s2$s) } } if(plot){ par(mfrow=c(2,2)) if (fit1$family$family == "gaussian"){ scatter.smooth(score2,score1,ylab="Full model",xlab="Null model") } else { plot(score2,score1,ylab="Full model",xlab="Null model") } abline(0,1,lty=2) plot(density(score1,na.rm=T),xlab="CRPS scores",ylab="",main="") lines(density(score2,na.rm=T),lty=2,col=2) plot(density((score2-score1)/score2,to=1,na.rm=T),xlab="Skill",ylab="",main="") abline(v=0,lty=2) y = fit1$y; if(n2>n1) y = fit2$y if (fit1$family$family == "gaussian"){ scatter.smooth(y,(score2-score1)/score2,ylab="Skill", xlab=names(fit1$model)[1]) } else { plot(y,(score2-score1)/score2,ylab="Skill", xlab=names(fit1$model)[1]) } abline(h=0,lty=2) if (existsFunction("X11")) X11() if (existsFunction("windows")) windows() mfrow <- set.mfrow(0, dim(fit$model)[2], 1) par(mfrow = mfrow) for (j in 1:dim(fit1$model)[2]){ y = (score2-score1)/score2 if(length(b)>1) y = (score2[-b]-score1[-b])/score2[-b] scatter.smooth(fit1$model[, j],y,axes=F,ylab="Skill",xlab= names(fit1$model)[j]) abline(h=0,lty=2) axis(1);axis(2) } } score1 = mean(score1,na.rm=T) score2 = mean(score2,na.rm=T) #if (n1>0) score1 = score1/n1 #if (n2>0) score2 = score2/n2 if (score2>0) crps = (score2 - score1)/score2 a1 = summary(fit1) a2 = summary(fit2) aic.skill = (a2$aic-a1$aic)/a2$aic r2.1 = (a1$null.deviance-a1$deviance)/a1$null.deviance r2.2 = (a2$null.deviance-a2$deviance)/a2$null.deviance r2.skill = r2.1 cat("Score 1 = ",score1, ", Score 2 = ", score2, ", Skill =", crps, "\n", "AIC skill =", aic.skill, ", R^2 skill =", r2.skill,"\n") ans = list(score1 = score1, score2 = score2, skill = crps, aic.skill = aic.skill, r.squared = r2.skill) } # source('../Rcode.R') # skill(fit1,fit2) obs.glm.prob = function(s1,s2, pic = TRUE, labels = list(xlab="New Observables", leg1 = "S1" , leg2 = "S2", main="") ){ if (s1$family == "gaussian"){ # normal approx to difference in t m = s1$central-s2$central s = sqrt(s1$scale^2+s2$scale^2) cat("Posterior probability that s1 < s2 = ", pnorm(0,m,s),"\n") fd = NULL if (pic){ x = seq(min(s1$central,s2$central)-3*s, max(s1$central,s2$central)+3*s, length.out = 100) plot(x,dt((x-s1$central)/s1$scale,s1$df)/s1$scale, type='l', xlab=labels$xlab,ylab="Density",main=labels$main) lines(x,dt((x-s2$central)/s2$scale,s2$df)/s2$scale,lty=2) legend(max(s1$central,s2$central)+s,dnorm(max(s1$central,s2$central),max(s1$central,s2$central),s),c(labels$leg1,labels$leg2),lty=1:2 ,bty="n") # for comparison # lines(x,dnorm(x, s1$central, s1$scale),lty=3,lwd=3,col=3) # pt((x-s2$central)/s2$scale,s2$df) } } else { # most likely a binomial # crude approximation, waiting for a better man # s1 first n = s1$newdata$n probs = s1$s s1d = 0 for (i in 0:n){ s1d[i+1] = mean(dbinom(i,n,probs)) } n = s2$newdata$n probs = s2$s s2d = 0 for (i in 0:n){ s2d[i+1] = mean(dbinom(i,n,probs)) } p01 = 0 p10 = 0 p0011 = 0 if (s1$newdata$n == s2$newdata$n){ for (i in 0:n){ for (j in 0:n){ if (j>i) {p01 = p01 + s1d[i+1]*s2d[j+1]} } } p0011 = sum(s1d*s2d) p10 = 1-p01-p0011 cat("Posterior probability that s1 < s2 = ", p01,"\n") cat("Posterior probability that s1 > s2 = ", p10,"\n") cat("Posterior probability that s1 = s2 = ", p0011,"\n") } else { cat("Future data sizes (Ns) have to be the same to get Pr (s1 most likely s1 value: ", sum(s1d[which.max(s1d):n]),"\n") cat("Probability of s2 > most likely s1 value: ", sum(s2d[which.max(s1d):n]),"\n") fd = list(s1d=s1d, s2d=s2d) if (pic){ nn = max(length(s1d),length(s1d)) if (nn>49) {th = 'l'; ; putpoints = FALSE} else { th = 'h'; putpoints = TRUE}; plot(0:(length(s1d)-1), s1d, xlab=labels$xlab, main=labels$main, ylab="Probability",type=th,xlim=c(0,max((length(s1d)-1), (length(s2d)-1))) , ylim=c(0,max(s1d,s2d))) if (putpoints) points(0:(length(s1d)-1), s1d) lines(0:(length(s2d)-1), s2d,type=th,lty=2) if (putpoints) points(0:(length(s2d)-1), s2d,pch=2) if (putpoints){ legend(nn/4, max(s1d,s2d),c(labels$leg1,labels$leg2),lty=1:2 ,bty="n", pch=c(1,2)) } else { legend(nn/4, max(s1d,s2d),c(labels$leg1,labels$leg2),lty=1:2 ,bty="n") } } } fout = fd } # how to run # obs.glm.prob(s1,s2) ## skill plot skill.plot <- function(fit,fitnull){ # same as in skill but for just the individual variables n = length(fit$y) # null is always the same newdata = data.frame(fitnull$model[1,]) snull = obs.glm(fitnull,newdata) score = 0 scorenull = 0 for(i in 1:n){ newdata = data.frame(fit$model[i,]) s = obs.glm(fit,newdata) # crps if (!is.na(fit$y[i])) score[i] = crps(fit$y[i], c(s$central,s$scale)) if (!is.na(fit$y[i])) scorenull[i] = crps(fitnull$y[i], c(snull$central,snull$scale)) } mfrow <- set.mfrow(0, dim(fit$model)[2], 1) par(mfrow = mfrow) for (j in 1:dim(fit$model)[2]){ scatter.smooth(fit$model[, j],scorenull-score,axes=F,ylab="NULL - Full CRPS",xlab= names(fit$model)[j]) abline(h=0,lty=2) axis(1);axis(2) } } ## predictive plot predictive.plot <- function(fit,fitnull,which.var = 1, jitter.p=FALSE){ n = length(fit$model[, which.var]) plot(fit$model[, 1],fit$model[, which.var],type="n",axes=F,ylab=names(fit$model)[which.var],xlab= names(fit$model)[1]) lim = c(0,0) for(i in 1:n){ newdata = data.frame(fit$model[i,]) s = obs.glm(fit,newdata) #lim = qt(c(0.025,0.975),s$df,s$central/s$scale) lim = qnorm(c(0.025,0.975),s$central,s$scale) a = fit$model[i, which.var] if(jitter.p) a = jitter(a) lines(lim,c(a,a),col="#888888") points(fit$model[i, 1],a,pch=20,col="#148833") } axis(1);axis(2) # null model red line newdata = data.frame(fitnull$model[i,]) s = obs.glm(fitnull,newdata) lim = qnorm(c(0.025,0.975),s$central,s$scale) a = max(fit$model[, which.var]) a = a*1.01 lines(lim,c(a,a),col="#FF6666",lwd=2) } ## calibration plot calibration.plot <- function(fit,fitnull){ n = length(fit$model[, 1]) x = seq(.9*min(fit$y),1.1*max(fit$y),length.out = n*10) p = 0 s.m = 0 s.s = 0 for(i in 1:n){ newdata = data.frame(fit$model[i,]) s = obs.glm(fit,newdata) if (fit$family$family == "gaussian") { p[i] = pnorm(fit$y[i],s$central,s$scale) s.m = s.m + s$central s.s = s.s + s$scale } else { p[i] = mean(s$s) } } par(mfrow=c(2,2)) if (fit$family$family == "gaussian"){ newdata = data.frame(fitnull$model[i,]) s = obs.glm(fitnull,newdata) hist(p,breaks=seq(0,1,.05),probability=T) lines(density(p,from=0,to=1,adjust=.5),xlab="p",ylab="",main="",col="#118844") s.m=s.m/n s.s=s.s/n a = ecdf(fit$y) plot(pnorm(a$vals,s.m,s.s),a$ecdf,type='b', ylab="eCDF(y)",xlab="P(y)") abline(0,1,lty=2) lines(pnorm(a$vals,s$central,s$scale),a$ecdf,col=2,lty=2,type='b') plot(a$vals,pnorm(a$vals,s.m,s.s)-a$ecdf,type='l', ylab="mean(P(y))-eCDF(y)",xlab=names(fit$model)[1]) abline(h=0,lty=2) lines(a$vals,pnorm(a$vals,s$central,s$scale)-a$ecdf,col=2,lty=2) } else { hist(p,breaks=seq(0,1,.05),xlab=paste("Pr(",names(fit$model)[1],")=1"),main="") abline(v=mean(fit$y,na.rm=T),lty=3) plot(calibration(p,fit$y),type='b',xlab=paste("Pr(",names(fit$model)[1],")=1"),ylab="Observed relative freq.") abline(0,1) a = mean(p,na.rm=T)-mean(fit$y,na.rm=T) plot(0,a,ylab="mean(P(y))-eCDF(y)",xlab="", ylim=c(min(a,0), max(a,0)) , axes=F) axis(2) abline(h=0,lty=3) z = sort(p, index.return=T) howdy<-skill.curve(z$x,fit$y[z$ix],CI=TRUE) plot(z$x,howdy$k,type="b",ylab=expression( hat(K) [group( "","1/2" ,"")] ), xlab=paste("Pr(",names(fit$model)[1],")=1"), ylim=c(-1.2,.3)) #lines(p,howdy$ciHI,lty=2,col="red") #lines(p,howdy$ciLO,lty=2,col="red") abline(h=0,col="blue",lty=3) w<-sort(howdy$k,index.return=TRUE, decreasing=TRUE) abline(v=z$x[w$ix[1]],lty=2) abline(v=mean(p,na.rm=T),lty=3) } } ## calibration code for prob predictions of binary events calibration <- function(x,y,binned=TRUE){ if(binned){ x<-round(x*10)/10; } u<-sort(unique(x)) p<-0; for (i in 1:length(u)){ j<-(x==u[i]) p[i]<-sum(y[j]==1,na.rm=T)/length(is.na(y[j])) } answer<-list(x=u,y=p) } ## basic package modified ecdf <- function (x) { x <- sort(x) n <- length(x) vals <- unique(x) rval <- cumsum(tabulate(match(x, vals)))/n return(list(ecdf=rval,vals=vals)) } ## modified from the verification package from Pocernich ## complete rank probability score from Gneiting et al. crps <-function (obs, pred) { if (is.null(dim(pred)) & length(pred) == 2) { mu <- pred[1] sigma <- pred[2] } else { mu <- as.numeric(pred[, 1]) sigma <- as.numeric(pred[, 2]) } z <- (obs - mu)/sigma crps <- sigma * (z * (2 * pnorm(z) - 1) + 2 * dnorm(z) - 1/sqrt(pi)) } ## stolen from the MSBVAR package ## code by Patrick T. Brandt rmultnorm <- function (n, mu, vmat, tol = 1e-10) { p <- ncol(vmat) if (length(mu) != p) stop(paste("mu vector is the wrong length:", length(mu))) vs <- La.svd(vmat) vsqrt <- t(t(vs$vt) %*% (t(vs$u) * sqrt(vs$d))) ans <- matrix(rnorm(n * p), nrow = n) %*% vsqrt ans <- sweep(ans, 2, mu, "+") dimnames(ans) <- list(NULL, dimnames(vmat)[[2]]) ans } ## set.mfrow from the coda package by Martyn Plummer et al. #mfrow <- set.mfrow(Nchains = nchain(x), Nparms = nvar(x), nplots = trace + density) set.mfrow <- function (Nchains = 1, Nparms = 1, nplots = 1, sepplot = FALSE) { mfrow <- if (sepplot && Nchains > 1 && nplots == 1) { if (Nchains == 2) { switch(min(Nparms, 5), c(1, 2), c(2, 2), c(3, 2), c(4, 2), c(3, 2)) } else if (Nchains == 3) { switch(min(Nparms, 5), c(2, 2), c(2, 3), c(3, 3), c(2, 3), c(3, 3)) } else if (Nchains == 4) { if (Nparms == 1) c(2, 2) else c(4, 2) } else if (any(Nchains == c(5, 6, 10, 11, 12))) c(3, 2) else if (any(Nchains == c(7, 8, 9)) || Nchains >= 13) c(3, 3) } else { if (nplots == 1) { mfrow <- switch(min(Nparms, 13), c(1, 1), c(1, 2), c(2, 2), c(2, 2), c(3, 2), c(3, 2), c(3, 3), c(3, 3), c(3, 3), c(3, 2), c(3, 2), c(3, 2), c(3, 3)) } else { mfrow <- switch(min(Nparms, 13), c(1, 2), c(2, 2), c(3, 2), c(4, 2), c(3, 2), c(3, 2), c(4, 2), c(4, 2), c(4, 2), c(3, 2), c(3, 2), c(3, 2), c(4, 2)) } } return(mfrow) } skillTABLE <- function(n11,n01,n10,n00,theta=0.5,CI=FALSE,t=1,u=0) { # No error checking, but either n_ij can be a vector or theta can be; # both can NOT be vectors. # briggs matt@wmbriggs.com n<-n11+n01+n10+n00; # to determine which of RARE of COMM(ON) is used for each set z<-((n10+n11)/n); z<-(z<=theta)*1; q11<-n11/(n11+n01); p11<-(q11-u)/(t-u); q00<-n00/(n10+n00); p00<-(q00-(1-t))/(t-u); px<-(n11+n01)/n py<-(n11+n10-n*u)/(n*(t-u)) #k_RARE <- (n11*(1-u-theta*(t-u))-n01*(u+theta*(t-u)))/((n11+n10-n*u)*(1-theta)); k_RARE <- (p11-theta)*px/((1-theta)*py); G_RARE <- 2*n11*log(p11/theta) + 2*n01*log((1-p11)/(1-theta)); #k_COMM <- (n00*(t-(1-theta)*(t-u))-n10*(1-t+(1-theta)*(t-u)))/((n00+n01-n*(1-t))*theta); k_COMM <- (p00-1+theta)*(1-px)/(theta*(1-py)); G_COMM <- 2*n00*log(p00/(1-theta)) +2*n10*log((1-p00)/(theta)); k<-k_RARE*z+k_COMM*(1-z); G<-G_RARE*z+G_COMM*(1-z); # if k<0 then G=0 and p=.5 z<-(k>0)*1; G<-G*z; p<-pchisq(G,1,lower.tail=FALSE)/2; p<-p*z+0.5*(1-z); bigP<-p11*z+p00*(1-z) if(CI){ ciLO<-0;ciHI<-0; rootLO<-0;rootHI<-0; f_RARE<-function(p11,n11,n01,theta) 2*n11*log(p11/theta) + 2*n01*log((1-p11)/(1-theta))- (2*n11*log((n11/(n11+n01))/theta) + 2*n01*log((n01/(n11+n01))/(1-theta))-3.84); # 3.84 for climate # 4.25 for markov f_COMM<-function(p00,n00,n10,theta) 2*n00*log(p00/(1-theta)) + 2*n10*log((1-p00)/theta)- (2*n00*log((n00/(n00+n10))/(1-theta)) + 2*n10*log((n10/(n00+n10))/theta)-3.84); if(length(theta)>1){ for (i in 1:length(theta)){ zz<-((n10[i]+n11[i])/n[i]); zz<-(zz<=theta[i]); tol<-0.0001 if(zz){ hi<-p11[i]-tol lo<-p11[i]+tol; # print(c(0,hi,lo,n11[i],n01[i],theta[i])) if (is.na(hi) | is.na(lo)){ ciLO[i] <- NA; ciHI[i] <- NA; } else { rootLO[i]<-uniroot(f_RARE,lower=0.001,upper=hi,n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta[i])$root ciLO[i] <- (rootLO[i]-theta[i])*px[i]/((1-theta[i])*py[i]); if(n01[i]==0) { rootHI[i]<-1 } else { rootHI[i]<-uniroot(f_RARE,lower=lo,upper=(1-tol),n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta[i])$root } ciHI[i] <- (rootHI[i]-theta[i])*px[i]/((1-theta[i])*py[i]); } } else { hi<-p00[i]-tol lo<-p00[i]+tol; #print(c(1,hi,lo,n00[i],n10[i],theta[i])) if (is.na(hi) | is.na(lo)){ ciLO[i] <- NA; ciHI[i] <- NA; } else { rootLO[i]<-uniroot(f_COMM,lower=0.001,upper=hi,n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta[i])$root ciLO[i] <- (rootLO[i]-1+theta[i])*(1-px[i])/(theta[i]*(1-py[i])); rootHI[i]<-uniroot(f_COMM,lower=lo,upper=(1-tol),n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta[i])$root ciHI[i] <- (rootHI[i]-1+theta[i])*(1-px[i])/(theta[i]*(1-py[i])); } } } } else { for (i in 1:length(n11)){ zz<-((n10[i]+n11[i])/n[i]); zz<-(zz<=theta); tol<-0.0001 if(zz){ hi<-p11[i]-tol lo<-p11[i]+tol; # print(c(0,hi,lo,n11[i],n01[i],theta)) if ((is.na(hi) | is.na(lo)) | (hi<0 | lo<0)){ ciLO[i] <- NA; ciHI[i] <- NA; } else { #print(hi); rootLO[i]<-uniroot(f_RARE,lower=0.001,upper=hi,n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta)$root ciLO[i] <- (rootLO[i]-theta)*px[i]/((1-theta)*py[i]); if(n01[i]==0) { rootHI[i]<-1 } else { rootHI[i]<-uniroot(f_RARE,lower=lo,upper=(1-tol),n11=(n11[i]+tol),n01=(n01[i]+tol),theta=theta)$root } ciHI[i] <- (rootHI[i]-theta)*px[i]/((1-theta)*py[i]); } } else { hi<-p00[i]-tol lo<-p00[i]+tol; #print(c(1,hi,lo,n00[i],n10[i],p00[i],theta)) if ((is.na(hi) | is.na(lo)) | (hi<0 | lo<0)){ ciLO[i] <- NA; ciHI[i] <- NA; } else { rootLO[i]<-uniroot(f_COMM,lower=0.001,upper=hi,n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta)$root ciLO[i] <- (rootLO[i]-1+theta)*(1-px[i])/(theta*(1-py[i])); rootHI[i]<-uniroot(f_COMM,lower=lo,upper=(1-tol),n00=(n00[i]+tol),n10=(n10[i]+tol),theta=theta)$root ciHI[i] <- (rootHI[i]-1+theta)*(1-px[i])/(theta*(1-py[i])); } } } } } else { ciLO <- CI ciHI <- CI } answer<-list(z=z,n=n,k=k,G=G,p=p,theta=theta,ciLO=ciLO,ciHI=ciHI); } skill.curve <- function(x,y,theta=0.5,CI=FALSE,t=1,u=0) { if(length(theta)<2){ n<-length(x); x<-sort(x); n11<-0;n01<-0;n10<-0;n00<-0; for (i in 1:n){ n11[i]<-0;n01[i]<-0;n10[i]<-0;n00[i]<-0; n11[i]<-sum(((x>x[i])*1==1 & y==1),na.rm=T); n01[i]<-sum(((x>x[i])*1==1 & y==0),na.rm=T); n10[i]<-sum(((x>x[i])*1==0 & y==1),na.rm=T); n00[i]<-sum(((x>x[i])*1==0 & y==0),na.rm=T); } answer<-skillTABLE(n11,n01,n10,n00,theta,CI,t=1,u=0); } else { n<-length(theta); n11<-0;n01<-0;n10<-0;n00<-0; for (i in 1:n){ fit<-table((x>theta[i])*1,y) n11[i]<-0;n01[i]<-0;n10[i]<-0;n00[i]<-0; if(fit[2,2]) n11[i]<-fit[2,2]; if(fit[2,1]) n01[i]<-fit[2,1]; if(fit[1,2]) n10[i]<-fit[1,2]; if(fit[1,1]) n00[i]<-fit[1,1]; } answer<-skillTABLE(n11,n01,n10,n00,theta,CI,t=1,u=0); } } skill.curve.3d <- function(y,x1,x2){ require(lattice) mf <- match.call() if(is.factor(y)){ y = as.numeric(y) y = y-1 } w1 = sort(unique(x1)) w2 = sort(unique(x2)) n1 = length(w1) n2 = length(w2) skill = matrix(0,n1,n2) for (i in 1:n1){ for (j in 1:n2){ n11<-sum(((x1>w1[i] & x2>w2[j])*1==1 & y==1),na.rm=T); n01<-sum(((x1>w1[i] & x2>w2[j])*1==1 & y==0),na.rm=T); n10<-sum(((x1>w1[i] & x2>w2[j])*1==0 & y==1),na.rm=T); n00<-sum(((x1>w1[i] & x2>w2[j])*1==0 & y==0),na.rm=T); answer<-skillTABLE(n11,n01,n10,n00); skill[i,j] = answer$k } } #contour(w1,w2,skill,nlevels=20) # row.values = seq_len(nrow(x)) at = c(-2,-1.5,-1,-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1) print(levelplot(skill, at = at, col.regions=terrain.colors(n=length(at)), contour=TRUE, xlab=as.character(mf[3]),ylab=as.character(mf[4]), row.values=w1, col.values=w2, aspect="fill")) return(skill) }