############################################## # Beispiel 3.59, nach Chap T. Le (S. 357ff.) # ############################################## beschwerden <- 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); visiten <- 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); facharzt <- c(1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1); geschlecht <- c(0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1); verguetung <- c(263.03, 334.94, 206.42, 226.32, 288.91, 275.94, 295.71, 224.91, 249.32, 269.00, 225.61, 212.43, 211.05, 213.23, 257.30, 326.49, 290.53, 268.73, 231.61, 241.04, 238.65, 287.76, 280.52, 237.31, 218.70, 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.40, 289.56, 305.67, 252.35, 276.86, 308.84); stunden <- c(1287.25, 1588.00, 705.25, 1005.50, 1667.25, 1517.75, 967.00, 1609.25, 1747.75, 906.25, 1787.75, 1480.50, 1733.50, 1847.25, 1433.00, 1520.00, 1404.75, 1608.50, 1220.00, 1917.25, 1506.25, 589.00, 1552.75, 1518.00, 1793.75, 1548.00, 1446.00, 1858.25, 1486.25, 933.75, 1168.25, 877.25, 1387.25, 1312.00, 1551.50, 973.75, 1638.25, 1668.25, 1652.75, 1029.75, 1127.00, 1547.25, 1499.25, 1747.50); n <- length(stunden); einsen <- rep(1.0, n); # Designmatrix erstellen design <- cbind(einsen, facharzt, geschlecht, verguetung, stunden); # negative Log-Likelihoodfunktion (ohne invariante Teile in beta) minus_loglikelihood_poisson_reduziert <- function(beta0, beta1, beta2, beta3, beta4) { beta <- c(beta0, beta1, beta2, beta3, beta4); dim(beta) <- c(5, 1); eta <- design %*% beta; offset <- log(visiten); return(-sum(beschwerden*(eta + offset)-exp(eta+offset))); } library("bbmle"); res_hand <- mle2(minus_loglikelihood_poisson_reduziert, start=list(beta0=0, beta1=0, beta2=0, beta3=0, beta4=0), method="BFGS"); summary(res_hand); #Maximaler Wert der Log-Likelihoodfunktion mit den optimierten Parametern -minus_loglikelihood_poisson_reduziert(res_hand@coef[1], res_hand@coef[2], res_hand@coef[3], res_hand@coef[4], res_hand@coef[5]); #Maximaler Wert der Log-Likelihoodfunktion mit den Buch-Parametern -minus_loglikelihood_poisson_reduziert(-8.1338, 0.2090, -0.1954, 0.0016, 0.0007) #Maximaler Wert der Log-Likelihoodfunktion mit den Buch-Parametern, #Vorzeichen von beta_1 und beta_2 geaendert -minus_loglikelihood_poisson_reduziert(-8.1338, -0.2090, 0.1954, 0.0016, 0.0007) # Eingebaute Routine glm() res_R <- glm(beschwerden ~ facharzt + geschlecht + verguetung + stunden, family=poisson(link = "log"), offset=log(visiten)); summary(res_R); -minus_loglikelihood_poisson_reduziert(res_R$coefficients[1], res_R$coefficients[2], res_R$coefficients[3], res_R$coefficients[4], res_R$coefficients[5]); # negative Log-Likelihoodfunktion (mit invarianten Teilen in beta) minus_loglikelihood_poisson <- function(beta0, beta1, beta2, beta3, beta4) { beta <- c(beta0, beta1, beta2, beta3, beta4); dim(beta) <- c(5, 1); eta <- design %*% beta; offset <- log(visiten); return(-sum(beschwerden*(eta + offset)-exp(eta+offset) -log(factorial(beschwerden)))); } -minus_loglikelihood_poisson(res_R$coefficients[1], res_R$coefficients[2], res_R$coefficients[3], res_R$coefficients[4], res_R$coefficients[5]); model_loglikelihood <- logLik(res_R); model_loglikelihood; # Stimmt ueberein #################################### # Likelihood-Quotienten-Globaltest # #################################### #Null-Modell res0<-glm(beschwerden~1, family=poisson(link = "log"), offset=log(visiten)); summary(res0); null_loglikelihood <- logLik(res0); null_loglikelihood; Devianz <- 2*(model_loglikelihood - null_loglikelihood); Devianz; p_wert <- 1.0 - pchisq(Devianz, df=4); p_wert; #################################################### # Alternative Bestimmung der Devianz (wie im Buch) # #################################################### redu_loglike_full <- -minus_loglikelihood_poisson_reduziert(res_R$coefficients[1], res_R$coefficients[2], res_R$coefficients[3], res_R$coefficients[4], res_R$coefficients[5]); redu_loglike_full; redu_loglike_null <- -minus_loglikelihood_poisson_reduziert(res0$coefficients[1], 0, 0, 0, 0); redu_loglike_null; Devianz_neu <- 2*(redu_loglike_full - redu_loglike_null); Devianz_neu; Devianz_neu - Devianz; ###################################################### # Devianz-Bestimmung ohne Anpassung des Null-Modells # ###################################################### summary(res_R); Devianz_ohne <- res_R$null.deviance - res_R$deviance; Devianz_ohne; abs(Devianz_ohne - Devianz_neu); ############################################### # Berechnung der geschaetzten Kovarianzmatrix # ############################################### kovarianzmatrix <- vcov(res_R); kovarianzmatrix; beta_hat <- res_R$coefficients; dim(beta_hat) <- c(5, 1); eta <- design %*% beta_hat; var_y <- exp(eta + log(visiten)); Cov_Y <- diag(as.vector(var_y)); F <- t(design) %*% Cov_Y %*% design; kovarianzmatrix2 <- solve(F); kovarianzmatrix2; abs(max(kovarianzmatrix-kovarianzmatrix2)); ############################################### # Pruefen einzelner Regressionskoeffizienten # ############################################### se <- sqrt(diag(kovarianzmatrix)); se; z_scores <- beta_hat / se; z_scores; p_werte <- 2.0 * (1.0 - pnorm(abs(z_scores))); p_werte; summary(res_R); ######################################### # Saturiertes Modell, Bestimmtheitsmass # ######################################### likeli_mu <- function(mu) { return(prod(mu^beschwerden / factorial(beschwerden) * exp(-mu))) } saturierte_loglikelihood <- log(likeli_mu(beschwerden)); saturierte_loglikelihood; #BESSER: y_groesser_null <- beschwerden[beschwerden > 0] sum(y_groesser_null * log(y_groesser_null) - log(factorial(y_groesser_null)) - y_groesser_null); # Nachvollziehen des Outputs summary(res_R) 2.0 * (saturierte_loglikelihood - null_loglikelihood); 2.0 * (saturierte_loglikelihood - model_loglikelihood); #Bestimmtheitsmass D_null <- 2.0 * (saturierte_loglikelihood - null_loglikelihood); D_mu <- 2.0 * (saturierte_loglikelihood - model_loglikelihood); D_null; D_mu; R_square <- 1.0 - D_mu / D_null; R_square;