title: Nearest neighbor information for recommendation engines

text: My friend and I are working on a project on distributed datastructures. We were wondering how much is nearest neighbor information used in modern recommendation systems and whether it would be worthwhile to work on a distributed datastructure (say a kd-tree) for that purpose. Thanks


title: How to find outliers in a data series?

text: I have a series of 100 points My dataset can be found here . Each row is a data series. The plot for 90th row is ![enter image description here](http://i.stack.imgur.com/yLkq5.png) It's easy to detect outliers visually by plotting example. I tried using hampel to find outliers assuming it as time series.               x <- read.table("anomaly_s57.dat")     data <- as.matrix(x)          plot_hampel = function(row, k = 2, t = 3) {              plot.ts(data[row,])         hp <- hampel(data[row,], k , t)         y <- data[row,hp$ind]             x <- hp$ind         points(x,y,col="red")          }       But it is not good enough. It misses some small peaks. Is the time series assumption not correct? Any statistical fitting is possible? How to detect outliers in this data series considering each row of the data set as independent data series? It is known that total number of outliers is around 500.


title: I want to build a crime index and political instability index based in news stories

text: I have this side project where I crawl the local news websites in my country and want to build a crime index and political instability index. I have already covered the information retrieval part of the project. My plan is to do:   * Unsupervised topic extraction.   * Near duplicates detection.   * Supervised classification and incident level (crime/political - high/medium/low). I will use python and sklearn and have already research the algorithms that I can use for those tasks. I think 2. could give me a relevancy factor of a story: the more news papers publish about an story or topic the more relevant for that day. My next step is to build the monthly, weekly and daily index (nation-wide and per cities) based on the features that I have, and I'm a little lost here as the "instability sensitivity" might increase to the time. I mean, the index from the major instability incident of the last year could be less than the index for this year. Also if to use fixed scale 0-100 or not. Later I would like to be able to predict incidents based on this, e.g. whether the succession of events in the last weeks are leading to a major incident. But for now I will be happy with getting the classification working and building the index model. I would appreciate any pointer to a paper, relevant readings or thoughts. Thanks. PD: Sorry if the question does not belong here. **UPDATE** : I haven't yet "make it", but recently there was a news about a group of scientists that are working in a system to predict the events using news archives and released a relevant paper Mining the Web to Predict Future Events (PDF).


title: Logarithmic scale on a plot with negative values

text: I would like to plot two time-series on a same graph. One series takes much larger values than the other, so I thought a semilog scale might be appropriate (i.e. linear X (dates) and log Y). However, both series take on negative and positive values. Does it still make sense to use a log scale? If so, should I transform both series as follows?                   if observation > 0             log_observation = log(observation)         elseif observation < 0             log_observation = -log(-observation)         else             log_observation = 0


title: Nonlinear estimates of regression coefficients?

text: About simple regression: It's well known that usual OLS estimators of $β_{0}$ and $β_{1}$ have minimum variance of all unbiased linear estimators. I wonder if there are popular biased linear estimators which have smaller variance or unbiased nonlinear estimators or biased nonlinear estimators which have better characteristics such as smaller variance. Recenlty, I found these!!   1. related question : OLS is BLUE. But what if I don't care about unbiasedness and linearity?   2. related comment : " It is probably the same narrow place of understanding that OLS is just a best linear unbiased estimator (BLUE), and there exist better biased or nonlinear examples like James-Stein estimator or LASSO estimator that minimizes MSE further." from Intuition behind why Stein's paradox only applies in dimensions $\ge 3$


title: Using Python for building machine learning application

text: I'm currently using R to find the best approach to solving a machine learning problem. Once I've got the approach sorted, I will need to build this into an application which can be used by end users. My background is as a .NET developer. I see there are a few questions related to this, but my question is more about what I should use to build an end user product which incorporates machine learning. From what I've seen so far, R is very powerful but does not integrate very well with other programming languages (and even less well with .NET). So I'm trying to figure out the best approach for building the app. I see that Python is widely used by the ML community. Is this a good choice for building an app which will be delivered to users, or is it better as a scripting tool for prototyping etc? One benefit I can see is the range of machine learning libraries available, whereas .NET does not have a large range of libraries available. Performance concerns me given that it is interpreted. Is Python my best choice, or would it be better to build the algorithms I need from scratch in C++, C# etc.?


title: Effect of corrections for multiple testing on sensitivity and specificity

text: Let's imagine we have one variable (factor.to.explain) that we want to explain by 10 other variables using 10 linear model (no interactions computed). We should correct for multiple testing. It is possible that one, some or several of them affects the factor.to.explain. Note: TP = True Positive, FN = False Negative, FP = False Positive, TN = True Negative   * How does a Bonferroni correction affects sensitivity (= TP / TP + FN) and specificity (= TN / TN + FP) ?   * How does other corrections for multiple comparisons affects sensitivity and specificity ?   * Do sensitivity and specificity change if I do only 5 or all (10) my linear models ?   * Is there a objective way of deciding what is the best sensitivity to specificity ratio ?* Specific question (not important if I don't get an answer):   * What sensitivity and specificity ratio is usually aimed in ecology ?


title: Clustering time series based on correlation

text: I want to cluster a set of time series regarding their pairwise correlation. If I normalize the series by subtracting their average value and then scaling to a standard deviation of one, the correlation coefficient (Pearson's r value) between original series is the same as the dot product of the normalized counterparts. Which clustering algorithms would be apropriate in this situation ?


title: Ways to find a confidence interval for robust and quantile regressions

text: I'm trying to compare a few regression models for my data. For linear regression everything is quite understandable, but robust and quantile regressions are not so obvious. I could not find almost anything about calculating confidence interval for these regression models unless I looked for something wrong. This is my code in R               # Robust linear modeling     library(MASS)     library(robustbase)     library(robust)     set.seed(343);      x <- rnorm(1000)     y <- x + 2*rnorm(1000)          lm1<-lm(y~x); rlm1<-rlm(y~x); rlm2 <- lmRob(y~x); rlm3 <- lmrob(y~x)     cbind(summary(lm1)$coeff,  confint(lm1))      cbind(summary(rlm1)$coeff, confint(rlm1))     cbind(summary(rlm2)$coeff, confint(rlm2))     cbind(summary(rlm3)$coeff, confint(rlm3))      This code produces the following result:               > cbind(summary(lm1)$coeff,  confint(lm1))                        Estimate Std. Error   t value     Pr(>|t|)      2.5 %     97.5 %     (Intercept) -0.06973191 0.06408983 -1.088034 2.768429e-01 -0.1954982 0.05603438     x            0.97647196 0.06619635 14.751145 1.071805e-44  0.8465720 1.10637196     > cbind(summary(rlm1)$coeff, confint(rlm1))                       Value Std. Error    t value 2.5 % 97.5 %     (Intercept) -0.06131788 0.06714405 -0.9132288    NA     NA     x            0.96016596 0.06935096 13.8450275    NA     NA     > cbind(summary(rlm2)$coeff, confint(rlm2))         Error in UseMethod("vcov") :            no applicable method for 'vcov' applied to an object of class "lmRob"     > cbind(summary(rlm3)$coeff, confint(rlm3))                   Estimate Std. Error    t value     Pr(>|t|)      2.5 %     97.5 %     (Intercept) -0.0568964 0.06608987 -0.8608945 3.895029e-01 -0.1865874 0.07279464     x            0.9612520 0.06821558 14.0913850 2.921913e-41  0.8273896 1.09511448      It's easy to spot that linear model works OK and only one robust regression gives a sensible result. Another observation is that lmrob(), which produces some actual confidence interval, calculates it in the same manner as lm(), with using 1.96 as the student coefficient. Is it a correct way to produce a confidence interval for the robust regression model? May the same method be used for the quantile regression model?


title: Hilbert curves: bounds / probabilities of preserving neighbors

text: Hilbert Curves (Wikipedia) are space-filling curves said to "fairly well preserve locality". Do you know any _theoretical_ results here, such as bounds that neighbors within a radius of $\varepsilon$ are preserved with probability $p$ or anything like that? Wikipedia also only has a "citation needed" flag there. I'd like to be able to give probabilistic bounds on the error when using Hilbert curves, for nearest neighbor search in $d$ dimensions on clustered data (i.e. not uniform on the data space). I believe I saw some results for $d=2$, but I'm interested in the more general case. Thank you for references to theoretical work on the quality of the hilbert curve.


title: Which model for panel data with dependent variables from [0,1]?

text: I'm stuck with a regression modeling problem. I have panel data where the dependent variable is a probability. Below is an excerpt from my data. The complete panel covers more countries and years, however it is unbalanced. What I can observe is the number of events and the number of trials. The event probability was derived from those values (estimation of this probability should be quite good, given the large number of trials). All independent variables are county-year specific.                    country  year  event_prob  events trials    x    x_lag2 ... more variables       1   Cyprus  2008  0.03902140  11342  290661   4.60   4.13  ...       2   Cyprus  2009  0.04586650  13482  293940   4.60   4.48  ...       3   Cyprus  2010  0.05188398  15206  293077   4.60   4.60  ...       4   Cyprus  2011  0.06433411  18505  287639   5.79   4.60  ...       5  Estonia  2008  0.07872978  21686  275449   6.02   4.11  ...       6  Estonia  2009  0.09516270  33599  353069  13.18   4.91  ...       7  Estonia  2010  0.08645905  36180  418464   7.95   6.03  ...       8  Estonia  2011  0.07731997  31590  408562   5.53  13.18  ...       ...     165  USA  2011  0.06100000  9192822  150702000   2.73  3.27  ...      My goal is to use regression analysis to find out which variables are significant for the event probability. In R-terminology, I'm looking for a model of the form `event_prob ~ x + x_lag2 + ...` . The problem is as follows: `event_prob` has to be between 0 and 1, hence using `event_prob ~ x + x_lag2 + ...` might not be the best idea. So I was thinking of using the logit transform of `event_prob` such that `logit(event_prob)` ranges from $-\infty$ to $\infty$. The first idea was to use the R's `plm` package, i.e. `plm(logit(event_prob)~x+x_lag2,data,index=c("country","year"),model="random")` or `model="within"` (see below). Is that a reasonable approach or am I violating some essential assumptions? I was also thinking of using panel generalized linear models from the package `pglm` (with the logit link function), however since I don't know the outcome of the binary events (only the total number of events and trials) is known, I got stuck there. Maybe someone can help me how to proceed here. Since I have panel data, I'd like to compute both fixed-effects models and random-effects model and then apply the Hausman (1978) test to decide which model is more appropriate. Do my first attempts at modeling make sense? I'm really not sure how to correctly address this problem. I hope the description of my problem is detailed enough. If not, I'm happy to provide more details In terms of software, I'd prefer R. SAS and SPSS are also ok since my university has licences for them. I just don't have much experience with them.


title: Implementing Logistic Regression in MATLAB

text: I have a data set of 13 attributes where some are categorical and some are continuous (can be converted to categorical). I need to use logistic regression to create a model that predicts the responses of a row and find the prediction's accuracy, sensitivity, and specificity.   * Can/Should I use cross validation to divide my data set and get the results?   * Is there any code sample on how to go about doing this? (I'm new to all of this)   * Should I be using mnrfit/mnrval or glmfit/glmval? What's the difference and how do I choose?


title: What is the difference between frequent pattern mining, subspace clustering, and biclustering?

text: I have seen the terms "frequent pattern mining", "subspace clustering", and "biclustering". They all pertain to finding clusters using subsets of the data attributes. What's the difference?


title: Estimating an overall summary statistic (mean) for a value reported by ten published studies

text: ### Background I am writing a systematic review and meta-analysis of the association of exposure to X with the outcome Y. I have identified ten studies that report on this subject. This will be a meta-analysis of the published literature, not of individual participant data. I want to report the overall mean age and standard deviation of all patients in all of the studies together. The studies all have different sample sizes. How should I do this? ### Possibilities (1) I imagine that I cannot just add up the reported mean values and divide by ten. This would give smaller studies an undue weight in contributing to the overall mean. (2) I thought perhaps that I could just do a simple weighted average, using the sample size of each study as the weight. Is this appropriate? Can I also just do a weighted average of the reported standard deviations and ranges? If not, how would I compute those? (3) Another thought I had was actually taking each of the ten reported mean weights along with their standard deviations and actually plugging them into a meta-analysis program to get a "pooled weight" with a fixed effects model. But perhaps this would be a strange approach? Any advice on which technique to use? Is there another possibility that I'm not thinking of? Thanks!


title: Sampling - fitting of sample

text: Does sample size influence the underlying distribution of population, i.e., the bigger the sample you generate the more it resembles whatever distribution of the population it's drawn from? What would you say are the key factors for measuring the quality of a given sample with regards to its fitting to the distribution of its population? I'm currently studying discrete random variables, eg. the geometric distribution. Thank you.


title: Using the Cox axioms to derive unknown probabilities from known probabilities

text: To strengthen my understanding of fundamental probability theory, I am working my way through Professor Aaron Hertzmann's _Introduction to Bayesian Learning_ course notes. Section 3.8 of these course notes includes the following exercise. ### Original Problem Derive a formula for $P(\mathbf{A})$, assuming you know $P(\mathbf{A}|\mathbf{B}_1,\mathbf{C})$ and $P(\mathbf{A}|\mathbf{B}_2,\mathbf{C})$ and $P(\mathbf{B}_1|\mathbf{C})+P(\mathbf{B}_2|\mathbf{C})=1$. ### My Attempt So, my goal is to express $P(\mathbf{A})$ in terms of known quantities. $$ \begin{align} P(\mathbf{B}_1|\mathbf{C}) &+ P(\mathbf{B}_2|\mathbf{C}) &=& 1 && \text{}\\ P(\mathbf{B}_1|\mathbf{C})P(\mathbf{C}) &+ P(\mathbf{B}_2|\mathbf{C})P(\mathbf{C}) &=& P(\mathbf{C}) && \text{Algebra}\\ P(\mathbf{B}_1,\mathbf{C}) &+ P(\mathbf{B}_2,\mathbf{C}) &=& P(\mathbf{C}) && \text{Product Rule}\\ P(\mathbf{B}_1,\mathbf{C}|\mathbf{A}) &+ P(\mathbf{B}_2,\mathbf{C}|\mathbf{A}) &=& P(\mathbf{C}|\mathbf{A}) && \text{Condition on } \mathbf{A}\\ \frac{P(\mathbf{A}|\mathbf{B}_1,\mathbf{C})P(\mathbf{B}_1,\mathbf{C})}{P(\mathbf{A})} &+ \frac{P(\mathbf{A}|\mathbf{B}_2,\mathbf{C})P(\mathbf{B}_2,\mathbf{C})}{P(\mathbf{A})} &=& P(\mathbf{C}|\mathbf{A}) && \text{Bayes Theorem} \\ \end{align} $$ $$ \begin{align} &\frac{P(\mathbf{A}|\mathbf{B}_1,\mathbf{C})P(\mathbf{B}_1,\mathbf{C}) + P(\mathbf{A}|\mathbf{B}_2,\mathbf{C})P(\mathbf{B}_2,\mathbf{C})}{P(\mathbf{A})} &=& P(\mathbf{C}|\mathbf{A}) && \text{Algebra}~~~~~~~~~~~~~~\\ &\frac{P(\mathbf{A}|\mathbf{B}_1,\mathbf{C})P(\mathbf{B}_1,\mathbf{C}) + P(\mathbf{A}|\mathbf{B}_2,\mathbf{C})P(\mathbf{B}_2,\mathbf{C})}{P(\mathbf{C}|\mathbf{A})} &=& P(\mathbf{A}) && \text{Algebra}~~~~~~~~~~~~~~\\ \end{align} $$ Here is where I get stuck, because I don't know $P(\mathbf{C}|\mathbf{A})$, $P(\mathbf{B}_1,\mathbf{C})$, or $P(\mathbf{B}_2,\mathbf{C})$. It is not clear to me how I can express these probabilities in terms of the quantities I know. These course notes are very much a work in progress. The author frequently writes _TO-DO_ notes to himself about things to add and remove. Is it possible that this exercise is actually under-determined? Or am I simply going about this derivation the wrong way? ### Similar Questions When searching StackExchange, I found this similar question. However, it was still not directly applicable to my problem. ### Tagging as Homework Although this question is not homework for me, I am tagging it as homework anyway because: (1) it is very introductory in nature; (2) a similar question could conceivably be homework for someone else; and (3) I would be happy with hints.


title: What is the empirical Bayes estimator for a gamma-Poisson with more than 1 observation for each Poisson parameter?

text: I am looking at the Wikipedia entry for empirical Bayes, but it's a bit confusing - it seems to me the solution must apply only to the case in which there's only $n=1$ sample $y$ for each $\theta$ and the "sample mean" that's referred to is really just a single value. If this were not the case, I would think that the shrinkage factors would depend on the sample size for each $\theta$. What if I want to generalize this to where there are multiple samples of $y$ for each Poisson mean $\theta$? Is this the formula for the posterior on $\theta_k$ (the Poisson mean for the $k^{\rm th}$ group): $$ \theta_k = \frac n {(n + \beta)} \times \text{sample_mean}(y_k) + \frac \beta {(n + \beta)} \times (\alpha \times \beta) $$ if so what's the derivation for this?


title: Taylor approximation to penalty function

text: I'm reading through Fan's SCAD paper and I feel like I'm not getting a simple step. On page 1354 where he is talking about quadratic approximations to a penalty function he has $$\left[\rho_\lambda(|\beta_j|)\right]' = \rho'_\lambda(|\beta_j|)\text{sgn}(\beta_j) \approx \left\{\rho'_\lambda(|\beta_{j0}|)/\beta_{j0}\right\}\beta_j$$ which I am comfortable with, but then he says $$\rho_\lambda(|\beta_j|) \approx \rho_\lambda(|\beta_{j0}|) + \dfrac{1}{2}\left\{\rho'_\lambda(|\beta_{j0}|)/\beta_{j0}\right\}(\beta_j^2-\beta_{j0}^2)$$ for $\beta_j \approx \beta_{j0}$. This looks pretty close to a 1st order taylor expansion about $\beta_{j0}$, but I'm not quite sure how the usual quadratic term (e.g. $(\beta_j - \beta_{j0})^2$) becomes the term above.


title: How to find implementation details for a SAS procedure?

text: I was was wondering about where to track down implementation details for SAS: I have certain outputs (depending on small amounts of data) which I would like to understand better by trying to check against a theoretical solution. I'm looking at something pretty specific (how standard errors on predicted probabilities from a regression are found), so I have no idea how to track this sort of thing down.


title: Identifying if text from 2 sources is about the same topic

text: I have many articles and they have headlines and body text. I am trying to deduce a way to identify with some (ideally)high degree of probability if 2 articles are about the same thing. One example of a probable match here would be two news articles from two separate publications discussing a natural disaster. Anyway, I currently am not sure where to start so any guidance would be greatly appreciated.


title: Guidance on classification analysis

text: I am doing a project where I would like to predict some characteristics from large data set and I am expecting that I would should use some Machine Learning techniques, but not sure how to proceed. I have a background in Econometrics, but I have not really done projects like this before. I have a large data set consisting of hour-based electricity consumption for around 200,000 households over a couple of years. For about 1,000 of these households, I know that they have a specific heating device that shows a specific pattern in the consumption (i.e. large variation from hour to hour, and through the seasons etc.). I would like to predict which of the remaining household’s consumption shows similar behavior with at certain probability. In other words, I would like to know which households uses the specific heating device. Which techniques are recommended for this type of analysis? I know that there probably isn’t any straight forward solution, but I would really like to know where I should focus my research.


title: Obtaining random number from a mixture of two normal distributions

text: I want to sample from mixed normal distribution, first one is $N(1,2)$, second one is $N(5,4)$. I used `rnorm(100, c(mean=c(1,5), sd=c(2,4)))`. Is this correct? The problem I am trying to solve is sampling from the 2 distribution above, first one with 75%, second one with 25%. Am I on the right track? **Edit:** I will rewrite the problem for clearance, with easier numbers. :) I want to sample from $N(0,1)$ with 70% probability, and $N(100,10)$ with 30% probability. Of course, that's just for sake of discussion, the actually distribution I am working with is n(21, 3.3), n(26,4).


title: left-censored dependent variables and prediction

text: I'm coding up a monte-carlo analysis; I've got a deterministic model that depends on parameters that are uncertain. One of those uncertain parameters is a partially-observed vector of prices by country -- I want to predict missing values in a manner similar to what is used in multiple imputation via chained equations. I've got a few candidates for modeling those missing values:   1. A linear model with censored predictions -- I simply censor negative predictions at zero. Average predictions will probably be biased upwards because of the left-censoring.   2. A linear model fit to a log-transformed price. Average predictions will be biased upwards because of the unnaturally fat right tail, and the probability of responses at the right tail will be unnaturally fat.   3. A glm with gaussian errors and a logarithmic link. I don't really have the intuition for this, but examining a few predictions gives me results that are even crazier than for #2.   4. A tobit model. This is the best choice, but how does one even think of a prediction interval in that case? How does one code it? Edit: when I say "prediction" I'm talking about a random draw from the fitted model's predictive distribution.


title: Do I need to drop variables that are collinear before running kmeans?

text: I am running kmeans to identify clusters of customers. I have approximately 100 variables to identify clusters. Each of these variables represent the % of spend by a customer on a category. So, if I have 100 categories, I have 100 variables such that sum of these variables is 100% for each customer. Now, these variables are correlated with each other. Do I need to drop some of these to remove colinearity before I run kmeans? Here is the sample data. In reality I have 100 variables and 10 million customers.               Customer CatA CatB CatC        1         10%  70%  20%        2         15%  60%  25%


title: How to install fGarch package

text: Can someone give me detailed instructions how to install and use fGarch? Eg. I would like to find the MLE parameter estimates for a skew normal dist using it - could someone give me instructions up to that point?


title: Confused about the proper way of testing for an improvement or worsening

text: I want to check if a new equipment improves a score compared with the older equipment. My database contains the score with the old equipment, and the scores with each new equipment. I googled, userR'd, checked similar journal articles and saw that one of them used the sign test to assess the improvement in the grading. # EDIT To make it clear, the journal article which studied a very similar question used the sign test. They and I have data like this one:                              New equipment     Old equipment  1  2  3                 1 33  0  0                 2 37  8  1                 3 12  4  0                 4  2  0  0      This table should be read like this:   * Numbers 1-4 are scores for some condition (1 is better, 4 is worst).   * The first row starting with 1 means that the old equipment provided a score of 1 in 33 patients, and the situation did not change with the new equipment.   * The second row starting with 2 means that the old equipment provided a score of 2 in 46 patients (derived from 37+8+1). The new equipment provided a score of 3 in 1 patient, a score of 2 in 8 patients, and a score of 1 in 37 patients. This is a good improvement, but in ~2% of cases, things may get worse.   * The third row starting with 3 means that the old equipment provided a score of 3 in 16 patients (derived from 12+4). The new equipment provided a score of 2 in 4 patients, and a score of 1 in 12 patients. This is a major improvement.   * The fourth row starting with 4 means that the old equipment provided a score of 4 in 2 patients. The new equipment provided a score of 1 in both of these patients. This is a major improvement. However, only 2 out of 97 patients scored 4 with no instrument. Therefore we may conclude that only patients who score high (3 or 4) with no instrument may benefit from using an instrument.  **I think that the best option is to test each row separately with wilcoxon test, by comparing the medians of both conditions (old equipment and new equipment).** As an example:               score.2 <- subset(original.df, old-equipment %in% 2)     wilcox.test(score.2$old-equipment, score.2$new-equipment, alternative="great")      and repeat this for each score with the old equipment. I couldn't see a way to use the sign test for this. What is your suggestion to test this condition properly? # EDIT#2 I think the Wilcoxon signed rank test is appropriate to compare the differences of two new methods/instruments over a common old method/instrument. But when it comes to test the significance of the difference of one method over the other, it has a shortcoming. Let me explain with an example: Lets say, on 4 trials, I could eat 3,4,4,5 peas with a fork. With a spoon, I could eat 4,5,5,6 peas. Since the median of the spoon-trial is higher than the median of the fork-trial, we can conclude that spoon is better for eating peas. It won't differ if the spoon-trial ends up as 6,5,5,4 or 5,5,5,5 or even 1,1,10,10. The spoon-median is higher than the fork-median and the Wilcoxon test will yield a low p-value. But when it comes to grading or notes, or anything that requires consistency, it differs. In my case, I'm working with old and new scoring systems. When a new method/instrument consistently achieves higher scores compared with the old one, it is simply better than the old one. But when it achieves very good scores in half of the trials, and very bad scores in the other half (but still better sum or median than the old one), it is not better. It may score better in some circumstances but it is not reliable. That is the weakness of the Wilcoxon test. It does not take the paired nature of the data into account. So my exact question is, how can I make the Wilcoxon test see this problem from a "paired samples" aspect. So I have data, where I test each subject with three different instruments and would like to know if it is worth the effort to switch to any of the new instruments.               #dput     structure(list(Mallampati = structure(c(2L, 3L, 3L, 3L, 2L, 2L,      2L, 1L, 3L, 2L, 4L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 1L,      2L, 1L, 2L, 1L, 2L, 2L, 3L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L,      1L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), VL = c(1,      1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), DL = c(1,      1, 3, 2, 1, 2, 1, 1, 2, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 1,      1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2)), .Names = c("Mallampati",      "VL", "DL"), row.names = c(NA, 40L), class = "data.frame")


title: Selecting most realistic C and g params after gridsearch

text: I just ran an extended SVC gridsearch in libsvm on about 9000 multi- dimensional vectors representing a time series. Here are the highest scoring results:               [local] 3 -7 72.4729 (best c=0.5, g=0.5, rate=76.9618)     ..     [local] -1 -5 71.79 (best c=8.0, g=0.5, rate=77.4432)     ..     [local] 15 -11 73.0326 (best c=2048.0, g=0.03125, rate=77.5887)     ..     [local] -5 -3 66.0249 (best c=32768.0, g=0.03125, rate=77.723)     ..     [local] 7 -3 77.9245 (best c=128.0, g=0.125, rate=77.9245)      Here's the contour of the cross validation accuracy: ![enter image description here](http://i.stack.imgur.com/5JiOg.png) Given the above I'm very tempted to settle on either a C value of 0.5 or 8.0 - anything above 100 seems heavily over fitted to me. Especially given the rather minuscule success rate differences. Any thoughts, insights, or pointers would be appreciated. Thanks in advance.


title: ML estimation of parameters that do not completely specify the model

text: I was wondering how ML is defined when the parameter does not completely specify the model. More concretely, suppose $X_1, X_2, \cdots, X_n$ are drawn iid such that $P(X_1=i)=\theta_i$, $ 1 \leq i \leq k$. I want to find the ML estimate of $\phi= \max_{1 \leq i \leq k} \theta_i$. To me it is not even clear if the ML estimate of $\phi$ is well defined in this case. Am I missing something here?


title: How is the "classification error rate" of an artificial neural network calculated?

text: Frequently I see artificial neural networks compared by their "classification error rates" or "error rates", particularly for multi-class problems like CIFAR-10. What does this error rate actually refer to? Hamming loss? How is it calculated?


title: Variance components estimation in Poisson mixed models

text: In the mixed Poisson regression model, with vector of random effects $w \sim N(0, \Sigma)$, how are the parameters in $\Sigma$ estimated? Is it the same REML method as in the linear mixed model?


title: Connection between power law and Zipf's law

text: I am trying to better understand the connection between the power law distribution and Zipf's distribution (law). There is a neat explanation in [1]. The article suggests that as we can derivate the power law function from Pareto's law, combined with the relationship between Pareto's law and Zipf's law, the power law parameter alpha is 1 + 1/b. From my understand, this would mean that we can directly determine Zipf's law's b parameter by simply having the power law alpha parameter. So e.g., an alpha of 2 would lead to a b of 1. Is this true? So can I calculate alpha by e.g., using the methods by Clauset et al. [2] on my data and then directly determine the Zipf parameter b by the definition? This would allow me to use the exact methods by Clauset instead of the non-exact methods like fitting a straight line on a log-log plot. So I would also overcome the necessity of producing rankings etc. [1] http://www.hpl.hp.com/research/idl/papers/ranking/ranking.html [2] http://arxiv.org/abs/0706.1062


title: Statistical significance for comparison of linear regression models

text: I have two linear regression models (with the same predictors) that try to estimate two different (although related) features of the same population. I am analyzing the hypothesis that these predictors are not as good for the second feature as they are for the first. Indeed, the RMSE of the second model is 7% worse than the RMSE of the first: it is not a huge difference, but for my purposes it is fine. The problem is: I have been asked the statistical significance of this result (the 7% worse). I am a computer scientist, not a statistics expert, and I am at a loss here on which test I should use, given that I am trying to compare the squared errors of the two (paired) sets of predictions (I am using SPSS software).


title: Non parametric one sample t-test alternative with a binary variable

text: I am doing an analysis of items (I1, I2, I3, etc.). The items could be either correctly answered (1) or incorrectly answered (0). Visually, most of the participants answered the items correctly. I want to know whether each of the items was different from 1. I thought to run a Wilcoxon signed rank test for each of the items, but then I realized that I only can do this with scale variables... What can I do instead?   Any help would be appreciated. Thanks


title: Want to make a function which allows for recursive window forecasting

text: I have been looking for a function that can make recursive window out-of- sample forecasts, but seems there is none. So I'm thinking about about making a function that can be used for recursive window forecasting in an ARIMA model. However I know little about programming, so I'm seeking for help. What I want to do is use the function `forecast.Arima` ( **forecast** package) to predict future values in a expanding window. Suppose 20 years is the initial window, and I expand the window by 1 year on each iteration until it is of size 30 years. More specifically, use 20 years data to predict one value, use 21 years data to predict the next value, etc.


title: Can MCMC iterations after burn in be used for density estimation?

text: After burn-in, can we directly use the MCMC iterations for density estimation, such as by plotting a histogram, or kernel density estimation? My concern is that the MCMC iterations are not necessarily independent, although they are at most identically distributed. What if we further apply thinning to the MCMC iterations? My concern is that the MCMC iterations are at most uncorrelated, and not yet independent. The ground I learned for using an empirical distribution function as an estimation of the true distribution function is based on  Glivenko–Cantelli theorem, where the empirical distribution function is calculated based on an iid sample. I seemed to see some grounds (asymptotic results?) for using histograms, or kernel density estimates as density estimations, but I can't recall them.


title: Regression techniques similar to Kriging/Gaussian process regression

text: I am looking for regression techniques which are similar to Kriging/Gaussian process regression, in that no explicit model needs to be specified. (Discounting the prior over functions) I have three independent variables and one dependent variable to which I want to apply such a procedure. The independent variables specify coordinates (locations in 3D), while the dependent variable specifies Wi-Fi signal strength at the given coordinates. Since it is hard to appropriately visualize such high dimensional data, techniques without explicit model dependence are of primary interest. The only similar technique I found was the somewhat unsophisticated Nearest-neighbour interpolation. Are the above-mentioned techniques the only choices for such a problem?


title: Two different methods....different results?

text: I have a response variable and 10 predictor variables (all ordinal). I wanted to see if there was any evidence of a relationship between the response and predictors. I used a two proportion z-test to see if there was a significant difference in my predictor at the two different levels of my response. I did this for each predictor variable. At the end, my results were significant for each comparison. I also used Binary logistic regression where my response is dependent and I put in all my predictor variables under "covariates" (because that's the only option in SPSS). Based on the results, only 6 of my 10 variables were significant. I know that logistic regression is for modeling and prediction, but do these results contradict each other? Is one method "stronger" than the other? Or maybe I'm not looking at this from a correct perspective. Any input is appreciated, thank you.


title: Equi-probable sampling in R using the prob argument in sample()?

text: If I have data of length `n` and want to generate a random sample of length `N`, does the following use of the `prob` argument make each observation equi- probable?               random.sample = sample(mydata, N, replace=TRUE, prob=rep(1/n, times=n))        As an example, the following density plot shows the density of a sample generated with (red) and without (black) the above `prob` argument: ![enter image description here](http://i.stack.imgur.com/duWRP.png) The overall shape of the curve and the position of peaks and troughs has not changed so how does `prob` influence the result? Following from the above train of thought:   1. Is it even possible, statistically speaking, to make such an "equi-probable" sample? If yes, how?   2. More generally, how does the `prob` argument work and in which cases is it used in random sampling?


title: Suggestions for clustering ordinal, non-normal data (unsupervised)

text: I have data from about 250 self-injurers and would like to cluster analyze their reported motivations for self injury (and then explore cluster membership vis-a-vis various psychological measures). I have scores on 9 motivation scales which are likert-type but non-normal ($\log_{10}$-transformations yield <1 skewness for maybe 5 of the scales but 3 remain greater than 2 --as many participants only endorse one or two motivation types and give 0's on the other scales). I understand that some algorithms are more robust against normality violations and I have been reading about possible routines to use in `R` including DBSCAN, MClust, and EMMIX but am struggling to wrap my mind around a best approach. Density vs. similarity approaches etc. I understand that DBSCAN struggles as dimensions increase: are 9 variables a concern for this approach?. I can live with losing cases as 'noise' and even non-exclusive cluster membership (though am less excited about that outcome). Any thoughts/references are most appreciated!


title: Calculating covariance for a non-strictly-stationary white noise process

text: I'm reading Hayashi's "Econometrics", and in Chapter 2, page 101 he discusses the following white noise process: ![enter image description here](http://i.stack.imgur.com/BncMb.png) I understand the calculations for the expected mean and variance, but I can't understand why the covariance would be 0 for $i\neq j$. Since $E(z_i)=0$, the expression for the covariance would be reduced to $E(\cos(iw)\cdotp \cos(jw))$. Would someone mind showing me why this goes to zero? Assuming for the moment that $cos(iw)$ is positive, then $\cos((i+1)w)$ would be negative and so on. But why would there be 0 covariance for any two $i$ and $j$, $i\neq j$? I know this is a simple question, but can't seem to wrap my head around it.


title: Feature Selection Limitations - Machine Learning FPR Test

text: I am trying to improve accuracy of my logistic regression model by selecting the best features. I did an FPR test and ranked the features based on their F-score. The problem is that selecting the best say 3 features doesn't perform as well due to what I suspect is a cross correlation between these features. I noticed this after I saw better performance by picking the best 1 or 2 in combination with ones ranked about 9 or 10. I am a bit of a newbie, but is there a way I can systematically find what the best combination is?


title: Propensity score matching with multiple treatments

text: Is anyone aware of propensity score matching methods for when there are more than 2 treatment groups? I am working on a project with 4 treatment groups:   1. A   2. B   3. A and B   4. Neither A nor B Calculating propensity scores using multinomial logistic regression might work, but then I'd get multiple scores for each observation so I'm not sure how I'd match/analyze the matched data. Thanks!


title: Binomial maximum likelihood expressed as a sum

text: If I have a binomial maximum likelihood function $$ L(\theta | k) = \Pi_i^N p(\theta)^{k_i}(1-p(\theta))^{(n_i-k_i)}$$ I know that the index terms can be made into a sum, simplifying the expression for calculation purposes. I have a model which predicts the probability of success as a function of a variable (x) and a parameter. I have measured k (knowing n) for a number of values of x and now I want to find the best fit value for my parameter which I am trying to determine experimentally. Am I right in thinking my likelihood function is in this case $$ L(\theta | x, k) = \Pi_i^N p(\theta, x_i)^{k_i}(1-p(\theta, x_i))^{(n_i- k_i)}$$ Assuming this is the case (if not can someone correct me), is there anyway I can convert my product to a sum, or, do I have to keep it in product form?


title: likelihood of a model

text: the likelihood of a model is defined as the probability of data given model: > Likelihood(Model) = p(DATApoints | Model) which is equivalent to the product of all p(datapoint | Model) for each datapoint in DATApoints. so, to get the likelihood of a model given normally distributed data, i run the following [r] function: > prod(dnorm(DATApoints, mean=ModelMeans, sd=ModelStDevs)) my question is in regards to the DATApoints: are they central tendencies for each dependent variable (DV) of interest? or are there _n_ datapoints for each DV given _n_ participants in a study? To make this concrete: Assuming that i ran an experiment and found the following Dependent Variable (DV) means across my results: DV 1, 2, 3, 4, 5 = .12, .23, .34, .45, .56 and then i ran a model to try and account for these results, finding the following central tendencies (averaged across all model runs): DV 1, 2, 3, 4, 5 = .1, .2, .3, .4, .5 with standard deviations for each DV being .1 across all model runs, would the likelihood of my model be 644.344? calculated as follows: > Model=c(.1, .2, .3, .4, .5) > > DATA=c(.12, .23, .34, .45, .56) > > prod(dnorm(DATA, mean=Model, sd=.1)) > >> [1] 644.344 does this seem right?


title: Link function for log-logistic shared gamma frailty model

text: I've been asked to replicate a study that models an accelerated failure time survival model with a log-logistic distribution and gamma distributed frailty (a 'log-logistic shared gamma frailty model') estimated with the streg command in Stata [1]. We have made some changes to the original model and now want to estimate a quantity of interest, namely expected values, to see if the changes made result in substantially different results. I would do this using the Clarify command in Stata, but unfortunately it does not support this specific model (which the authors of the study I am replicating also acknowledge). Since I am more familiar with R, my idea was first to replicate the findings there. However, the parfm package, which is roughly analogue to Stata's streg, does not estimate the accelerated failure time model with a log-logistic survival function and a gamma frailty. I have therefore, somewhat inefficiently, extracted the information I need from the estimated model in Stata (the VCV-matrix and the estimated parameters), to estimate expected values in R according to Gary King's simulation approach [2]. To the specific issue; I believe I have all components in place except for the link function. The Stata manual for streg provides the stochastic distribution, but not the link function needed to calculate the expected survival time [3]. According to the streg manual, the log-logistic survival function has the following form: $$ S(t) = \{1 + (\lambda_{j} t_{j})^{1/\gamma}\}^{-1} $$ where $ \lambda_{j}=exp(-\textbf{XB})$ and $\log t = \textbf{XB}$, which implies $t=exp(\textbf{XB})$. I might be mistaken, but I believe that the frailty is multiplicative and it should therefore be possible to simply multiply the the survival function with the frailty parameter $\alpha_i$. Thus, my question simply concerns how I should specify the link function to get interpretable values from $y$ generated in the simulation. I apologize if my question is unclear, but I have no prior experience of working with survival models before being asked to replicate this. Any input is greatly appreciated, [1] The Stata syntax reads: `stset time.var, failure(fail.var) streg covariates, dist(llo) frailty(gamma) shared(cluster.var)` [2] King, G., Tomz, M., & Wittenberg, J. (2000). Making the most of statistical analyses: Improving interpretation and presentation. American journal of political science, 44(2), 347-361. [3] [http://www.stata.com/manuals13/ststreg.pdf] See page 9 and 25 respectively.


title: Finding interaction between several factors and multiple dependent variables

text: Asking a question essentially for a colleague, since it's far enough afield of what I do that I want to make sure I don't lead her astray: Assume a study with 5 factors and 3 dependent variables - we'd like to test the effect of interaction between any two factors on a given dependent variable (for Factor A at some level and Factor B at some level increases the mean of dependent variable A). Should this be handled with two-way ANOVA or MANOVA? My instinct is to say MANOVA unless there's reason to believe there's no correlation between any of the dependent variables, but I figured it's best to run it by folks who do that more than I do.


title: Comparing two F-test statistics

text: I have two groups A and B, each of which consists of 5 samples. Each sample is described in a vector of length (>1000) of continuous numeric values (characteristics). I want to test if the sample in the first group varies between each other more than the sample in the second group. I tried one-way ANOVA for each group independently and calculated the F-test statistic for each group. Now, can I compare the two obtained statistics directly (so if statistic 1 > statistic 2 I say the first group has more variance between its samples than the second group)? or do I need to perform statistical significance to compare the two statistics?


title: Confidence interval for Uniform($\theta$, $\theta + a$)

text: I am encountering a difficulty with the following task. Have I made a mistake, or is this an inherent flaw in the notion of confidence intervals? (Other such flaws exist.) Consider a random sample $X_1,\ldots , X_n$ from a Uniform($\theta$, $\theta + a$) distribution, where $\theta$ is unknown and $a$ is known. We wish to determine a confidence interval for $\theta$. The reader may verify the following details: The statistics $Y=\text{min}_i X_i$ and $Z=\text{max}_i X_i$ are jointly sufficient for $\theta$. For $\theta \le c_1 \le c_2 \le \theta + a$, $P\{c_1 \le Y \le Z \le c_2\} = [(c_2 - c_1)/a]^n$. For $0 < \gamma < 1$, set $d_1 =(1-\sqrt[n]\gamma)/2$ and $d_2 =(1+\sqrt[n]\gamma)/2$. Then $\gamma = P\{\theta + ad_1 \le Y \le Z \le \theta + ad_2\} = P\{Z -ad_2 \le \theta \le Y-ad_1\}$. Thus, $[Z -ad_2, Y -ad_1]$ is a $\gamma$ confidence interval for $\theta$. Now here's the difficulty: If we observe $Z - Y > a\sqrt[n]\gamma$, then $Z -ad_2 > Y -ad_1$, so our formula yields a nonsensical answer. Have I made an error in my calculations? Or is this one of those problems with confidence intervals? (Homework? I guess so - but a homework problem that I _wrote_ for my students. Inspired by another problem in DeGroot & Schervish.) [Also posted at math.stackexchange. I didn't know about this group, and I received no satisfactory answer there.]


title: Issues with fitting a variogram

text: I am trying to fit a spherical variogram to some synthetic data using the code available at http://www.mathworks.com/matlabcentral/fileexchange/25948-variogramfit. However, I have some doubts. I simulated some synthetic data using GMRF with the defined graph structure. I mean, I generated some samples where each sample is a multivariate gaussian distribution over space. I then took pairwise semi variances and got the overall semivariogram. I then fit the model variogram to this empirical one. This is what I got. However, I have some doubts, when I generated the synthetic data it had no noise. So I was assuming no nugget effect. However, I am seeing significant nugget effect why is it so. I just drew samples from a multivariate gaussian distribution over space. I have a spatial grid of size 10x10, the variables in this form a multivariate gaussian distribution. Any suggestions? I am not sure but the more samples I add, and calculate the pair wise semi variance, the nugget starts increasing. Even though the empirical semi variogram is smoother, the nugget starts appearing more. What is the reason behind this? ![enter image description here](http://i.stack.imgur.com/RbiNh.jpg)


title: What is the name of this chart, and is there an online tool to create one?

text: I would like to create a chart like the one on the left here in order to visualize maximum probabilities: ![enter image description here](http://i.stack.imgur.com/0GUgv.png) Is there a name for this type of chart, and is there a tool to create a similar one with my own data?


title: Is it important to scale data before clustering?

text: I found this tutorial, which suggests that you should run the scale function on features before clustering (I believe that it converts data to z-scores). I'm wondering whether that is necessary. I'm asking mostly because there's a nice elbow point when I don't scale the data, but it disappears when it's scaled. :)


title: Parameter estimation from a Normal distribution

text: Please can you check if am I correct? I have a random variable $X$ normally distributed with mean $\mu$ and variance $\sigma^2$. I generate two independent sample $T_1$ and $T_2$ with $T_1 < T_2$ where $\bar{X}_1$ and $\bar{X}_2$ are the sample means respectively. **1)The proposed estimators below for $\mu^2$ are biased or unbiased?**   * $\bar{X}_1^2$   * $\bar{X}_2^2$   * $(\bar{X}_1^2+\bar{X}_2^2)/2$   * $\bar{X}_1\times\bar{X}_2$ I think that in order to solve this question I have to use the property $Var(X) = E(X^2)-(E(X))^2$ and the fact that the sample mean $\bar{X} \sim (\mu,\sigma^2/T)$. Therefore, _First estimator $\bar{X}_1^2$,_ $E(X^2)= (E(X))^2 + Var(X) \rightarrow E(\bar{X_1}^2) = (E(\bar{X_1}))^2 + Var(\bar{X_1})$ which implies that $E(\bar{X_1}^2) = \mu^2 +\sigma^2/T_1\neq\mu^2 \Rightarrow$ **BIASED** _Second estimator $\bar{X}_2^2$, the same as the first one. **BIASED**._ _Third estimator $\frac{\bar{X}_1^2+\bar{X}_2^2}{2}$,_ using also information above $E\left(\frac{\bar{X}_1^2+\bar{X}_2^2}{2}\right)=E\left(\frac{\bar{X}_1^2}{2}+\frac{\bar{X}_1^2}{2}\right)=E\left(\frac{\bar{X}_1^2}{2}\right)+E\left(\frac{\bar{X}_1^2}{2}\right)=\frac{\mu^2+\sigma^2/T_1}{2}+\frac{\mu^2+\sigma^2/T_2}{2}= \mu^2 + \frac{\sigma^2}{2}\times\left({\frac{1}{T_1}+\frac{1}{T_2}}\right)$ **BIASED** _Forth estimator $\bar{X}_1\times\bar{X}_2$ ,_ $E(\bar{X}_1\times\bar{X}_2) = E(\bar{X}_1)\times E(\bar{X}_2)$ because the two samples are independent. Therefore, $E(\bar{X}_1\times\bar{X}_2) = \mu\times\mu=\mu^2$ **UNBIASED** **2) Now show that both the sample mean are unbiased estimator for $\mu$, that is, $E(\bar{X}_1) = \mu$ and $E(\bar{X}_2) = \mu$.** It's straight forward because the sample mean of a normal is normal distributed, $\bar{X} \sim (\mu,\sigma^2/T)$. **3) Now Define the follow estimator of $\mu$ and check if it biased or unbiased estimator of $\mu$, $\bar{X}_3 = \frac{(T_1\bar{X}_1+T_2\bar{X_2})}{T_1+T_2}$** So, $E\left(\frac{(T_1\bar{X}_1+T_2\bar{X_2})}{T_1+T_2}\right)$=$\frac{1}{T_1+T_2}\times E\left(T_1\bar{X}_1+T_2\bar{X}_2\right)=\frac{1}{T_1+T_2}\times T_1 E(\bar{X}_1)+T_2\times E(\bar{X}_2)=$ $\frac{1}{T_1+T_2}\times (T_1 +T_2)\times\mu=\mu$ **4) You known that an unbiased estimator is more efficient compared to an other unbiased if its variance is smaller. What of $\bar{X_1}$, $\bar{X_2}$, $\bar{X_3}$ has the smallest variance?** Of course it depends by sample size T, $\bar{X}\sim\left(\mu,\frac{\sigma^2}{T}\right)$ and in this case we have $T_3>T_2>T_1$ therefore $\bar{X}_3$ is the estimator with the smallest variance.


title: Confidence interval for values for a fitted line

text: I'm using JMP to analyze some sample data to make predictions about the population. My sample is from a destructive QC test, so I obviously want to minimize my sample. I have a response (my Y) and a known factor (a very strong and consistent correlation that is measurable by non-destructive means) but the exact relationship between them varies from lot to lot (the slope and y offset vary). So, in JMP, I am fitting a line and then showing the "confidence limits for an individual predicted value" which I believe gives me an indicator of how the population is likely to behave. So I'm using that plot to make disposition decisions. I want to automate this process, perhaps using R, but I'm a total novice at R. I could do the math if I was just dealing with a mean and standard deviation, but I don't know how to do it with a fit line and a known factor. Can someone please give me either the general information on how to get the confidence limits around the line, or else tell me how to do the whole thing in R? Thankss much.


title: Should I use decision trees to predict user preferences?

text: I'm designing a web service that will predict and recommend new items a user might like based on their expressed preferences on previous items (simple thumbs up/down interface). I was told to look into decision trees, since they're a simple way to group things based on their characteristics, and they can be easily used to classify new items. But after having done some reading and looked at a number of examples, I'm wondering if they're the best solution in this instance. Some problems I see:   * There are thousands of items in the database, which means there are potentially thousands of items on which to split the tree. Running all of these splits and then calculating which ones provide the most gain (loss of entropy) seems inefficient to me.   * The tree would be changing constantly. Most of the examples I've looked at build a tree with an initial sample set of data, and then run new items through it to classify them. While I might be able to keep a stable tree in memory after awhile, it would need to be updated pretty frequently for the first few users.   * The tree is stored in memory, not a database. This could lead to issues if a server dies or goes offline. Also, I'm planning to user multiple machines to host the project, so there'd need to be some sort of shared memory space they could all access.   * I don't have any categorical or descriptive data on each of the items to help with the decisions. The entire tree would have to be "did they like this other item or not?" I'm not sure if that poses any problems, but I haven't seen any other examples that use that sort of data for building the tree. This is my first attempt at machine learning, so it could be that these are all trivial problems compared to using other solutions, and I should just stick with decision trees. But if there's another solution that seems more appropriate, or if you can address some of these concerns, I'd love to hear it!


title: Choosing alternatives to Poisson regression for overdispersed count data

text: I'm currently analyzing data from a series of behavioral experiments that all use the following measure. The participants in this experiment are asked to select clues that (fictitious) other people could use to help solve a series of 10 anagrams. The participants are led to believe that these other people will either gain or lose money, depending on their performance in solving the anagrams. The clues vary in how helpful they are. For example, for the anagram NUNGRIN, an anagram of RUNNING, three clues might be:   1. Moving quickly (unhelpful)   2. What you do in a marathon race (helpful)   3. Not always a healthy hobby (unhelpful) To form the measure, I count the number of times (out of 10) a participant chooses an unhelpful clue for the other person. In the experiments, I'm using a variety of different manipulations to affect the helpfulness of the clues that people select. Because the helpfulness / unhelpfulness measure is fairly strongly positively skewed (a large proportion of people always choose the 10 most helpful clues), and because the measure is a count variable, I've been using a Poisson Generalized Linear Model to analyze these data. However, when I did some more reading on Poisson regression, I discovered that because Poisson regression does not independently estimate the mean and variance of a distribution, it often underestimates the variance in a set of data. I started to investigate alternatives to Poisson regression, such as quasipoisson regression or negative binomial regression. However, I admit that I'm rather new to these kinds of models, so I'm coming here for advice. Does anybody have any recommendations about which model to use for this kind of data? Are there any other considerations that I should be aware of (for example, is one particular model more powerful than another?)? What sort of diagnostics should I look at to determine if the model I select is handling my data appropriately? Thanks in advance for your advice.


title: How does glm.nb work?

text: I have been working with `glm.nb` from `MASS` package for quite a while now. However, there are somethings I seem to not quite able to get my head around. Suppose I have a data that looks like this:               Expression  Species  timePoint  Replicate     40            A          T1       R1     60            A          T1       R2     48            A          T1       R3     52            A          T2       R1     58            A          T2       R2     64            A          T2       R3     39            B          T1       R1     48            B          T1       R2     54            B          T1       R3     448           B          T2       R1     490           B          T2       R2     378           B          T2       R3      Now, if I would like to check if there is expression difference between `speciesA` and `speciesB` between time points `T1` and `T2`, then, I do:               require(MASS)     df <- data.frame( Expression=c(40,60,48,52,58,64,39,48,54,448,490,378), Species=c(rep("A",6), rep("B",6)), timePoint=rep(c(rep("T1",3), rep("T2",3)), 2), Replicate=rep(c("R1","R2","R3"),4), stringsAsFactors=T)     nb.fit <- glm.nb( Expression ~ Species * timePoint, data=df, control=glm.control(maxit=25, trace=T) )     summary(nb.fit)          Call:       glm.nb(formula = Expression ~ Species * timePoint, data = df,      control = glm.control(maxit = 25, trace = T), init.theta = 163.3237449,      link = log)            Deviance Residuals:       Min        1Q    Median        3Q       Max       -1.57348  -0.78584   0.06399   0.71550   1.27660            Coefficients:                               Estimate Std. Error z value Pr(>|z|)         (Intercept)           3.89860    0.09380  41.565   <2e-16 ***     SpeciesB             -0.04845    0.13391  -0.362    0.717          timePointT2           0.16184    0.12879   1.257    0.209          SpeciesB:timePointT2  2.07175    0.16888  12.268   <2e-16 ***           (Dispersion parameter for Negative Binomial(163.3237) family taken to be 1)         Null deviance: 947.708  on 11  degrees of freedom     Residual deviance:  10.024  on  8  degrees of freedom     AIC: 102.06          Number of Fisher Scoring iterations: 1                   Theta:  163                Std. Err.:  138       2 x log-likelihood:  -92.06       Now, the `estimate` obtained can be computed by log( T2/T1 of B) - log( T2/T1 of A) as follows:               > meanVal <- c( t( sapply( split(df, df[,2:3] ), function(x) mean(x[,1] ) ) ) )     > estimate <- log( meanVal[4]/meanVal[2] ) - log( meanVal[3]/meanVal[1] )     > estimate     > [1] 2.071749      Until this I follow. However, from here, I would like to know these:   1) How is the standard error estimated?   3) And how is the fitting of negative binomial distribution influence the std. error, z-value or the p-value? I mean, where does the `dispersion parameter` calculated used? I have read and tried to understand from quite a few tutorials and books. But I don't seem to understand. I would be really grateful if any of you could boil it down for me. Thank you!


title: Why does my research findings contradict my hypothesis result?

text: From my research findings/results it was clear that lecturers and students use different web 2.0 applications. But my null hypothesis result is contradicting this one of my null hypotheses is 'there is no significant difference between the web 2.0 application commonly used by students and those used by lecturers'. Please help, I am confused


title: SPSS:While performing ANCOVA, I get a warning: large number of columns may cause integer overflow. How do I resolve this?

text: SPSS Warning The large number of columns in the design matrix will cause integer overflow. The large number of columns may be due to too many levels in one or more factors, or to higher-order interactions among factors with many levels. It may also be due to too many factors This command is not executed. How do I resolve this?


title: Weighted Variance, one more time

text: Unbiased weighted variance was already addressed here and elsewhere but there still seems to be a surprising amount of confusion. There appears to be a consensus toward the formula presented in the first link as well as in the Wikipedia article. This also looks like the formula used by R, Mathematica, and GSL (but not MATLAB). However, the Wikipedia article also contains the following line which looks like a great sanity check for a weighted variance implementation: > For example, if values {2,2,4,5,5,5} are drawn from the same distribution, > then we can treat this set as an unweighted sample, or we can treat it as > the weighted sample {2,4,5} with corresponding weights {2,1,3}, and we > should get the same results. My calculations give the value of 2.1667 for variance of the original values and 2.9545 for the weighted variance. Should I really expect them to be the same? Why or why not?


title: Use cross validation to calculate a ridge regression parameter

text: I have a dataset of text messages of which I'm trying to filter out the spam from the legitimate ones. I have roughly 4600 pieces of data spread among 57 features and then their classification as spam or not. I have four 'versions' of the data, one the regular data and the other three I have applied various types of preprocessing. I'm supposed to be fitting each of these to a regression model using ridge regressing, and must use cross validation to calculate the ridge regression parameter. I have a loose understanding of what cross validation is but I'm confused as to how I apply it to this situation, particularly how to split my data. Could I get some guidelines/pointers for this?


title: On combining SVMs

text: Suppose we have a supervised training set $T=\{ (x_1, y_1),\dots, (x_n,y_n)\}$ where $x_i$ is an example and $y_i \in \{-1,+1\}$ is its label. Further suppose that examples are only observable through a feature extraction function $f(x;s)$ where $x$ is an example and $s \in \{s_1,\dots,s_m\}$ is an argument for feature extraction. For each possible value of $s$, we train a linear SVM (on the set $\{ (f(x_1;s), y_1),\dots, (f(x_n;s),y_n)\}$). Let $w_i$ be learned weights of the SVM for $s=s_i$. My question is on combining subsets of these SVMs for improved classification. Specifically, for a test example $x$, suppose that we have the scores of only the first two SVMs (feature extraction is costly): $w_1^T f(x;s_1)$ and $w_2^Tf(x;s_2)$. How can we combine these scores (optimally) to obtain a final decision? A trivial answer would be to train a SVM for each subset of $s$ values but this is not tractable. Ideally, I'm interested in a probabilistic interpretation. Assuming each SVM models $P(y|f(x;s_i))$, I want to express $P(y|f(x;s_1), f(x;s_2))$ using $P(y|f(x;s_1))$ and $P(y|f(x;s_2))$.


title: Why is my idea about correlation wrong here?

text: Autocorrelation is defined as $\rho(\tau)=E((x_i-\mu)(x_{i+\tau}-\mu))/\sigma^2$ which should have a value between $\pm1$. However, what if both $x_i$ and $x_{i+\tau}$ are significantly larger than $\mu$? Would it not be possible to have correlation greater than one? I'm particularly thinking about large lags where you only have a few sets of points before you run out of data. I'm sure this shouldn't be the case, but I can't see why.


title: Minimize a function with respect to a matrix

text: I have two sets of vectors, `A` and `B`. Vectors from set `A` live in an m-dimensional space, while those from set `B` are n-dimensional. I also have a mapping `f` from `A` to `B`. I would like to learn an `n x m` matrix Q such that the following sum is minimized: $$ \sum_{v \in A}{\|Qv - f(v)\|^2}$$ Are there any standard approaches to tackle this, for example in Matlab?


title: Variable selection / Dataset reduction for large datasets (in R)

text: I'm working on a behavoural scorecard modelling exercise, and many of the decisions taken to date have been based on the experience of a consulting credit analyst (whose experience software-wise is SAS) as I am primarily in BI. So far I have:   * a linux pc with 32gb of ram and an i7 processor   * an observation window    * ~90 potential characteristics    * a binary outcome In **R** , I have   1. loaded the dataset (225k obs of 88 vars, 1 outcome)   2. split the dataset up based on the recommendations/examples in the package **caret** i.e. predictors and outcomes split up (150k obs in training sample)   3. removed any variables showing a high degree of correlation (caret::findCorrelation)   4. cut all continuous variables into categorical intervals    5. reduced the number of variables based on near zero variance, missing values, and low information value (IV) (150k obs of 48 vars)   6. tried bestglm::bestglm, caret::train (with glm and glmnet), FWDselect::selection, FWDselect::qselection but eventually had to interrupt each of these due to not completing after 4 hours of 100% CPU usage   7. used FactoMineR:MCA to perform a multiple correspondence analysis (on predictors only) What I would like to do is have a selection of logistic regression models for say 4, 8, 12, and 16 variables that are the most predictive models at each point. I'm not sure if I'm going in the correct direction here with MCA as I've mainly been simply trying to find something that works in a timely fashion for reducing my variables further or going directly to variable selection steps. I would appreciate any advice on how to do any of step 6 better, whether 7 makes sense and what step 8 should be. Thanks, Steph PS Design decisions up to 6&7 can't be revised so please, no telling me off for them!


title: Does false discovery rate estimate some population quantity?

text: False discovery rate (FDR) is defined as > FDR = FP / (TP + FP). Does it estimate some population quantity, independent of sample, when the sample satisfies some condition? To make it clearer, True positive rate is defined as > TPR = TP/ (TP + FN) If all the points in the sample come from a sample distribution $F$ in the alternative hypothesis, true positive rate estimates $P_{X \sim F}($ reject null $)$. Thanks and regards!


title: R module for creating plots of prototypical individuals from fitted models?

text: In order to help interpret fitted models — especially those with interaction terms and non-linear components — I've found it useful to plot predicted values of a dependent variables for what we might think of as _prototypical individuals_. I would like to know if there are existing R modules that help do this. The idea behind these plots of prototypical individuals is to use parameter estimates from a model to plot the predicted values of some dependent variable across a range of an independent variables (often the x-asis is time in the context of longitudinal models) for different values of a (usually categorical) variable while holding everything at the sample mean or median. This is a technique discussed in depth by Judith Singer and John Willett in their book on Applied Longitudinal Data Analysis (see slide 20 in this presentation (PPT) by the authors) Another example might be this figure (Figure 2) from this paper (PubMedCentral) on development trajectories or this plot (PNG) from a recent paper of my own. I have written R code that takes a series of different model objects along with a long list of optional arguments that let you specify a wide variety of options for creating these plots. My module then gives you back a data frame that you can use with ggplot2 or your plotting software of choice. As I invest more time into my own code, I'm wondering if its worth continuing development of this and the release of a package in CRAN or if there is an existing system I might use and build upon. The only similar module I know about is visreg which is great but which, as far as I can tell, is designed to visualize the fit of model rather than creating plots of prototypical individuals from fitted models that one might use in publications.


title: Bayes factor for selecting between two beta-distributions

text: I have two beta-distributions: $H_1 = Beta(\alpha_1, \beta_1) $ and $H_2 = Beta(\alpha_2, \beta_2) $ (parameters are known), and I'd like to estimate whether a new sample $D$ rather comes from $H_1$ or $H_2$. It seems that Bayes factor is the solution $K = \frac{\int \Pr(\theta_1|H_1)\Pr(D|\theta_1,H_1)\,d\theta_1} {\int \Pr(\theta_2|H_2)\Pr(D|\theta_2,H_2)\,d\theta_2}$ but I don't really see how to compute it. How do I get and integrate $Pr(\theta)$ and $Pr(D|\theta)$? Also, how does it generalize to Dirichlets? Thanks, JP.


title: Statistics, Expected values and Philosophy

text: First a note: I am not a statistician. I studied maths at university (but opted out of every single stats class), and now find myself in a job where I'm doing stats. My question is a little bit philosophical, and I'm sure that I have my brain in a twist by thinking too hard about the wrong things, but I'm struggling with the concept of a random variable. Suppose a typist can type on average $\mu$ words per minute. If we take a sample of his (or her) typing we could approximate a normal distribution with mean $\mu $. I think that is pretty standard. So we can say we expect the typist to type $\mu$ words per minute. OK, what if we replace the typist with another person? Are we to expect the second person to type at $\mu$ words per minute? Of course not, because we only have a sample for the first person and not the second, right? The second typist's mean could be completely different. That makes sense, but then what if I put the original typist back in and tie one hand behind his back. Obviously the expected value is different, because the typist clearly can't type as quickly. But then what conditions can we expect the typist's average to be $\mu$? It seems any perturbation in the typist's state changes the probability. But we never really defined the typist's state to begin with, the expected value was calculated from a sample. And obviously every time the typist types he or she is in a slightly different state, i.e. maybe has a headache, maybe fingers ache, maybe daydreaming, but we still count those in our sample. So it seems to me, that technically, if we can accept some small changes in state, then why can't we accept big changes in state? Or where is the cutoff? i.e. If I replace the typist, I have just as much right to say that the expected words per minute is still $\mu$. I would like someone to explain why the above is wrong. Maybe it's an obvious answer.


title: GLMM in SPSS: what does significance of the fixed coefficients mean?

text: I conducted a GLMM analysis in SPSS with `y~a+b+a*b`. The output result contains two tables:   1. the fixed effects, for which the main effects and the interaction are significant;    2. the fixed coefficients, for which the significance changes according to the order of `a`, `b` and `a*b` in the model. My questions are: How do I interpret the fixed coefficients? (Are they important for inferring conclusions?) What does their significance mean?


title: Variance and bias

text: Consider the mean squared error $$\text{MSE}(x_0) = E_{T}(f(x_0)-\hat{y}_0)^{2}$$ $$ = E_{T}(\hat{y}_0-E_{T}(\hat{y}_0))^{2}+(f(x_0)-E_{T}(\hat{y}_0)^{2})$$ Is the first term after the equal sign in the second line basically saying the following: Get an estimate and find its deviation from the average estimates over all training samples. The estimate $\hat{y}_0$ will depend on the training set which is why we take the expected value. Also is $E_{T}(\hat{y}_0)$ constant?


title: Population model to model year to year dynamics

text: My task is to assess how various environmental variables affect annual population fluctuations. For this, I would use a model like: $$ \mbox{log} ( \mu_{i,j+1} ) = \mbox{log} ( \mu_{i,j} ) + R_{j} + \sum\limits_{k} \alpha_k x_{k,j} \\ N_{i,j} \sim \mbox{Poiss} ( \mu_{i,j} ) $$ Where $N_{i,j}$ is number of observed individuals at site `i` in year `j`, $\mu_{i,j}$ is the expected number of individuals at site `i` in year `j`, $x_{k,j}$ accross all $k$ is vector of environmental variables in year $j$, $\alpha_k$ are coefficients and $R_{j}$ is model coefficient meant to handle "background" population growth (it is not measured, this is just a model coefficient). However I think I will remove the year index $j$ and use it as intercept only, so that it doesn't hide possible global effect of environmental variables. This was the simplest version of the model - in the next stage I would like to handle overdispersion (not sure how yet) and maybe add some per-site random effect. **Questions:**   1. **How can I fit this model in R?** I can write model in WinBUGS but prefer to have "frequentist" solution in R, because it is much faster and inference is easier (one has p-values, t-tests, F-tests...). But which function or package use to fit it? I don't think this can be implemented using GLM! I spotted that my equation can be converted to: $$ \log\left({\mu_{i,j+1} \over \mu_{i,j}}\right) = R_{j} + \sum\limits_{k} \alpha_k x_k$$ Which resembles logistic regression: $$\log\left({\mu_{i,j+1} \over \mu_{i,j}}\right) = \text{logit}\left({\mu_{i,j+1} \over { \mu_{i,j}} + \mu_{i,j+1}}\right) = R_{j} + \sum\limits_{k} \alpha_k x_k$$ $$N_{i,j+1} \sim \text{Binom}\left(N_{i,j} + N_{i,j+1}, p = {\mu_{i,j+1} \over { \mu_{i,j}} + \mu_{i,j+1}}\right)$$ However, I am not sure this gives equivalent result; this converts the Poisson counts to Binomial; and possibly, it would not be quite straightforward to handle overdispersion (the Poisson overdispersion for animal counts is well covered and published; it is not clear how would it work in the binomial version). For this reason, **I prefer to compute the original model** as it is (Poisson).   2. How to incorporate a negative **density dependence** (i.e. the population growth is lower where there is a lot of individuals)? Add something like $\beta * \ln(\mu_{i,j})$ to the right side? Seems little strange to me...


title: What to do with interactions that improve the model fit but are not significant?

text: I'm using Poisson loglinear regression (from Generalized Linear models, SPSS) to analyze different patient variables that influence LNT-lymph node yield-(poisson distribution). After removing the main effects from the model that did not influence the outcome I tried looking for possible interactions. I used the Likelihood ratio test to identify those that significantly improved the model. I now have 9 main effects (8 of them are categorical variables) with 22 interactions. Now the problem is: Only 10 have significant levels according to the Wald test. My question: 1) Is it ok to put the main effects and all the interactions in 1 model? I know that you can't interpret the main effects properly like this, but does it influence the interactions aswell? 2)I also added interactions with variables whos main effect was not significant. Do they need to be added back to the model as a main effect to properly interpret the interactions? 3)How would I go about reporting all this? Do I need to mention the interactions that are not significant in any dummy category? I guess I should report the main effects in a model without interactions? I'm working with SPSS version 22. There are 4000+ patients in the database. To complicate things a bit, I'm workiing on an imputed database. I'm sorry if I'm totally wrong, I'm quite new at this and, as a medicine student, have little statistical knowledge.. Thanks in advance!


title: Marginal parameter estimation in copula with copula (dependence) parameter known

text: Suppose we have data $x_i, i=1,2,3,...n$ that are _dependent_ and identically distributed with marginal $f(\cdot|\alpha)$. If we model this with the likelihood $ L = c(F(x_1|\alpha),F(x_2|\alpha),...F(x_n|\alpha)|\theta)\prod_{i=1}^n f(x_i|\alpha) $ and the dependence parameter $\theta$ is known, can we apply some variant of the Expectation Maximization algorithm to estimate $\alpha$ using an iterative procedure with relatively simple steps? For instance, I considered a simple problem with exponential marginals and Gaussian copula (with known correlation), and did something procedural (and hokey). I introduced the unknown independent samples $\tilde{x}_i, i=1,2,...,n$ which you would compute by knowing the correct value of $\alpha$, mapping the $x_i$ to correlated Gaussians $y_i = \Phi^{-1}(F(x_i|\alpha))$ and then "undoing" the correlations $z = B^{-1}y$ and mapping the $z$ forward again to produce the independent $\tilde{x_i}=F^{-1}(\Phi(z_i)|\alpha)$. Here $C=BB^T$ is the correlation matrix. Here $\Phi$ is the (0,1)-normal cdf. If you turn this into an iterative procedure, using $x$ as the initial guess for the independent data, it seems to produce a series that (at least in my trials) converged. However, the whole thing is doubtful since it depends entirely on what you choose for $B$ (only defined up to a unitary matrix). I think it's the unitary invariance of the Gaussian hitting you when you try to basically do an inversion. Is there an obvious way to turn this kind of problem into a sane iterative procedure using simple steps like the EM? I feel like I'm missing a simple trick.


title: Is the following a valid MCMC sampler?

text: I have a model with two variables $X_1$ and $X_2$. In general, I would like to create a Gibbs sampler based on the conditionals $p(X_1 | X_2)$ and $p(X_2 | X_1)$. However, computing either $p(X_2 | X_1)$ or $p(X_1 | X_2)$ is intractable. Though, they can be computed up to a normalization constant. Fortunately, there is a function $f$ such that $p(X_2 | f(X_1))$ and $p(X_1| f(X_2))$ is tractable. I wanted to know if the following is a valid sampler that has the guaranteed MCMC convergence.   1. randomly initialize $x_1$ and $x_2$   2. until convergence is reached, perform the following: 2a. perform _one_ Hastings step to get a new value for $x_1$ (i.e. one accept/reject step) for $p(X_1 | X_2)$ based on the proposal $p(X_1 | f(X_2))$. 2b. perform _one_ Hastings step analogously for $x_2$. I am just not sure whether in order to get guaranteed convergence, 2a and 2b have to be run until convergence too (instead of taking one accept/reject step).


title: Distance correlation and prediction

text: If the distance correlation (ref. Gabor J. Szekely) $R_n(X,Y)>R_n(Z,Y)$ would the expected generalization error of a prediction model over $(Z,Y)$ be lower than $(X,Y)$ in predicting $Y$, where $Z$, $X$ are features?


title: Calculate standard errors: interaction between 2 factors, one of which has 3 levels in a regression model

text: I have the following model: $y \sim b_0 + b_1x_1 + b_2x_2 + b_3x_1x_2$. $x_1$ is a factor with 2 levels (0 and 1), and $x_2$ is a factor with 3 levels. I know that to calculate standard errors for the interaction term I should use $\sqrt{\text{var}(b_1) + \text{var}(b_2) + 2\text{cov}(b_1,b_2)}$. When I use `vcov()` in R to get the variance/covariance matrix I get a 6x6 matrix with the following column/row names:   * (Intercept)   * factor(x1)level1   * factor(x2)level1   * factor(x2)level2   * factor(x1)level1:factor(x2)level1   * factor(x1)level1:factor(x2)level2.  Where on this matrix is the `cov(b1,b2)`? Is it in the cell [factor(x1)level1, factor(x2)level1] or in the cell [factor(x1)level1,factor(x2)level2] or neither?   Or to put it better, if I want the standard errors for the marginal effects of x1 and x2 on y, how do I calculate them? When we have more than 2 levels in a factor does the above mentioned standard error equation change? Many thanks in advance! EDIT (migrated from OP answer): Maybe I should rephrase: I need the standard error for the marginal effect of x1 on y and the standard error for the marginal effect of x2 on y. the equation I put there supposedly gives that standard error (please see in this website How to calculate the interaction standard error of a linear regression model in R? and in Brambor T. et al Understanding Interaction Models: Improving Empirical Analyses. Political Analysis (2006) 14:63–82. doi:10.1093/pan/mpi014, equation 8). But maybe I am not using the correct terminology or references here? any ideas on how to calculate that standard error would be very wellcome.


title: Given tolerance intervals for multivariate hypergeometric, estimate a mean

text: Suppose we have a multivariate hypergeometric distribution on balls in $k$ colors. Let $n_i$ be the number of balls for the $i$-th color, so that we are drawing from $N = \sum_i n_i$ balls. Each sample for this experiment consists of drawing 30 balls from the urn. Finally, suppose we are given a function which maps each color to a range of percentages representing a tolerance interval. In particular, we might have: $$r_t(i) = [0.35, 0.5]$$ representing $$\mathbf{P}\left(0.35 N < X_i < 0.5 N\right) = t$$ (As a matter of fact, $t$ is not given in the data, but it seems reasonable to assume it's fairly high) I am trying to get a tight estimate on the expected number of blue balls in terms of $N$, the $n_i$, given that we know $r_t(i)$. Since the distribution is multivariate hypergeometric, we have: $$ \mathbf{E}(X_i) = 30 \frac{n_i}{N} $$ and $$ \mathbf{V}(X_i) = \mathbf{E}(X_i)\left(1 - \frac{n_i}{N}\right) \frac{N-30}{N-1}$$ My idea is to use something like Chebyshev's inequality to simultaneously estimate the variance and the mean, since they are related by a scalar multiple. $$ \mathbf{P}\left[\left|X_i - \mathbf{E}(X_i)\right| \leq \frac{\sigma_i}{\sqrt{t}} \right] > t$$ This _seems_ to encode the same (or, at least very similar) kind of information as $r_t(i)$, but I'm not sure how to plug $r_t(i)$ in to Chebyshev without making assumptions.


title: Maximum Likelihood for shifted Geometric Distribution

text: Really struggling with this please help. Find MLE for p and c \begin{equation} \ {f}(x,p,c) = (1-p)^{x-c}p \end{equation} x=c,c+1,c+2,..... p is between 0 and 1 c is element of the integers I am more interested in the mle for c For the mle of p is it fine to just take the derivative of the log likelihood function.


title: Estimating the variance of a bootstrap aggregator performance?

text: When performing cross-validation or bootstrap re sampling to estimate the performance of some machine learning algorithm, one commonly records the mean and variance of the errors obtained in of all the trials. This is commonly used in model selection, such as choosing the simplest (or quickest to run) model whose mean is within 1 sigma of the best mean (or some other rule of thumb). Typically I have seen it stated that when using bootstrap aggregators (such as random forest) the out of bag error should be an unbiased estimate of the generalisation error, and therefore we don't need to bother doing cross- validation to estimate the error, we get that for free. In my experience this is usually true, however I don't see how I can get an estimate of the variance of that estimate, comparable to that which you get from CV or bootstrap resampling. Is there a reasonable way to estimate this variance? We can't simply take the variance of the individual ensemble members (e.g. the trees in the case of random forest) since they will all have a much worse error than the ensemble (that being the point of bagging!). Is there some other reasonable approach? Please note that statistical/machine learning is not my main field and I tend to get a bit muddled with the terminology at times. Please correct any misuse or confusing terminology.


title: Where can I find online notes about statistical process control?

text: Please could you someone help me to find online resources related to the topic of quality control (statistical process control)


title: Interpreting definition of stable distributions

text: I am trying to interpret the following definition:               A non-degenerate distribution is a stable distribution if it satisfies the            following property:     Let X1 and X2 be independent copies of a random variable X. Then X is said to           be stable if for any constants a>0 and b>0 the random variable aX1 + bX2 has           the same distribution as cX + d for some constants c>0 and d. The distribution                is said to be strictly stable if this holds with d = 0 (Nolan 2009).      I have searched the internet for more intuitive explanations but couldn't find anything. Can someone help? I couple of questions:   1. what does independent copies of a random variable mean?   2. is the random variable X already a distribution? So does this say basically that if we add for example two normal distributions then it will also be normal?


title: Can someone explain Gibbs sampling in very simple words?

text: I'm doing some reading on topic modeling (with Latent Dirichlet Allocation) which makes use of Gibbs sampling. As a newbie in statistics -- well, I know things like binomials, multinomials, priors etc -- I find it difficult to grasp how Gibbs sampling works. Can someone please explain it in simple English and/or using simple examples? (If you are not familiar with topic modeling, any examples will do.) Thank you all in advance.


title: HMM ever better than CRF?

text: For classifying a sequence of instances, are there any specific circumstances that make Hidden Markov Models (HMMs) more accurate than Conditional Random Fields (CRFs)? I have seen several papers that show CRFs outperforming HMMs, but none showing the reverse (granted, my sample size is fairly small so far). Any real world examples would be great, since I'm still working to understand the differences conceptually.


title: How to use global weights when normalizing a mean value?

text: This is a very silly question, however, it confuses me a lot! Let's suppose that each observation consists of the prob of $A$ and the prob of $B$, with values $1/2$ and $1/2$ respectively. Also, let's assume that the global weight for $A$ (i.e. the mean of all observations for $A$) is $3/4$ and for $B$ is $1/4$. My problem is how I should weight $A$ up and $B$ down by considering the global weights. In other words are $0,5625$ and $0,4375$ the normalised weights for $A$ and $B$ respectively or I am getting something wrong in here? **EDIT:** Is this the correct process? $$ 0,5 * (1 + (0,59193 - 0,5)) = 0,5 * (1 + (+0,09193)) = 0,5 * 1,09193 = 0,545965$$ $$0,5 * (1 + (0,40807 - 0,5)) = 0,5 * (1 + (-0,09193)) = 0,5 * 0,90807 = 0,454035$$


title: An unfair coin with different face value problem

text: Could any one give me some hints for the following problem: A coin has two face: 1 and -1/2, and the probability for 1 is $P(x=1)=1/3$ and for $-1/2$ is $P(x=-1/2)=2/3$. The expected value is therefore $E(x)=1\cdot(1/3)+(-1/2)\cdot(2/3) = 0$.   Supposing for an infinite tossing, what's the chance that the average value would be bigger than $k$? (supposing the sample number is sufficiently large) Since if the face value are 1 and 0, it is simply a binomial distribution and can be approximate by normal distribution, so I also want to know if this problem can be transformed into binomial distribution $P(x>k)$ ? This is not any homework problem, but a problem I just think of and try to solve.


title: Joint probability of two correlated RVs

text: I am trying to get the joint PDF of two RVs $X$ and $Y$ where $aX<Y<bX$, so I am stuck in calculating the probability of $\mathbb{P}(X<x,Y<y|aX<Y<bX)$ any idea?


title: PCA with all categorical factors prior a regression with a continuous response

text: I realize that similar questions have already been asked and answered, but I am in need of a bit more detail and specific advice as I am new to PCA and statistical methods in general. My question is also a bit broader because I will be putting it in context and I need to know if I'm even headed in the right direction. I have a great deal of data. For each datapoint, there is one continuous response variable that I am interested in examining. Let's call it X. There are also five or six categorical variables, most of which have between three and ten possible values. One of them, however (let's call it A), has 154 possible values, and to complicate things further, each datapoint can fall in 1-4 of those 154 categories. For the vast majority of them, they just take one of the 154 values, but about 10% of them take two or three values, and maybe 0.5% of them take four values. (I am actually considering including a discrete but quantitative variable that will be equal to the number of values taken by S, as I think it might also be a relevant factor affecting X.) My ultimate goal here is twofold: to create a predictive model with multiple regression, and to use ANOVA to determine how much each of my variables' variance explains the variance in X. Someone more familiar with statistics than I suggested that I start with PCA because both multiple regression and ANOVA assume that all factors are independent. I'm pretty sure there are some correlations between a few of my factors (though I have no idea what they are) so I figured PCA would be a good way to begin disentangling. My questions are:   1. Can I perform PCA given the categorical nature of my data? If so, what method should I use to "dummy code" the variables? If not, what method would be more effective?   2. Will including a single discrete, quantitative variable (the number of values taken by A) complicate matters?   3. Will PCA even do what I want? (namely, disentangling the variables so I can then use multiple regression and ANOVA)   4. Whatever you recommend, is it possible in R, and if so, how? (I haven't even downloaded R yet but it's been recommended to me and it's free so I'm inclined to give it a swing. I have some programming experience in Python and C++ so in theory I can learn it without too much difficulty.) Thanks very much in advance.


title: Visualizing change in risk ratio along with confidence limits

text: I want to graphically show the change in relative risk along with confidence limits over 17 years. Can I use a forest plot without meta-analysis because it will have the added advantage to tabulate Z-values and p-values also? Here I want to show the progressive change in the relative risk and NOT meta-analysis. Or is there a better alternative for the graphic presentation?


title: What test to use for qPCR time course data with multiple genes and multiple strains?

text: My knowledge is severely lacking when it comes to statistics.. I have qPCR data at different times after induction and I don't know what kind of statistics to use on it. I have values for 9 genes at 6 different time points from 2 different mutants and I want to show for every gene if there is a change over time and a differences between the mutants. I just don't know what is possible with my data. **Example for 1 gene both mutants:** ![enter image description here](http://i.stack.imgur.com/1kWjx.png) I only have access to `R`.


title: How to produce a CI for a value predicted in CART?

text: **My goal is to create CI for the CART prediction of new_x** Consider the following code:               require(rpart)         set.seed(147830)     n <- 100     x1 <- runif(n)     x2 <- rnorm(n)     y <- x1 + 2*x2 + rnorm(n, 0, .5)     DAT <- data.frame(y,x1,x2)     fit <- rpart(y ~ x1 + x2, data = DAT)     new_x <- data.frame(x1 = .5 , x2 = .25)     predict(fit, newdata = new_x) # 1.353142 # "should" be ~1      The only way I can think of at the moment is to bootstrap on x1/x2 from DAT, build many prediction models with them, for each predict new_x, and then use the lower/upper percentiles. Is there some other way to do it? What are the advantages/disadvantages of it? Thanks.


title: should the sample be homogenous for developing a tool supposed to be generalized across 4 countries??

text: I am planning to develop a new tool and planning to use it across 4 countries. For the tool validation, is it right to include sample from very large population of 4 countries? Will it affect the normality of the distribution of data and tamper the data analysis process? Is it important that the sample need to be homogeneous for applying factor analysis?


title: Learning curves and SPSS

text: I followed a course on machine learning. The analysis of learning curves on the cross-validation set was recommended, as a way to see in what respect your model can be improved (more data, more features etc.). Does anybody know if SPSS supports the analysis of learning curves?


title: Is the MCMC simply a probabilistic gradient descent?

text: I'm learning about Markov Chain Monte Carlo methods, and to my undifferentiated mind, they basically resemble gradient descent with a stochastic component replacing the gradient computation. Is this a correct understanding. If not what key difference am I missing?


title: Clustering Techniques

text: I'm a little new to data mining and would definitely appreciate some tips. I'm using clustering algorithms looking for possible grouping in some variables described below. I've been using the Excel data mining add-in which connects to SSAS and uses the EM algorithm by default. I'm also using R, so far with the Kmeans algorithm. I have two independent variables X, and G which are integers, and three dependent variables A, B, C which are related to X and G by the equations below: $$A=100\left(\frac{2G\cdot.06+K}{2G\cdot.4+2\cdot.06+K+X}\right)$$ $$B=100\left(\frac X{2G\cdot.4+2G\cdot0.06+K+X}\right)$$ C=100-((((2xG)x(0.06)+K)/((2xG)x(.4)+(2xG)x(0.06)+K+X)))x100 + (X/((2xG)x(.4)+(2xG)x(0.06)+K+X))X100) where K is an arbitrary constant which varies. Most of the data follows the formulas above with some variation. I'm using clustering algorithms to look for groupings in the A,B, C variables. Any advice on strategies for looking for clusters, and how to tell when I've been successful would be greatly appreciated.


What is relationship between "significant correlation" and "significant difference" between two variables?


title: How to conduct regression analysis on ordinal variables in SPSS?

text: My data basically includes a number of ordinal (and some nominal) variables on perceptions of volunteers. I essentially want to see whether perceived personal benefit from volunteering (based on a Likert scale) is dependent upon certain factors, such as whether they were a good teacher (also an ordinal Likert scale). I tried multiple linear regression and it came out with some pretty good results, and I'm largely statistically illiterate so it took me some time to even do that. But reading a bit more into it, it seems that conducting ordinal regression specifically is the way to go. I don't quite understand what the difference is, I looked at some youtube vids and still can't figure it out. How do I interpret the table? Will multiple linear regression do? Its for my undergraduate dissertation and the article I'm basing this off is published and it just says 'regression analysis', which i assumed was just your average linear multiple regression. Thanks in advance. Help!


title: When optimizing a logistic regression model, sometimes more data makes things go *faster*. Any idea why?

text: I've been toying around with logistic regression with various batch optimization algorithms (conjugate gradient, newton-raphson, and various quasinewton methods). One thing I've noticed is that sometimes, adding more data to a model can actually make training the model take much less time. Each iteration requires looking at more data points, but the total number of iterations required can drop significantly when adding more data. Of course, this only happens on certain data sets, and at some point adding more data will cause the optimization to slow back down. Is this a well studied phenomenon? Where can I find more information about why/when this might happen?


title: statistical handling of lab values below limit of quantitation (BLQ)

text: There were several samples BLQ because of the lower limit of quantitation (LLQ) of the method, e.g. 5 ng/ml or less. Using the statistical program PRISM6 I marked these values together with the outliers (determined with the Rout method, 1% rule, before doing the calculations). The question was raised by the laboratory if excluding BLQ values from the analysis would cause a bias on the data.