R notes

Smooth spline fitting and interpolation

The R script below reads user-defined x,y values from an Excel file, fits smooth spline and generates Y values for the interpolated X values supplied by the user (e.g. X = seq(0,4,0.2) below)

Output shows various fits achieved by varying df parameter in smooth.spline()

write.table() command at the end writes and saves all the Y values in a tab-delimited text file

# This script needs xlsx package to read the data from an Excel file.

# Run install.packages("xlsx") and library("xlsx") before running this script

# for the first time to install and load the package in the workspace.

# Author: Ved P. Sharma, May 18, 2018

# Albert Einstein College of Medicine, New York

rm(list = ls())

if(dev.cur()!=1)

    dev.off()

dir = "D:/example/"

filename = "dataFile.xlsx"

file = paste(dir, filename, sep = "")

txtfile = paste(dir, "interpolatedData.txt", sep = "") # text file where interpolated X,Y values will be saved

par(mfrow=c(1,4))

X = seq(0,4,0.2)

x_DF = read.xlsx(file, "sheetName", rowIndex=56:77, colIndex=10, header=F)

y_DF = read.xlsx(file, "sheetName", rowIndex=56:77, colIndex=11, header=F)

x = x_DF$X10

y = y_DF$X11

plot(x, y, xlim = c(0, 4), main ="default")

f1 = smooth.spline(x, y)

lines(predict(f1, X), lwd=2)

plot(x, y, xlim = c(0, 4), main ="df = 20")

f2 = smooth.spline(x, y, df=20)

lines(predict(f2, X), lwd=2)

plot(x, y, xlim = c(0, 4), main ="df = 10")

f3 = smooth.spline(x, y, df=10)

lines(predict(f3, X), lwd=2)

plot(x, y, xlim = c(0, 4), main ="df = 5")

f4 = smooth.spline(x, y, df=5)

lines(predict(f4, X), lwd=2)

write.table(cbind(X, predict(f1, X)$y, predict(f2, X)$y, predict(f3, X)$y, predict(f4, X)$y), file = txtfile, sep="\t")

Output:

Line plot with error bars as shaded area

The example code below generates 101 x data points(from 0 to 10, 0.1 spaced), the corresponding y values from the sin curve (in red) and the inverted sin curve (in green). It also generates random error values using sample() function.

# Author: Ved P. Sharma

# Albert Einstein College of Medicine, New York

# generate first graph and the error values

x=seq(0,10,0.1)

y = sin(x*pi/10)

errorValues = c(0.05, 0.07, 0.09, 0.11)

error = sample(errorValues, 101, replace=TRUE)

y_plus_error = y+error

y_minus_error = y-error

dev.off()

dev.new(width=4, height=3, res=300)

par(mar=c(3.1,3.8,0.4,0.4), cex.main=1, font.main=1, las=1)

plot(x, y, type = "l", xlim = c(0, 10), ylim = c(-0.1, 1.1), lwd=2, xlab="")

polygon(c(x, rev(x)), c(y_plus_error, rev(y_minus_error)), col = adjustcolor("red", alpha.f=0.2), border=NA)

title(xlab="x", line=2)

lines(x, y, lwd = 2, col="red3")

# generate next graph and the error values

y2 = 1-sin(x*pi/10)

errorValues = c(0.05, 0.07, 0.09, 0.11)

error = sample(errorValues, 101, replace=TRUE)

y2_plus_error = y2+error

y2_minus_error = y2-error

par(new=TRUE) # to add new plot to the existing graphics window

plot(x, y2, type = "l", xlim = c(0, 10), ylim = c(-0.1, 1.1), lwd=2, xlab="", ylab="")

polygon(c(x, rev(x)), c(y2_plus_error, rev(y2_minus_error)), col = adjustcolor("green", alpha.f=0.2), border=NA)

lines(x, y2, lwd = 2, col="green4")

Output:

Scatter plot with vertical and horizontal error bars

We will use the plotCI function from the package plotrix. To install and load the package:

> install.packages("plotrix")

> library("plotrix")

set.seed(2)

x = 1:10

y = runif(10)

err.y = (1/10)*runif(10)

err.x = runif(10)

dev.off()

dev.new(width=4, height=3, res=300)

par(mar=c(3.1,3.8,0.4,0.4)) #default is c(5, 4, 4, 2) + 0.1

plotCI(x,y,err.y,xlim=c(0,10),ylim=c(0,1),lwd=2,las=1,ylab="",mgp=c(2,1,0),col="blue",scol="red", gap=T) # margin 2 lines for x title

plotCI(x,y,err.x, lwd=2,col="blue",scol="red", err="x",add=TRUE, gap=T)

title(ylab="y", line=3) # add y title at 3 lines of margin

Output:

Multiple plots in the same window

par(mfrow=c(rows,columns)) option can be used to plot multiple plots in the same window

dev.off()

dev.new(width=9, height=3, res=300)

par(mar=c(3.1,3.8,1,0.4), mfrow=c(1,3), mgp=c(2,1,0), font.main=1,las=1, cex.axis=1.2, cex.lab=1.5)

x = 1:10

for(i in 1:3) {

    set.seed(i)

    y = runif(10)

    err.y = (1/10)*runif(10)

    err.x = runif(10)

    title = paste("seed =", i)

    plotCI(x,y,err.y,xlim = c(0,10),ylim=c(0,1),main=title,ylab="",lwd=2,col="blue", scol="red", gap=T)

    plotCI(x,y,err.x, lwd=2,col="blue",scol="red", err="x",add=TRUE, gap=T)

    title(ylab="y", line=3) # add y title at 3 lines of margin

Output:

Useful plot commands

if(dev.cur() != 1)

dev.off() # closes the current plot window, if it exists

dev.new(width=4, height=3, res=300) # to create a new device (plot window) of a custom size and resolution (default resolution is 72)

par(new=TRUE) # insert this before the plot() command to add the plot to the existing plot window

par(mar=c(3.1,3.8,0.4,0.4)) # plot margins at (bottom, left, top, right); default is c(5, 4, 4, 2) + 0.1

las = 1 # to make axis titles horizontal

cex.axis=n and cex.lab=n in the plot() or par() to increase/decrease sizes for axes annotations and labels respectively

mpg: margins for the axis title (both x and y), axis labels and axis line. The default is c(3, 1, 0).

ylab="" # to turn y-axis label off

use title() function to add labels to the plot function. For example:

title(ylab="y axis label", line=3) # will print "y axis label" at 3 lines margins