Grupo R : RMyP

library(rms)
 
b3$BroodC <- NULL ##Eliminar una columna (datadist da problemas si en alguna columna hay NAs
dd <- datadist(b3); options(datadist='dd')  #opción de longitud necesaria para el Gls(  Gls {rms})
Form <- formula(log(Media) ~ logLc+fBrood2)
Ident<-varIdent( ~1 | fBrood2) 
M2<- Gls(Form,data=b3,weights = Ident)
 
###Ploteamos los residuales estandarizados
 
E <- resid(M, type = "normalized")
Fit <- fitted(M)
op <- par(mfrow = c(1, 3))
plot(x = Fit, y = E,
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals versus fitted values")
text(Fit,E,b2$Cod,cex=0.6)
hist(E, nclass = 15)
plot(E ~fBrood2, data = b3, ylab = "Normalised residuals")
text(as.numeric(b3$fBrood2),E,b3$Cod,cex=0.6) ##Para saber identificar las observaciones
par(op) 
 
####Plot de los predict values and 95% confidence intervals in Gls
 
plot(Predict(M2))

boxplot.with.outlier.label <- function(y, label_name, ..., spread_text = T, data, plot = T, range = 1.5, label.col = "blue", 
										jitter_if_duplicate = T, jitter_only_positive_duplicates = F)
{	
	# jitter_if_duplicate - will jitter (Actually just add a bit of numbers) so to be able to decide on which location to plot the label when having identical variables...
	require(plyr) # for is.formula and ddply
 
	# a function to jitter data in case of ties in Y's
	jitter.duplicate <- function(x, only_positive = F)
	{
		if(only_positive) {
			ss <- x > 0
		} else {
			ss <- T
		}	
		ss_dup <- duplicated(x[ss])
		# ss <- ss & ss_dup
		temp_length <- length(x[ss][ss_dup])	
		x[ss][ss_dup] <- x[ss][ss_dup] + seq(from = 0.00001, to = 0.00002, length.out = temp_length)
		x
	}
	# jitter.duplicate(c(1:5))
	# jitter.duplicate(c(1:5,5,2))
	# duplicated(jitter.duplicate(c(1:5,5,2)))
	# jitter.duplicate(c(0,0,1:5,5,2))
	# duplicated(jitter.duplicate(c(0,0,1:5,5,2)))
 
 
 
	# handle cases where 
	if(jitter_if_duplicate) {
		# warning("duplicate jutter of values in y is ON")
		if(!missing(data)) {	#e.g: we DO have data
			# if(exists("y") && is.formula(y)) {		# F && NULL # F & NULL
			y_name <- as.character(substitute(y))	# I could have also used as.list(match.call())
												# credit to Uwe Ligges and Marc Schwartz for the help
												# https://mail.google.com/mail/?shva=1#inbox/12dd7ca2f9bfbc39
			if(length(y_name) > 1) {	# then it is a formula (for example: "~", "y", "x"
				model_frame_y <- model.frame(y, data = data)
				temp_y <- model_frame_y[,1]
				temp_y  <- jitter.duplicate(temp_y, jitter_only_positive_duplicates)	# notice that the default of the function is to work only with positive values...
				# the_txt <- paste(names(model_frame_y)[1], "temp_y", sep = "<<-") # wrong...
				the_txt <- paste("data['",names(model_frame_y)[1],"'] <- temp_y", sep = "")				
				eval(parse(text = the_txt))	# jutter out y var so to be able to handle identical values.
			} else {	# this isn't a formula
				data[,y_name] <- jitter.duplicate(data[,y_name], jitter_only_positive_duplicates)
				y <- data[,y_name]	# this will make it possible for boxplot(y, data) to work later (since it is not supposed to work with data when it's not a formula, but now it does :))
			}		
		} else {	# there is no "data"		 
			if(is.formula(y)) { # if(exists("y") && is.formula(y)) {		# F && NULL # F & NULL
				temp_y <- model.frame(y)[,1]
				temp_y  <- jitter.duplicate(temp_y, jitter_only_positive_duplicates)	# notice that the default of the function is to work only with positive values...
				environment(y) <- new.env()
				assign(names(model.frame(y))[1], temp_y, environment(y))
					# Credit and thanks for doing this goes to Niels Richard Hansen (2 Jan 30, 2011)
					# http://r.789695.n4.nabble.com/environment-question-changing-variables-from-a-formula-through-model-frame-td3246608.html
				# warning("Your original variable (in the global environemnt) was just jittered.")	# maybe I should add a user input before doing this....
				# the_txt <- paste(names(model_frame_y)[1], "temp_y", sep = "<<-")
				# eval(parse(text = the_txt))	# jutter out y var so to be able to handle identical values.
			} else {
				y <- jitter.duplicate(y, jitter_only_positive_duplicates)
			}		
		}
	}
	# the_txt <- paste("print(",names(model_frame_y)[1], ")")
	# eval(parse(text = the_txt))	# jutter out y var so to be able to handle identical values.
	# print(ls())
 
 
	# y should be a formula of the type: y~x, y~a*b
	# or it could be simply y
	if(missing(data)) {
			boxdata <- boxplot(y, plot = plot,range = range ,...)
		} else {
			boxdata <- boxplot(y, plot = plot,data = data, range = range ,...)
		}
	if(length(boxdata$names) == 1 && boxdata$names =="") boxdata$names <- 1	# this is for cases of type: boxplot(y) (when there is no dependent group)
	if(length(boxdata$out) == 0 ) stop("No outliers detected for this boxplot")
 
	if(!missing(data)) attach(data)	# this might lead to problams I should check out for alternatives for using attach here...
 
 
	# creating a data.frame with information from the boxplot output about the outliers (location and group)
	boxdata_group_name <- factor(boxdata$group)
	levels(boxdata_group_name) <- boxdata$names[as.numeric(levels(boxdata_group_name))]	# the subseting is for cases where we have some sub groups with no outliers
	boxdata_outlier_df <- data.frame(group = boxdata_group_name, y = boxdata$out, x = boxdata$group)
 
 
	# Let's extract the x,y variables from the formula:
	if(is.formula(y))
	{
		model_frame_y <- model.frame(y)
		y <- model_frame_y[,1]
		x <- model_frame_y[,-1]
		if(!is.null(dim(x))) {	# then x is a matrix/data.frame of the type x1*x2*..and so on - and we should merge all the variations...
			x <- apply(x,1, paste, collapse = ".")
		}
	} else {
		# if(missing(x)) x <- rep(1, length(y))
		x <- rep(1, length(y))	# we do this in case y comes as a vector and without x
	}	
 
	# and put all the variables (x, y, and outlier label name) into one data.frame
	DATA <- data.frame(label_name, x ,y)
	if(!missing(data)) detach(data)	# we don't need to have "data" attached anymore.
 
	# let's only keep the rows with our outliers 
	boxplot.outlier.data <- function(xx, y_name = "y")
	{
		y <- xx[,y_name]
		boxplot_range <- range(boxplot.stats(y, coef = range )$stats)
		ss <- (y < boxplot_range[1]) | (y > boxplot_range[2])
		return(xx[ss,])	
	}
	outlier_df <-ddply(DATA, .(x), boxplot.outlier.data)
 
 
	# create propor x/y locations to handle over-laping dots...
	if(spread_text) {
		# credit: Greg Snow
		require(TeachingDemos)		
		temp_x <- boxdata_outlier_df[,"x"]
		temp_y1 <- boxdata_outlier_df[,"y"]
		temp_y2 <- temp_y1
		for(i in unique(temp_x))
		{
			tmp <- temp_x == i
			temp_y2[ tmp ] <- spread.labs( temp_y2[ tmp ], 1.3*strheight('A'), maxiter=6000, stepsize = 0.05) #, min=0 )
		}
 
	}
 
 
	# max(strwidth(c("asa", "a"))
	# move_text_right <- max(strwidth(outlier_df[,"label_name"]))	
 
	# plotting the outlier labels :)  (I wish there was a non-loop wise way for doing this)
	for(i in seq_len(dim(boxdata_outlier_df)[1]))
	{
		# ss <- (outlier_df[,"x"]  %in% boxdata_outlier_df[i,]$group) & (outlier_df[,"y"] %in% boxdata_outlier_df[i,]$y)
 
		# if(jitter_if_duplicate) {
			# ss <- (outlier_df[,"x"]  %in% boxdata_outlier_df[i,]$group) & closest.number(outlier_df[,"y"]  boxdata_outlier_df[i,]$y)
		# } else {
		ss <- (outlier_df[,"x"]  %in% boxdata_outlier_df[i,]$group) & (outlier_df[,"y"] %in% boxdata_outlier_df[i,]$y)
		# }
 
		current_label <- outlier_df[ss,"label_name"]
		temp_x <- boxdata_outlier_df[i,"x"]
		temp_y <- boxdata_outlier_df[i,"y"]		
		# cbind(boxdata_outlier_df,		temp_y2)
		# outlier_df
 
 
		if(spread_text) {
			temp_y_new <- temp_y2[i] # not ss			
			move_text_right <- strwidth(current_label)
			text( temp_x+move_text_right, temp_y_new, current_label, col = label.col)			
			# strwidth
			segments( temp_x+(move_text_right/6), temp_y, temp_x+(move_text_right*.47), temp_y_new )
		} else {
			text(temp_x, temp_y, current_label, pos = 4, col = label.col)
		}		
	}
 
	# outputing some of the information we collected
	list(boxdata = boxdata, boxdata_outlier_df = boxdata_outlier_df, outlier_df=outlier_df)
}

boxplot.with.outlier.label (Lip~Egon,Cod, varwidth=T, ylab="Lipidos(ug/mg)", xlab="Estado gonada", data=P1_98gon, ylim=c(50,350))
sample.size <- tapply(P1_98gon$Lip, P1_98gon$Egon, length) ####Codigo  para determinar el tamaño muestral
ss.ch <- paste("N=", sample.size, sep="")            ####Codigo  pegar N= a los datos de tamaño muestral
mtext(ss.ch, at=1:length(unique(P1_98gon$Lip)), line=2, side=3)   ####Codigo  poner una linea de texto en el eje x en cada boxplot

boxplot.with.outlier.label (Lip~Egon,Cod, varwidth=T, ylab="Lipidos(ug/mg)", xlab="Estado gonada", data=A99gon, ylim=c(50,350), spread_text = F) #añadimos spread_text=F

 
#introducir tabla de datos 
 
 
   z<- structure(list(Intermoult = structure (c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 5L, 6L, 7L,
8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L),.Label = c("Instar1","Instar2","Instar3 ","Instar4","Instar5 ","Instar6","Instar7","Instar8",
"Instar9","Instar10","Instar11","Instar12","Instar13","Instar14" ,"Instar15 ","Instar16"),class = "factor"), 
Way = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L , 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), 
.Label = c(" Initial 200 larvae", " 80 larvae transferred from the mass culture"), class = "factor"),
Mean = c(8.4,7.6,9.4,11.3,16.7,10.3,17.3,14.2,13.4,14.4,14.2,16.4,19.0,22.3,27.0,NA,16.7,17.0,18.3,15.0,15.6,15.1,14.7,17.6,21.1,18.8,27.0),
SD = c(0.5,1.3,1.0,1.6,2.4,1.7,4.0,3.5,2.1,4.1,2.5,2.2,3.5,3.9,0.0,NA,5.0,6.7,5.8,4.7,4.0,2.6,2.1,3.2,4.5,1.3,0.0),
upp = c(11,19,13,15,20,15,22,22,17,25,19,20,25,28,27,NA,24,33,30,28,26,20,19,23,30,20,27),
low = c(8,6,8,9,13,9,12,10,11,12,10,13,14,19,27,NA,12,10,11,8,11,11,11,13,16,17,27)) ,
class = "data.frame",  ,row.names = c(NA, -27L) )
colnames(z)<- c("Intermoult", "Way", "Mean", "SD", "upper", "lower")
 
###Crear plot
 
library(ggplot2)
 
pdf(file = "intermuda2.pdf", width = 6, height = 6, dpi= 600)    #Guardarlo como pdf
 
 
 
p<-ggplot(z,aes(x=Intermoult,y=Mean,ymin=Mean-SD,ymax=Mean+SD,
             groups=Way,colour=Way))+  
  geom_point(position=position_dodge(width=0.5))+ 
  geom_pointrange(width=0.5,position=position_dodge(width=0.5))+
   labs(x="Intermoult",y="Intermoult Period (Days)")+ theme_bw()+ opts(legend.position = c(0.18,0.6))+ #colocar la leyenda 
 scale_colour_manual("", c(" Initial 200 larvae" = "#5CB0E2", " 80 larvae transferred from the mass culture" = "darkblue"))
 
   p + ylim(c(-1, 30))
 
 
 
 
 #Características del plot#########################
 
last_plot() + opts(panel.grid.minor = theme_line(colour = NA),
panel.grid.major = theme_line(colour = NA),
plot.background =  theme_rect(colour = NA, fill = NA),
axis.title.x = theme_text( face = "bold"),
axis.title.y = theme_text( face = "bold", angle = 90),
legend.key = theme_rect ( colour = NA, fill = NA))
 
 
dev.off
 
#colocar graficos en una página (creamos un nuevo grid y seleccionamos la forma)
library(gridExtra)
grid.newpage()
print(grid.arrange(plot1, plot2, plot3,plot4, plot5, nrow=1, ncol=5))   #ponerlos en columnas y filas

pm.stg<-function (response, factor1, factor2, error.bars = c("se",
"sd", "conf.int", "none"), level = 0.95, xlab =
deparse(substitute(factor1)),
    ylab = paste("mean of", deparse(substitute(response))), legend.lab
= deparse(substitute(factor2)),
    main = "Plot of Means", pch = 1:n.levs.2, lty = 1:n.levs.2,
    col = palette(), ylimits=range(response))
{
    if (!is.numeric(response))
        stop(gettextRcmdr("Argument response must be numeric."))
    xlab
    ylab
    legend.lab
    error.bars <- match.arg(error.bars)
    if (missing(factor2)) {
        if (!is.factor(factor1))
            stop(gettextRcmdr("Argument factor1 must be a factor."))
        valid <- complete.cases(factor1, response)
        factor1 <- factor1[valid]
        response <- response[valid]
        means <- tapply(response, factor1, mean)
        sds <- tapply(response, factor1, sd)
        ns <- tapply(response, factor1, length)
        if (error.bars == "se")
            sds <- sds/sqrt(ns)
        if (error.bars == "conf.int")
            sds <- qt((1 - level)/2, df = ns - 1, lower.tail = FALSE) *
                sds/sqrt(ns)
        sds[is.na(sds)] <- 0
        yrange <- if (error.bars != "none")
            c(min(means - sds, na.rm = TRUE), max(means + sds,
                na.rm = TRUE))
        else range(means, na.rm = TRUE)
        levs <- levels(factor1)
        n.levs <- length(levs)
        plot(c(1, n.levs), yrange, type = "n", xlab = xlab, ylab = ylab,
            axes = FALSE, main = main, ylim=ylimits)
        points(1:n.levs, means, type = "b", pch = 16, cex = 2)
        box()
        axis(2)
        axis(1, at = 1:n.levs, labels = levs)
        if (error.bars != "none")
            arrows(1:n.levs, means - sds, 1:n.levs, means + sds,
                angle = 90, lty = 2, code = 3, length = 0.125)
    }
    else {
        if (!(is.factor(factor1) | is.factor(factor2)))
            stop(gettextRcmdr("Arguments factor1 and factor2 must be factors."))
        valid <- complete.cases(factor1, factor2, response)
        factor1 <- factor1[valid]
        factor2 <- factor2[valid]
        response <- response[valid]
        means <- tapply(response, list(factor1, factor2), mean)
        sds <- tapply(response, list(factor1, factor2), sd)
        ns <- tapply(response, list(factor1, factor2), length)
        if (error.bars == "se")
            sds <- sds/sqrt(ns)
        if (error.bars == "conf.int")
            sds <- qt((1 - level)/2, df = ns - 1, lower.tail = FALSE) *
                sds/sqrt(ns)
        sds[is.na(sds)] <- 0
        yrange <- if (error.bars != "none")
            c(min(means - sds, na.rm = TRUE), max(means + sds,
                na.rm = TRUE))
        else range(means, na.rm = TRUE)
        levs.1 <- levels(factor1)
        levs.2 <- levels(factor2)
        n.levs.1 <- length(levs.1)
        n.levs.2 <- length(levs.2)
        if (length(pch) == 1)
            pch <- rep(pch, n.levs.2)
        if (length(col) == 1)
            col <- rep(col, n.levs.2)
        if (length(lty) == 1)
            lty <- rep(lty, n.levs.2)
        if (n.levs.2 > length(col))
            stop(sprintf(gettextRcmdr("Number of groups for factor2,
%d, exceeds number of distinct colours, %d."),
                n.levs.2, length(col)))
        plot(c(1, n.levs.1 * 1.4), yrange, type = "n", xlab = xlab,
            ylab = ylab, axes = FALSE, main = main, ylim=ylimits)
        box()
        axis(2)
        axis(1, at = 1:n.levs.1, labels = levs.1)
        for (i in 1:n.levs.2) {
            points(1:n.levs.1, means[, i], type = "b", pch = pch[i],
                cex = 2, col = col[i], lty = lty[i])
            if (error.bars != "none")
                arrows(1:n.levs.1, means[, i] - sds[, i], 1:n.levs.1,
                  means[, i] + sds[, i], angle = 90, code = 3,
                  col = col[i], lty = lty[i], length = 0.125)
        }
        x.posn <- n.levs.1 * 1.1
        y.posn <- sum(c(0.1, 0.9) * par("usr")[c(3, 4)])
        text(x.posn, y.posn, legend.lab, adj = c(0, -0.5))
        legend(x.posn, y.posn, levs.2, pch = pch, col = col,
            lty = lty)
    }
    invisible(NULL)
}
 
####################################
 
Esta función, pm.stg, permite ya definir los limites de y:
set.seed(1)
dv <- rnorm(100)
iv1 <- factor(sample(letters[1:2], 100, replace=T))
iv2 <- factor(sample(letters[14:17], 100, replace=T))
 
plotMeans(dv, iv1, iv2, error.bars="se") # los limites de Y los impone R
pm.stg(dv, iv1, iv2, error.bars="se", ylimits=c(-1,1)) # los limites de Y los imponemos nosotros

figure4 <- read.csv("l8_k_estimates.txt",header=T,  sep=" ")
 
pdf("figure4b.pdf", width=8, family="Times") 
op <- par(mfrow=c(1,3), mai=c(0.5,0,0.5,0.2))
 
# Se añaden los labels 
plot.new()
axis(2,at=seq(0.01,0.99, length.out=19), labels=figure4$name, tick= F, las=2, pos=c(1.1,1))
# Se añaden los estimadores de l8
dotchart(figure4$l8.value, pch=20, xlab="L8", color="grey20", xlim=c(min(figure4$l8.value)-max(figure4$l8.std),max(figure4$l8.value)+max(figure4$l8.std)))
segments(figure4$l8.value,1:19,figure4$l8.value+figure4$l8.std,1:19,col="grey20")
segments(figure4$l8.value,1:19,figure4$l8.value-figure4$l8.std,1:19,col="grey20")
# Se añaden los estimadores de k 
dotchart(figure4$k.value, pch=20, xlab="K", color="grey20", xlim=c(min(figure4$k.value)-max(figure4$k.std),max(figure4$k.value)+max(figure4$k.std)))
segments(figure4$k.value,1:19,figure4$k.value+figure4$k.std,1:19,col="grey20")
segments(figure4$k.value,1:19,figure4$k.value-figure4$k.std,1:19,col="grey20")
 
par(op)
dev.off()

library(nlme) 
fm1 <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject, data = Orthodont) 
 
plot(Orthodont)
newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female")) 
 
newdat$pred <- predict(fm1, newdat, level = 0)
 
Designmat <- model.matrix(eval(eval(fm1$call$fixed)[-2]), newdat[-3]) 
predvar <- diag(Designmat %*% fm1$varFix %*% t(Designmat)) 
newdat$SE <- sqrt(predvar) 
 
newdat$SE2 <- sqrt(predvar+fm1$sigma^2)
 
library(ggplot2) 
pd <- position_dodge(width=0.4) 
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+ 
   geom_point(position=pd)+ 
  geom_linerange(aes(ymin=pred-2*SE,ymax=pred+2*SE), position=pd)
 
 
## prediction intervals 
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+ 
   geom_point(position=pd)+ 
   geom_linerange(aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)

Galicia<- read.table("Galicia.txt")
names(Galicia) <- c("Longitud","Latitud")
# Usando plot() la costa sale desproporcionada
plot(Galicia$Long, Galicia$Lat, type="l", xlab="Longitud",ylab="Latitud")
title(main="Galicia")
# Usar mejor "eqscplot()"
eqscplot(Galicia$Long,Galicia$Lat,type="l",xlab="Longitud",ylab="Latitud",tol=-0.10)
#tol premite ajustar los margenes para encajar el mapa
title(main="Galicia")

# Colores en R
pdf("ColoursInR.pdf", width=11.7, height=8.5)
plot.new()
legend(x="topleft",colors(), pch=15, col=colors(), bty="n", ncol=10, cex=0.5, pt.cex=1, inset=0)
dev.off()

# Símbolos en R
pdf("SymbolsInR.pdf", width=5.5, height=4.5)
plot.new()
legend(x="topleft", as.character(0:25), pch=c(0:25), bty="p", ncol=4, cex=1.5, pt.cex=1, inset=0)
dev.off()

opar <-par(mar = c(20, 4, 4, 4)) #Con par guardamos las opciones gráficas que queremos
#en este caso  aumentamos los  margenes
plot(1:10)
lines(1:10)
par(xpd=TRUE)  #xpd para poder usar los margenes
legend(locator (1),lty=1, col="black", legend="straigh line")
#locator nos  permite colocar una o varias leyendas donde le  señalemos con el ratón
par(opar) #restauramos las opciones de graficos iniciales