title: Nested outcome categories in a multinomial logistic

text: I’m modeling a target-shooting video game. A player can shoot and hit the target, shoot and miss, or switch weapons. For simplicity, the outcomes are `HIT`, `MISS`, `NEW_WEAPON`. I have covariates on player and weapon properties. The data are unbalanced: some players play a lot, some a little, some on one visit, some on multiple, and the targets are randomly generated with varying difficulty. **Some weapons are better for some targets, but switching has a penalty in the game, so players don’t like to switch. What I want to ask is whether switching is a good strategy on average, and whether it’s more effective for some players than others.** I’m thinking I’ll use a two-level logistic regression, but not a baseline- categories multinomial model.   1. A logistic regression of the outcomes `HIT`+`MISS` as 0 vs. `NEW_WEAPON` as 1. Including an intercept for each player will tell me a player’s tendency to fire weapons vs. switch weapons.   2. Then I can include this intercept as a covariate in a second logistic regression, this time of the outcomes `MISS` as 0 vs. `HIT` as 1. This covariate would have slope $\beta$, which tells me the effect of the tendency to fire a weapon vs. switch weapons on average.   3. As an extension to (2), I can build a logistic regression of `MISS` as 0 vs. `HIT` as 1 with slopes $\beta_p$ per player, which tells me whether switching is more effective for some players than others. **Does this sound reasonable? Any suggestions? I'm concerned that reusing the observations of`HIT` and `MISS` in (1) and in (2) is somehow not right.**


title: PCA or no PCA. Which results are true

text: I am trying to correlate hundreds of geoclimatic variables with genetic data. When I do individual tests, I find several significant correlations. Then I went ahead and performed PCA on the variables and found that first 5 PCs explain most variance. The correlations between the PCs and the genetic data reveals very few (as compared to when using individual variables) significant correlations. So I am wondering which results are true.


title: High heteroscedasticity level

text: My dependent variable - logincome - and one independent variable - age - are continuous. All other explanatory variables are categorical including BA_degree, race, occupation, region, homeownership. When I run OLS in stata, it reveals very high levels of heteroscedasticity. Do you have any suggestions how to fix this problem?


title: How to predict missing values in time series?

text: I have the following time series `1.3578511 0.5119648 1.3189847 0.9214787 1.2272616 4.9167998 1.2272616 1.2272616 0.8854192 2.3386331 1.6132899 0.2030302 0.8426226 1.2277843 NA 1.3189847 1.3578511 0.8530141 2.3386331 1.0541099 0.7747481 0.5764672 1.3189847 1.2160533 1.2272616 0.6715839 0.9651803 1.6132899 1.2006974 0.6875047 1.3245534 1.2006974 0.8221709 1.3101684 1.6132899 1.6132899 1.2006974 1.3189847 1.0018480 1.2277843 1.4424190 1.6132899 1.2277843 1.2006974 0.7779642 0.9381081 0.8854192 NA NA 1.3189847 1.1070461 0.8221709 4.9167998 0.9214787 1.3189847 1.3189847 1.2277843 1.4424190 1.6132899 1.6132899 4.9167998 0.8235792 0.9708839 1.1070461 1.2160533 0.8354292 1.4424190 1.1958634 0.5119648 1.4424190 1.4424190 1.6132899 1.6132899 0.6710844 1.2272616 0.9708839 0.8890464 1.4424190 0.8890464 0.8221709 1.1958634 0.8132233 0.4630722 4.9167998 0.8890464 1.3189847 0.7373181 1.1070461 1.2279813 0.8890464 0.3588158 1.4424190 0.8132233 0.4297043 1.3578511 4.9167998 1.2272616 0.8426226 1.4424190 1.6132899 NA ` in which `NA` are missing values, and I want to predict/forecast it. I searched over internet, but I haven't found that the Amelia package can impute missing values. I used it as follows:               library(Amelia)     t <- read.table("C:\\Users\\exam\\Desktop\\missing_ts.txt")     a.out <- amelia(t)      but I got the following error:               Amelia Error Code:  42      There is only 1 column of data. Cannot impute      Likewise,               amelia(x=as.matrix(1:101,t$V1))      resulted in               Amelia Error Code:  39      Your data has no missing values.  Make sure the code for      missing data is set to the code for R, which is NA.      Is there something wrong in the way I'm trying to forecast this time series? If yes, then what method should I use?


title: Multinomial probit for varying choice set

text: I apologize in advance for the "how do I run this model in R" question, but it turns out that's the shape my problem has assumed :-). Hopefully I have enough substantive questions surrounding it to be interesting, and the question will come out more like, "Does this command in R correspond to this statistical model?" I have collected data in a conjoint analysis (design matrix created using the Federov method from the `AlgDesign` package)--since this term seems to have a few meanings, I'll clarify. I asked individuals to choose among four offered options, which varied (orthogonally) on several plan dimensions `c("cm","pio","prev","price")` . The fourth option was always the base option and did not vary between choice sets. I also collected characteristics about each individual `c(chi, hexp, sex, age)`. Each individual was presented with a choice set, made their choice, then was presented with another choice set for another choice, and so on. They were instructed to ignore their previous choices when making the next choice. I can run a logit on this fairly trivially:               library(mlogit)     form <- choice~ cm*pio*prev+price | chi+hexp+sex+age     choices.logit <- mlogit(form, data=choices.mldata,reflevel=4)      However the logit as run misses a lot of the complexities of the data: the repeated samples from the same individual, the fact that individuals are presented with choice sets that do not always have the same options, etc. My model is similar to this one, but without the quantity discounting. `MNP:mnp` seems to be able to handle this case just fine, in that it notes: > If the choice set varies from one observation to another, use the syntax, > cbind(y1, y2, y3) ~ x1 + x2, in the case of a three choice problem, and > indicate unavailable alternatives by NA. If only the most preferred choice > is observed, y1, y2, and y3 are indicator variables that take on the value > one for individuals who prefer that choice and zero otherwise. The last > column of the response matrix, y3 in this particular example syntax, is used > as the base category. I could dummy out each of the plans (where plan is defined as a unique combination of `c("cm","pio","prev","price")`), include the plan characteristics as covariates, and run the probit. However, what makes a choice distinct in my survey is not the ID of the choice itself but the combination of characteristics which are offered. My sense is that what `MNP` would estimate is the probability of choosing e.g. plan 3, when what I care about is not plan 3 but the willingness to pay for pio, cm, and prev. Am I just misinterpreting the details of what `mnp` would do, and I should go ahead and run it as described? Or is my concern legitimate and I should look for another way to run this model? Thanks! Here's a sample of my data:               choices.mldata <- structure(list(PersonId = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 12L,      12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L,      22L, 22L, 22L, 22L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 31L, 31L, 31L, 31L, 31L, 31L,      31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 33L, 33L, 33L,      33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L,      36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L, 36L,      36L, 36L, 36L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L, 37L,      37L, 37L, 40L, 40L, 40L, 40L, 40L, 40L, 40L, 40L, 40L, 40L, 40L,      40L, 40L, 40L, 40L, 40L), option = structure(c(1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,      1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("1",      "2", "3", "4"), class = "factor"), cm = c(FALSE, TRUE, FALSE,      FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE,      TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE,      FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE,      FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,      TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,      FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE,      FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE,      TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,      FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,      FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,      TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE,      FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE,      FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,      TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE), pio = c(TRUE,      FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE,      FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE,      FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE,      TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,      FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE,      TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE,      FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE,      FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE,      FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE,      TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE,      FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE,      FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE,      FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), prev = c(TRUE, TRUE,      TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE,      FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,      TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE,      FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE,      TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,      FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE,      FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE,      FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE,      TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE,      TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE,      FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,      FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE,      FALSE, FALSE, TRUE, TRUE, TRUE, FALSE), price = c(2500, 2987,      2714, 2500, 2500, 2604, 2501, 2500, 2884, 2500, 2584, 2500, 2747,      2536, 2970, 2500, 2883, 2503, 2700, 2500, 2512, 2578, 2500, 2500,      2500, 2544, 2685, 2500, 2537, 2967, 2500, 2500, 2586, 2798, 2552,      2500, 2851, 2532, 2643, 2500, 2529, 2644, 2581, 2500, 2542, 2501,      3012, 2500, 2549, 2500, 2927, 2500, 2515, 2923, 2592, 2500, 2515,      2761, 2500, 2500, 2553, 2714, 2671, 2500, 2500, 2788, 2516, 2500,      2679, 2848, 2518, 2500, 2593, 2500, 2731, 2500, 2840, 2500, 2502,      2500, 2726, 2583, 2924, 2500, 2670, 2535, 2580, 2500, 2505, 2624,      2781, 2500, 2512, 2573, 2500, 2500, 2888, 2500, 2500, 2500, 2500,      2633, 2507, 2500, 2500, 2520, 2500, 2500, 2509, 2586, 2502, 2500,      2745, 2529, 2581, 2500, 2672, 2517, 2513, 2500, 2560, 3021, 2519,      2500, 2511, 2516, 2543, 2500, 2681, 2501, 3024, 2500, 2574, 2652,      2510, 2500, 2933, 2500, 2592, 2500, 2835, 2505, 2755, 2500),          choice = c(FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE,          FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,          FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,          TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE,          TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE,          FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,          TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,          FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,          FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE,          FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE,          FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,          FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,          TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE,          FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,          FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE,          FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE,          FALSE, FALSE), id = c("2_1", "2_1", "2_1", "2_1", "2_3",          "2_3", "2_3", "2_3", "12_1", "12_1", "12_1", "12_1", "12_2",          "12_2", "12_2", "12_2", "12_3", "12_3", "12_3", "12_3", "19_1",          "19_1", "19_1", "19_1", "19_2", "19_2", "19_2", "19_2", "19_3",          "19_3", "19_3", "19_3", "19_4", "19_4", "19_4", "19_4", "22_1",          "22_1", "22_1", "22_1", "22_2", "22_2", "22_2", "22_2", "22_3",          "22_3", "22_3", "22_3", "22_4", "22_4", "22_4", "22_4", "26_1",          "26_1", "26_1", "26_1", "26_2", "26_2", "26_2", "26_2", "26_3",          "26_3", "26_3", "26_3", "26_4", "26_4", "26_4", "26_4", "31_1",          "31_1", "31_1", "31_1", "31_2", "31_2", "31_2", "31_2", "31_3",          "31_3", "31_3", "31_3", "31_4", "31_4", "31_4", "31_4", "33_1",          "33_1", "33_1", "33_1", "33_2", "33_2", "33_2", "33_2", "33_3",          "33_3", "33_3", "33_3", "33_4", "33_4", "33_4", "33_4", "36_1",          "36_1", "36_1", "36_1", "36_2", "36_2", "36_2", "36_2", "36_3",          "36_3", "36_3", "36_3", "36_4", "36_4", "36_4", "36_4", "37_1",          "37_1", "37_1", "37_1", "37_3", "37_3", "37_3", "37_3", "37_4",          "37_4", "37_4", "37_4", "40_1", "40_1", "40_1", "40_1", "40_2",          "40_2", "40_2", "40_2", "40_3", "40_3", "40_3", "40_3", "40_4",          "40_4", "40_4", "40_4"), chi = structure(c(2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),          hexp = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 5L, 5L,          5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,          5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,          5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,          4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,          4L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("$1,000 - $2,499",          "$2,500 - $4,999", "$5,000 or More", "$500 - $999", "Less than $500"         ), class = "factor"), sex = structure(c(1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L), .Label = c("Female", "Male"), class = "factor"),          age = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L,          4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 4L, 4L,          4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,          4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,          3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L,          2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("18-29",          "30-39", "40-49", "50-64"), class = "factor")), .Names = c("PersonId",      "option", "cm", "pio", "prev", "price", "choice", "id", "chi",      "hexp", "sex", "age"), row.names = c("1.1", "1.2", "1.3", "1.4",      "2.1", "2.2", "2.3", "2.4", "3.1", "3.2", "3.3", "3.4", "4.1",      "4.2", "4.3", "4.4", "5.1", "5.2", "5.3", "5.4", "6.1", "6.2",      "6.3", "6.4", "7.1", "7.2", "7.3", "7.4", "8.1", "8.2", "8.3",      "8.4", "9.1", "9.2", "9.3", "9.4", "10.1", "10.2", "10.3", "10.4",      "11.1", "11.2", "11.3", "11.4", "12.1", "12.2", "12.3", "12.4",      "13.1", "13.2", "13.3", "13.4", "14.1", "14.2", "14.3", "14.4",      "15.1", "15.2", "15.3", "15.4", "16.1", "16.2", "16.3", "16.4",      "17.1", "17.2", "17.3", "17.4", "18.1", "18.2", "18.3", "18.4",      "19.1", "19.2", "19.3", "19.4", "20.1", "20.2", "20.3", "20.4",      "21.1", "21.2", "21.3", "21.4", "22.1", "22.2", "22.3", "22.4",      "23.1", "23.2", "23.3", "23.4", "24.1", "24.2", "24.3", "24.4",      "25.1", "25.2", "25.3", "25.4", "26.1", "26.2", "26.3", "26.4",      "27.1", "27.2", "27.3", "27.4", "28.1", "28.2", "28.3", "28.4",      "29.1", "29.2", "29.3", "29.4", "30.1", "30.2", "30.3", "30.4",      "31.1", "31.2", "31.3", "31.4", "32.1", "32.2", "32.3", "32.4",      "33.1", "33.2", "33.3", "33.4", "34.1", "34.2", "34.3", "34.4",      "35.1", "35.2", "35.3", "35.4", "36.1", "36.2", "36.3", "36.4"     ), class = c("mlogit.data", "data.frame"), index = structure(list(         chid = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L,          3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L,          7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L,          10L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L,          13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L,          16L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 19L, 19L, 19L,          19L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L, 22L, 22L, 22L,          22L, 23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 25L, 25L, 25L,          25L, 26L, 26L, 26L, 26L, 27L, 27L, 27L, 27L, 28L, 28L, 28L,          28L, 29L, 29L, 29L, 29L, 30L, 30L, 30L, 30L, 31L, 31L, 31L,          31L, 32L, 32L, 32L, 32L, 33L, 33L, 33L, 33L, 34L, 34L, 34L,          34L, 35L, 35L, 35L, 35L, 36L, 36L, 36L, 36L), .Label = 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"), class = "factor"), alt = structure(c(1L,          2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,          1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L,          4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L,          3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L,          2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,          1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L,          4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L,          3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L,          2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L,          1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("1", "2", "3",          "4"), class = "factor"), id = structure(c(1L, 1L, 1L, 1L,          1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,          2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,          3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,          4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,          5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,          6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,          7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,          8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,          9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,          10L, 10L, 10L, 10L, 10L, 10L, 10L), .Label = c("2", "12",          "19", "22", "26", "31", "33", "36", "37", "40"), class = "factor")), .Names = c("chid",      "alt", "id"), row.names = c("1.1", "1.2", "1.3", "1.4", "2.1",      "2.2", "2.3", "2.4", "3.1", "3.2", "3.3", "3.4", "4.1", "4.2",      "4.3", "4.4", "5.1", "5.2", "5.3", "5.4", "6.1", "6.2", "6.3",      "6.4", "7.1", "7.2", "7.3", "7.4", "8.1", "8.2", "8.3", "8.4",      "9.1", "9.2", "9.3", "9.4", "10.1", "10.2", "10.3", "10.4", "11.1",      "11.2", "11.3", "11.4", "12.1", "12.2", "12.3", "12.4", "13.1",      "13.2", "13.3", "13.4", "14.1", "14.2", "14.3", "14.4", "15.1",      "15.2", "15.3", "15.4", "16.1", "16.2", "16.3", "16.4", "17.1",      "17.2", "17.3", "17.4", "18.1", "18.2", "18.3", "18.4", "19.1",      "19.2", "19.3", "19.4", "20.1", "20.2", "20.3", "20.4", "21.1",      "21.2", "21.3", "21.4", "22.1", "22.2", "22.3", "22.4", "23.1",      "23.2", "23.3", "23.4", "24.1", "24.2", "24.3", "24.4", "25.1",      "25.2", "25.3", "25.4", "26.1", "26.2", "26.3", "26.4", "27.1",      "27.2", "27.3", "27.4", "28.1", "28.2", "28.3", "28.4", "29.1",      "29.2", "29.3", "29.4", "30.1", "30.2", "30.3", "30.4", "31.1",      "31.2", "31.3", "31.4", "32.1", "32.2", "32.3", "32.4", "33.1",      "33.2", "33.3", "33.4", "34.1", "34.2", "34.3", "34.4", "35.1",      "35.2", "35.3", "35.4", "36.1", "36.2", "36.3", "36.4"), class = "data.frame"), choice = "choice")


title: How to learn similarity of typed/attributed graphs?

text: I have a question for graph machine learning gurus :). For this project I'm working on, I need to be able to learn similarity between typed graphs. By typed I mean that every vertex and every edge of the graph has it's type. Also, both nodes and vertices can be attributed (assigned with some "key- value" attributes). Now, I could employ "standard" feature-based supervised ML to perform this task, but then I would have to come up with good features that represent similarity between two typed graphs. In this case I would operate with the computed features and not at the level of graphs themselves directly. I've read about graph kernels (and their corresponding kernel machines learning algorithms) that sort of deal precisely with this issue - they handle graphs directly, instead of inventing numeric features that reflect similarity. What I also figured out is that these graph kernel methods are used to classify graphs. I couldn't fugre out how to use them to learn a good "similarity function" between a pair of graphs? Are kernel methods even the right way to go for this type of problem? I'd really appreciate any advice. Are there perhaps any good/well-known "unsupervised" algorithm for measuring the similarity between two typed graphs? The similarity measure should consider both structure and semantics of the graphs (by semantics I mean types and attributes assigned to nodes/edges). I would also be very interested in some algorithms that identify inexact matches between such graphs (supervised, unsupervised, any). Any recommendations for those? Any advice would be really appreciated! Thanks in advance!


What are the pros and cons of segmented regression and regression with ARIMA errors for interrupted time-series analysis? In what situation would I choose one method over the other? Are there any papers comparing the two techniques?


title: How to assess performance of an algorithm delivering solutions of varying quality?

text: I lack formal knowledge about statistics and got lost in the Wikipedia articles about the subject. I have an algorithm that produces a solution for a problem. There might be several solutions possible. Not all solutions are equally good, and I am able to give score for each solution. I want to be able to tell if the algorithm is good or not. For each problem I can calculate bounds on the best and worst possible solution (only bounds, not exact values). I can also generate a random solutions for a problem, but not every random solution is a valid one (high chances for getting a valid one, though). The data set I am going to use have problems of approximately the same difficulty. Clearly I should use some statistical test, but which one? It seems that they require knowledge about the distribution and are more suitable for comparing algorithms. **Edit:** I do not want to add too much irrelevant specifics, so I will try to simplify. I have an algorithm that solves a problem. for example, it finds a path between two nodes in a graph. The input is a graph. The output is a path. Obviously, shorter paths are better than longer ones. I want to know if my algorithm is good. I do not know exactly how to define "good", but that is where statistics can help me, is not it? What I though about so far was to do the following:   1. Run the algorithm on a set of problems (problem is a graph). The sizes of all these graphs are of the same order of magnitude.   2. Find a random path (in my specific problem it is possible to do) for each graph, and do this many times.   3. Compare the results of 1 and 2 to see if my algorithm is better than choosing a path at random. Am I doing it right? How should I perform step 3?


title: Calculating energy output of a pellet stove using time series

text: I want to know which brand of wood pellets (i.e., the fuel) is best suited for my stove, to know this I had 2 temperature loggers installed at the air intake and hot air vent. Essentially, if you know how cold the intake air was and how hot it came out, it should be possible to calculate the energy that was achieved during the given time interval. Now, I don't need to know how much kW was developed during this time, a simple 'points' system would suffice, e.g., I could award 0.1 points for each minute for each degree the outgoing air is hotter than the ingoing air. For example, if the ingoing air was 20°C and the outgoing air was 60°C for one minute that would 'earn' this brand of pellets 4 points. For each type of pellet brand, I've burned about 30kg of it and had my 2 temperature loggers on the whole time, which is approximately 27 hours for each brand of pellets. If I could calculate the points each brand has 'earned', it would allow me to see what's the best pellet brand for me. The problem I'm facing, however, is that the temperature loggers log data samples each 30-60 seconds and not each time with the same interval, and also not at the same timestamps with respect to each other. So my data looks kinda like this:               In  - 10/03/2012 - 00.01.15h : 20°C     Out - 10/03/2012 - 00.00.55h : 60°C     In  - 10/03/2012 - 00.01.50h : 21°C     Out - 10/03/2012 - 00.01.35h : 63°C      Now I've hit a wall in trying to apply my points system over this data, since the time stamps do not match up. How can I process data like the samples above to achieve my points value for each burning session?


title: Traffic visualization in Excel: five sources, two metrics, month-by-month

text: I hope you'll suffer a question from a novice :) I'm trying to find the best way to graph (in Excel) the month-by-month share of site-wide visits and pageviews in five referral sources, grouping the remaining share as one: "other sources." My thought was to do a clustered, stacked column chart, where each month has a cluster of two columns—visits and pageviews. Each column in a cluster would be stacked by source. The visits column would use one y-axis, and pageviews would use a secondary y-axis—pageviews have much higher values than visits. Please forgive me if any of this terminology is wrong. Am I on the right track or is there a better way to visualize this data? If this is the way to go, any direction on how to arrange my data to coax Excel into getting all these parameters into one chart? Or does that question belong on a different StackExchange site? Thanks!


title: How to recognize repeated movement patterns of mouse behavior using OpenCV?

text: I have a situation where I have collected mouse movement points from a website. I have a series of (x, y)-points, and I need to detect different repeated patterns of mouse movement from this data. For example, mouse moving very slowly, mouse moving very fast toward a direction and then stopping for a while, mouse scrolling, etc. Is there a way to do that with OpenCV, or maybe some other Library? P.S. Please keep in mind, that I am a beginner in this kind of stuff.


title: Converting Bs parameter estimates into Betas

text: I have unstandadized B coeff.s in my parameter estimates for ANOVA. How can I use Beta=B(SDx/SDy) to convert them into standardized B's (Beta's). Why doesn't SPSS do this for me, like it does in regression? THanks.


title: Sampling for quality control

text: I would like to do some quality control study for a set of manufactured components. For a given set, I can draw several samples. Each drawn sample can be either good or bad. My plan is to create a probability index that can indicate whether this set is good or bad. For instance, 9 of 10 samples are evaluated as good. Then a naïve index for the original set can be 0.9. Are there any statistically sound methods to build this kind of index? I think “sample theory” maybe related, but I am not sure which method will be the best fit so that I can explore without reading whole book without any target. The method should also be based on the size of the set. In practice, the size of the set can vary between 500 to 10000. Certainly, we can have some prior knowledge on the quality of set. For instance, it is very likely that the set can be good or bad.


title: Standard deviation of a coin flip bet

text: What is the standard deviation for the pay-offs of a coin flipping game where you get 2 dollars if heads and lose 1 dollar if tails. How does the standard deviation change if the probabilities remain constant but the pay-off changes to 3 dollars for heads and lose 1 for tails. What is the generalized formula for standard deviation in a coin flip with given probabilities and pay-offs?


title: Can regression to the mean cause this?

text: Let's say I have two independent samples (DS and VS). My data is also highly multivariate (e.g. one dependent variable, 300+ potential predictors, N<300). To reduce the dimensionality of this problem I decided to aggregate the effect of my 300+ potential predictors into one index parameter. Therefore, I regressed each predictor on my dependent variable, which provides me with an effect score (beta value) for each predictor. By multiplying these beta values with the predictors, I get the accumulated net effect (risk index) of all 300+ predictors for a particular observation. Finally, I tested if this risk index is associated with my dependent variable in the Validation Sample (VS). I expected that the beta value for the risk index would fall between 0 (no association) and positive values (predictive power). However, to my surprise I consistently got negative beta values for some sets of predictor variables, especially those that I didn't actually expect to be good predictors. I'm wondering if regression to the mean might be a sufficient explanation for this?


title: Statistical workload estimation

text: I am a software developer and am frequently asked by prospects to provide a workload estimate for a potential project. I would like to implement a statistical approach to estimation that works as follows:   1. Create a list of tasks that make up the project.   2. Have N colleagues provide a low and high estimate (in man/days) for each task.   3. Average the low estimates and high estimates for each task.   4. Assume that the distribution of real possible workloads for each task (again in man/days) is normally distributed with the low and high estimates each one standard deviation from the mean.   5. Calculate workload estimates for the entire project (in man/days) with 10%, 50% and 90% confidence. In other words, based on the assumptions in step 4, what is the total workload estimate where I can be 10/50/90% confident that the real workload will be equal to or less than the estimate? Where I need guidance is step 5. Is there a standard statistical approach that can be used? (Side question: I'd also be interested about how to refine the estimate once we have historical data to feed into the model. Obviously it would be ideal to have the real workload data for each task, but this will be hard to know with certainty since programmers usually work on several tasks in parallel. Realistically the most reliable historical data will be the actual workload required to complete the entire project.)


title: Measuring multiple samples from a specific point vs multiple samples each from a random point

text: I have been reading about sampling and I am wondering about the following: If your interested in measuring the rate of change of your population, ideally I assume you would choose a point at random in the series and then measure the rate of change x minutes later (where x is a randomly generated number between your minimum and maximum value of the period of interest). The aforementioned process would then be repeated a certain number of time to reach the required number of samples. Is it possible to reduce the number of randomly selected start points and instead measure the rate of change, from that start point, over multiple x minute periods? For example, measuring the rate of change experienced after 2,10,17 minutes (all randomly generated) from a single start point before selecting a new start point and calculating the rate of change for 3 more random time durations. I expect the number of start points would still have to be relatively high to avoid any issues due to the fact that each group of 3 rates would be correlated? Are there any rules of thumb regarding the total number of samples vs how many could be calculated from each start point? Walter


title: Problem with character vectors and linear regression in R

text: Quick question. I want to perform a linear regression that looks like this:               lm(y ~ x1 + x2 + x3 + x4 +x5, mydata)      This works fine if I manually write out this code. However, the independent variables that I want to use are stored as a character, like this:               > vars     [1] "x1 + x2 + x3 + x4 +x5"      I tried typing this:               lm(y ~ vars, mydata)     Error in model.frame.default...      But it gives an error! So then I tried this:               lm(y ~ noquote(vars), mydata)     Error in model.frame.default...      And then this               lm(y ~ print(vars, quote = FALSE), mydata)     Error in model.frame.default...      Anyone have a clue how I can get around this problem? The character string in "vars" is being provided to me by a program upstream, so I can't work around it at that level. Thanks!


title: unbalanced response variable

text: I have some general questions. Suppose I have a dataset with binary response Y, and some predictors. If the Y is not quite balanced, with say, many more 1's than 0's.   1. Will it bring any potential problem for further analysis(for instance, logistic regression or classification)?   2. If so, how to define the "much more" and what kind of impact might be introduced?   3. Stratified sampling seems to be an option with sufficient data. Any good reference talking about the empirical advantages and disadvantages?   4. Any other treatment/discussion to this problem? Thank you in advance.


title: Interpret Eviews Output: EGARCH - ARCH and GARCH term

text: I am having some difficultires, figuering out what (and why) the ARCH term in the following output is:   ![enter image description here](http://i.stack.imgur.com/umHVj.png)   _Please note that above's output is fromIntroductory Econometrics for Finance from p. 407_     1. I figuered C(5) out to be the GARCH term, but I am uncertain about the ARCH term. Do you know what it is? And if so: How did you know it?   2. If the sum of the ARCH and GARCH term is =>1, does that mean it is non-stationary (i.e. has a unit root)?   3. Do I interpret the ARCH and GARCH term? I think the ARCH term can be interpreted as the effect of the previous error term on the current error term? **Thank you** so much in advance!


title: Step-wise feature selection with caret

text: can anyone direct me to a package/commands in R for performing step-wise feature selection, preferably using the `caret` package. I have already used linear discriminant analysis (LDA), Random forest, PCA and a wrapper using a support vector machine. I was thinking of including a partial least sqaures or a gradient boosting method, but while trying to use them on multi-class data, they cause R to crash. People have reported similar experiences on multi-class data using `caret` when attempting to use `gbm`. I realize that I haven't used a step-wise approach and I was searching for one that can be implemented on highly correlated, dependent variables for selecting the best performing 20 variables (for example) to create a parsimonious model. Any suggestions would be welcomed


title: From Marginal Exp-Norm Distributions to What Conditionals and Joint?

text: I have (trivariate: multivariate with three variables) data that appears to be good empirical and reasonable theoretical fit for a (univariate) convolution of an exponential and a normal distribution (some times called exp-norm or exGauss distributions). My data are samples from the joint distribution: J(R,G,B) It appears that the marginals of R,G,B follow: $R=R_N+R_E$ with $R_N \sim N(\mu_R,\sigma^2_R)$ and $R_E \sim Exp(\lambda_R)$ (and likewise for G,B). I would like to effectively summarize the marginal, conditional, and joint distributions. The main purpose of the summary is to compare these distributions (marginal, conditional, and joint) with other distributions generated in a like manner. My difficulty is that (1) I don't know a form for the joint (or the conditional) distribution and (2) following from that, I have no parameters to estimate. Thus, I need either a distribution free (non-parametric) approach to working with this data or I need to figure out a multivariate distribution that matches the data and the form of the marginals. Or should I think about other options?


title: How do you abbreviate "Millions"?

text: M or MM? I prefer "M" -- it's shorter, but someone here was pushing for "MM".


title: Need a suitable Metric to indicate closeness between a set of expected vs observed values

text: I have the expected ratio of four parameters A, B, C and D ` A : B : C : D 3.9 : 12.0 : 24.3 : 59.8 ` I perform three experiments which provides me ratios of A,B,C and D. As the results show, the ratios obtained in Experiment 1 are more closer to the Expected ratio Experiment 3 seems to be the worst ` A : B : C : D Experiment 1 3.9 : 12.0 : 24.3 : 59.8   Experiment 2 0.7 : 2.0 : 25.6 : 71.7   Experiment 3 3.9 : 29.9 : 53.3 : 7.2   ` Please suggest a single metric that can clearly indicate the closeness of values obtained in various experiments (as compared to the Expected ratio) Need to do this for a lot of experiments


title: Determining the posterior distribution for an Autoregressive or order 1 model

text: **Question:** For this question, note that the notation $y_{1:T} = (y_1, y_2, \cdots, y_T)$, ie, a vector of random variables. Consider the following AR(1) model: \begin{align*} y_{t+1} = \phi y_t + \sigma \eta_t \ \ \cdots (a) \end{align*} for $t = 1, 2, \cdots, T$ where \begin{align*} \eta_t \stackrel{iid}{\sim} N(0,1) \ \ \cdots (b) \end{align*} with $\eta_1$ independent of $y_k$ for $k \le t$, and where \begin{align*} y_1 \sim N\left(0, \frac{\sigma^2}{1-\phi^2}\right) \ \ \cdots (c) \end{align*} Define $\theta = (\phi, \sigma)$ and consider a prior distribution given by: \begin{align*} p(\theta) \propto \frac{1}{\sigma} \ \text{for} \ \infty < \phi < \infty \ \text{and} \ 0<\sigma<\infty \ \ \cdots (1) \end{align*} Also define the conditional likelihood function to be: \begin{align*} p(y_{2:T} \mid y_1, \theta) = \left(2\pi \sigma^2 \right)^{-\frac{T-1}{2}}\exp\left[- \frac{1}{2\sigma^2}\sum_{t=2}^T \left(y_t - \phi y_{t-1} \right)^2\right] \ \ \cdots (2) \end{align*} Show that the conditional posterior distribution, corresponding to the conditional likelihood function in $(2)$ and the prior density in $(1)$, may be obtained analytically, whereas the full posterior distribution, corresponding to the model in $(a)-(c)$ and the prior density in $(1)$, requires an alternative computational approach. * * * **My Working:** So firstly, I worked out the likelihood function corresponding to the AR(1) model described by $(a)-(c)$. My working is as follows: We are given $y_1 \sim N\left(0, \frac{\sigma^2}{1-\phi^2}\right)$, so the pdf is given by: \begin{align*} p(y_1 \mid \theta) = \left(2\pi\left( \frac{\sigma^2}{1-\phi^2}\right) \right)^{-\frac{1}{2}}\exp\left[ - \frac{y_1^2}{2\left(\frac{\sigma^2}{1-\phi^2}\right)} \right] \end{align*} Conditionally, we have $y_2 \mid y_1 \sim N\left( \phi y_1, \sigma^2\right)$, so the pdf is given by: \begin{align*} p(y_2 \mid y_1, \theta) = \left(2\pi \sigma^2\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2\sigma^2} \left(y_2 - \phi y_1\right)^2 \right] \end{align*} Similarly, $y_3 \mid y_2, y_1 \equiv y_3 \mid y_2 \sim N\left(\phi y_2, \sigma^2\right)$, so the pdf is given by: \begin{align*} p(y_3 \mid y_2, y_1, \theta) = \left(2\pi \sigma^2\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2\sigma^2} \left(y_3 - \phi y_2\right)^2 \right] \end{align*} So in general, $y_t \mid y_{t-1}, y_{t-2}, \cdots, y_1 \equiv y_t \mid y_{t-1} \sim N\left(\phi y_{t-1}, \sigma^2\right)$, with pdf: \begin{align*} p(y_t \mid y_{t-1}, y_{t-2}, \cdots, y_1, \theta) = \left(2\pi \sigma^2\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2\sigma^2} \left(y_t - \phi y_{t-1}\right)^2 \right] \end{align*} Using the method of composition, we have: \begin{align*} p(y_{1:T} \mid \theta) & = p(y_1 \mid \theta)p(y_2 \mid y_1, \theta)p(y_3 \mid y_1, y_2, \theta) \cdots p(y_T \mid y_1, y_2, \cdots, y_{T-1}, \theta) \\\ & = p(y_1 \mid \theta) \prod_{t=2}^T p(y_t \mid y_{t-1}, \theta) \end{align*} Thus the likelihood function computed for a given value of $\theta = \left(\phi, \sigma^2\right)$, is given by: \begin{align} L(\theta) & = \left\\{ \left(2\pi\left( \frac{\sigma^2}{1-\phi^2}\right) \right)^{-\frac{1}{2}}\exp\left[ - \frac{y_1^2}{2\left(\frac{\sigma^2}{1-\phi^2}\right)} \right]\right\\} \prod_{t=2}^T \left(2\pi \sigma^2\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2\sigma^2} \left(y_t - \phi y_{t-1}\right)^2 \right] \\\ & \propto (1-\phi^2)^{\frac{1}{2}}\sigma^{-T} \exp\left[-\frac{\sum_{t=2}^T\left(y_t - \phi y_{t-1}\right)^2+y_1^2\left(1-\phi^2\right)}{2\sigma^2} \right] \ \ \cdots (W1) \end{align} Deriving the conditional likelihood given in $(2)$ can be done as follows: \begin{align*} p(y_{2:T} \mid y_1, \theta) & = \frac{p(y_{1:T} \mid \theta)}{p(y_1 \mid \theta)} \\\ & = \prod_{t=2}^T p(y_t \mid y_{t-1}, \theta) \\\ & = \left(2\pi \sigma^2 \right)^{-\frac{T-1}{2}}\exp\left[- \frac{1}{2\sigma^2}\sum_{t=2}^T \left(y_t - \phi y_{t-1} \right)^2 \right] \end{align*} * * * **My Query:** I don't really get what the question is trying to ask me to do? What does it mean by conditional posterior? Full posterior? Any assistance will be appreciated! * * * **EDIT 1 PROGRESS:** Okay, so I've played around a bit more and made a bit of progress. I interpret 'conditional posterior' as follows: Notice that when $T$ is large, then the factor $(1-\phi^2)$ in Eqn. $(W1)$ is small, so we can approximate the full likelihood with the conditional likelihood: \begin{align*} p(y_{2:T} \mid y_1, \theta) \propto \sigma^{-(T-1)}\exp\left[- \frac{1}{2\sigma^2}\sum_{t=2}^T \left(y_t - \phi y_{t-1} \right)^2 \right] \end{align*} So under the prior $p(\theta) \propto \frac{1}{\sigma}$, we have: \begin{align*} p(\theta \mid y_{1:T}) & \propto \sigma^{-T}\exp\left[- \frac{1}{2\sigma^2}\sum_{t=2}^T \left(y_t - \phi y_{t-1} \right)^2 \right] \\\ & = \sigma^{-T}\exp\left[-\frac{1}{2\sigma^2}\sum_{t=2}^T \left(y_t^2 -2y_t \phi y_{t-1} + \phi^2 y_{t-1}^2 \right) \right] \\\ & = \sigma^{-T}\exp\left[-\frac{1}{2\sigma^2}\left(\underbrace{\sum_{t=2}^Ty_t^2}_{C} - 2\phi\underbrace{\sum_{t=2}^Ty_t y_{t-1}}_{B} + \phi^2 \underbrace{\sum_{t=2}^T y_{t-1}^2}_{A} \right) \right] \end{align*} First, note that: \begin{align*} p(\phi \mid \sigma, y_{1:T}) & \propto \exp\left[-\frac{1}{2\sigma^2}\left(C-2\phi B + \phi^2A \right) \right] \\\ & = \exp\left[-\frac{A}{2\sigma^2} \left(\phi^2 - 2\phi \frac{B}{A} + \frac{C}{A} \right) \right] \\\ & = \exp\left[-\frac{A}{2\sigma^2}\left(\left(\phi - \frac{B}{A} \right)^2-\left(\frac{B}{A}\right)^2 + \frac{C}{A} \right) \right] \\\ & \propto \exp\left[-\frac{A}{2\sigma^2}\left(\phi-\frac{B}{A}\right)^2 \right] \\\ & = \exp\left[-\frac{1}{2\left(\frac{\sigma^2}{A} \right)}\left(\phi-\frac{B}{A}\right)^2 \right] \end{align*} So the distribution of $\phi \mid \sigma, y_{1:T}$ is given by: \begin{gather} \phi \mid \sigma, y_{1:T} \sim N\left(\frac{B}{A}, \frac{\sigma^2}{A} \right) \\\ \implies \phi \mid \sigma, y_{1:T} \sim N\left(\frac{\sum_{t=2}^T y_t y_{t-1}}{\sum_{t=2}^T y_{t-1}^2}, \frac{\sigma^2}{\sum_{t=2}^T y_{t-1}^2} \right) \end{gather} Next, note that: \begin{align*} p(\sigma \mid y_{1:T}) & \propto \int_{\phi} \sigma^{-T} \exp\left[-\frac{1}{2\sigma^2} \sum_{t=2}^T \left(y_t - \phi y_{t-1}\right)^2 \right]d\phi \\\ & = \sigma^{-T} \int_{\phi} \exp\left[-\frac{A}{2\sigma^2} \left(\left(\phi - \frac{B}{A} \right)^2 - \left(\frac{B}{A}\right)^2 + \frac{C}{A} \right) \right]d\phi \\\ & = \sigma^{-T} \int_{\phi} \exp\left[-\frac{A}{2\sigma^2} \left(\phi - \frac{B}{A} \right)^2 + \left(\frac{A}{2\sigma^2}\right) \left(\frac{B}{A}\right)^2 - \left(\frac{A}{2\sigma^2}\right)\left(\frac{C}{A}\right) \right]d\phi \\\ & = \sigma^{-T} \exp\left[\frac{B^2/A - C}{2\sigma^2} \right] \int_{\phi} \exp\left[-\frac{A}{2\sigma^2} \left(\phi - \frac{B}{A} \right)^2\right] d\phi \\\ & = \sigma^{-T} \exp\left[\frac{B^2/A - C}{2\sigma^2} \right] \left(2\pi\left(\frac{\sigma^2}{A} \right) \right)^{\frac{1}{2}} \\\ & \propto \frac{1}{\sigma^{(T-2)+1}}\exp\left[\frac{B^2/A - C}{2\sigma^2} \right] \end{align*} Now define $b = \frac{B}{A}$, notice: \begin{align*} Q(y_{2:T}, b) & = \sum_{t=2}^T \left(y_t - b y_{t-1}\right)^2 \\\ & = \sum_{t=2}^T \left(y_t^2 - 2y_t b y_{t-1} + b^2 y_{t-1}^2 \right) \\\ & = \sum_{t=2}^T y_t^2 - 2b\sum_{t=2}^T y_t y_{t-1} + b^2 \sum_{t=2}^T y_{t-1}^2 \\\ & = C - 2bB+b^2A \end{align*} Then clearly, \begin{align*} -Q(y_{2:T}, b) = B^2/A-C \end{align*} So the distribution of $\sigma \mid y_{1:T}$ is given by: \begin{gather} \sigma \mid y_{1:T} \sim IG\left(v = T-2, \widehat{\sigma^2} = \frac{1}{T-2} \left(B^2/A - C \right) \right) \\\ \implies \sigma \mid y_{1:T} \sim IG\left(v = T-2, \widehat{\sigma}^2 = \frac{1}{T-2} \left(\frac{\left(\sum_{t=2}^T y_t y_{t-1} \right)^2}{\sum_{t=2}^T y_{t-1}^2} - \sum_{t=2}^T y_t^2 \right) \right) \end{gather} Or equivalently: \begin{gather*} \sigma \mid y_{1:T} \sim IG\left(v = T-2, \widehat{\sigma^2} = -\frac{1}{T-2}Q(y_{2:T}, b) \right) \end{gather*} Thus, we can derive the conditional posterior distribution analytically as: \begin{align*} p(\phi, \sigma \mid y_{1:T}) = p(\phi \mid \sigma, y_{1:T})p(\sigma \mid y_{1:T}) \end{align*} where $p(\phi \mid \sigma, y_{1:T})$ and $p(\sigma \mid y_{1:T})$ are derived above. * * * However, if we use the full likelihood (which is what the second part of the question is asking), how can we derive the joint posterior? I do not see any obvious ways to find the appropriate integrating constants. I'm assuming I need to use Gibbs/M-H/or some other kind of MCMC sampling scheme?


title: Genetic algorithms, genetic programming or machine learning algorithms for solving this problem

text: I have a problem that consists of finding the optimal solution based on the following criteria:   1. Logic for identifying that event A has occurred (i.e. "find" logic that most accurately categorises an A type event)   2. Logic for identifying that event B has occurred (i.e. "find" logic that most accurately categorises a B type event)   3. Transition states (finding a parsimonious number of states in which the system can be said to occupy)   4. State transition logic (Which state to transition to based on an event) I have been wrestling with this problem for over a year now, and I have managed to identify that the problem decomposes into combining the four sub problems above into one that produces the optimal solution. For the purpose of brevity, assume that I have an objective function (or fitness evaluation function). The fitness function will return a score based on how a proposed solution performs with new data. My question has been stated with an unintentional bias to GA, but thats simply because that is the approach that (to my perhaps naive mind), appears to be the most promising line of enquiry. I think GP may be useful here (though I have no experience with it), because as requirement 1,2 and 4 suggest, the "optimal" logic to use is not yet known. As an aside, I have _some_ ideas about what possible solutions using this logic are likely to entail. It's just that I'm not sure which approach would produce the "fittest" result (pardon the GA terminology, I don't mean to bias the answers). My question is that given such an "optimisation" problem, what is the best approach for solving it - GA, GP, ML (or some combination of all three)? More (equally?) importantly, what are the practical steps I need to take to implement a solution that can help me solve this problem? Last but not the least (it may be too early to start talking about the implementation details), but I would like to implement the solution using Python, if possible. However, I can also use R, Octave etc if there are existing packages/toolkits etc to prevent me from reinventing the wheel.


title: Mixed model vs. Pooling Standard Errors for Multi-site Studies - Why is a Mixed Model So Much More Efficient?

text: I've got a data set consisting of a series of "broken stick" monthly case counts from a handful of sites. I'm trying to get a single summary estimate from two different techniques: Technique 1: Fit a "broken stick" with a Poisson GLM with a 0/1 indicator variable, and using a time and time^2 variable to control for trends in time. That 0/1 indicator variable's estimate and SE are pooled using a pretty straight up and down method of moments technique, or using the tlnise package in R to get a "Bayesian" estimate. This is similar to what Peng and Dominici do with air pollution data, but with fewer sites (~a dozen). Technique 2: Abandon some of the site-specific control for trends in time and use a linear mixed model. Particularly:               lmer(cases ~ indicator + (1+month+I(month^2) + offset(log(p)), family="poisson", data=data)      My question involves the standard errors that come out of these estimates. Technique 1's standard error, which is actually using a weekly rather than monthly time set and thus should have _more_ precision, has a standard error on the estimate of ~0.206 for the Method of Moments approach and ~0.306 for the tlnise. The lmer method gives a standard error of ~0.09. The effect estimates are reasonably close, so it doesn't seem to be that they're just zeroing in on different summary estimates as much as the mixed model is vastly more efficient. Is that something that's reasonable to expect? If so, why are mixed models so much more efficient? Is this a general phenomena, or a specific result of this model?


title: Gibbs sampler for local linear trend model

text: **Question:** Consider the local linear trend model given by: \begin{align*} y_t = \mu_t + \tau \varepsilon_t \ \cdots \ \text{Observation equation} \\\ \mu_{t+1} = \phi \mu_t + \eta_t \ \cdots \ \text{State equation} \end{align*} for $t = 1, 2, \cdots, T$, where $(\varepsilon_t, \eta_t)'$ is independent of $\mu_k$ for $k \le t$ and where: \begin{align*} \begin{bmatrix} \varepsilon_t \\\ \eta_t \end{bmatrix} \stackrel{i.i.d}{\sim} N\left(\begin{bmatrix} 0 \\\ 0 \end{bmatrix}, \begin{bmatrix} 1 & 0 \\\ 0 & 1 \end{bmatrix}\right) \end{align*} and \begin{align*} \mu_1 \sim N\left(0, \frac{1}{1-\phi^2} \right) \end{align*} (Note everything is a scalar quantity in this question). Consider a Bayesian analysis of this model under the prior distribution given by: \begin{align*} p(\theta) \propto \frac{1}{\tau} \ \text{for} \ -\infty < \phi < \infty \ \text{and} \ 0<\tau<\infty \end{align*} where $\theta = (\phi, \tau)$. Devise a Gibbs sampler to sample from the joint posterior distribution, $p(\mu_{1:T}, \theta \mid y_{1:T})$ where the notation $\mu_{1:T}$ denotes $(\mu_1, \mu_2, \cdots, \mu_T)$ and similarly, $y_{1:T}$ denotes $(y_1, y_2, \cdots, y_T)$. * * * **My Working So Far:** The joint posterior distribution is given by: \begin{align*} p(\mu_{1:T}, \theta \mid y_{1:T}) & \propto p(\mu_{1:T}, \theta, y_{1:T}) \\\ & = \underbrace{p(y_{1:T} \mid \mu_{1:T}, \theta)}_{\text{'likelihood'}}\underbrace{p(\mu_{1:T} \mid \theta)p(\theta)}_{\text{prior}} \\\ & = \left[\prod_{t=1}^{T} p(y_t \mid \mu_t, \theta)\right]\left[\prod_{t=1}^{T-1}p(\mu_{t+1} \mid \mu_t, \theta) \right]p(\mu_1 \mid \theta)p(\theta) \ \ \cdots \ \ (1) \end{align*} Since $y_t \mid \mu_t, \theta \sim N(\mu_t, \tau^2)$ for $t=1, 2, \cdots, T$, the pdf is given by: \begin{gather*} p(y_t \mid \mu_t, \theta) = \left(2\pi \tau^2\right)^{-\frac{1}{2}}\exp\left[-\frac{1}{2\tau^2}\left(y_t - \mu_t\right)^2 \right] \end{gather*} Similarly, $\mu_{t+1} \mid \mu_t, \theta \sim N\left(\phi \mu_t, 1 \right)$, so the pdf is given by: \begin{gather*} p(\mu_{t+1} \mid \mu_t, \theta) = \left(2\pi\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2}\left(\mu_{t+1} - \phi \mu_t \right)^2 \right] \end{gather*} We know that $\mu_1 \mid \theta \sim N\left(0, \frac{1}{1-\phi^2} \right)$, so the pdf is given by: \begin{gather*} p(\mu_1 \mid \theta) = \left(2\pi\left(\frac{1}{1-\phi^2} \right) \right)^{-\frac{1}{2}} \exp\left[-\frac{\mu_1^2}{2\left(\frac{1}{1-\phi^2}\right)} \right] \end{gather*} Finally, we are given that $p\left(\theta\right) \propto \frac{1}{\tau}$. Substituting all of the above into Eqn. $(1)$, yields the joint posterior distribution: \begin{align*} p(\mu_{1:T}, \theta \mid y_{1:T}) & \propto \left[\prod_{t=1}^{T} \left(2\pi \tau^2\right)^{-\frac{1}{2}}\exp\left[-\frac{1}{2\tau^2}\left(y_t - \mu_t\right)^2 \right]\right]\left[\prod_{t=1}^{T-1}\left(2\pi\right)^{-\frac{1}{2}} \exp\left[-\frac{1}{2}\left(\mu_{t+1} - \phi \mu_t \right)^2 \right] \right] \\\ & \times \left(2\pi\left(\frac{1}{1-\phi^2} \right) \right)^{-\frac{1}{2}} \exp\left[-\frac{\mu_1^2}{2\left(\frac{1}{1-\phi^2}\right)} \right] \left( \frac{1}{\tau}\right) \end{align*} I will implement a "blocked" Gibbs sampler sampling $\mu_{1:T}^{(i)}$ together, as follows: For $i = 1, 2, \cdots, M$, sample: \begin{align*} \mu_{1:T}^{(i)} & \sim \mu_{1:T} \mid \phi^{(i-1)}, \tau^{(i-1)}, y_{1:T} \ \ \cdots \ \ (2)\\\ \phi^{(i)} & \sim \phi \mid \mu_{1:T}^{(i)}, \tau^{(i-1)}, y_{1:T} \ \ \cdots \ \ (3) \\\ \tau^{(i)} & \sim \tau \mid \mu_{1:T}^{(i)}, \phi^{(i)}, y_{1:T} \ \ \cdots \ \ (4) \end{align*} Sampling from $(2)$ is straightforward by using the Forward Filter Backwards Sampling (FFBS) **My Query:** I am stuck on how to sample from $(3)$ and $(4)$, in order to use a Gibbs sampler on $(3)$ and $(4)$, we need to find the full conditional of $\phi \mid \mu_{1:T}, \tau, y_{1:T}$ and $\tau \mid \mu_{1:T}, \phi, y_{1:T}$, but how do you find these full conditionals? I do not see any obvious way by examining the joint posterior distribution. Perhaps a Metropolis Hastings subchain could work? But then what should I pick for my candidate density for the $\phi$ and $\tau$? Thanks in advance.


title: Trace(AB)=Trace(BA)? even if A, B are vectors?

text: I have read this vector manipulation in standard books: $$E[XX^T] =E[\mathrm{trace}(XX^T)]$$ where $X^T$ is the transpose of $X$, $X$ has a normal distribution and has dimension $n\times 1$, $XX^T$ has dimension $n\times n$. How can they introduce trace into expectation?


title: Residuals from model missing interaction

text: In a plot of residuals against fitted values from a generalised linear model, I'm wondering what the plot would look like if an interaction was missing from the model. Can anyone simulate a model that is missing an interaction. Without the interaction included in the model, the residuals against fitted values plot shows clear patterns. With the interaction, the residuals against fitted values does not show any patterns.


title: energy adjusted nutrient intake

text: I have a data on 7 nutrients (intake of 150 people - 75 men and 75 women). Can i do factor analysis to derive the nutrient intake pattern? If so shouldn't each nutrient be adjusted with the total energy intake? I would like to do a nutrient intake pattern analysis with the data that I have. Intakes are available for protein, fat, carbohydrates, calcium, iron, vitamin A and vitamin C. Only protein and fat seems to be normally distributed. What should I do?


title: $\nu$ SVM in terms of C-SVM

text: I have an implementation that solves the $C$-SVM optimization problem. Is it possible to use this algorithm to create a $\nu$-SVM? I know there is a connection saying > $\nu$-SVM leads to $\rho > 0$, then $C$-SVM with $C$ set a > priori to > $1/\rho$ leads to the same decision function. (New Support Vector > Algorithms) But can I use a given $\nu$ to derive the parameter $C$ and train the $C$-SVM?


I am developing a questionnaire to measure four factors which constitute spirituality, and I would like to ask the following question: Are data transformations on non-normal data necessary for an exploratory factor analysis when using the principal axis factoring extraction method? I finished screening my data yesterday, and I found that 3 out of 20 questions are positively skewed while 1 out of 20 is negatively skewed (Question 6 = 4.88, Question 9 = 7.22, Question 12 = 11.11, Question 16 = -6.26). I also found that 1 of the questions (out of 20) is leptokurtic (Question 12 = 12.21). I chose the principal axis factoring extraction method because I read that it is used on "severely non-normal data" while maximum-likelihood is used on normal data, but:   1. How would I know if my data is "severely" non-normal?    2. If my data is "severely non-normal", does this mean that I can leave the data as it is now (not transform it) and analyze it using the principal axis factoring extraction method? Or do I need to transform the data before proceeding with the EFA?    3. If I do need to transform the data, what transformations would I use for positively skewed, negatively skewed, and leptokurtic items?


title: "Convergence for 1st lambda value not reached"-error using GLMNET package and specifying lambda parameter

text: I get a weird problem when I specify lambda in the function glmnet, that does not appear if I let the function go through all the lambdas. When I run the following lines, it works great:               glmnet.out = glmnet(y = yy[train_id], x=xx[train_id,], family = "binomial", alpha = 1)     Warning message:     from glmnet Fortran code (error code -82); Convergence for 82th lambda value not reached after maxit=100000 iterations; solutions for larger lambdas returned           glmnet.out$lambda     [1] 2.126060e-02 1.937187e-02 1.765092e-02 1.608287e-02 1.465411e-02 1.335228e-02 1.216610e-02 1.108530e-02 1.010051e-02 9.203207e-03     [11] 8.385619e-03 7.640664e-03 6.961888e-03 6.343413e-03 5.779882e-03 5.266413e-03 4.798560e-03 4.372269e-03 3.983849e-03 3.629934e-03      and more rows, up to the 82th value. The warning is ok, as I don't use that small lambdas. The cv.glmnet gives a lambda around 0.001. When I specify a lambda, e.g. the last in the output above,               glmnet.out = glmnet(y = yy[train_id], x=xx[train_id,], family = "binomial", alpha = 1, lambda = 3.629934e-03)      then I get the following warning message:               Warning messages:     1: from glmnet Fortran code (error code -1); Convergence for 1th lambda value not reached after maxit=100000 iterations; solutions for larger lambdas returned      2: In getcoef(fit, nvars, nx, vnames) :an empty model has been returned; probably a convergence issue      Note that the 1st lambda value is infinity, so nothing is returned. How is this possible? I can get around the problem by not specifying lambda, and then manually pick the estimates for prediction with the cross validation dataset like this:               coefficients = as.matrix(coef(glmnet.out)[,selected_lamdda])     xBeta = cbind(rep(1,length(cv_id)),xx[cv_id,])%*%as.matrix(coefficients)     prediction = 1/(1+exp(-xBeta))      but it is quite inefficient, as the glmnet does a lot of unnecessary calculations. Thank you,


title: Multilevel Modeling and multiple testing

text: I am interested in working out the correct way to correct for multiple testing in multilevel models for longitudinal data, where I am investing a potential interaction between two predictors. My data consists of 7 outcome variables, 2 predictor variables with a interaction, and 5 covariates. For the model building strategy I am using a set up approach. Model 1: Null model   Model 2: Time is introduced   Model 3: Covariates are introduced   Model 4: Predictors are introduced   Model 5: Interaction between Predictors are introduced. This strategy is then repeated for each of the outcome variables. Currently I am using bonferroni correction to account for multiple testing, .05/35 = .001, were 35 is 7 outcome variables multiplied by 5 models. Is this the correct way to perform the calculation for the bonferroni correction? Alternatively would you divide .05 by the number of final models in the analysis (7) or the number of parameters in the models been analyzed?


title: A way to use GPU hardware in matlab

text: I know that theano is a python library for using gpu hardware and make effective implementations. Is there any such library or a way to do the same in matlab?


title: Metropolis independence sampler

text: I need to implement an algorithm of Metropolis Independence Sampler, where the proposal distribution is a Multivariate Normal with 3 parameters. Since my function $f$ is really difficult to analyize,I was thinking if is reasonable to choose as $mean$ of my normal the point that maximizes $f$ (this I can find!) and only adjust the covariance matrix as desired. Thanks in advance :)


title: p-values for p-values?

text: For the purpose of this post, accept that there is some critical threshold below which a finding is declared to be statistically significant. The _p_ -value is a statistic because it is a function of the data. Statistics are random. We can put confidence intervals around statistics. Why don't we put confidence intervals around _p_ -values? Or, equivalently, why don't we test the hypothesis that the observed _p_ -value is below 0.05, say? One answer could be that it would be hard to calculate analytically, but bootstrapping could solve that issue.


title: Expressing a time series in more detailed time interval using mathematical approximations

text: I have a time series with a dependent variable (electricity comsumption) measured on a daily basis. I have also an independent variable (unemployment) measures in monthly intervlas. I want to use an ARIMAX type scenario analysis to see the impact of changes in future unemployment to the electricity consumption. First oi want to buid a model of electricity consumption with the unemployment as the independent variabe but i want to use the time series on the daily basis. What should i do to transform the unemployment on a daily basis? I cannot find data of unemployment on a daily basios since it is reported by the statistical service only in monthly ibtervals. Idea 1) Divide the monthly unemployment rate with the number of days in a speoicif month. It does not sound good to me. Idea 2) Produce a polynomial inteprolation (spline) in the montly series of unemployment and create continuois smooth approximation of the unemployment each day of the month. What do you think is the best idea? DO you have any other ideas to propose on how to tackle the problem? Thnakls in advance, Andreas


title: Estimating parameters of Log Pearson III distribution in R

text: How to estimate the Parameters of "Log Pearson III" distribution and how to generate the random nos. using these estimated parameters in R? Kindly guide Akshata


title: Do you need to stratify / match sample for ANCOVA comparisons?

text: I have data about body image and weight change behaviours from $N = 397$ adolescent boys that we collected in 2000. Age ranged from $13.5$ years to $18.67$ years and the mean age was $15.92~ [SD= 1.02]$. We recently conducted another round of data collection in the same schools (plus additional schools that have a similar socioeconomic profile), using the same questionnaire, and we wish to compare body image and weight change behaviours between 2000 and 2012. As there were no year seven students in the initial cohort, the year 7 students were removed from the 2012 cohort. We therefore have $n = 892$ participants in the 2012 cohort. Age in this group ranges from $12.33-21.08$ years and the mean age was $15.09~ [SD=1.26]$. ANOVA indicated that there was still a significant difference between the age of the participants in the two waves of data collection $F(1,1281)= 135.91,~ p <.000$. Body image and weight change behaviours are closely linked to age. My questions:   * Is there a problem comparing data from groups with these sample sizes ($397$ vs. $892$) using ANCOVA?    * Is ANCOVA enough to limit the impact of age on analyses? or should I match the sample for age by selecting at random from each age group in the recent data?


title: Pre-Post analysis with repeated measure

text: I've 50 subjects and ran an experiement 4 times for each subject and collected the pre and post vlaues for each run. I need to do pre post analysis. Here a **sample** of my data:               **ID  experiement   PrePost Value**       1         1        pre    1.2       1         2        post   5.8       1         3        post   4.6       1         4        pre    2.7     ....       2         1        pre    3.2       2         2       post    7.9       2         3       pre     3.5       2         4       post    6.6     ....       3         1       pre     5.1       3         2       post    4.8       3         3       pre     5.7       3         4       post    2.5      How can I do this in SAS and check for the normality assumptions?


title: Coefficients of Linear Discriminants in R

text: I've read the answers in What are "coefficients of linear discriminants" in LDA?, but I still don't understand what **coefficients of linear discriminants** on output of R means. What is it? (How) Is it related to the decision boundary? nb: my knowledge about LDA can be summed up in this slide.               library(ISLR, MASS)     train <- (Smarket$Year < 2005)     lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)         Call:     lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)          Prior probabilities of groups:         Down       Up      0.491984 0.508016           Group means:                 Lag1        Lag2     Down  0.04279022  0.03389409     Up   -0.03954635 -0.03132544          Coefficients of linear discriminants:                 LD1     Lag1 -0.6420190     Lag2 -0.5135293


title: How to do meta-analysis and create forest plots using Beta correlation coefficients

text: Does anyone knows how to do meta-analysis and create forest plots using beta- correlation coefficients?


title: Chi-square using factors with multiple levels in R

text: I'm not sure if I have the right concept of how to perform the chi-square test. I have a variable called `race` which is a factor with multiple levels for different races, and I would like to see whether there is a correlation between each race and the response variable (`mortality`). If I were to perform chi-square like:               chi.race <- chisq.test(mortality, race)      then it gives me an output like:               Pearson's Chi-squared test     data:  mortality and race     X-squared = 4.9626, df = 9, p-value = 0.8376      but it doesn't tell me anything about the correlation between each level in `race` (e.g. White, Black, Asian, Hispanic, etc.) and `mortality`. What am I doing wrong and/or what is wrong with my concept of chi-square? Thanks for the help.


title: How to interpret coefficients in a Poisson regression with interaction terms?

text: This question is a prolongation of this question: How to interpret coefficients in a Poisson regression? If we follow the (almost) exact same routine, but we add correlation between the variablese treatment and improved (just for the sake of my question, which is interpreting the output), we get:               treatment     <- factor(rep(c(1, 2), c(43, 41)),                              levels = c(1, 2),                             labels = c("placebo", "treated"))     improved      <- factor(rep(c(1, 2, 3, 1, 2, 3), c(29, 7, 7, 13, 7, 21)),                             levels = c(1, 2, 3),                             labels = c("none", "some", "marked"))         numberofdrugs <- rpois(84, 10) + 1         healthvalue   <- rpois(84, 5)        y             <- data.frame(healthvalue, numberofdrugs, treatment, improved)     test          <- glm(healthvalue~numberofdrugs+treatment+improved + treatment:improved, y, family=poisson)     summary(test)      Note the $\textbf{ treatment:improved}$ term I added inside the glm function. Now, we get the following output:                   Call:     glm(formula = healthvalue ~ numberofdrugs + treatment + improved +          treatment:improved, family = poisson, data = y)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.9261  -0.8733  -0.0296   0.5473   2.3358            Coefficients:                                      Estimate Std. Error z value Pr(>|z|)         (Intercept)                      1.553051   0.184229   8.430   <2e-16 ***     numberofdrugs                    0.004298   0.014242   0.302   0.7628         treatmenttreated                 0.007399   0.149440   0.050   0.9605         improvedsome                     0.358897   0.164891   2.177   0.0295 *       improvedmarked                  -0.178360   0.203756  -0.875   0.3814         treatmenttreated:improvedsome   -0.330336   0.265310  -1.245   0.2131         treatmenttreated:improvedmarked  0.050617   0.260203   0.195   0.8458         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 97.805  on 83  degrees of freedom     Residual deviance: 89.276  on 77  degrees of freedom     AIC: 383.29          Number of Fisher Scoring iterations: 5      If we ignore what seems to be insignificant coefficients, I can ask my question: I understand that, as in the original post, treatment=placebo and improved=none is the base level for those variables, and thus are set to zero. My question is, why does it not exist any interaction terms with the base lavels for treatment=placebo and improved=none? I thought setting the base levels to zero was just a construct, and in my mind there should still exist correlation between them...(?)


title: Unknown name of this test

text: I would like to know the exact name of a test for measuring the similarity of 2 noise systems. Assuming m1 being the measurement for system 1 and m2 the measurement for system 2 this test computes S1 = (m1 + m2)/2 and S2 = (m1 - m2)/2. Then it plots S1 vs S2 for each measurement point. Thanks a lot ps no one of the 2 systems is the ground truth


title: how can the probabilities of sampling subset proportions be tested?

text: Assume we have 10 different jigsaw puzzles. Each puzzle has a different number of pieces, e.g.,               1st puzzle = 100 pieces     2nd puzzle = 120 pieces     3rd puzzle = 240 pieces     .     .     .      Further, assume we take all the pieces from all the different puzzles and put them together in one pile (n = 100 + 120 + 240 + ...). Now, assuming we take a sample of X pieces from the pile and find that we have picked a certain number of pieces from each puzzle such that:               X = 100     3/100 cover 1st puzzle 3/100     23/100 cover 2nd puzzle 23/120     12/100 cover 3rd puzzle 12/240      I am interested in calculating the probability of getting a certain coverage (or proportion) of each puzzle (separately) from sample X. I am not interested in the probability of the proportion distribution (as can be found using chi squared goodness of fit). Essentially, I would like to show that a very high coverage of a small puzzle has a very low probability in a small sample. My questions are:   1. Assuming the sampling is random is it valid to use a one-tailed z-statistic to find the p-value.   2. Assuming the sampling is initially random, but from sample X only those pieces that answer a certain criteria (e.g., colour or wight) are considered in the puzzle coverage calculation, is the z-statistic test still valid?   3. Finally, what test would be appropriate if some pieces fit more than one puzzle? (Just to clarify, this is not homework!)


title: Adding noise to a column of data

text: I am working with this optdigits data set from UCI machine learning repository and want to create a new training dataset with noise. How to randomly add noise to a vector in R? say corrupt 10% of the values of the vector.


title: How big of a dataset can R or SAS handle for regression?

text: On a standard computer (~3-6gb RAM) with 2 - 4 processors. What are the size restrictions on a matrix of data for R or SAS?


title: Using Davies-Bouldin index in clustering

text: I am clustering data using k-medoid. I used **Davies–Bouldin** index for $2$ to $n-1$ clusters. Here $n = 100$ (using smaller test case). I find minimal value of the index for 98 clusters. But the overall accuracy rate for 98 cluster is very small (smaller than 1). Here accuracy rate is how accurately test data is matched to training data. What should I do in that situation. If dataset is larger then finding Davies–Bouldin value from $2$ to $n-1$ is large task. What should I do for larger dataset? Here is my plot of Davies–Bouldin index value for cluster solutions [X axis "index" is actually the number of clusters in a solution]. ![enter image description here](http://i.stack.imgur.com/Lakqz.png)


title: Association analysis returns 0 useful rules

text: Since new to association rules need help for identifying the most frequent extra service ordered together with the products. And to derive the association rules, I have used `arules` R package on past data. I have collected and shaped the dataset as like you can download sample dataset from link (https://www.dropbox.com/s/4g5sqag9afsttmo/finaldf2.csv). When I tried it with the following set of R codes, it can't calculate the association rules between product and services. But it calculates the rules between all services. Can any one guide me why I am getting set of 0 rules for calculation of association rules of product. Am I doing something wrong to calculate the associations. Or My dataset should be different. Here are the R codes generated for `arules`:               # Loading the library     library(arules)              # Loading the dataset      product.transaction.set <- read.csv("finaldf2.csv")              # calculating the association with apriori algo     rules <- apriori(product.transaction.set, parameter = list(support = 0.4,conf = 0.6,target = "rules"))          ## 1 - for product     # Deriving the association rules for productname     rules.sub.product <- subset(rules, subset = lhs %pin% "LD") # output -> set of 0 rules               # for print the calculated association rules     inspect(rules.sub.product)           ## 2- for add on services ()     # Deriving the association rules for add on services     rules.sub.services <- subset(rules, subset = lhs %pin% "AO1") # output -> set of 1 rules           # for print the calculated association rules     inspect(rules.sub.services)


title: Variance of average of estimate computed by importance sampling

text: I have a random variable X of which I sample N values [$x_{1}$...$x_{N}$]. From these values I calculate the estimate P of function H(x) using Importance Sampling, i.e. $P = \sum_{i=1}^{N} w_{i}H(x_{i})$   ( $H(x_{i})$ gives as output either 0 or 1) .   This is done for T trials and for each trial ' $t$ ', I have two output vectors -     1. [$P_{0}$, $P_{1}$, ..., $P_{t}$] , and      2. [$Pavg_{0}$, $Pavg_{1}$, ..., $Pavg_{t}$] where $Pavg_{i}$ is the average of all P's upto trial t i.e. $Pavg_{i}= 1/t*\sum_{k=1}^{t}P_{k}$ Question is what is the variance of $Pavg_{t}$. I want this value so that I can know how close my simulation results are to the actual case. Should I just calculate the variance($Pavg_{t}$) from the vector(2)?   I also came across the formula to calculate variance of estimate in Monte Carlo, given as $Var(P_{MC}) = P_{MC}(1-P_{MC})/N$, where $P_{MC}$ is the monte-carlo estimate. Do I use this to calculate $Var(Pavg_{t})$?


title: Statistics regarding different data sets having a common response variable

text: Lets assume one has 2 datasets: with different number of rows (samples) and columns (features). Each of these 2 datasets have a column as a binary response variable. Lets say healthy or not. What sort of statistical methods can be used to help us for better feature selection results, or better performance in classification models, utilizing both datasets?


title: Which numbers are least likely to be selected by people in a lottery?

text: The Mega Millions is over $500 million today. I remember reading a JSTOR paper about some numbers that are most unlikely to be chosen. For example lots of people choose 7 because it's their lucky number, and I want the opposite of that. However my JSTOR membership has ran out. Which numbers are least likely to be chosen by people as their lottery selections in a lottery of numbers between 1 and 80? **Note:** Each number has an equal chance of being selected; I want to choose numbers no one else does so I don't have to share the prize with anyone if I win.


title: Small and unequal sample sizes

text: I have a data set consisting of a small number of participants: n=5, n=5, n=2. I want to compare the means on a particular variable. I've done a 1-way anova and the p-value was .07, the Welch test gives me p = .005. Which is correct?


title: How to apply clustering analysis to help identify criminal entities out from credit card usage data?

text: Let's say I have a ton of credit card usage data and have also some means to predict if a given transaction is fraudulent. Now I want to know what kind of criminal entities are behind these frauds. Maybe some are big international/domestic crime organizations, some are small groups and some are individuals etc. I'd assume that clustering analysis would help making the first step in this situation. What combination of techniques are useful in this scenario, and how do you apply them? _EDIT_ : To clarify, what I want to do is clustering "similar" fraudulent transactions so that I can make some guess on from whom/what they originated. Maybe certain types of frauds tend to come from a certain region. Some other types target certain combinations of flaw in the system (e.g. lot of small transactions vs. bigger transactions) etc. In principle, I want to extract some pattern that will help the investigators.


title: Calculating quantiles for chi squared distribution

text: For chi squared distribution, how would I find quantiles for the following three cases: A) $P(X^2> X^2_{\alpha})=0.01$ when $v = 21$ B) $P(X^2 < X^2_{\alpha})=0.95$ when $v =6$ C) $P(X^2_{\alpha} < X^2 <23.209) = 0.015$ when $v = 10$ Here $X^2_{\alpha}$ is the $\alpha$-quantile of the $\chi^2_{v}$ distribution.


title: How to analyze GEE with unevenly spaced observations?

text: I am interested in using Generalized Estimating Equations (GEE) to model longitudinal count data. I recorded animal count observations on the same sites on many days but the spacing of the observations was uneven. I was wondering if there is a way to specify this uneven structure when using an autoregressive correlation structure? for example, when using lme from the nlme package to analyze linear mixed models it is possible to specify the correlation structure as:               correlation = corAR1(form = ~DAY | Subject)      That way the spacing of the observations for the autoregressive structure is adjusted based on DAY. Is there a way to do this in R in a GEE model using a package such as geepack, gee, or any others? Thanks, Dan


title: Validating a linear model with R, lm()

text: I've created a model (cue ominous music) in R based on previous months data using lm(). Now, I would like to see how well it predicts the current months data. For example, my model predicts sales figures. I have figures for January, February and March. My model is based on figures from January and February and I would like to estimate how good of a predictor the model is of March's sales. Is there a method for doing this without having to manually apply all the coeficients and weights?


title: How to input self-defined distance function in R?

text: I want to know how to to input a self-defined distance in R, in hierarchical clustering analysis. R implements only some default distance metrics, for example "Euclidean", "Manhattan" etc. Suppose I want to input a self-defined distance '1-cos(x-y)'. Then what should I do? Writing a function is obviously a solution. But, it will be quite complicated, and also difficult to write. Please help me. I am unable to write the code.


title: Finding similarities using Wavelet transform

text: I have a time serie and I want to find similarities in it. For the first step I have calculated Haar-wavelet coefficients for this time serie, and now I don't know exactly how should I continue should I extract features from this transformed data, to find similarities? how can I do that?


title: How to analyze this incomplete block design in R?

text: I was just wondering if you could point me in the right direction. I have a dataset with 5 tree clones planted at 10 different sites, i.e. same clones are replicated twice at different sites.                 Site Clone           1     A             2     A           3     B           4     B           5     C           6     C           7     D           8     D           9     E          10     E      At each site each clone is replicated multiple times. Ideally, I would want to know what is the effect of clone on my dependent variable y and whether a site effect is present. To me it looks like an **incomplete** block design with clone as a fixed effect and site as a random effect (and block). Using lmer from the lme4 package in R, I would specify the model as follows:               lmer(y~clone + (1|site), data=mydata)          Is this a correct way of analyzing this dataset? I could also average by clone over sites and eliminate sites. But this way I will lose potential important information as to whether a site effect is present. Any pointers are appreciated!


title: Econometrics - choosing the best model when removing variables

text: So I am looking into a regression model that is supposed to predict the value of a house based on numerous independent variables. What I don't quite understand is how to select the "best" model when eliminating insignificant variables. The original model contained four independant variables of which three turned out to be insignificant (p-value > 0.05). So I removed, say X1, which had the largest p-value (small t-stat). The new model containing three variables clearly appeared to be better, as Adjusted R2 rose (I understand Adjusted R Square is suitable for comparing models with different number of parameters). The F-tests also seemed to confirm that removing X1 was the right thing to do due to its insignificance. However, in the three variable model there was still one insignificant variable, X3. So I removed it as well. Now, in the two variable model, both variables were reliable (small p-values), but Adjusted R Square was actually somewhat lower than that of the previous model. F-tests seemed to indicate that X3 had a small enough impact on the model to be eliminated. Which is actually the better model statistically and economically? The three variable model with the highest Adjusted R Square or the two variable model with no insignificant variables but lower Adjusted R2? I couldn't quite figure out the logic behind this.


title: Expectation of Quotient of Sums of IID Random Variables (Cambridge University worksheet)

text: I'm preparing for an interview which requires a decent knowledge of basic probability (at least to get through the interview itself). I'm working through the sheet below from my student days as revision. It's mostly been fairly straightforward, but I am completely stumped on question 12. http://www.trin.cam.ac.uk/dpk10/IA/exsheet2.pdf Any help would be appreciated. Edit: the question is: Suppose that $X_1, X_2, ... $ are independent identically distributed positive random variables with $\mathbb{E}(X_1) = \mu < \infty$ and $\mathbb{E}(X_1^{-1}) < \infty$. Let $S_n = \sum_{i=1}^n X_i$. Show that $\mathbb{E}(S_m/S_n) = m/n$ when $m<=n$, and $\mathbb{E}(S_m/S_n) = 1 + (m-n)\mu\mathbb{E}(S_n^{-1}))$ when $m>=n$. In fact, in the process of typing this up, I've solved the second part. For $m>=n$, $\mathbb{E}(S_m/S_n) = \mathbb{E}(X_1+ . . . +X_m)/\mathbb{E}(X_1+ . . . +X_n)$ $=\mathbb{E}(1 + (X_{n+1} + ... + X_m)/(X_1 + ... + X_n)) $ and the numerator and denominator of the ratio above are clearly independent, so: $ = 1 + \mathbb{E}(X_{n+1} + ... + X_m)\mathbb{E}(S_n^{-1})$ and we obtain the desired result. I'm still stuck on the first part though.


title: Multivariate normal density function

text: I am trying to compute multivariate normal distributions at some points. I am using Matlab's `mvnpdf` function: `y = mvnpdf(X,MU,SIGMA)` The first argument is the point where I compute the density, MU is the mean and SIGMA the covariance. I am puzzled by the following result:                mvnpdf([0 0 0],[0 0 0],0.001*eye(3))          ans =             2.0078e+03          >> mvnpdf([0 0 0.002],[0 0 0],0.001*eye(3))          ans =             2.0038e+03      I am going at $2\sigma^2$ from the mean and the density is almost the same? Shouldn't the result be close to zero?


title: Order Statistics-Expected Value of Random Length

text: Let $Y_1<Y_2 $ denote the order statistics of a random sample of size 2 from a distribution that is $N\left( \mu,\sigma^2 \right) $, where $\sigma^2$ is known. Compute the expected value of the random length $Y_2-Y_1$. I can see that the answer is $\frac{2\sigma}{\sqrt{\pi}}$ but I do not know how to get there since I cannot evaluate the double integral: $$ \int_{-\infty}^{\infty} \int_{-\infty}^{y_2} \left( y_2-y_1 \right) \frac{1}{2\pi \sigma^2} exp \left\\{ -\frac{1}{2\sigma^2}\left[ \left( y_1-\mu \right)^2 +\left( y_2-\mu \right)^2 \right]\right\\} \mathrm {dy_1 dy_2}$$ Any ideas on how to compute this are greatly appreciated, thank you!


title: Similarity distance score to remove outliers for survey data

text: I'm still a beginner at data mining. I'm working on finding the association rules from hypothesis X to conclusion Y. To this end, I've conducted a survey with questions that go something like this:               Q1: Do you have any relatives with ABC property? (Ans: yes or no)     Q2: Do you have an interest in XYZ field? (Ans: yes or no)     Q3: Which institute are you from? (Ans: Option 1, 2, 3, 4 or 5)      and so on.. So there are lots of "parameters" or "features" or "dimensions" to my data. My data is formatted similar to this: http://www.hakank.org/weka/weather.arff, and I'll be using WEKA. However I'm still currently in the data pre-processing stage. Removing duplicate entries and dealing with missing values is no issue. What I'm worried about is removing outliers. Firstly, how can I represent this type of record data in such a way that similarity measures like Euclidean or Minkowski (or perhaps any!) distance can even be applied to it? And secondly, what's the most reasonable similarity measure to use in this type of case? I've looked at the Mahalanobi distance and it seems useless for my project because I have no "ideal" set of features against which I could compare other sections of data. Is it usual to even need to detect outliers before finding association rules? Or are outliers usually detected after the rules have been learned? I've been thinking about this for a while but can't seem to reach a sensible conclusion. Could any of the more experienced data miners help please?


title: Which toolbox for Belief Propagation and other inference methods in graphical models?

text: Which open-source software (toolbox) do you think is the best for modelling graphical models (e.g. factor graphs), and doing inference on them? (the language doesn't matter)