library(yags) #help(yags) ### first create data frame pr with yerr=response that does not converge "pr"<- structure(.Data = list(ROWNAMES = structure(.Data = c(1, 50, 61, 72, 83, 94, 105, 116, 127, 2, 13, 24, 35, 44, 45, 46, 47, 48, 49, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 36, 37, 38, 39, 40, 41, 42, 43), .Label = c("1", "10", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "11", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "12", "120", "121", "122", "123", "124", "125", "126", "127", "128", "129", "13", "130", "131", "132", "133", "134", "135", "136", "137", "14", "15", "16", "17", "18", "19", "2", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "3", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "4", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "5", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "6", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "7", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "8", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "9", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99"), class = c("AsIs", "factor")), YERR = c(1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1), YORIG = c(1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1), FEMALE = c(1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1), AGE = c(0.1, 0.6, 0.2, 1.1, 0.2, 0.3, 1.4, 0.1, 0.8, 0.1, 1.3, 1.2, 0.7, 0.3, 0.3, 0.1, 0.1, 0.3, 0.7, 0.3, 0.2, 0.7, 0.3, 0.6, 0.4, 0.1, 0.3, 0.3, 0.2, 0.6, 0.7, 0, 0.1, 0.2, 1.3, 0.2, 1, 0, 1, 0.2, 0, 0.2, 0.5, 0.1, 0.8, 0.2, 1.5, 0.8, 0.4, 0.2, 0.6, 0.6, 0.7, 1.3, 0.4, 0.2, 0.1, 0.2, 0.4, 0.8, 0.7, 0.7, 0.2, 0, 0.1, 0.2, 0.1, 0.8, 0.5, 0.2, 0.3, 0.4, 0.7, 1, 0.6, 0, 0.2, 0.2, 0.3, 1.2, 0.1, 0.7, 0.8, 0.4, 0.6, 0.6, 0.9, 0.8, 0.9, 0.1, 0.2, 0.9, 0, 0.7, 0.9, 0.4, 0.9, 0.4, 0.2, 0.1, 0.2, 0.3, 0.2, 0.7, 0.4, 0.6, 1, 0.1, 1, 0.5, 1.1, 1, 0.5, 0.5, 0.3, 1.3, 1.1, 1.1, 0.9, 1.2, 1, 0.3, 0, 0.9, 1.1, 0.5, 1.1, 0.2, 1, 0.5, 0.6, 1.1, 0.3, 0.9, 0.8, 0, 0), DAYACC = c(7, 1, 7, 0.286, 2, 1, 15, 9.286, 3, 14.857, 6, 6, 2, 0.714, 0.143, 1, 1.429, 0.571, 0.143, 0.571, 3, 2, 1, 0.571, 2, 6, 9, 1, 0.286, 1, 0.143, 0.571, 0.429, 2, 2, 3, 0.143, 3, 2, 1, 1, 6, 2, 3, 13, 6, 6, 1, 1, 3, 4, 2, 6, 9, 1.143, 10, 6, 0.143, 15, 2, 2, 0.286, 0.429, 0.286, 3, 3, 0.286, 5, 10, 3, 1, 3.571, 2, 2, 4.286, 3, 0.429, 0.429, 0.429, 0.143, 2, 0.286, 0.286, 5, 1.571, 9, 1, 0.714, 5.571, 4, 2, 0.143, 4, 0.286, 4, 6, 16.714, 11.143, 3, 1.714, 0.143, 1, 0.143, 0.143, 5, 0.286, 3, 0.286, 4, 4, 1, 2, 4, 6, 4, 3, 6, 0.143, 2.857, 3, 0.143, 0.143, 7, 0.714, 4, 3, 1, 8, 0.286, 0.429, 10, 7, 4, 0.143, 5, 0.143, 5), SEVERE = c(3, 1, 3, 2, 2, 2, 4, 1, 2, 2, 2, 2, 3, 2, 1, 2, 3, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 1, 2, 2, 2, 2, 2, 2, 2, 3, 2, 4, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 4, 2, 2, 3, 2, 2, 3, 3, 2, 4, 2, 1, 3, 1, 2, 2, 1, 3, 2, 2, 1, 1, 3, 2, 2, 1, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 3, 2, 2, 2, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 3, 3, 1, 2, 2, 2, 3, 2, 4, 3, 2, 2, 2, 2, 2, 1, 3, 3, 3, 2, 3, 1, 2, 2), TOILET = c(8, 3, 6, 6, 4, 4, 20, 10, 4, 15, 5, 3, 10, 6, 5, 5, 7, 7, 3, 4, 4, 4, 6, 8, 6, 5, 12, 5, 4, 10, 5, 5, 4, 4, 3, 5, 8, 12, 4, 4, 6, 10, 10, 20, 6, 5, 8, 10, 4, 15, 5, 8, 6, 10, 10, 10, 7, 5, 6, 5, 5, 2, 5, 6, 4, 6, 4, 6, 7, 5, 10, 8, 5, 10, 6, 10, 3, 4, 8, 4, 5, 10, 4, 8, 11, 5, 5, 5, 6, 6, 6, 5, 7, 4, 8, 8, 8, 5, 5, 5, 4, 6, 6, 12, 3, 7, 2, 8, 6, 5, 6, 3, 6, 6, 8, 6, 8, 8, 4, 7, 4, 4, 6, 4, 3, 6, 6, 4, 3, 6, 9, 7, 5, 3, 3, 10, 6), PRACTID = c(8, 8, 8, 24, 24, 27, 27, 27, 27, 27, 41, 41, 41, 41, 45, 45, 45, 45, 55, 55, 56, 56, 56, 60, 60, 60, 60, 60, 60, 60, 65, 65, 65, 89, 89, 89, 89, 102, 102, 102, 102, 102, 107, 107, 107, 108, 108, 108, 108, 111, 111, 113, 113, 118, 118, 118, 118, 124, 124, 124, 124, 125, 125, 125, 125, 125, 125, 125, 125, 127, 127, 127, 127, 130, 132, 132, 132, 137, 137, 137, 137, 146, 146, 146, 153, 156, 156, 156, 182, 182, 182, 182, 182, 182, 185, 185, 185, 185, 185, 195, 195, 195, 201, 201, 206, 206, 206, 206, 207, 207, 207, 208, 208, 211, 211, 211, 211, 216, 216, 216, 216, 220, 220, 220, 220, 220, 228, 228, 228, 232, 232, 232, 232, 235, 235, 235, 235)), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "120", "121", "122", "123", "124", "125", "126", "127", "128", "129", "130", "131", "132", "133", "134", "135", "136", "137"), class = c("data.sheet", "data.frame")) y<-pr$YERR female<-pr$FEMALE age<-pr$AGE dayacc<-pr$DAYACC severe<-pr$SEVERE toilet<-pr$TOILET pract.id<-pr$PRACTID ### VJ Carey's LZ.exchalp function LZ.exchalp.original<- function(R, ID, cor.met = NULL, scalefun, p) { # example 3, section 4 of LZ 1986 residlist <- split.preserveord(R, ID) A0 <- 0 den <- 0 n <- unlist(lapply(residlist, length)) for(i in 1:length(residlist)) { if(n[i] > 1) { # # A0 <- A0 + (sum(residlist[[i]] %*% t(residlist[[i]])) - t(residlist[[i]]) %*% residlist[[i]])/2 # den <- den + 0.5 * n[i] * (n[i] - 1) } } # P <- scalefun(R, ID, cor.met, p) # A0/(P * (den - p)) } ## Modified LZ.exchalp function modLZ.exchalp<- function(R, ID, cor.met = NULL, scalefun, p) { # example 3, section 4 of LZ 1986 # modified with a minimum (and maximum?) bound similar to SAS version 8 residlist <- split.preserveord(R, ID) A0 <- 0 den <- 0 n <- unlist(lapply(residlist, length)) for(i in 1:length(residlist)) { if(n[i] > 1) { # # A0 <- A0 + (sum(residlist[[i]] %*% t(residlist[[i]])) - t(residlist[[i]]) %*% residlist[[i]])/2 # den <- den + 0.5 * n[i] * (n[i] - 1) } } # P <- scalefun(R, ID, cor.met, p) # fuzz <- 1e-005 min(max(A0/(P * (den - p)), -1/(max(n) - 1) + fuzz), 1 - fuzz) } ## The following produces an error LZ.exchalp<-LZ.exchalp.original #yags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr="exchangeable") ## But this converges LZ.exchalp<-modLZ.exchalp yags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr="exchangeable")