see also: http://cran.r-project.org/doc/contrib/Baggott-refcard-v2.pdf
x <- c(1, 2, 3, 4, 5, 6, 7)
x
## [1] 1 2 3 4 5 6 7
x <- c(1:7)
x
## [1] 1 2 3 4 5 6 7
x <- 1:7
x
## [1] 1 2 3 4 5 6 7
5 + (6 + 7) * pi^2
## [1] 133.3
log(exp(1))
## [1] 1
log(1000, 10)
## [1] 3
sin(pi/3)^2 + cos(pi/3)^2
## [1] 1
log2(32)
## [1] 5
sqrt(2)
## [1] 1.414
seq(0, 5, length = 6)
## [1] 0 1 2 3 4 5
2^seq(0, 5, length = 6)
## [1] 1 2 4 8 16 32
Logical
x <- T
y <- F
x
## [1] TRUE
y
## [1] FALSE
logical operators
x & y
## [1] FALSE
x | y
## [1] TRUE
!x
## [1] FALSE
c(T, F, T, F) && c(T, T, F, F)
## [1] TRUE
c(T, F, T, F) || c(T, T, F, F)
## [1] TRUE
c(T, F, T, F) | c(T, T, F, F)
## [1] TRUE TRUE TRUE FALSE
c(T, F, T, F) & c(T, T, F, F)
## [1] TRUE FALSE FALSE FALSE
Numerical
a <- 5
b <- sqrt(2)
a
## [1] 5
b
## [1] 1.414
a == 5
## [1] TRUE
a < 5
## [1] FALSE
a <= 5
## [1] TRUE
Character
a <- "1"
b <- 1
a
## [1] "1"
b
## [1] 1
a <- "character"
b <- "a"
c <- a
a
## [1] "character"
b
## [1] "a"
c
## [1] "character"
letters
## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q"
## [18] "r" "s" "t" "u" "v" "w" "x" "y" "z"
a %in% letters
## [1] FALSE
"a" %in% letters
## [1] TRUE
factor
fac <- c("small", "small", "medium", "tall", "huge", "huge", "small", "tiny")
fac
## [1] "small" "small" "medium" "tall" "huge" "huge" "small" "tiny"
fac <- as.factor(fac)
levels <- c("small", "medium", "huge")
ordFac <- factor(sample(levels, 10, replace = T), levels = levels, ordered = T)
ordFac
## [1] small medium huge small medium medium small small small medium
## Levels: small < medium < huge
Vector
x <- c(5.2, 1.7, 6.3)
log(x)
## [1] 1.6487 0.5306 1.8405
y <- 1:5
z <- seq(from = 1, to = 1.4, by = 0.1)
y + z
## [1] 2.0 3.1 4.2 5.3 6.4
length(y)
## [1] 5
mean(y + z)
## [1] 4.2
Matrix
m <- matrix(1:12, 4, byrow = T)
m
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9
## [4,] 10 11 12
y <- -1:2
m.new <- m + y
m.new
## [,1] [,2] [,3]
## [1,] 0 1 2
## [2,] 4 5 6
## [3,] 8 9 10
## [4,] 12 13 14
t(m.new)
## [,1] [,2] [,3] [,4]
## [1,] 0 4 8 12
## [2,] 1 5 9 13
## [3,] 2 6 10 14
m.new
## [,1] [,2] [,3]
## [1,] 0 1 2
## [2,] 4 5 6
## [3,] 8 9 10
## [4,] 12 13 14
m/m
## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 1 1
## [3,] 1 1 1
## [4,] 1 1 1
dim(m)
## [1] 4 3
dim(t(m.new))
## [1] 3 4
missing values
x <- c(1, 2, 3, NA)
x + 3
## [1] 4 5 6 NA
“NaN”, +/- Inf
log(c(0, 1, 2))
## [1] -Inf 0.0000 0.6931
0/0
## [1] NaN
x <- c("a", "b", "c", "d", "e", "f", "g", "h")
x <- letters[1:8]
x[1]
## [1] "a"
x[3:5]
## [1] "c" "d" "e"
x[-(3:5)]
## [1] "a" "b" "f" "g" "h"
x[c(T, F, T, F, T, F, T, F)]
## [1] "a" "c" "e" "g"
x[x <= "d"]
## [1] "a" "b" "c" "d"
m[, 2]
## [1] 2 5 8 11
m[3, ]
## [1] 7 8 9
Lists are very flexible
Empl <- list(employee = "Anna", husband = "Fred", children = 3, child.ages = c(4,
7, 9))
str(Empl)
## List of 4
## $ employee : chr "Anna"
## $ husband : chr "Fred"
## $ children : num 3
## $ child.ages: num [1:3] 4 7 9
Empl[[4]]
## [1] 4 7 9
Empl$child.a
## [1] 4 7 9
Empl[4] #' a sublist consisting of the 4th component of Empl
## $child.ages
## [1] 4 7 9
names(Empl) <- letters[1:4]
Empl <- c(Empl, service = 8)
unlist(Empl) #' converts it to a vector. Mixed types will be converted to character, giving a character vector.
## a b c d1 d2 d3 service
## "Anna" "Fred" "3" "4" "7" "9" "8"
with matrices
x.mat <- matrix(c(3, -1, 2, 2, 0, 3), ncol = 2)
x.mat
## [,1] [,2]
## [1,] 3 2
## [2,] -1 0
## [3,] 2 3
dimnames(x.mat) <- list(c("L1", "L2", "L3"), c("R1", "R2"))
x.mat
## R1 R2
## L1 3 2
## L2 -1 0
## L3 2 3
subsetting
cw <- chickwts
`?`(chickwts)
## starting httpd help server ...
## done
Let's look at it:
str(cw)
## 'data.frame': 71 obs. of 2 variables:
## $ weight: num 179 160 136 227 217 168 108 124 143 140 ...
## $ feed : Factor w/ 6 levels "casein","horsebean",..: 2 2 2 2 2 2 2 2 2 2 ...
summary(cw)
## weight feed
## Min. :108 casein :12
## 1st Qu.:204 horsebean:10
## Median :258 linseed :12
## Mean :261 meatmeal :11
## 3rd Qu.:324 soybean :14
## Max. :423 sunflower:12
summary(cw$weight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 108 204 258 261 324 423
cw[3, 2]
## [1] horsebean
## Levels: casein horsebean linseed meatmeal soybean sunflower
View(cw)
Labels in data frames
labels(cw)
## [[1]]
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"
## [29] "29" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42"
## [43] "43" "44" "45" "46" "47" "48" "49" "50" "51" "52" "53" "54" "55" "56"
## [57] "57" "58" "59" "60" "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
## [71] "71"
##
## [[2]]
## [1] "weight" "feed"
you can edit the data using fix:
fix(cw)
but this is not a good idea! better edit manually:
x <- which(cw$weight == 136)
cw[x, "weight"] <- 140
plot(fac)
plot(ordFac)
pie(table(fac))
boxplot(weight ~ feed, cw)
pie(table(cw$feed))
plot(sin(seq(0, 2 * pi, length = 100)))
plot(sin(seq(0, 2 * pi, length = 100)), type = "l")
Histogram:
par(mfrow = c(1, 2)) #' set up multiple plots
simdata <- rchisq(100, 8)
hist(simdata) #' default number of bins
hist(simdata, breaks = 2) #' etc,4,20
Scatterplot
par(mfrow = c(1, 1))
z1 <- rnorm(50)
z2 <- rnorm(50)
rho <- 0.75 #' (or any number between -1 and 1)
x2 <- rho * z1 + sqrt(1 - rho^2) * z2
plot(z1, x2)
Grouped expressions in R
x = 1:9
if (length(x) <= 10) {
x <- c(x, 10:20)
print(x)
} else {
print(x[1])
}
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Loops in R
for (i in 1:length(x)) {
x[i] <- rnorm(1)
}
x
## [1] 0.56435 -0.09729 1.23762 -0.03614 0.43765 -1.12782 1.17394
## [8] 0.51190 0.40218 -1.36156 -1.95546 -0.03511 0.27251 0.46406
## [15] -0.80665 0.03622 -0.46300 0.84632 2.00269 0.36275
j = 1
while (j < 10) {
print(j)
j <- j + 2
}
## [1] 1
## [1] 3
## [1] 5
## [1] 7
## [1] 9
start with the letters “d” for density, “p” for cumulative dsitribution function, “q” for quantiles or “r” fro random numbers
curve(dnorm, -3, 3)
curve(dunif, -1, 2)
plot(dpois(1:10, 4), type = "h", ylab = "density", xlab = "value")
points(dpois(1:10, 4), type = "p")
hist(rexp(100))
hist(rexp(10000))
boxplot(rexp(100))
hist(rnorm(100))
hist(rnorm(10000))
plot.ecdf(rnorm(100))
curve(pnorm, col = "red", add = T)
add <- function(a, b) {
result = a + b
return(result)
}
add(2, 5)
## [1] 7
General Form of Functions function(arguments) { expression }
larger <- function(x, y) {
if (any(x < 0))
return(NA)
y.is.bigger <- y > x
x[y.is.bigger] <- y[y.is.bigger]
x
}
larger(2, 5)
## [1] 5
If you are in doubt…
help(predict)
`?`(predict)
or press “F1” when the cursor is “in” the function's name get some delivered data
data(cars)
plot(cars)
str(cars)
## 'data.frame': 50 obs. of 2 variables:
## $ speed: num 4 4 7 7 8 9 10 10 10 11 ...
## $ dist : num 2 10 4 22 16 10 18 26 34 17 ...
carLm <- lm(dist ~ speed, cars)
cars$fit <- predict(carLm, cars)
signif(cars$fit, 2)
## [1] -1.8 -1.8 9.9 9.9 14.0 18.0 22.0 22.0 22.0 26.0 26.0 30.0 30.0 30.0
## [15] 30.0 34.0 34.0 34.0 34.0 37.0 37.0 37.0 37.0 41.0 41.0 41.0 45.0 45.0
## [29] 49.0 49.0 49.0 53.0 53.0 53.0 53.0 57.0 57.0 57.0 61.0 61.0 61.0 61.0
## [43] 61.0 69.0 73.0 77.0 77.0 77.0 77.0 81.0
plot(dist ~ speed, cars)
plot(cars$speed, cars$dist, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
las = 1)
points(cars$speed, cars$fit, type = "l", col = "blue")
legend("topleft", c("obs.", "modelled"), pch = c(1, NA), lty = c(NA, 1), col = c("black",
"blue"))
x1 <- rnorm(100)
y1 <- rnorm(100)
t.test(x1, y1, var.equal = F, conf.level = 0.99)
##
## Welch Two Sample t-test
##
## data: x1 and y1
## t = 0.206, df = 194.3, p-value = 0.837
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
## -0.3306 0.3874
## sample estimates:
## mean of x mean of y
## 0.10441 0.07598
t.test(var.equal = F, conf.level = 0.99, x1, y1)
##
## Welch Two Sample t-test
##
## data: x1 and y1
## t = 0.206, df = 194.3, p-value = 0.837
## alternative hypothesis: true difference in means is not equal to 0
## 99 percent confidence interval:
## -0.3306 0.3874
## sample estimates:
## mean of x mean of y
## 0.10441 0.07598
x <- 1:10
y <- 2 * x + rnorm(10, 0, 1)
plot(x, y, type = "p") #' Try also l,b,s,o,h,n
#' axes=T, F
#' xlab='age', ylab='weight'
#' sub='sub title', main='main title'
#' xlim=c(0,12), ylim=c(-1,12)
loc <- locator(x, type = "p")
text(loc$x, loc$y, "Hey")
## Error: no coordinates were supplied
identify(x, y, labels = LETTERS[x])
## integer(0)
Plots for Multivariate Data
x <- 1:20/20
y <- 1:20/20
z <- outer(x, y, function(a, b) {
cos(10 * a * b)/(1 + a * b^2)
})
contour(x, y, z)
persp(x, y, z)
image(x, y, z)
install.packages(“rgl”) install.packages(“copula”) install.packages(“evd”) install.packages(“spcopula”, repos=“http://R-Forge.R-project.org”) peak, duration and volume triples library(spcopula) data(“simulatedTriples”) write.csv(triples, file=“simulatedTriples.csv”, row.names=F)
triples <- read.csv("simulatedTriples.csv")
library(rgl)
find demos:
demo(package = "rgl")
demo(abundance)
##
##
## demo(abundance)
## ---- ~~~~~~~~~
##
## > # RGL-Demo: animal abundance
## > # Authors: Oleg Nenadic, Daniel Adler
## > # $Id: abundance.r 564 2007-02-22 09:56:01Z dmurdoch $
## >
## > rgl.demo.abundance <- function()
## + {
## + open3d()
## + clear3d("all") # remove all shapes, lights, bounding-box, and restore viewpoint
## +
## + # Setup environment:
## + bg3d(col="#cccccc") # setup background
## + light3d() # setup head-light
## +
## + # Importing animal data (created with wisp)
## + terrain<-dget(system.file("demodata/region.dat",package="rgl"))
## + pop<-dget(system.file("demodata/population.dat",package="rgl"))
## +
## + # Define colors for terrain
## + zlim <- range(terrain)
## + zlen <- zlim[2] - zlim[1] + 1
## + colorlut <- terrain.colors(82)
## + col1 <- colorlut[9*sqrt(3.6*(terrain-zlim[1])+2)]
## +
## + # Set color to (water-)blue for regions with zero 'altitude'
## + col1[terrain==0]<-"#0000FF"
## +
## + # Add terrain surface shape (i.e. population density):
## + surface3d(
## + 1:100,seq(1,60,length=100),terrain,
## + col=col1,spec="#000000", ambient="#333333", back="lines"
## + )
## +
## + # Define colors for simulated populations (males:blue, females:red):
## + col2<-pop[,4]
## + col2[col2==0]<-"#3333ff"
## + col2[col2==1]<-"#ff3333"
## +
## + # Add simulated populations as sphere-set shape
## + spheres3d(
## + pop[,1],
## + pop[,2],
## + terrain[cbind( ceiling(pop[,1]),ceiling(pop[,2]*10/6) )]+0.5,
## + radius=0.2*pop[,3], col=col2, alpha=(1-(pop[,5])/10 )
## + )
## +
## + }
##
## > rgl.demo.abundance()
demo(envmap)
##
##
## demo(envmap)
## ---- ~~~~~~
##
## > # RGL-Demo: environment mapping
## > # Author: Daniel Adler
## > # $Id$
## >
## > rgl.demo.envmap <- function()
## + {
## + open3d()
## + # Clear scene:
## + clear3d("all")
## + light3d()
## + bg3d(sphere=T, color="white", back="filled"
## + , texture=system.file("textures/refmap.png",package="rgl")
## + )
## + data(volcano)
## +
## + surface3d( 10*1:nrow(volcano),10*1:ncol(volcano),5*volcano
## + , texture=system.file("textures/refmap.png",package="rgl")
## + , texenvmap=TRUE
## + , color = "white"
## + )
## + }
##
## > rgl.demo.envmap()
## Warning: RGL: PNG Pixmap Loader Warning: iCCP: known incorrect sRGB
## profile
## Warning: RGL: PNG Pixmap Loader Warning: iCCP: gamma value does not match
## sRGB
## Warning: RGL: PNG Pixmap Loader Warning: iCCP: known incorrect sRGB
## profile
## Warning: RGL: PNG Pixmap Loader Warning: iCCP: gamma value does not match
## sRGB
demo(lollipop3d)
##
##
## demo(lollipop3d)
## ---- ~~~~~~~~~~
##
## > cone3d <- function(base,tip,rad,n=30,...) {
## + degvec <- seq(0,2*pi,length=n)
## + ax <- tip-base
## + ## what do if ax[1]==0?
## + if (ax[1]!=0) {
## + p1 <- c(-ax[2]/ax[1],1,0)
## + p1 <- p1/sqrt(sum(p1^2))
## + if (p1[1]!=0) {
## + p2 <- c(-p1[2]/p1[1],1,0)
## + p2[3] <- -sum(p2*ax)
## + p2 <- p2/sqrt(sum(p2^2))
## + } else {
## + p2 <- c(0,0,1)
## + }
## + } else if (ax[2]!=0) {
## + p1 <- c(0,-ax[3]/ax[2],1)
## + p1 <- p1/sqrt(sum(p1^2))
## + if (p1[1]!=0) {
## + p2 <- c(0,-p1[3]/p1[2],1)
## + p2[3] <- -sum(p2*ax)
## + p2 <- p2/sqrt(sum(p2^2))
## + } else {
## + p2 <- c(1,0,0)
## + }
## + } else {
## + p1 <- c(0,1,0); p2 <- c(1,0,0)
## + }
## + ecoord2 <- function(theta) {
## + base+rad*(cos(theta)*p1+sin(theta)*p2)
## + }
## + for (i in 1:(n-1)) {
## + li <- ecoord2(degvec[i])
## + lj <- ecoord2(degvec[i+1])
## + triangles3d(c(li[1],lj[1],tip[1]),c(li[2],lj[2],tip[2]),c(li[3],lj[3],tip[3]),...)
## + }
## + }
##
## > lollipop3d <- function(data.x,data.y,data.z,surf.fun,surf.n=50,
## + xlim=range(data.x),
## + ylim=range(data.y),
## + zlim=range(data.z),
## + asp=c(y=1,z=1),
## + xlab=deparse(substitute(x)),
## + ylab=deparse(substitute(y)),
## + zlab=deparse(substitute(z)),
## + alpha.surf=0.4,
## + col.surf=fg,col.stem=c(fg,fg),
## + col.pt="gray",type.surf="line",ptsize,
## + lwd.stem=2,lit=TRUE,bg="white",fg="black",
## + col.axes=fg,col.axlabs=fg,
## + axis.arrow=TRUE,axis.labels=TRUE,
## + box.col=bg,
## + axes=c("lines","box")) {
## + axes <- match.arg(axes)
## + col.stem <- rep(col.stem,length=2)
## + x.ticks <- pretty(xlim)
## + x.ticks <- x.ticks[x.ticks>=min(xlim) & x.ticks<=max(xlim)]
## + x.ticklabs <- if (axis.labels) as.character(x.ticks) else NULL
## + y.ticks <- pretty(ylim)
## + y.ticks <- y.ticks[y.ticks>=min(ylim) & y.ticks<=max(ylim)]
## + y.ticklabs <- if (axis.labels) as.character(y.ticks) else NULL
## + z.ticks <- pretty(zlim)
## + z.ticks <- z.ticks[z.ticks>=min(zlim) & z.ticks<=max(zlim)]
## + z.ticklabs <- if (axis.labels) as.character(z.ticks) else NULL
## + if (!missing(surf.fun)) {
## + surf.x <- seq(xlim[1],xlim[2],length=surf.n)
## + surf.y <- seq(ylim[1],ylim[2],length=surf.n)
## + surf.z <- outer(surf.x,surf.y,surf.fun) ## requires surf.fun be vectorized
## + z.interc <- surf.fun(data.x,data.y)
## + zdiff <- diff(range(c(surf.z,data.z)))
## + } else {
## + z.interc <- rep(min(data.z),length(data.x))
## + zdiff <- diff(range(data.z))
## + }
## + xdiff <- diff(xlim)
## + ydiff <- diff(ylim)
## + y.adj <- if (asp[1]<=0) 1 else asp[1]*xdiff/ydiff
## + data.y <- y.adj*data.y
## + y.ticks <- y.adj*y.ticks
## + ylim <- ylim*y.adj
## + ydiff <- diff(ylim)
## + z.adj <- if (asp[2]<=0) 1 else asp[2]*xdiff/zdiff
## + data.z <- z.adj*data.z
## + if (!missing(surf.fun)) {
## + surf.y <- y.adj*surf.y
## + surf.z <- z.adj*surf.z
## + }
## + z.interc <- z.adj*z.interc
## + z.ticks <- z.adj*z.ticks
## + zlim <- z.adj*zlim
## + open3d()
## + clear3d("all")
## + light3d()
## + bg3d(color=c(bg,fg))
## + if (!missing(surf.fun))
## + surface3d(surf.x,surf.y,surf.z,alpha=alpha.surf,
## + front=type.surf,back=type.surf,
## + col=col.surf,lit=lit)
## + if (missing(ptsize)) ptsize <- 0.02*xdiff
## + ## draw points
## + spheres3d(data.x,data.y,data.z,r=ptsize,lit=lit,color=col.pt)
## + ## draw lollipops
## + apply(cbind(data.x,data.y,data.z,z.interc),1,
## + function(X) {
## + lines3d(x=rep(X[1],2),
## + y=rep(X[2],2),
## + z=c(X[3],X[4]),
## + col=ifelse(X[3]>X[4],col.stem[1],
## + col.stem[2]),lwd=lwd.stem)
## + })
## + bbox <- par3d("bbox")
## + if (axes=="box") {
## + bbox3d(xat=x.ticks,xlab=x.ticklabs,
## + yat=y.ticks,ylab=y.ticklabs,
## + zat=z.ticks,zlab=z.ticklabs,lit=lit)
## + } else if (axes=="lines") { ## set up axis lines
## + bbox <- par3d("bbox")
## + axis3d(edge="x",at=x.ticks,labels=x.ticklabs,
## + col=col.axes,arrow=axis.arrow)
## + axis3d(edge="y",at=y.ticks,labels=y.ticklabs,
## + col=col.axes,arrow=axis.arrow)
## + axis3d(edge="z",at=z.ticks,labels=z.ticklabs,
## + col=col.axes,arrow=axis.arrow)
## + box3d(col=col.axes)
## + }
## + decorate3d(xlab=xlab, ylab=ylab, zlab=zlab, box=FALSE, axes=FALSE, col=col.axlabs)
## + }
##
## > x <- 1:5
##
## > y <- x*10
##
## > z <- (x+y)/20
##
## > open3d()
## wgl
## 3
##
## > spheres3d(x,y,z)
##
## > axes3d()
##
## > set.seed(1001)
##
## > x <- runif(30)
##
## > y <- runif(30,max=2)
##
## > dfun <- function(x,y) { 2*x+3*y+2*x*y+3*y^2 }
##
## > z <- dfun(x,y)+rnorm(30,sd=0.5)
##
## > ## lollipops only
## > lollipop3d(x,y,z)
##
## > ## lollipops plus theoretical surface
## >
## > lollipop3d(x,y,z,dfun,col.pt="red",col.stem=c("red","blue"))
##
## > ## lollipops plus regression fit
## >
## > linmodel <- lm(z~x+y)
##
## > dfun <- function(x,y) {predict(linmodel,newdata=data.frame(x=x,y=y))}
##
## > lollipop3d(x,y,z,dfun,col.pt="red",col.stem=c("red","blue"))
##
## > ####
plot(triples)
plot3d(triples)
fit margins
library(evd)
gevFit <- fgev(triples[, 1])$estimate
dgevFun <- function(x) dgev(x, gevFit[1], gevFit[2], gevFit[3])
loglik
sum(log(dgevFun(triples[, 1]))) #' -2682
## [1] -2682
fit a log-normal
optFunLogNorm <- function(param) {
-sum(log(dlnorm(triples[, 1], param[1], param[2])))
}
lnormFit <- optim(c(1, 1), optFunLogNorm)$par
dlnormFun <- function(x) dlnorm(x, lnormFit[1], lnormFit[2])
loglik
sum(log(dlnormFun(triples[, 1]))) #' -2730
## [1] -2730
hist(triples[, 1], freq = F, n = 20, ylim = c(0, 0.012), xlim = c(-50, 500))
curve(dgevFun, add = T, col = "red")
curve(dlnormFun, add = T, col = "blue")
#' qq-plots
qqplot(qgev(ppoints(500), gevFit[1], gevFit[2], gevFit[3]), triples[, 1])
qqline(triples[, 1], distribution = function(x) qgev(x, gevFit[1], gevFit[2],
gevFit[3]))
qqplot(qlnorm(ppoints(500), lnormFit[1], lnormFit[2]), triples[, 1])
qqline(triples[, 1], distribution = function(x) qlnorm(x, lnormFit[1], lnormFit[2]))
assessing dependence
cor(triples)
## peak duration volume
## peak 1.0000 0.4958 0.8688
## duration 0.4958 1.0000 0.7481
## volume 0.8688 0.7481 1.0000
cor(triples, method = "kendall")
## peak duration volume
## peak 1.0000 0.4196 0.8479
## duration 0.4196 1.0000 0.5384
## volume 0.8479 0.5384 1.0000
cor(triples, method = "spearman")
## peak duration volume
## peak 1.0000 0.5843 0.9611
## duration 0.5843 1.0000 0.7230
## volume 0.9611 0.7230 1.0000
striping margins
rtTriples <- apply(triples[, c(1, 3)], 2, function(x) rank(x)/(nrow(triples) +
1))
plot(rtTriples, asp = 1)
library(copula)
copFit <- fitCopula(claytonCopula(1), rtTriples, method = "ml")
persp(copFit@copula, dCopula, zlim = c(0, 20))
## Warning: surface extends beyond the box
library(sp)
data(meuse)
str(meuse)
## 'data.frame': 155 obs. of 14 variables:
## $ x : num 181072 181025 181165 181298 181307 ...
## $ y : num 333611 333558 333537 333484 333330 ...
## $ cadmium: num 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num 85 81 68 81 48 61 31 29 37 24 ...
## $ lead : num 299 277 199 116 117 137 132 150 133 80 ...
## $ zinc : num 1022 1141 640 257 269 ...
## $ elev : num 7.91 6.98 7.8 7.66 7.48 ...
## $ dist : num 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om : num 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num 50 30 150 270 380 470 240 120 240 420 ...
coordinates(meuse) <- ~x + y
str(meuse)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 155 obs. of 12 variables:
## .. ..$ cadmium: num [1:155] 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## .. ..$ copper : num [1:155] 85 81 68 81 48 61 31 29 37 24 ...
## .. ..$ lead : num [1:155] 299 277 199 116 117 137 132 150 133 80 ...
## .. ..$ zinc : num [1:155] 1022 1141 640 257 269 ...
## .. ..$ elev : num [1:155] 7.91 6.98 7.8 7.66 7.48 ...
## .. ..$ dist : num [1:155] 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## .. ..$ om : num [1:155] 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## .. ..$ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## .. ..$ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## .. ..$ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## .. ..$ dist.m : num [1:155] 50 30 150 270 380 470 240 120 240 420 ...
## ..@ coords.nrs : int [1:2] 1 2
## ..@ coords : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] 178605 329714 181390 333611
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
## .. .. ..@ projargs: chr NA
meuse[1, ]
## coordinates cadmium copper lead zinc elev dist om ffreq soil
## 1 (181072, 333611) 11.7 85 299 1022 7.909 0.001358 13.6 1 1
## lime landuse dist.m
## 1 1 Ah 50
meuse[, "zinc"]
## coordinates zinc
## 1 (181072, 333611) 1022
## 2 (181025, 333558) 1141
## 3 (181165, 333537) 640
## 4 (181298, 333484) 257
## 5 (181307, 333330) 269
## 6 (181390, 333260) 281
## 7 (181165, 333370) 346
## 8 (181027, 333363) 406
## 9 (181060, 333231) 347
## 10 (181232, 333168) 183
## 11 (181191, 333115) 189
## 12 (181032, 333031) 251
## 13 (180874, 333339) 1096
## 14 (180969, 333252) 504
## 15 (181011, 333161) 326
## 16 (180830, 333246) 1032
## 17 (180763, 333104) 606
## 18 (180694, 332972) 711
## 19 (180625, 332847) 735
## 20 (180555, 332707) 1052
## 21 (180642, 332708) 673
## 22 (180704, 332717) 402
## 23 (180704, 332664) 343
## 24 (181153, 332925) 218
## 25 (181147, 332823) 200
## 26 (181167, 332778) 194
## 27 (181008, 332777) 207
## 28 (180973, 332687) 180
## 29 (180916, 332753) 240
## 30 (181352, 332946) 180
## 31 (181133, 332570) 208
## 32 (180878, 332489) 198
## 33 (180829, 332450) 250
## 34 (180954, 332399) 192
## 35 (180956, 332318) 213
## 37 (180710, 332330) 321
## 38 (180632, 332445) 569
## 39 (180530, 332538) 833
## 40 (180478, 332578) 906
## 41 (180383, 332476) 1454
## 42 (180494, 332330) 298
## 43 (180561, 332193) 167
## 44 (180451, 332175) 176
## 45 (180410, 332031) 258
## 46 (180355, 332299) 746
## 47 (180292, 332157) 746
## 48 (180283, 332014) 464
## 49 (180282, 331861) 365
## 50 (180270, 331707) 282
## 51 (180199, 331591) 375
## 52 (180135, 331552) 222
## 53 (180237, 332351) 812
## 54 (180103, 332297) 1548
## 55 (179973, 332255) 1839
## 56 (179826, 332217) 1528
## 57 (179687, 332161) 933
## 58 (179792, 332035) 432
## 59 (179902, 332113) 550
## 60 (180100, 332213) 1571
## 61 (179604, 332059) 1190
## 62 (179526, 331936) 907
## 63 (179495, 331770) 761
## 64 (179489, 331633) 659
## 65 (179414, 331494) 643
## 66 (179334, 331366) 801
## 67 (179255, 331264) 784
## 69 (179470, 331125) 1060
## 75 (179692, 330933) 119
## 76 (179852, 330801) 778
## 79 (179140, 330955) 703
## 80 (179128, 330867) 676
## 81 (179065, 330864) 793
## 82 (179007, 330727) 685
## 83 (179110, 330758) 593
## 84 (179032, 330645) 549
## 85 (179095, 330636) 680
## 86 (179058, 330510) 539
## 87 (178810, 330666) 560
## 88 (178912, 330779) 1136
## 89 (178981, 330924) 1383
## 90 (179076, 331005) 1161
## 123 (180151, 330353) 1672
## 160 (179211, 331175) 765
## 163 (181118, 333214) 279
## 70 (179474, 331304) 241
## 71 (179559, 331423) 317
## 91 (179022, 330873) 545
## 92 (178953, 330742) 505
## 93 (178875, 330516) 420
## 94 (178803, 330349) 332
## 95 (179029, 330394) 400
## 96 (178605, 330406) 553
## 97 (178701, 330557) 577
## 98 (179547, 330245) 155
## 99 (179301, 330179) 224
## 100 (179405, 330567) 180
## 101 (179462, 330766) 226
## 102 (179293, 330797) 186
## 103 (179180, 330710) 198
## 104 (179206, 330398) 187
## 105 (179618, 330458) 199
## 106 (179782, 330540) 157
## 108 (179980, 330773) 203
## 109 (180067, 331185) 143
## 110 (180162, 331387) 136
## 111 (180451, 331473) 117
## 112 (180328, 331158) 113
## 113 (180276, 330963) 130
## 114 (180114, 330803) 192
## 115 (179881, 330912) 240
## 116 (179774, 330921) 221
## 117 (179657, 331150) 140
## 118 (179731, 331245) 128
## 119 (179717, 331441) 166
## 120 (179446, 331422) 191
## 121 (179524, 331565) 232
## 122 (179644, 331730) 203
## 124 (180321, 330366) 722
## 125 (180162, 331837) 210
## 126 (180029, 331720) 198
## 127 (179797, 331919) 139
## 128 (179642, 331955) 253
## 129 (179849, 332142) 703
## 130 (180265, 332297) 832
## 131 (180107, 332101) 262
## 132 (180462, 331947) 142
## 133 (180478, 331822) 119
## 134 (180347, 331700) 152
## 135 (180862, 333116) 415
## 136 (180700, 332882) 474
## 161 (180201, 331160) 126
## 162 (180173, 331923) 210
## 137 (180923, 332874) 220
## 138 (180467, 331694) 133
## 140 (179917, 331325) 141
## 141 (179822, 331242) 158
## 142 (179991, 331069) 129
## 143 (179120, 330578) 206
## 144 (179034, 330561) 451
## 145 (179085, 330433) 296
## 146 (179236, 330046) 189
## 147 (179456, 330072) 154
## 148 (179550, 329940) 169
## 149 (179445, 329807) 403
## 150 (179337, 329870) 471
## 151 (179245, 329714) 612
## 152 (179024, 329733) 601
## 153 (178786, 329822) 783
## 154 (179135, 329890) 258
## 155 (179030, 330082) 214
## 156 (179184, 330182) 166
## 157 (179085, 330292) 496
## 158 (178875, 330311) 342
## 159 (179466, 330381) 162
## 164 (180627, 330190) 375
meuse@data$zinc
## [1] 1022 1141 640 257 269 281 346 406 347 183 189 251 1096 504
## [15] 326 1032 606 711 735 1052 673 402 343 218 200 194 207 180
## [29] 240 180 208 198 250 192 213 321 569 833 906 1454 298 167
## [43] 176 258 746 746 464 365 282 375 222 812 1548 1839 1528 933
## [57] 432 550 1571 1190 907 761 659 643 801 784 1060 119 778 703
## [71] 676 793 685 593 549 680 539 560 1136 1383 1161 1672 765 279
## [85] 241 317 545 505 420 332 400 553 577 155 224 180 226 186
## [99] 198 187 199 157 203 143 136 117 113 130 192 240 221 140
## [113] 128 166 191 232 203 722 210 198 139 253 703 832 262 142
## [127] 119 152 415 474 126 210 220 133 141 158 129 206 451 296
## [141] 189 154 169 403 471 612 601 783 258 214 166 496 342 162
## [155] 375
meuse@data$soil
## [1] 1 1 1 2 2 2 2 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
## [36] 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 2
## [106] 3 3 3 3 3 3 3 3 2 1 1 1 2 2 2 1 1 1 1 1 2 2 2 1 1 3 2 1 2 2 2 3 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3
## Levels: 1 2 3
plot this data:
plot(meuse)
spplot(meuse)
## Warning: non-vector columns will be ignored
spplot(meuse, "zinc", col.regions = bpy.colors())
spDists(meuse[1:10, ])
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.00 70.84 118.8 259.2 366.3 473.6 258.3 252.0 380.2 471.0
## [2,] 70.84 0.00 141.6 282.9 362.6 471.2 234.4 195.0 328.9 441.5
## [3,] 118.85 141.57 0.0 143.2 251.0 356.9 167.0 222.1 323.5 375.0
## [4,] 259.24 282.85 143.2 0.0 154.3 242.2 175.2 296.8 347.4 322.8
## [5,] 366.31 362.64 251.0 154.3 0.0 108.6 147.5 281.9 266.1 178.5
## [6,] 473.63 471.20 356.9 242.2 108.6 0.0 250.4 377.3 331.3 182.8
## [7,] 258.32 234.40 167.0 175.2 147.5 250.4 0.0 138.2 174.2 212.8
## [8,] 252.05 195.01 222.1 296.8 281.9 377.3 138.2 0.0 136.1 282.9
## [9,] 380.19 328.87 323.5 347.4 266.1 331.3 174.2 136.1 0.0 183.2
## [10,] 471.01 441.53 375.0 322.8 178.5 182.8 212.8 282.9 183.2 0.0