Estadística y Elementos de Quimiometría - 2021

Curvas de calibrado

Ejercicio 1

Enunciado

Los siguientes valores se obtuvieron para la concentración de \(SO_2\) (dióxido de azufre) en el aire próximo a una fábrica de papel

1.96,1.91,1.88,1.94,1.89,1.95 ppm

1. Calcular (a) la media,(b) la desviación estándar absoluta, (c) el coeficiente de variación,(d) el intervalo de confianza al 99 % de confianza y (e) expresar el resultado final de modo correcto.

Solución

SO2 <- c(1.96,1.91,1.88,1.94,1.89,1.95)
n <- length(SO2)
mean(SO2)

## [1] 1.921667

sd(SO2)

## [1] 0.03311596

SE <- sd(SO2)/sqrt(n) 
# Asumimos normalidad?
qqnorm(SO2)
qqline(SO2)

shapiro.test(SO2)

## 
##  Shapiro-Wilk normality test
## 
## data:  SO2
## W = 0.91747, p-value = 0.4873

alpha = 0.01
tcrit <- qt(alpha/2,n-1,lower.tail = FALSE)
(IC <- mean(SO2)+c(-1,1)*tcrit*SE)

## [1] 1.867154 1.976179

1. Supongamos que la legistación vigente establece que la concentración tiene que ser inferior a 1.7 ppm.

# Test unilateral 1
xbarra <- mean(SO2)
mu0 <- 1.7     #H0: mu = 1.7 vs H1: mu > 1.7 (el problema es si la concentración de SO2 es mayor al límite permitido)
Tobs <- abs((xbarra-mu0)/(SE))
p.valor <- pt(Tobs,n-1,lower.tail = FALSE)
t.test(SO2,mu=mu0,alternative = 'greater',conf.level = (1-alpha))

## 
##  One Sample t-test
## 
## data:  SO2
## t = 16.396, df = 5, p-value = 7.699e-06
## alternative hypothesis: true mean is greater than 1.7
## 99 percent confidence interval:
##  1.876174      Inf
## sample estimates:
## mean of x 
##  1.921667

Pues entonces sabemos, a un 99% de confianza, que esta fábrica no está cumpliendo con la legislación y por tanto se le puede sancionar.

2. Imaginemos ahora que el límite que nos a la legislación es 1.9 ppm. En tal caso no podemos ganatizar si la empresa está cumpliendo o no con el límite establecido.

# Test unilateral 2
xbarra <- mean(SO2)
mu0 <- 1.9     #H0: mu = 1.9 vs H1: mu > 1.9 (el problema es si la concentración de SO2 es mayor al límite permitido)
Tobs <- abs((xbarra-mu0)/(SE))
p.valor <- pt(Tobs,n-1,lower.tail = FALSE)
t.test(SO2,mu=mu0,alternative = 'greater',conf.level = (1-alpha))

## 
##  One Sample t-test
## 
## data:  SO2
## t = 1.6026, df = 5, p-value = 0.08496
## alternative hypothesis: true mean is greater than 1.9
## 99 percent confidence interval:
##  1.876174      Inf
## sample estimates:
## mean of x 
##  1.921667

3. Finalmente si el límite que establece la legislación fuese 2.1 ppm, entonces si puede decir que cumple con la legislación (No existe evidencia para decir que la concentración de \(SO_2\) es superior a 2.1 ppm.

# Test unilateral 3
xbarra <- mean(SO2)
mu0 <- 2.1     #H0: mu = 2.1 vs H1: mu > 2.1 (el problema es si la concentración de SO2 es mayor al límite permitido)
Tobs <- abs((xbarra-mu0)/(SE))
p.valor <- pt(Tobs,n-1,lower.tail = FALSE)
t.test(SO2,mu=mu0,alternative = 'greater',conf.level = (1-alpha))

## 
##  One Sample t-test
## 
## data:  SO2
## t = -13.191, df = 5, p-value = 1
## alternative hypothesis: true mean is greater than 2.1
## 99 percent confidence interval:
##  1.876174      Inf
## sample estimates:
## mean of x 
##  1.921667

Ejercicio 2

Enunciado

1. La Tabla 1 proporciona un ejemplo de datos de respuesta-concentración para su análisis, incluyendo respuestas medidas por triplicado. Grafique los datos de respuesta en función de la concentración y compruebe un forma visual que se desvían de la linealidad. Establezca un límite superior del rango lineal en forma cualitativa para luego compararlo con el calculado mediante una prueba estadística apropiada.

2. Con los datos restringidos al rango lineal, estimar la recta de mínimos cuadrados.

3. Estimar las incertidumbres asociadas a las estimaciones de la pendiente y la ordenada al orien de la recta de calibración ajustada e informe los los valores de \(\beta_1\) y \(\beta_0\) con el número correcto de cifras significativas.

4. La Tabla 2 muestra los valores de la señal para cuatro muestras incógnita, todos por triplicado. Estime la concentación de analito en las cuatro muestras problema, calcule sus desvíos estándar e informe el resultado con el numero apropiado de cifras significativas.

Solución

Comenzamos graficando los puntos para chequear la linealidad de manera exploratoria.

concentracion <- 0:10
resp1 <- c(0.06,1.44,2.82,4.15,5.29,6.61,7.79,8.89,10.03,10.84,11.87)
resp2 <- c(0.08, 1.56,2.76,4.20,5.46,6.54,7.70,8.97,9.88,10.91,11.81)
resp3 <- c(-0.06,1.41,2.90,4.08,5.52,6.69,7.69,8.83,9.77,10.65,11.90)
mydata <- data.frame(resp1,resp2,resp3,concentracion)
mydata = tidyr::gather(mydata,"id", "value", 1:3) 
ggplot(mydata,aes(x=concentracion, y=value)) +
  geom_point(size=3,color="blue3") +
  labs(x="Concentración", y="Señal") +
  scale_x_continuous(expand = c(0, 0),limits = c(-0.5, 10.5),breaks = seq(0,10,1),labels = seq(0,10,1)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-1, 12.5)) +
  geom_segment(aes(x = 0, y = -1, xend = 0, yend=12.5))+
  geom_segment(aes(x = -0.5, y = 0, xend = 10.5, yend=0))+
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")

Si ajustamos la recta y la graficamos sobre la nube de puntos…

ajuste <- lm(value~concentracion,data=mydata)
summary(ajuste)

## 
## Call:
## lm(formula = value ~ concentracion, data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50833 -0.20352  0.09718  0.20097  0.36270 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.44833    0.08562   5.236 1.09e-05 ***
## concentracion  1.17724    0.01447  81.340  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2629 on 31 degrees of freedom
## Multiple R-squared:  0.9953, Adjusted R-squared:  0.9952 
## F-statistic:  6616 on 1 and 31 DF,  p-value: < 2.2e-16

ggplot(mydata, aes(x=concentracion, y=value))+ 
  geom_abline(intercept=coefficients(ajuste)[1],slope=coefficients(ajuste)[2], color="#e67e1b")+
  geom_point(size=3,color="blue3") +
  labs(x="Concentración", y="Señal") +
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")

Si miramos los residuos, notamos también la falta de ajuste del modelo (la curvatura nos está indicando que pueder ser que agregar un término cuadrático podría mejorar el ajsute).

par(mfrow=c(1,2))
plot(ajuste,c(1,2))

De hecho, si ajustamos un modelo cuadrático, todo mejora…

ajuste2 <- lm(value~concentracion+I(concentracion^2),data=mydata)
summary(ajuste2)

## 
## Call:
## lm(formula = value ~ concentracion + I(concentracion^2), data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22902 -0.05854  0.00129  0.05595  0.13631 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.037086   0.036637   1.012     0.32    
## concentracion       1.451407   0.017046  85.149   <2e-16 ***
## I(concentracion^2) -0.027416   0.001642 -16.700   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08329 on 30 degrees of freedom
## Multiple R-squared:  0.9995, Adjusted R-squared:  0.9995 
## F-statistic: 3.31e+04 on 2 and 30 DF,  p-value: < 2.2e-16

plot(ajuste2,1)

Pero lo que estamos estudiando es calibración mediante un modelo lineal, entonces veamos el rango donde tenemos linealidad.

concentracion <- 0:10
resp1 <- c(0.06,1.44,2.82,4.15,5.29,6.61,7.79,8.89,10.03,10.84,11.87)
resp2 <- c(0.08, 1.56,2.76,4.20,5.46,6.54,7.70,8.97,9.88,10.91,11.81)
resp3 <- c(-0.06,1.41,2.90,4.08,5.52,6.69,7.69,8.83,9.77,10.65,11.90)
x <- rep(concentracion,3)
y = c(resp1,resp2,resp3)
N <- length(y)
k <- length(concentracion)
mod1 <- lm(y~x)
medias <- tapply(y,x,mean)
var <- tapply(y,x,var)
ns <- tapply(y,x,length) 
sy <- sum((ns-1)*var)/(N-k)  # esta es la pooled
syx <- sigma(mod1)

Fobs <- syx^2/sy

pf(Fobs,N-2,N-k,lower.tail = FALSE)

## [1] 7.814483e-07

# Rechazo linealidad. Quitamos la concentración más grande (10 en este ejemplo) y repetimos la prueba.
N <- N-3
k <- k-1
x <- rep(concentracion[1:k],3)
y = c(resp1[1:k],resp2[1:k],resp3[1:k])
mod1 <- lm(y~x)
medias <- tapply(y,x,mean)
var <- tapply(y,x,var)
ns <- tapply(y,x,length) 
sy <- sum((ns-1)*var)/(N-k)  # esta es la pooled
syx <- sigma(mod1)

Fobs <- syx^2/sy

pf(Fobs,N-2,N-k,lower.tail = FALSE)

## [1] 2.428533e-05

Repetimos el procedimiento hasta no rechazar la hipótesis de linealidad.

# Lo hacemos con un lazo while???
concentracion <- 0:10
resp1 <- c(0.06,1.44,2.82,4.15,5.29,6.61,7.79,8.89,10.03,10.84,11.87)
resp2 <- c(0.08, 1.56,2.76,4.20,5.46,6.54,7.70,8.97,9.88,10.91,11.81)
resp3 <- c(-0.06,1.41,2.90,4.08,5.52,6.69,7.69,8.83,9.77,10.65,11.90)
x <- rep(concentracion,3)
y = c(resp1,resp2,resp3)
N <- length(y)
k <- length(concentracion)

mod1 <- lm(y~x)
medias <- tapply(y,x,mean)
var <- tapply(y,x,var)
ns <- tapply(y,x,length) 
sy <- sum((ns-1)*var)/(N-k)  # esta es la pooled
syx <- sigma(mod1)

Fobs <- syx^2/sy
print(c(k,N))

## [1] 11 33

p.valor <- pf(Fobs,N-2,N-k,lower.tail = FALSE)
# vamos a guardar lo Fobs y los p.value
Fobss <- NULL
p.values <- NULL
while(p.valor<0.05){
  k <- k-1
  N <- N-3
  print(c(k,N))
  x <- rep(concentracion[1:k],3)
  y = c(resp1[1:k],resp2[1:k],resp3[1:k])
  mod1 <- lm(y~x)
  medias <- tapply(y,x,mean)
  var <- tapply(y,x,var)
  ns <- tapply(y,x,length) 
  sy <- sum((ns-1)*var)/(N-k)  # esta es la pooled
  syx <- sigma(mod1)
  
  Fobs <- syx^2/sy
  Fobss <- c(Fobss,Fobs)
  p.valor <-pf(Fobs,N-2,N-k,lower.tail = FALSE)
  p.values <- c(p.values,p.valor)
}

## [1] 10 30
## [1]  9 27
## [1]  8 24
## [1]  7 21
## [1]  6 18

k; N; p.values; Fobss  

## [1] 6

## [1] 18

## [1] 2.428533e-05 7.496356e-04 6.495406e-03 3.016222e-02 2.139387e-01

## [1] 6.692972 4.622001 3.502072 2.732476 1.579207

Por lo tanto debemos quedarnos con las primeras 6 concentraciones (y las primeras 6 señales de cada repetición)

Volvemos a graficar la recta y los puntos experimentales:

concentracion <- 0:10
resp1 <- c(0.06,1.44,2.82,4.15,5.29,6.61,7.79,8.89,10.03,10.84,11.87)
resp2 <- c(0.08, 1.56,2.76,4.20,5.46,6.54,7.70,8.97,9.88,10.91,11.81)
resp3 <- c(-0.06,1.41,2.90,4.08,5.52,6.69,7.69,8.83,9.77,10.65,11.90)
mydata <- data.frame(resp1,resp2,resp3,concentracion)
mydata <- mydata[1:k,]
mydata = tidyr::gather(mydata,"id", "value", 1:3) 


ajuste <- lm(value~concentracion,data=mydata)
ggplot(mydata, aes(x=concentracion, y=value))+ 
  geom_abline(intercept=coefficients(ajuste)[1],slope=coefficients(ajuste)[2], color="#e67e1b")+
  geom_point(size=3,color="blue3") +
  labs(x="Concentración", y="Señal") +
   scale_x_continuous(expand = c(0, 0),limits = c(-0.5, 5.5),breaks = seq(0,5,1),labels = seq(0,5,1)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-0.5, 7.5)) +
  geom_segment(aes(x = 0, y = -0.5, xend = 0, yend=7))+
  geom_segment(aes(x = -0.5, y = 0, xend = 5.5, yend=0))+
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")

Ejercicio 3

Enunciado

Los errores aleatorios que se cometen en las pasadas de una cierta balanza siguen una distribución normal de media 0 y desviación estándar 1 mg. Después de n=9 pesadas se han registrado los errores (en mg): -0.07, 0.3, 1.8, -0.1, 2.0, 2.3, 0.62, 0.12 y 1.4. Se pide contrastar la hipótesis de que la balanza funciona correctamente frente a la alternativa de que tiene un error sistemático positivo, con un nivel de significación de 0.1.

Solución

Debemos testear las hipótesis \[\begin{cases} H_0: \mu=0\\ H_1: \mu>0 \end{cases}\]

error <- c(-0.07,0.3,1.8,-0.1,2.0,2.3,0.62,0.12,1.4) 
qqnorm(error)
qqline(error)

shapiro.test(error)

## 
##  Shapiro-Wilk normality test
## 
## data:  error
## W = 0.88194, p-value = 0.1646

alpha = 0.1
n <- length(error)
sigma <- 1
mu0 = 0
xbarra <- mean(error)
zcrit <- qnorm(1-alpha,0,1/sqrt(n))
zcrit<xbarra

## [1] TRUE

p.valor <- pnorm(xbarra,0,1/sqrt(n),lower.tail = FALSE)

Ejercicio 4

Enunciado

Diez mediciones de recuperación de bromuro potásico por cromatografía de gas-líquido en muestras de tomates de cierta partida arrojaron una media muestral de \(782.0\) \(\mu\)g/g y un desvío \(s = 16.2\) \(\mu\)g/g. Suponga que las mediciones tienen distribución normal.

1. Halle un intervalo de confianza del 95% para la media (\(\mu\)) de las mediciones en esta partida de tomates.

2. Suponiendo que el error de medición debido al método es despreciable respecto a la variabilidad entre los tomates y que las 10 mediciones se realizan para cada uno de 10 tomates, ¿cómo debe interpretarse la media \(\mu\)? ¿Cómo debe interpretarse \(\sigma^2\)?

3. En cambio, si las 10 mediciones se realizaron sobre el mismo tomate, ¿cómo puede interpretarse \(\mu\)? ¿Cómo debe interpretarse \(\sigma^2\)?

4. Halle un intervalo de confianza del 95% para la varianza \(\sigma^2\) de las mediciones.

5. Halle un intervalo de confianza del 95% para la desviación estándar (\(\sigma\)) de las mediciones.

Solución

n <- 10
s <- 16.2
media <- 782
alpha <- 0.05
chi1 <- qchisq(alpha/2,n-1)
chi2 <- qchisq(1-alpha/2,n-1)
IC <- c((n-1)*s^2/chi2,(n-1)*s^2/chi1)
ICs <- sqrt(IC)
ICs/media*100

## [1] 1.424927 3.781955

Ejercicio 5

Enunciado

En un proceso químico es necesario que una solución que se usa como reactivo tenga pH 8.21. Se dispone de un método para determinar el pH que produce mediciones normalmente distribuidas con media igual al verdadero valor del pH y desvío 0.02. Diseñar un test de hipótesis de manera que: si el pH es 8.21, se obtenga esa conclusión con probabilidad 0.95, y si el pH difiere de 8.21 en 0.03 (en cualquiera de las direcciones), se descarta la solución con probabilidad no inferior a 0.95.

1. ¿Qué se puede concluir si la media de las mediciones observadas fuese 8.32?

2. Si el pH es 8.33, hallar la probabilidad de concluir que el pH no es 8.21.

Solución

Ejercicio 7

Enunciado

Sospechamos que nuestro cromatógrafo está estropeado, y queremos determinar si los resultados que nos proporciona son lo suficientemente precisos. Para ello, realizamos una serie de 8 mediciones del contenido de una solución de referencia que, sabemos, contiene 90% de un determinado compuesto. Los resultados que obtenemos son: 93.3, 86.8, 90.4, 90.1, 94.9, 91.6, 92.3, 96.5 Construir un intervalo de confianza al nivel de 95% para la varianza poblacional. ¿Que conclusiones podemos realizar?

Solución

med <- c(93.3, 86.8, 90.4, 90.1, 94.9, 91.6, 92.3, 96.5)
mean(med)

## [1] 91.9875

mu <- 90
n <- length(med)
var <- var(med)
alpha <- 0.05
crit1 <- qchisq(alpha/2,n-1)
crit2 <- qchisq(1-alpha/2,n-1)
IC <- c((n-1)*var/crit2,(n-1)*var/crit1)

# Pero como conozco mu en la construcción del intervalo
# lo más adecuado sería usar el estadístico
X <- sum((med-mu)^2)

crit1 <- qchisq(alpha/2,n)
crit2 <- qchisq(1-alpha/2,n)
c(X/crit2,X/crit1)

## [1]  5.441259 43.771461

sqrt(c(X/crit2,X/crit1))/mu*100

## [1] 2.591834 7.351112

Si consideramos que nosotros conocemos el valor de \(\mu\), podemos usar el estadístico \(\frac{\sum_{i=1}^n(x_i-\mu)^2}{\sigma^2}\) que tiene distribución \(\chi^2_n\) (siempre que \(x_1,\ldots, x_n\) sea una muestra aleatoria de una distribución normal con media \(\mu\) y varianza \(\sigma^2\)) (ver pág. 222 y 244 del Walpole)

Ejercicio 8

Enunciado

Para determinar la constante de un resorte, se puede medir la oscilación de una masa \(m\) fijada en uno de sus extremos mediante la fórmula \[\begin{align} k &= \frac{4\pi^2 m}{T^2} \end{align}\] donde \(k\) es la constante del resorte y \(T\) es el período de oscilación de dicha masa. Se realizan mediciones de la masa y del período. Asumamos que cada medición de la masa y del período siguen distribuciones normales, con varianza \(0.1\) y \(0.05\), respectivamente. De 5 mediciones de la masa y 4 del período se obtuvieron medias muestrales \(0.52\;\) y \(1.3\). Asumiendo que todas las mediciones son independientes, estime la constante del resorte e indique el valor del error de estimación asociado a la estimación.

Solución

Hallemos primero la estimación puntial de la contante \(k\) del resorte:

m <- 0.52
T <- 1.3
Var_m <- 0.1
Var_T <- 0.05
k <- 4*pi^2*m/(T^2)

Ejercicio 9

Enunciado

Se quiere determinar el volumen, en litros, de un tanque cilíndrico. Para ello se realizan cuatro mediciones del diámetro interno del tanque, obteniéndose los siguientes resultados: 5.2, 5.7, 5.3 y 5.5 m. Además se realizaron tres mediciones de la altura del tanque, y los resultados fueron: 7.9, 7.8 y 7.6 m. Calcule el volumen del tanque indicando su correspondiente incertidumbre.

Solución

Sabemos que para calcular el volumen del tanque tenemos que usar la fórmula \[\begin{align} V &= \hbox{área de la base}\times\hbox{altura} = \pi r^2 h \end{align}\]

Como el diámetro y la altura vienen dados en metros,

diam <- c(5.2, 5.7, 5.3, 5.5)
mean(diam)

## [1] 5.425

# La media de las alturas
h <-c(7.9, 7.8, 7.6)
mean(h)

## [1] 7.766667

(V = pi*(mean(diam)/2)^2*mean(h))

## [1] 179.5246

# Calculamos ahora la incertidumbre de esta estimación.
var_d <- var(diam)/4
var_h <- var(h)/3

var_V <- pi^2*mean(diam/2)^2*mean(h)^2*var_d+pi^2*mean(diam/2)^4*var_h
sqrt(var_V)

## [1] 7.615601

Por lo tanto tenemos que \(V = 180\;(\pm \; 8)\)

Ejercicio 10

Enunciado

En una determinación volumétrica de un analito A, los datos obtenidos y sus incertidumbres asociadas fueron como sigue:

• Lectura inicial de la bureta 0.19 mL, 0.02 mL

• Lectura final de la bureta 9.26 mL, 0.03 mL

• Masa de la muestra 45.0 mg, 0.2 mg

A partir de los datos, encuentre el coeficiente de variación del resultado final para el %A que se obtiene al utilizar la siguiente ecuación y considerando que no hay incertidumbre en la masa equivalente.

\[\begin{align} \%A &= \frac{\hbox{volumen del titulante}\times \hbox{masa equivalente}}{\hbox{masa de la muestra}}\times 100\% \end{align}\]

Solución

Dado que el volumen del valorante (o titulante?) es igual a la lectura final de la bureta menos la lectura inicial de la bureta, podemos introducir los valores dados en la ecuación para \(\%A\).

\[\begin{align} \%A &= \frac{(9.26-0.19)\times \hbox{masa equivalente}}{45.0}\times 100\% \end{align}\]

Entonces el valor y la incerteza asociada al volumen del titulante son: \[\begin{align} V &= 9.26-0.19\;\; \hbox{mL}\\ s_V & = \sqrt{(0.03)^2+(0.02)^2} \;\; \hbox{mL} \end{align}\]

li <- 0.19
lf <- 9.26
m <- 45.0
sli <- 0.02
slf <- 0.03
sm <- 0.2
V <- lf-li
sV <- sqrt(slf^2+sli^2)

Ahora podemos obtener el error realtivo del resultado final para el %A:

sArel <- sqrt((sV/V)^2+(sm/m)^2)

\[\begin{align} \frac{s_{\%A}}{\%A} = \sqrt{\Big(\frac{s_V}{V}\Big)^2+\Big(\frac{s_m}{m}\Big)^2} =0.00596... \end{align}\]

Por lo tanto, el \(CV\) será

\[\begin{align} CV(\%A) = \frac{s_{\%A}}{\%A}\times 100\%= 0.6\% \end{align}\]

Ejercicio 11

Enunciado

La masa de una probeta vacía de 10 mL es de 25.442(2) g. Después de agregar 8.5(1) mL de líquido, la masa aumenta a 32.402(2) g. Calcule la densidad del líquido en g/mL y kg/L.

Solución

Por dato del enunciado, se pesa la probeta vacía y seca (\(w_{i}\)), enseguida se llena con \(V = 8.5\) mL del líquido problema y luego se pesa todo el conjunto (\(w_f\)). La diferencia \(w_f - w_{i}\) corresponde a la masa del líquido. Por lo tanto: \[\begin{align} m_l =32.402-25.442 \; \hbox{g} \end{align}\]

wi <- 25.442
wf <- 32.402
Vl = 8.5
ml = (wf-wi)
rho <- ml/Vl

Dado que \(\rho =\frac{w_f-w_i}{V}\). Luego, resulta

\[\begin{align} \rho &= \frac{m_l}{V_l}=\frac{6.96\;\hbox{g}}{8.5\;\hbox{mL}}= 0.8188235 \;\frac{\hbox{g}}{\hbox{mL}} \end{align}\]

Determinemos ahora la incerteza de esta estimación:

\[\begin{align} s_{m_l} & = \sqrt{(0.002)^2+(0.002)^2}\;\hbox{g}=0.002828427 \; \hbox{g}\\ s_{\rho} & = \hat{\rho}\sqrt{\Big(\frac{s_{m_l}}{m_l}\Big)^2+\Big(\frac{s_V}{V}\Big)^2}\;\frac{\hbox{g}}{\hbox{mL}} \end{align}\]

sp <- 0.002
sml <- sqrt(2*sp^2)
sVl <- 0.1
srho_rel <- sqrt((sm/ml)^2+(sVl/Vl)^2)

Ejercicio 12

Enunciado

Para determinar el volumen de un cilindro se utiliza la fórmula \(V = \pi r^2h\;\) donde \(V\) , \(r\) y \(h\) son el volumen, radio de la base y altura del cilindro respectivamnete. Se asume que las mediciones son veraces, con distribución normal y varianza 0.1 y 0.5 para el radio y la altura, respectivamente. De 5 mediciones del radio y 3 de la altura, las medias muestrales obtenidas resultaron ser 1 y 3, respectivamente. Asumamos también que todas las mediciones son independientes.

1. Estimar el volumen del cilindro y dar la incertidumbre asociada a la estimación.

2. Si el metrólogo quisiera mejorar la precisión de su estimación pero sólo dispone de tiempo para volver a medir solo una vez el radio o la altura. ¿Cuál le recomendaría que repita? ¿Por qué?

Solución

r <- 1
nr <- 5
h <- 3
nh <- 3
sdr <- sqrt(0.1/nr)
sdh <- sqrt(0.5/nh)
V <- pi*r^2*h
# Calculo de la incertidumbre
var_V <- (2*pi*r*h)^2*sdr^2+(pi*r^2)^2*sdh^2
sqrt(var_V)

## [1] 2.958217

#  2.958217


###### Si aumento una medición del radio
r <- 1
nr <- 6
h <- 3
nh <- 3
sdr <- sqrt(0.1/nr)
sdh <- sqrt(0.5/nh)
V <- pi*r^2*h
# Calculo de la incertidumbre
var_V <- (2*pi*r*h)^2*sdr^2+(pi*r^2)^2*sdh^2
sqrt(var_V)

## [1] 2.750763

# 2.750763

###### Si aumento una medición de la altura
r <- 1
nr <- 5
h <- 3
nh <- 4
sdr <- sqrt(0.1/nr)
sdh <- sqrt(0.5/nh)
V <- pi*r^2*h
# Calculo de la incertidumbre
var_V <- (2*pi*r*h)^2*sdr^2+(pi*r^2)^2*sdh^2
sqrt(var_V)

## [1] 2.887874

# 2.887874

Ejercicio 13

Enunciado

Se han determinado cuantitativamente sulfatos en agua mediante una técnica instrumental adecuada y calibrando por el método de adiciones de estándar (adición de patrones). Para ello se preparó en un conjunto de matraces A, B, C, D, E y F la serie de disoluciones siguiente:

Todas las disoluciones obtenidas se homogeneizan y se aforan a 50 mL.

1. Representar los valores de absorbancia frente al aumento en la concentración de analito en en el volumen de muestra original. Comprobar que la linealidad de los datos se verifica para todo el intervalo estudiado.

2. Encontrar la recta de calibración a la que se ajustan los datos obtenidos.

3. Determinar la concentración de la muestra problema sabiendo que la concentración de la disolución patrón es de 100.0 ppm. Determinar también la incertidumbre asociada a la concentración estimada.

Solución

El método de calibración por adición de estándar requiere representar gráficamente la señal medida frente al aumento en la concentración de analito en el volumen de muestra original (de cada disolución de una serie de disoluciones de patrón de concentraciones conocidas).

En este caso la señal adecuada es es la absorbancia pues la transmitancia no tiene una relación lineal con la concentración. Por ello debemos transformar los valores de transmitancia en absorbancia.

Ley de Lambert-Beer.

\[\begin{align} A &= a b c \end{align}\] donde \(a\) es una constante de proporcionalidad llamada absortividad. Debido a que la absorbancia es una cantidad sin unidades, la absortividad debe tener unidades que cancelen las unidades de \(b\) y de \(c\). Si, por ejemplo, \(c\) tiene las unidades de g L\(^{-1}\) y \(b\) tiene las unidades de cm, la absortividad tiene las unidades de L g\(^{-1}\) cm\(^{-1}\).

Cuando expresamos la concentración en moles por litro y \(b\) en cm, la constante de proporcionalidad se llama absortividad molar y se le da el símbolo \(\varepsilon\). Por lo tanto, \[\begin{align} A &= \varepsilon b c \end{align}\] donde \(\varepsilon\) tiene unidades mol\(^{-1}\)cm\(^{-1}\).

1 Conversión de las T en tantos por 1:

Para calcular la absorbancia tendremos que transformar la transmitancia expresada en % en transmitancia en tanto por 1. Para esto, basta con dividir por 100:

2 Cálculo de las absorbancias:

Para obtener la absorbancia necesitamos calcular el logaritmo (en base 10) de la transmitancia (expresada en tanto por 1)

Disoluciones para calibrar por el método de adiciones de estándar.

Según dice el enunciado, tomamos 6 alícuotas de la muestra, cada un ade volumen \(V_x=25\;\) mL de volumen (una alícuotas es una parte representativa de la muestra, es decir, en la alícuota el analito tiene la misma concentración que en la muestra global, la cual llamaremos \(c_x\))

1. A la primera alícuota se le añade 1 mL de reactivo precipitante (para que el analito reaccione con él y de un producto que es el que vamos a determinar mediante la técnica instrumental adecuada). Después añadimos disolvente (agua destilada, por ejemplo) de modo que el volumen total, \(V_t\), sea 50 mL.

2. A la segunda alicuota le añadimos 1 mL del precipatante más 1 mL de patrón que contiene al analito en concentración \(c_s=100\;\) ppm, más agua destilada hasta que el volumen final sea 50 mL.

3. Así sucesivamente con las demás alícuotas, añadiendo cada vez más volumen del patrón.

Independientemente de lo que nos dice el enunciado, vamos a resolverlo de dos maneras diferentes.

A Método 1: usando los aumentos de las concentraciones del analito.

Representamos gráficamente las absorbancias medidas (eje Y) frente a los aumentos en la concentración de analito en el volumen de muestra original (\(V_x= 25\) mL) debido a los estándares añadidos, que como vimos antes, \(c'= \frac{c_s V_s}{V_x}\;\) (eje X) .

Figura 1: Gráfico de calibración lineal para el método de adiciones estándar.

Vs <- 0:5
Transmitancia <- c(59.0,47.0,37.5,30.0,24.0,20.0)
T1 <- Transmitancia/100
A <- -log(T1,10)
cs <- 100
Vx <- 25.00
aumento <- cs*Vs/Vx
ajuste <- lm(A~aumento)
summary(ajuste)

## 
## Call:
## lm(formula = A ~ aumento)
## 
## Residuals:
##          1          2          3          4          5          6 
## -0.0043652 -0.0005162  0.0026452  0.0046501  0.0066550 -0.0090689 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.2335132  0.0047800   48.85 1.05e-06 ***
## aumento     0.0237263  0.0003947   60.11 4.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006604 on 4 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9986 
## F-statistic:  3614 on 1 and 4 DF,  p-value: 4.586e-07

(c0 <- coefficients(ajuste)[1]/coefficients(ajuste)[2])

## (Intercept) 
##    9.841962

Al aplicar la regresión lineal al conjunto de datos \((A, c')\;\) se obtiene \(m = 0.0237\;(\pm 0.0004)\;\) , \(b = 0.234\;(\pm 0.005)\;\;\) y el cociente entre estas dos cantidades proporciona la concentración de sulfatos en la muestra de agua de \(c_x = 9.841962\;\) ppm.

B Método 2: usando los volúmenes de patrón.

Representamos gráficamente las absorbancias medidas (eje Y) frente a los volúmenes de patrón (eje X) añadidos a cada alícuota de 25 mL de muestra.

Figura 2: Gráfico de calibración lineal para el método de adiciones estándar.

Vs <- 0:5
Transmitancia <- c(59.0,47.0,37.5,30.0,24.0,20.0)
T1 <- Transmitancia/100
A <- -log(T1,10)
ajuste <- lm(A~Vs)
summary(ajuste)

## 
## Call:
## lm(formula = A ~ Vs)
## 
## Residuals:
##          1          2          3          4          5          6 
## -0.0043652 -0.0005162  0.0026452  0.0046501  0.0066550 -0.0090689 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.233513   0.004780   48.85 1.05e-06 ***
## Vs          0.094905   0.001579   60.11 4.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006604 on 4 degrees of freedom
## Multiple R-squared:  0.9989, Adjusted R-squared:  0.9986 
## F-statistic:  3614 on 1 and 4 DF,  p-value: 4.586e-07

# volumen cuando la absorbancia es 0
(Vs0 <- -coefficients(ajuste)[1]/coefficients(ajuste)[2])

## (Intercept) 
##   -2.460491

La ecuación de la recta de mínimos cuadrados resulta: \(\hat{A}=0.234+0.0949\;V_s\;\) .

Nos interesa determinar el valor de la abscisa (que llamaremos \(V_{s,0}\)) cuando la Absorbancia es 0. Reemplazando en la ecuaciónde la recta hallada, tenemos: \[\begin{align} 0 & = 0.234+0.0949\;V_{s,0}\\ V_{s,0} & = -\frac{0.234}{0.0949}=-2.460 \end{align}\]

\(V_{s,0}\) sería el volumen de patrón de 100 ppm en SO\(_4^{2-}\) que contendría igual masa de analito que la que existe en el primer matraz, la cual se debe solo al analito de la muestra. “Quitando” esa masa de analito, la señal caería a 0.

¿Cómo obtenemos a partir de este valor la concentración de analito en la muestra?

• En cada matraz en el que se hace la medida, ca masa total de analito será la suma de las masas de analito procedente del patrón y de la muestra. Tenemos entonces \(m_t=m_s+m_x\)

• Las masas se pueden relacionar con las concentraciones y los volumenes a través de la fórmula: \(c=\frac{m}{V}\), que es equivalente a \(m= cV\). En particular la masa de analito debida al patrón, \(m_s\), que habrá en cada matraz será igual al producto de la concentración original del patrón, \(c_s\) (valor constante) por el volumen de él que se toma en cada paso, \(V_s\) (valor variable). \[\begin{align} c_tV_t & = c_sV_s+c_xV_x\\ c_t \;50 \hbox{ mL } & = 100\hbox{ ppm }V_s+c_x\;25\hbox{ mL } \end{align}\]

• En particular, en el punto en que la prolongación de la recta corta al eje X, la señal es 0 y por lo tanto, en dicho punto, la concentración total de analito tiene que ser nula: \(c_t = 0\). Por ende también será \(C_tV_t=0\).

Esto convierte a la ecuación anterior en: \(0 = c_sV_{s,0}+c_xV_x\).

Despejando, podemos calcular la concentración de analito en la muestra problema a partir de valores conocidos: * \(V_{s,0}\) es el volumen (negativo) que se obtiene cuando la recta corta al eje X. * \(V_x\) es el volumen constante de muestra que se añade en cada matraz (25 mL en este ejemplo). * \(c_s\) es la concnetración de analito del patrón.

\[\begin{align} c_x &= -\frac{c_sV_{s,0}}{V_x}\\ c_x &= -\frac{100\hbox{ ppm}(-2.460) \hbox{ mL}}{25 \hbox{ mL}}\\ \end{align}\]

De donde resulta que \(c_x = 9.842\) ppm.

(cx <- -100*Vs0/25)

## (Intercept) 
##    9.841962

Error en la estimación

Para expresar correctamente el resultado de esta estimación, necesitamos calcular el error asociado a la misma (expresado por el desvío estándar)

Recordemos la fórmula de cálculo del mismo: \[\begin{align} S_{c_x} = \frac{S_{y|x}}{m}\sqrt{\frac{1}{N}+\frac{\bar{y}^2}{m^2\sum\limits_{i}(x_i-\bar{x})^2}} \end{align}\]

Vs <- 0:5
Transmitancia <- c(59.0,47.0,37.5,30.0,24.0,20.0)
T1 <- Transmitancia/100
A <- -log(T1,10)
cs <- 100
Vx <- 25.00
aumento <- cs*Vs/Vx
ajuste <- lm(A~aumento)
m <- coefficients(ajuste)[2]
b <- coefficients(ajuste)[1]
cx <- b/m
Sxx <- sum((aumento-mean(aumento))^2)
N <- length(Vs)
Se <- sigma(ajuste)
(desvio_cx <- Se/m*sqrt(1/N+mean(A)^2/(m^2*Sxx)))

##   aumento 
## 0.3490909

Luego, tenemos que sulfatos en la muestra de agua resulta \(9.8\pm 0.3\) ppm.

¿Cómo cambia la fórmula del calculo odel error de la estimación cuando utilizamos el método 2?

\[\begin{align} S_{c_x} = S_{V_s}\frac{c_s}{V_x} \end{align}\]

donde \(S_{V_s} = \frac{S_{y|x}}{m}\sqrt{\frac{1}{N}+\frac{\bar{y}^2}{m^2\sum\limits_{i}(x_i-\bar{x})^2}}\)

Vs <- 0:5
Transmitancia <- c(59.0,47.0,37.5,30.0,24.0,20.0)
T1 <- Transmitancia/100
A <- -log(T1,10)
cs <- 100
Vx <- 25.00
ajuste <- lm(A~Vs)
m <- coefficients(ajuste)[2]
b <- coefficients(ajuste)[1]
cx <- b/m*cs/Vx
Sxx <- sum((Vs-mean(Vs))^2)
N <- length(Vs)
Se <- sigma(ajuste)
(desvio_cx <- Se/m*sqrt(1/N+mean(A)^2/(m^2*Sxx))*cs/Vx)

##        Vs 
## 0.3490909

De paso observemos que, mediante este método, la concentración \(c_x\) se puede calcular con la fórmula

\[\begin{align} c_x = \frac{b}{m}\frac{c_s}{V_x} \end{align}\]

Ejercicio 14

Enunciado

Los valores en la Tabla 3 corresponden a la comparación entre las predicciones efectuadas para la determinación de teofilina en sangre mediante un método espectrofotométrico, comparado con un método de inmunofluorescencia polarizada (FPIA). No se determinaron las muestras por triplicado debido a la cantidad insuficiente de muestra (sueros de pacientes pediátricos). Sin embargo se estima que los desvíos estándar promedio para cada método son: 0.4 µg mL\(^{-1}\) para el método FPIA y 0.9 µg mL\(^{-1}\) para el método espectrofotométrico. Llevar a cabo la comparación de métodos mediante la construcción de la región de confianza apropiada, suponiendo que los desvíos estándar anteriores son constantes para todos los datos.

Solución

Como queremos comparar dos métodos analíticos, lo más adecuado es utilizar BLS. Pero dado que el desvío estándar para FPIA es menor que para el método espectrofotométrico, podríamos emplear un análisis de tipo WLS, con los valores de desvío estándar igual a 0.9 para todos los datos de la tabla anterior. Esto último, sin embargo, es idéntico al uso de un método OLS (varianza constante). Por lo tanto, podemos en este caso particular realizar una regresión lineal ordinaria empleando como variable Y los valores provistos por el método espectrofotométrico y como variable X los provistos por el método FPIA. Veamos las tres opciones.

Ajuste utilizando BLS

M1 <- c(0.0, 6.5, 33.2, 9.7, 12.2, 14.8, 20.1, 15.6, 19.3, 16.8, 24.2, 28.6, 0.0, 3.9, 8.0, 11.2, 11.4, 14.7, 16.5, 16.6, 19.8, 19.5, 23.0)
M2 <- c(1.4, 5.3, 30.6, 12.7, 14.9, 17.7, 19.9, 18.5, 20.4, 22.6, 27.1, 29.8, 0.0, 1.6, 5.7, 14.2, 15.3, 17.5, 17.6, 19.4, 18.7, 18.9, 21.2)
desvio1 <- rep(0.4,length(M1))
desvio2 <- rep(0.9,length(M2))

datos <- data.frame(M1=M1, V1=desvio1^2, M2=M2, V2=desvio2^2)

library(BivRegBLS)
res <- BLS(data = datos, xcol = 1, ycol = 3,var.x = 0.4, var.y = 0.9,npoints = 1000)
res$Estimate.BLS

##              H0          Estimate          Std Error       Lower 95%CI
## Intercept     0 0.993176862598204   1.03415999734129 -1.15747658553627
## Slope         1  1.00739853055625 0.0605378472025166 0.881503185383891
## Joint     (0,1)                                                       
##                Upper 95%CI    pvalue
## Intercept 3.14383031073268 0.3477968
## Slope     1.13329387572862 0.9038921
## Joint                      0.1140301

Tenemos entonces que las estimaciones de la pendiente y la ordenada al origen mediante BLS son: \[\begin{align} \hat{\beta}_1 &=1.01 \; \pm \; 0.06\\ \hat{\beta}_0 &= 1.0 \; \pm \; 1.0 \end{align}\]

y la región de confianza es

myFont1 <- "Montserrat"
puntos <- res$Ellipse.BLS
df <- data.frame(puntos)

beta0hat <- as.numeric(res$Estimate.BLS$Estimate[1])
beta1hat <- as.numeric(res$Estimate.BLS$Estimate[2])
ggplot(df, aes(x=x, y=y))+ 
  geom_point(size=0.5,col='darkorange')+
  annotate("point", x = beta1hat, y = beta0hat, colour = "orange", size=2)+
  annotate("point", x = 1, y = 0, colour = "darkred", size=2)+
  annotate("text", x=beta1hat, y=(beta0hat+0.3), parse = TRUE, label='paste(\'( \', hat(beta)[1], \',\',hat(beta)[0],\')\')',size=4,col='#536d4a')+
  labs(title='Región de confianza del 95%',y="Ordenada al Origen", x="Pendiente") +
  scale_x_continuous(expand = c(0, 0),limits = c(0.8, 1.25)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-2, 5)) +
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")+
  theme(plot.title = element_text(family=myFont1, face='bold', colour='black', size=12,hjust=0.5))

Ajuste utilizando WLS y OLS

M1 <- c(0.0, 6.5, 33.2, 9.7, 12.2, 14.8, 20.1, 15.6, 19.3, 16.8, 24.2, 28.6, 0.0, 3.9, 8.0, 11.2, 11.4, 14.7, 16.5, 16.6, 19.8, 19.5, 23.0)
M2 <- c(1.4, 5.3, 30.6, 12.7, 14.9, 17.7, 19.9, 18.5, 20.4, 22.6, 27.1, 29.8, 0.0, 1.6, 5.7, 14.2, 15.3, 17.5, 17.6, 19.4, 18.7, 18.9, 21.2)
desvio <- rep(0.9,length(M2))

q <- length(M1)
w <- (q/desvio^2)/(sum(1/desvio^2)) # pesos tal que suman q
ajuste.OLS <- lm(M2~M1)
ajuste.WLS <- lm(M2~M1,weights = w)

summary(ajuste.OLS)

## 
## Call:
## lm(formula = M2 ~ M1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5862 -1.8788  0.0655  1.7985  4.7247 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.35000    1.03036    1.31    0.204    
## M1           0.98365    0.06032   16.31 2.12e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.351 on 21 degrees of freedom
## Multiple R-squared:  0.9268, Adjusted R-squared:  0.9233 
## F-statistic:   266 on 1 and 21 DF,  p-value: 2.122e-13

summary(ajuste.WLS)

## 
## Call:
## lm(formula = M2 ~ M1, weights = w)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5862 -1.8788  0.0655  1.7985  4.7247 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.35000    1.03036    1.31    0.204    
## M1           0.98365    0.06032   16.31 2.12e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.351 on 21 degrees of freedom
## Multiple R-squared:  0.9268, Adjusted R-squared:  0.9233 
## F-statistic:   266 on 1 and 21 DF,  p-value: 2.122e-13

confidenceEllipse(ajuste.WLS,fill = T,lwd = 0,which.coef = c("(Intercept)", "M1"),main = "Región de confianza del 95%",xlab=expression(paste("Ordenanda al origen ", beta[0])), ylab=expression(paste("Pendiente  ", beta[1])),col='#f0c289')
text(coefficients(ajuste.WLS)[1],y=(coefficients(ajuste.WLS)[2]+0.005),labels=expression(paste('( ', hat(beta)[0], ',',hat(beta)[1],')')),col='#536d4a')
points(0,1,pch=19,col='darkred')

myFont1 <- "Montserrat"
puntos <- confidenceEllipse(ajuste.OLS,fill = T,lwd = 0,which.coef = c("(Intercept)", "M1"),segments=1000,draw=FALSE)
df <- data.frame(puntos)
colnames(df) <- c('ordenada','pendiente')
beta0hat <- as.numeric(ajuste.OLS$coefficients[1])
beta1hat <- as.numeric(ajuste.OLS$coefficients[2])
ggplot(df, aes(x=pendiente, y=ordenada))+ 
  geom_point(size=0.5,col='darkorange')+
  annotate("point", x = beta1hat, y = beta0hat, colour = "orange", size=2)+
  annotate("point", x = 1, y = 0, colour = "darkred", size=2)+
  annotate("text", x=beta1hat, y=(beta0hat+0.3), parse = TRUE, label='paste(\'( \', hat(beta)[1], \',\',hat(beta)[0],\')\')',size=4,col='#536d4a')+
  labs(title='Región de confianza del 95%',y="Ordenada al Origen", x="Pendiente") +
  scale_x_continuous(expand = c(0, 0),limits = c(0.8, 1.25)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-2, 5)) +
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")+
  theme(plot.title = element_text(family=myFont1, face='bold', colour='black', size=12,hjust=0.5))

Tenemos entonces que las estimaciones de la pendiente y la ordenada al origen mediante WLS son: \[\begin{align} \hat{\beta}_1 &=0.98 \; \pm \; 0.06\\ \hat{\beta}_0 &= 1.4 \; \pm \; 1.0 \end{align}\]

y vía OLS, las estimaciones son exactamente las mismas.

En cualquiera de los casos, como el punto ideal \((1,0)\) pertenece a la región de confianza, no rechazamos la hipótesis nula de equivalencia entre los dos métodos analíticos, es decir, no encontramos envidencia en los datos experimentales de que los métodos conduzcan a resultados diferentes.

plot(M1,M2,pch=19,col='blue')
abline(ajuste.OLS,col='red')
abline(ajuste.WLS,col='orange')
abline(a=as.numeric(res$Estimate.BLS$Estimate[1]), b=as.numeric(res$Estimate.BLS$Estimate[2]),col='green')
abline(a=0,b=1,col='magenta',lty=2)

Ejercicio 15

Enunciado

En la determinación del antibiótico ciprofloxacina en orina se emplean tres métodos multivariandos diferentes. La Tabla 4 proporciona datos para estudiar la exactitud de cada método, frente a una serie de muestras de referencia, cuya concentración de analito es conocida. Obtenga la región de confianza correspondiente a cada método y comente los resultados obtenidos.

Observación: como no hay datos disponibles acerca de los desvíos estándar de las predicciones, deberá realizarse un análisis OLS.

Solución

Nominal <- c(190,87,23,13,38,150,26,58,125,65,90,160,48,75,0,0,0,0)
M1 <- c(173,80,26,6,19,142,33,67,146,63,89,158,41,64,10,5,3,11)
M2 <- c(214,86,29,14,28,145,16,60,126,67,92,172,52,68,11,8,7,7)
M3 <- c(208,107,46,28,50,160,47,80,146,75,120,174,61,92,26,21,30,27)

ajuste.1 <- lm(M1~Nominal)
ajuste.2 <- lm(M2~Nominal)
ajuste.3 <- lm(M3~Nominal)

myFont1 <- "Montserrat"
puntos1 <- confidenceEllipse(ajuste.1,fill = T,lwd = 0,which.coef = c("(Intercept)", "nominal"),segments=1000,draw=FALSE)
puntos2 <- confidenceEllipse(ajuste.2,fill = T,lwd = 0,which.coef = c("(Intercept)", "nominal"),segments=1000,draw=FALSE)
puntos3 <- confidenceEllipse(ajuste.3,fill = T,lwd = 0,which.coef = c("(Intercept)", "nominal"),segments=1000,draw=FALSE)

df1 <- data.frame(puntos1)
colnames(df1) <- c('ordenada','pendiente')
beta0hat_1 <- as.numeric(ajuste.1$coefficients[1])
beta1hat_1 <- as.numeric(ajuste.1$coefficients[2])
df2 <- data.frame(puntos2)
colnames(df2) <- c('ordenada2','pendiente2')
beta0hat_2 <- as.numeric(ajuste.2$coefficients[1])
beta1hat_2 <- as.numeric(ajuste.2$coefficients[2])

df3 <- data.frame(puntos3)
colnames(df3) <- c('ordenada3','pendiente3')
beta0hat_3 <- as.numeric(ajuste.3$coefficients[1])
beta1hat_3 <- as.numeric(ajuste.3$coefficients[2])

df1 <- cbind(df1,df2,df3)

ggplot(df1, aes(x=pendiente, y=ordenada))+ 
  geom_point(size=0.5,col='darkorange')+
  geom_point(aes(x = pendiente2, y=ordenada2),color="#7a49a5",size=0.5)+
  geom_point(aes(x = pendiente3, y=ordenada3),color='red',size=0.5)+
  annotate("point", x = 1, y = 0, colour = "darkred", size=2)+
  annotate("text", x = 0.87, y = 15, colour = "#449406", label='Método 1',size=3)+
  annotate("text", x = 1.1, y = 5, colour = "#449406", label='Método 2',size=3)+
  annotate("text", x = 1.0, y = 15, colour = "#449406", label='Método 3',size=3)+
  labs(title='Región de confianza del 95%',y="Ordenada al Origen", x="Pendiente") +
  scale_x_continuous(expand = c(0, 0),limits = c(0.8, 1.15)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-10, 30)) +
  theme_ipsum(
    base_size = 11.5,
    axis_title_size = 12,
    axis_title_just = "rt")+
  theme(plot.title = element_text(family=myFont1, face='bold', colour='black', size=12,hjust=0.5))