File contents

# En quimica analítica se usa el modelo de regresión para calibrar 
# un método de medición. Para calibrar un fluorímetro se han 
# examinado 7 soluciones estándar de fluoresceína (de  las  que  se  
# conoce  la  concentración  medida con  mucha  precisión)  en  el  
# fluorímetro.  Los siguientes  datos  son  las concentraciones  
# y  la  intensidad  de  fluorescencia observada en el fluorímetro: 
# 
# Concentración, pg/ml:         0     2     4     6     8     10     12 
# Intensidad de fluorescencia: 2.1   5.0   9.0   12.6  17.3  21.0   24.7  

#######################################################################
# 1) Cargamos los datos a R
#######################################################################

conc<-scan()
0     2     4     6     8     10     12 

fluo<-scan()
2.1   5.0   9.0   12.6  17.3  21.0   24.7 

####################################################################
# 2) Grafico de dispersion
####################################################################

#scatter plot de los datos
plot(conc,fluo,xlab="concentracion",ylab = "fluorescencia",
     main="diagrama de dispersion 
          de fluorescencia vs. concentracion")
# ponemos la concentracion en las x, como variable independiente, 
# la fluorescencia como variable dependiente, en las y

####################################################################
# 3) Ajuste del modelo
####################################################################
#la instruccion lm (linear model) ajusta la recta de minimos cuadrados a los datos

salida<-lm(fluo~conc)
# el lm produce un objeto (una lista) que contiene muchas cosas
names(salida)

# para obtener los coeficientes, hay dos formas
salida$coefficients
coefficients(salida)

#superponemos la recta de minimos cuadrados al grafico
plot(conc,fluo,xlab="concentracion",ylab = "fluorescencia",
     main="diagrama de dispersion 
          de fluorescencia vs. concentracion
        con recta de minimos cuadrados")
abline(salida,col="red")

#esta tabla resume todo
summary(salida)

#tambien es un objeto de R
names(summary(salida))
######################################################################
# 4) estimamos el valor esperado de la fluorescencia cuando conc = 6 y 10
######################################################################

# a mano
coefficients(salida)[1]+coefficients(salida)[2]*6
coefficients(salida)[1]+coefficients(salida)[2]*10

# con la instruccion predict que calcula directamente el R.
predict(salida,newdata=data.frame(conc=6))
predict(salida,newdata=data.frame(conc=10))

# que calcula la instruccion fitted.values?
fitted.values(salida)

# que calcula la instruccion residuals?
residuals(salida)

####################################################################
# 5) intevalos de confianza para coeficientes
####################################################################
#el cuantil
qt(0.975,df=5) 

#los coeficientes estimados, y sus desvios estandares en las 
#primeras dos columnas
summary(salida)$coef

#los extremos superiores de los ic
summary(salida)$coef[,1]+ summary(salida)$coef[,2]*qt(0.975,df=5) 

#los extremos inferiores de los ic
summary(salida)$coef[,1]- summary(salida)$coef[,2]*qt(0.975,df=5) 

# el R calcula todo automaticamente
confint(salida)

###################################################
# 6) IC para E(Y) cuando X=6 y X=10
###################################################
# con el siguiente comando, que calcula el estimador e IC 
# correspondiente a X=6 y a X=10
predict(salida,newdata=data.frame(conc =  6), interval="confidence")
predict(salida,newdata=data.frame(conc = 10), interval="confidence")

###################################################
# 7) Intervalos de prediccion
###################################################

predict(salida,newdata=data.frame(conc =  6), interval="prediction")
predict(salida,newdata=data.frame(conc = 10), interval="prediction")

####################################################################
# 8) Validacion de supuestos
####################################################################
# Graficos para decidir si se satisfacen los supuestos de linealidad
# y homoscedasticidad 

# 1) el grafico que ya teniamos
plot(conc,fluo)
abline(salida,col=2)

mfrow=c(2,2)  #para hacer cuatro graficos juntos
# 2)
plot(conc,salida$residuals)
abline(a=0,b=0,col=4)

#3)
plot(salida$fitted.values,salida$residuals)
abline(a=0,b=0,col=4)

#4)
qqnorm(rstandard(salida),main="qqplot de residuos 
       estandarizados")


#5)
boxplot(rstandard(salida),main="boxplot de residuos 
        estandarizados")
mfrow=c(1,1)  #volver a poner un solo grafico
# tests
shapiro.test(rstandard(salida))