Modelos Estadísticos. Grado Biotecnología
Abstract
Guía para la resolución de ejercicios sobre lineales generalizados para respuesta Poisson. Se muestra el código necesario para resolver cada uno de los ejercicios planteados. Se recomienda ir respondiendo a las preguntas que se van planteando a lo largo de cada uno de los ejercicios.# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)
library(reshape2)
library(lmtest)
library(mgcv)
library(MASS)
library(modelr)
library(survival)
library(survminer)En el banco de datos siguiente aparecen los datos correspondientes a un pequeño experimento en el que 7 de un total de 16 ratones fueron aleatoriamente seleccionados para recibir un nuevo tratamiento médico. Los 9 ratones restantes fueron asignados a un grupo control en el que no se administró ningún tipo de tratamiento. El objetivo del tratamiento era prolongar el tiempo de supervivencia después de una operación quirúrgica. La tabla muestra los tiempos de supervivencia (en días) tras la operación para los 16 ratones. El tratamiento, ¿prolongaba la vida de los ratones tras la operación?
grupo <- c(rep("T",7),rep("C",9))
tiempo <- c(94,38,23,197,99,16,141,52,104,146,10,51,30,40,27,46)
caso1 <- data.frame(grupo,tiempo)Representación gráfica de la información. En este caso tenemos riesgo constante.
ggplot(caso1,aes(x = grupo, y = log(tiempo))) +
geom_boxplot() +
theme_bw()¿Cuáles son las principales consecuencias del modelo ajustado? ¿Como interpretamos el coeficiente asociado con si el sujeto ha sido o no tratado? ¿Cuál es la estimación del riesgo relativo para un sujeto tratado frente a uno que no ha sido tratado? Justifica tu respuesta.
ajuste <-glm(tiempo ~ grupo,family = poisson(), caso1)
# Resumen del modelo
summary(ajuste)# Estimación de riesgo entre controles y tratados
RR <- exp(-0.43495)
# Expresado en porcentaje
RRporcen <- 100*(1-RR)
RRporcen # Riesgo relativo entre controles y tratados¿Qué podemos concluir sobre el diagnóstico el modelo?
ajuste <-glm(tiempo ~ grupo,family = poisson(), caso1)
# Obtención de valores para el diagnóstico
fortify(ajuste)
# Gráficos
par(mfrow=c(2,2))
plot(ajuste)¿Qué podemos decir sobre la predicción del tiempo de supervivencia para cada grupo de sujetos?
newdata <- data.frame(grupo = c("C","T"))
prediccion <- data.frame(newdata,predict(ajuste, newdata,type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
prediccion <- prediccion %>% mutate(
lwr = fit - (1.96 * se.fit),
upr = fit + (1.96 * se.fit)
)
prediccionLos datos siguientes contienen el tiempo de supervivencia, que es el tiempo hasta la muerte (en semanas) desde el diagnóstico de leucemia, y el número inicial de células blancas en sangre (en escala \(log_{10}\)), para un grupo de 17 sujetos que fueron diagnosticados de leucemia. Cuáles son tus conclusiones sobre la utilidad de saber el número de células blancas en la sangre para predecir el tiempo de supervivencia.
lcelulas <- c(3.36,2.88,3.63,3.41,3.78,4.02,4.00,4.72,5.00,4.23,3.73,3.85,3.97,4.51,4.54,5.00,5.00)
tiempo <- c(65,156,100,134,16,108,121,5,65,4,39,143,56,26,22,1,1)
caso2 <- data.frame(lcelulas,tiempo)Representación gráfica de la información. En este caso tenemos riesgo constante.
ggplot(caso2,aes(x = lcelulas, y = log(tiempo))) +
geom_point() +
theme_bw()¿Cuáles son las principales consecuencias del modelo ajustado? ¿Cuál es la estimación del riesgo relativo en función del recuento del logaritmo de células? Justifica tu respuesta y representa la curva estimada.
ajuste <-glm(tiempo ~ lcelulas,family = poisson(), caso2)
# Resumen del modelo
summary(ajuste)# Riesgo relativo
RR <- exp(-0.98639)
RR
# Expresado en porcentaje
# Riesgo por cada aumneto de un unidiad de lcelulas
RRporcen <- 100*(1-RR)
RRporcen¿Qué podemos concluir sobre el diagnóstico el modelo?
# Obtención de valores para el diagnóstico
fortify(ajuste)
# Gráficos
par(mfrow=c(2,2))
plot(ajuste)¿Qué podemos decir sobre la predicción del tiempo de supervivencia en fución del recuento del logaritmo de células?
#newdata <- data.frame(lcelulas = seq(min(caso2$lcelulas),max(caso2$lcelulas),0.1))
#prediccion <- data.frame(newdata,predict(ajuste, newdata,type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
#prediccion <- prediccion %>% mutate(
# lwr = fit - (1.96 * se.fit),
# upr = fit + (1.96 * se.fit))
# Gráfico
# Gráfico de la predicción
#ggplot(prediccion, aes(x = lcelulas, y = fit)) +
# geom_line() +
# geom_point(data = caso2, aes(x = lcelulas, y = tiempo)) +
# geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
# labs(x = "log(células)", y = "Tiempo de superviv encia", title = "Bandas de predicción") +
# theme_bw()Los datos de este ejercicio provienen de un famoso estudio realizado por Sir Richard Doll y colegas. En 1951, a todos los médicos británicos se les envió un breve cuestionario sobre si fumaban. Desde entonces, la información sobre sus muertes han sido recogidas. El conjunto de datos muestra el número de muertes por enfermedad coronaria (deaths), si eran o no fumadores (smoking), el grupo de edad al fallecimiento (age), y el número total de personas-años de observación en el momento del análisis (person-years). Para aprovechar el carácter ordinal de la variable age se sugiere utilizar una codificación numérica que identifique con un 1 al grupo de menor edad y con un 5 al grupo de mayor edad (otra opción sería utilizar el punto medio del intervalo). Representa la tasa de muertes por 100000 habitantes-año para el grupo de edad identificando si eran o no fumadores y ajusta el modelo que consideres más oportuno atendiendo a ese comportamiento. Los datos aparecen en Breslow and Day (1987).
require(dobson)
caso3 <- doctors
caso3$persons <- caso3$`person-years`
# Construimos el vector numérico de edad
caso3$age_num <- c(1,2,3,4,5,1,2,3,4,5)
caso3$age_num2 <- c(40,50,60,70,80,40,50,60,70,80)Representamos la incidencia teniendo en cuenta que el riesgo no es contante ya que depende del número de personas años observadas. ¿Qué conclusiones obtenemos del gráfico realizado?
caso3$tasa <- caso3$deaths / caso3$persons
# Gráfico de tasas
ggplot(caso3, aes(x = age, y = log(tasa), group =smoking, color = smoking)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line") +
theme_bw()Utilizando la variable de edad con el código numérico de la marca de clase de cada grupo de edad ¿qué tipo de modelo podemos estimar? ¿Cómo interpretamos el modelo obtenido tras el proceso de selección de efectos? ¿Como varía el riesgo relativo entre fumar y no fumar en función de la edad?
ajuste <- glm(deaths ~ smoking * age_num2, family = poisson(), offset = log(persons), caso3)
# Resumen del modelo
summary(ajuste)
# Selección del mejor modelo
ajuste <- step(ajuste, test = "Chisq")
ajusteEvaluación del riesgo entre fumador y no fumador para las diferentes edades. ¿Cómo evaluamos la función de riesgo relativo estimado? ¿Cómo varaía el riesgo realtivo para un doctor de 40, 50, 60 0 70 años?
# grid de edad
edadsec <- seq(35,84,1)
# Coeficientes del modelo
coefs <- tidy(ajuste)$estimate
# predictor lineal
plfumador <- coefs[1] + coefs[2] + coefs[3]*edadsec + coefs[4]*edadsec
plnofumador <- coefs[1] + coefs[3]*edadsec
# Riesgo relativo
RR <- exp(plfumador - plnofumador)
# Gráfico
datos <-data.frame(edadsec,RR)
ggplot(datos,aes(edadsec,RR)) +
geom_line() +
theme_bw() +
scale_y_continuous(breaks = seq(0.9,3.5,0.1)) +
labs(x = "Edad", y ="Riesgo Relativo", title = "Fumador vs No fumador")¿Qué podemos concluir sobre el diagnóstico el modelo? ¿Consideras necesario modificar el modelo de partida? ¿Que modelo propones como alternativa?
# Obtención de valores para el diagnóstico
fortify(ajuste)
# Gráficos
par(mfrow=c(2,2))
plot(ajuste)¿Qué conclusiones extraemos del nuevo ajuste? ¿Como varía el riesgo relativo en el nuevo modelo?
ajuste <- glm(deaths ~ smoking * (age_num2 + I(age_num2^2) + I(age_num2^3) + I(age_num2^4)), family = poisson(), offset = log(persons), caso3)
# Resumen del modelo
summary(ajuste)
# Selección del mejor modelo
ajuste <- step(ajuste, test = "Chisq")
ajuste¿Qué podemos concluir sobre el diagnóstico el modelo?
# Obtención de valores para el diagnóstico
fortify(ajuste)
# Gráficos
par(mfrow=c(2,2))
plot(ajuste)¿Consideras que el modelo obtenido tiene buena capacidad predictiva?
# Predicción corregineod por offset para obtener la tasa predicha
datos <- caso3 %>% mutate(fit = predict(ajuste,type = "response", se.fit=FALSE)/persons)
#prediccion
# Gráfico comparativo tasa real con tasa predicha
ggplot(datos,aes(x = tasa,y = fit)) +
geom_point() +
geom_abline(intercept = 0, slope = 1) +
facet_grid(. ~ smoking) +
labs(x = "Tasas observadas" , y = "Tasas predichas", title = "") +
theme_bw()Los datos de este ejercicio contienen el número de pólizas de seguros de coches (n) y el número de reclamaciones o partes recibidos (y). Se registra además:
district), que toma valor 1 para las contratadas en la ciudad de Londres y 0 para el resto de ciudades,car), registrada con un valor numérico donde 1 representa la categoría más baja y 4 la más alta,age), registrada con un valor numérico donde 1 indica los conductores de menor edad y 4 los de mayor edad.Estas dos últimas variables son categóricas ordinales que se codifican numéricamente para indicar el carácter ordinal de dicho factor. Analiza los datos teniendo en cuenta:
Los datos provienen de Baxter, Coutts, and Ross (1980)
require(dobson)
caso4 <- insurance
caso4$car_f <- as.factor(caso4$car)
caso4$age_f <- as.factor(caso4$age)
# tasa observada
caso4$tasa <- caso4$y / caso4$nRepresentamos la incidencia teniendo en cuenta que el riesgo no es contante ya que depende del número de pólizas. ¿Qué conclusiones obtenemos del gráfico realizado?
# Gráfico de tasas
ggplot(caso4, aes(x = car_f, y = log(tasa), group =age_f, color = age_f)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line") +
theme_bw()Utilizando las variables en formato factor ¿qué tipo de modelo podemos estimar? ¿Cómo interpretamos el modelo obtenido tras el proceso de selección de efectos? ¿Como varía el riesgo relativo entre el grupo de edad más bajo y más alto para la póliza contratada más barata? ¿y para la más cara?. Interpreta los resultados obtenidos ¿por qué coinciden ambos riesgos? ¿cuando no serían iguales?
ajuste <- glm(y ~ age_f*car_f, family = poisson(), offset = log(n), caso4)
# Resumen del modelo
summary(ajuste)
# Selección del mejor modelo
ajuste <- step(ajuste, test = "Chisq")
ajuste# Edad más baja y póliza barata
tasa.bb <- -1.7953
# Edad más alta y póliza barata
tasa.ab <- -1.7953 -0.5269
# Edad más baja y póliza cara
tasa.bc <- -1.7953 +0.5696
# Edad más alta y póliza cara
tasa.ac <- -1.7953-0.5269 +0.5696
## Riesgo edad - poliza barata
RR1 <-exp(tasa.ab - tasa.bb)
RR1
RR2 <-exp(tasa.ac - tasa.bc)
RR2¿Qué podemos concluir sobre el diagnóstico el modelo?
# Obtención de valores para el diagnóstico
fortify(ajuste)
# Gráficos
par(mfrow=c(2,2))
plot(ajuste)¿Consideras que el modelo obtenido tiene buena capacidad predictiva? ¿Qué características tienen las observaciones con predicciones erróneas?
# Predicción corregineod por offset para obtener la tasa predicha
datos <- caso4 %>% mutate(fit = predict(ajuste,type = "response", se.fit=FALSE)/n)
#prediccion
# Gráfico comparativo tasa real con tasa predicha
ggplot(datos,aes(x = tasa,y = fit)) +
geom_point() +
geom_abline(intercept = 0, slope = 1) +
labs(x = "Tasas observadas" , y = "Tasas predichas", title = "") +
theme_bw()En este ejemplo están recogidos los datos de un estudio sobre 44 médicos que trabajaban en el servicio de urgencias de un hospital Le (1998). El objetivo del estudio era determinar qué variables estaban relacionadas con el número de quejas que recibían los médicos a lo largo de un año. Para ello se controló. para cada médico el número de consultas que realizó a lo largo del año (consultas) y el número de quejas recibidas (quejas). Se tuvieron en cuanta además otras posibles características que podrían influir como:
ingresos)horas)sexo)Representa la tasa de quejas para cada posible variable predictora y obtén el modelo correspondiente a estos datos que te permita responde sobre la influencia de cada una de las características recogidas en el número de quejas recibidas.
consultas <- c(2014, 3091, 879, 1780, 3646, 2690, 1864, 2782, 3071,
1502, 2438, 2278, 2458, 2269, 2431, 3010, 2234, 2906,
2043, 3022, 2123, 1029, 3003, 2178, 2504, 2211, 2338,
3060, 2302, 1486, 1863, 1661, 2008, 2138, 2556, 1451,
3328, 2927, 2701, 2046, 2548, 2592, 2741, 3763)
quejas <- c(2, 3, 1, 1, 11, 1, 2, 6, 9, 3, 2, 2, 5, 2, 7,
2, 5, 4, 2, 7, 5, 1, 3, 2, 1, 1, 6, 2, 1, 1, 1,
0, 2, 2, 5, 3, 3, 8, 8, 1, 2, 1, 1, 10)
residente <- c("Si", "No", "Si", "No", "No", "No", "Si", "No", "No", "Si", "No", "No", "No",
"No", "No", "Si", "Si", "No", "Si", "No", "No", "Si", "Si", "No", "Si", "No", "Si", "Si", "No",
"Si", "Si", "No", "No", "No", "No", "Si", "Si", "No", "No", "Si", "Si", "No", "Si", "Si")
sexo <- c("M", "H", "H", "H", "H", "H", "H", "H", "M", "H", "M", "H", "H",
"M", "H", "H", "H", "H", "H", "H", "M", "M", "M", "H", "M", "M", "H", "H", "H",
"M", "H", "H", "H", "H", "H", "M", "H", "H", "H", "H", "H", "H", "M", "H")
ingresos <-c(263.03, 334.94, 206.42, 226.32, 288.91, 275.94, 295.71, 224.91,
249.32, 269, 225.61, 212.43, 211.05, 213.23, 257.3, 326.49, 290.53,
268.73, 231.61, 241.04, 238.65, 287.76, 280.52, 237.31, 218.7,
250.01, 251.54, 270.52, 247.31, 277.78, 259.68, 260.92, 240.22,
217.49, 250.31, 229.43, 313.48, 293.47, 275.4, 289.56, 305.67,
252.35, 276.86, 308.84)
horas <- c(1287.25, 1588, 705.25, 1005.5, 1667.25, 1517.75, 967, 1609.25,
1747.75, 906.25, 1787.75, 1480.5, 1733.5, 1847.25, 1433, 1520,
1404.75, 1608.5, 1220, 1917.25, 1506.25, 589, 1552.75, 1518,
1793.75, 1548, 1446, 1858.25, 1486.25, 933.75, 1168.25, 877.25,
1387.25, 1312, 1551.5, 973.75, 1638.25, 1668.25, 1652.75, 1029.75,
1127, 1547.25, 1499.25, 1747.5)
caso5 <- data.frame(consultas,quejas,residente,sexo,ingresos,horas)Disponemos de los pesos (g) de dos pikas afganas preñadas, a lo largo de 14 periodos de tiempo igualmente espaciados entre sí desde la concepción al parto. Ajusta una curva de crecimiento para describir estos datos ¿Puedes detectar alguna diferencia entre los dos animales?
animal <- c(rep("Pika1",14),rep("Pika2",14))
peso <- c(251,254,267,267,274,286,298,295,307,318,341,342,367,370,258,263,269,266,282,289,295,308,338,350,359,382,390,400)
semana <- c(1:14,1:14)
caso6 <- data.frame(animal,peso,semana)Este problema se refiere a datos de un estudio de cangrejos herradura anidados. Cada cangrejo de herradura hembra en el estudio tenía un cangrejo macho unido a ella en su nido. El estudio investigó los factores que afectan si el cangrejo hembra tenía otros machos, llamados satélites, que residían cerca de ella. La característica que se cree que afecta es el ancho del caparazón de la hembra. Se recogió información sobre las hembras y los machos, y estas se agruparon en grupos según su anchura. Cada grupo se caracterizo con la anchura media. Además se registro el número de hembras correspondientes a cada anchura y el número de satélites alrededor de ellas. representa la tasa media de satélites por hembra en función de la anchura y construye un modelo que represente la información contenida en estos datos. Extrae todas las conclusiones derivadas de el.
anchura <- c(22.69, 23.84, 24.77, 25.84, 26.79, 27.74, 28.67, 30.41)
hembras <- c(14,14,28,39,22,24,18,14)
satel <- c(14,20,67,105,63,93,71,72)
caso7 <- data.frame(anchura,hembras,satel)Baxter, L.A., S.M Coutts, and G.A.F. Ross. 1980. Applications of Linear Models in Motor Insurance. Proceedings of the 21st International Congress of Actuaries.
Breslow, N.E., and N. E. Day. 1987. Statistical Methods in Cancer Research, Volume 2: The Design and Analysis of Cohort Studies. International Agency for Research on Cancer.
Le, C.T. 1998. Applied Categorical Data Analyis. Wiley.
Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.