title: p-value adjustment via permutation for multiple correlations

text: I have the following problem: I have performed 3200 correlations between a certain measure of connectivity and 1 behavioural outcome in a sample of subjects. As the single correlations (represented by R and p) are not independent from each other, Bonferroni or FDR correction (at given p-values) failed to leave any significant correlations. So I tried to follow this permuation approach: While permutating the available subject's connectivity data on the one hand, and also permutating the behavioural data on the other , I created a new pseudo-sample (in which subjects connectivity and behaviour are not matching anymore) for which I calculated all correlations again (3200). I did this 10000x and thus generated 10000x3200 random p-values based on the actual data set. Having ordered all p-values (low-to-high) I looked at the lowest 5% of those p-values (alpha): Random correlations would reach a p-value of 0.01 by 5% chance. In another group, the 5% lowest p-value would only reach 0.052. My question would be whether this would be a plausible way to correct for multiple correlations in a dataset with lots of non-independent single correlations. Particularly, is it reasonable to use an adjusted p-value that is even higher than the original one, 0.05. Maybe the chance to reach 0.05 on that data set is even lower than 5% due to variance etc. I am looking forward to your answers and comments. Thank you very much Best, Robert


title: Lift measure in data mining

text: I searched many websites to know what exactly lift will do? The results that I found all were about using it in applications not itself. I know about the support and confidence function. From Wikipedia, in data mining, lift is a measure of the performance of a model at predicting or classifying cases, measuring against a random choice model. But how? Confidence*support is the value of lift I searched another formulas too but I can't understand why the lift charts are important in accuracy of predicted values I mean I want to know what policy and reason is behind lift?


title: Are there other ways to represent time series part-to-whole data in a chart besides an area chart?

text: Are there other ways to represent time series part-to-whole data instead of using an area chart like the one in the example below? I want to highlight changes in population year to year. And I have some folks that find this a little confusing. ![enter image description here](http://i.stack.imgur.com/pTvHu.png)


title: Concordance and Discordance role in modelling

text: I am new to statistics and was asked to develop a statistical model, which I had started, they ask me to carry out concordance and discordance now, however I don't know anything about these terms except that the concordance is the probability that a pair of individuals will both have a certain characteristic, given that one of the pair has the characteristic and the opposite for discordance.   Still I don't know why I have to find them and what would be the appropriate value of both for a decent model.


title: Undefined link function in gamma distribution

text: The nonlinear model I am fitting gamma distribution with inverse or log is not converging. There is one observation having zero value in the response variable. Does this zero affects to model the data? Any answer warmly appreciated!


title: Is it better to assess the strength of regression predictors using semipartial correlations or standardized coefficients?

text: I have been told in the past that semipartial coefficients are better, but I do not remember why that is meant to be the case. A google search seems to reveal that many researchers use one or the other, and some report both. Assuming I'm going to only report one, is there a good reason to use one rather than the other, or to use one rather than the other in particular circumstances? Is there a good reason to report both?


title: R e1071 tune plot does not give me the best gamma?

text: I am using the R e1071 library for the SVM (Support Vector Machine) algorithm. And I used tune to find out the best Cost and gamma parameters. Though the plot doesn't seem to provide the actual best prediction. Here is some details:               gammalist <- c(0.005,0.01,0.015,0.02,0.025,0.03,0.035,0.04,0.045,0.05)     obj <- tune(svm, Class~., data = trainData, ranges = list(gamma = gammalist, cost = 2^(2:4)), tunecontrol = tune.control(sampling = "fix") )     plot(obj)      Here is the plot obtained: ![enter image description here](http://i.stack.imgur.com/Al0QG.png) The plot leads me to believe 0.02 is roundabout the best gamma. But I actually tested manually several others and find better results for 0.042. On a 200 sample, I get 23 errors with gamma=0.042, and 26 errors with gamma=0.02. How do you explain this?


title: Using the sample weights in regression

text: I have dataset A where households are interviewed each quarter for one year. I am focussing only on the households who have data for all quarters. With this imposition, observations are reduced by more than half. Next, I have dataset B where there is a variable "notinA" which is not in data A, but is in dataB. I am trying to impute this variable in data A using common demographic variables (available in both datasets) and nearest neighborhood matching. In this case all of the observations in data B are not used and some of the observations are used more than once when imputing the variable "notinA" in A from data B. Moreover, there are some observations in dataA for which the match is not possible, which are then dropped. So, my question is whether it is justifiable to use sample weights (pweights) of dataset A for the final data (with inputted variable from dataB). Can you please point to the literature that suggests how to deal with the sample weights for the example like this?


title: Multiple imputation - original model

text: I have a question. I have a dataset with some missing values that were not MCAR. I imputed them with fully conditional specification method iterations. I then executed my analysis on the basis of the imputed dataset. The results of the original model (with missing values: listwise deletion) did not change much in the eventual pooled model. My idea would be to go back the missing values dataset. What do you think?


title: Lisrel: Goodness-of-fit with FIML and missing data

text: When doing a missing-data analysis with FIML in Lisrel, I don't get a lot of goodness-of-fit statistics. I have discovered that you can calculate the fit- indices using the Chi-Square from an independence model. Now the question is first of all, how do I estimate an independence model in a missing data analysis with FIML in Lisrel (simplis syntax)? Is this possible? Thank you!


title: Classification performance and the feature set selection

text: I am now working on a classification problem. The generated feature set can be separated into two group. I did a comparison study: use all of the features; use the features of group 1 only; and use the features of group 2 only. It turns out that the performance of using all of the features are the worst. How to understand this kind of observation?


title: On the applicability of Benjamini-Hochberg

text: If I've understood correctly, the Benjamini-Hochberg (BH) correction is used to correct for the rate of false discoveries (FDR) when testing a collection of $m$ random variables, $\\{X_1, \ldots,X_m\\}$ against $m$ null hypotheses $\\{H_0^1, \ldots, H_0^m\\}$, of which $k\leq m$ can be true. Now consider a situation where you have a set of random variables that are _sorted_. For example, $\\{Y_1,\ldots,Y_m\\}$ with $Y_m\leq\ldots \leq Y_1$ and each $Y_i\sim F_i$ (i.e., distributed according to $F_i$). An example of such a scenario would be the eigenvalues of a random matrix. Suppose now, that given a single null hypothesis $H_0$, and test statistic $T_i$, the $Y_i$'s are tested in descending order — i.e., $Y_1$, then $Y_2$, and so on, until $H_0$ is false. Is the BH correction also applicable here or is it a fundamentally different scenario? Or does the question of controlling FDR not arise at all?


title: How is a uniform sample generated by the computer?

text: I know that uniform samples used on inverse cdf's can be used to generate normal samples. I was wondering how a technique would generate a uniform sample!?


title: review papers on P-values, adjusted P-values and common statistical tests in bioinformatics

text: I'm looking for review papers which give comprehensive reviews of P-values, adjusted P-values and common statistical tests in bioinformatics.


title: enrichment test and relations between variables

text: I'm in a situation in which I derived a list of 500 genes using a sophisticated method, of which I found a subgroup of 60 genes that are mutated according to a mutational annotation. To verify the enrichment was not due by chance, I generated 1000 random list of 500 genes (length) and I annotated the genes for mutations. Then I compared the number of mutated genes in random lists with the number of mutated genes in y reference list. However this comparison does not take into account the fact that in my reference list mutated genes may be related between them. How can I control this? I can not estimate the relationships in any way and I can not repeat the method by which I derived the reference list of genes on the list of 1000 randomly generated lists because of it is time consuming. Any idea or statistical test about this issue? Best, F.


title: Winsorizing data

text: I'm analysing data that relates to trade sizes over time and want to test that the changes are statistically significant. The following figure sums up what I want to do. ![enter image description here](http://i.stack.imgur.com/3Zxv3.jpg) For example, for node WW below (with a mean bet size of 60)and node W (mean bet size 55), I want to test whether the change in bet size from W to WW is significant. Up to now, I've just calculated the mean difference between the two samples using only the survivors who bet in WW (so matching 220 WW's with 220 W's and disregarding 5 of the customers that dropped out). Is this ok? I'm not taking attrition from the dataset into account in any way with that approach. . Secondly, the data is slightly skewed so I want to winsorized it. I had tried the approach of just winsorizing at each 'node' (10 in each tail, for example). So to get a winsorized mean is trivial, but the formula for the winsorized standard deviation takes account of the number of winsorized datapoints (STDEV*N-1/N*2K-1) in Excel. I was also winsorizing the mean difference, getting a standard error with the appropriate 'winsorized' formula and calculating a t-stat. That worked fine up to now. I no longer want to winsorize at each node. By winsorizing the top 10 and bottom 10 at each node, I'm trimming off values that may not be outliers in the dataset in terms of magnitude - it just happens that they're the top and bottom 10 in a particular node. I want to winsorize the entire dataset, so getting a vector of all the stakes placed between bet 1 and bet 4 and winsorizing the top and bottom 50, for example. Then go ahead and calculate the means at each node based on this data. The problem is in calculating the standard error for the t-test (or non- parametric test, if appropriate). Can I just go ahead and calculate a standard deviation as normal (even though the data HAS been winsorized)? I don't think I can even apply the formula for the winsorized standard deviation as I won't have a value for 'k' above (the number of winsorized datapoints). By winsorized at the start rather than at each node, I don't know which datapoints have been winsorized at each node. Am I on the wrong track completely with this?


title: How to derive the characteristic function of a polar coordinates representation of a bivariate normal

text: Suppose to have a bivariate normal variable $\mathbf{x}=(x_1,x_2)$ with mean $\mu$ and covariance matrix $\Sigma$. I move from $\mathbf{x}$ to $(\theta,r)$ where $x_1 = r \cos \theta$ and $x_2 = r \sin \theta$: the representation in polar coordinates of $\mathbf{x}$. My question is: is it possible to derive the characteristic function of $(\theta,r)$? and the characteristic function of $\theta$?


title: Full and fractional factorial designs

text: I have a data set consisting of six independent environmental variables (all binomial: present / absent) and one dependent variable (binomial: disease present / absent). In order to determine the combination of factors that have the highest probability of leading to disease, I first need to conduct an Expert Opinion poll where I will have several experts rank all possible combinations of variables according to their probability of leading to the occurrence of disease. Then, I will obtain regression parameters for each variable using a conjoint analysis approach where each expert conforms a level (hierarchical design), the six environmental variables are independent variables, and the rankings are the dependent variable. There being six factor variables, there exist a total of 64 possible orthogonal combinations. I reduced this overwhelming number of possible combinations (while retaining orthogonality) using the `AlgDesign` package of R. Here is the code followed only by relevant pieces of output:               levels.design = c(2,2,2,2,2,2)     full.design <- gen.factorial(levels.design)             X1 X2 X3 X4 X5 X6     1  -1 -1 -1 -1 -1 -1     2   1 -1 -1 -1 -1 -1     3  -1  1 -1 -1 -1 -1        .................     63 -1  1  1  1  1  1     64  1  1  1  1  1  1          set.seed(69)     fractional <- optFederov(~., data=full.design, approximate=FALSE, criterion="D")     fractional      The result is a subset of 12 combinations to be included in the conjoint analysis:               $design         X1 X2 X3 X4 X5 X6     4   1  1 -1 -1 -1 -1     5  -1 -1  1 -1 -1 -1         ................     57 -1 -1 -1  1  1  1       From what I understand, doing a regression analysis on all 64 combinations should lead to the same regression parameters as those obtained if I use only the reduced set (i.e. 12 combinations from the fractional factorial design). Questions:   1. Do the code and the resulting output make sense?   2. Could anyone point me to a good and simple reference on how this fractional design works? I am afraid I might be doing things wrong by selecting a subset that produces different results from those obtained if a full factorial design was employed.


title: Regression Trees / Boosted Regression Trees for Tweedie Distribution in R

text: I am currently working at work on a project that attempts to predict an environmental change variable. I am personally not a huge fan of the project, but I still want to do the best job possible. Anyhow, let me first describe the properties of the data, and then state my question. The environmental variable we are trying to model is continuous, and ranges from 0 to 25 (it can have values such as 0.1345 or 1.2335 or 5.674). The environmental variable also has a significant mass at zero (around 30% of the data are zero values), and most of the data is in the range between > 0 and < 1\. To complicate things further, around 3% of the data have extreme values of greater than 2. In my opinion, the distribution of the data resembles a Tweedie distribution. We have around 5 million observations in the dataset. We are going to predict the environmental variable using a set of eight explanatory variables. We have been modeling the prediction of the environmental for a week, and the predictive power of our results has been meager. Here are the different modeling approaches I have used:   1. GLM model with a tweedie distribution (glm function in r). This model approach overestimate the environmental change when the variables has low values. Most values between 0 - 1 are significantly overestimated, and the model does not predict high values very well either. The GLM approach produces a very bad fit for our data.   2. Generalized Additive Model with a tweedie distribution (bam function from mgcv package in r). This approaches fits the data slightly better than the model above, but still significantly overestimates values in the low range.   3. Regression tree model (rpart model in r). The regression tree model most accurately predicts values in the lower range of the distribution, but fails to predict zero values, and performs also poorly for values greater than 2.    4. Boosted regression model (dismo package in r using gbm.step). The model significantly overestimate values in the lower range of the distribution. For example, if values are 0.23 it predicts values to be 1.67 and so forth. I believe the gbm.step model to be inadequate for our data, since the family of distribution in the package only models the bernoulli (=binomial), poisson, laplace or gaussian family. None of which accurately describe our distribution. Since our dataset is very large, I split the data 50/50 into a test and training dataset, and evaluated model fit on a variety of test statistics for the test data set. Knowing the structure of our data, can anyone think of our modeling approaches that would possibly produce better predictive results? I am still new to machine learning techniques, and I hope that someone might know of other techniques that could be employed. My strongest program language is R, but I can also do this analysis in Python or Stata.


title: Naive Line Search for Gradient Descent

text: I'm trying to understand the Line Search approach to Gradient Descent (http://en.wikipedia.org/wiki/Line_search). It seems that a naive implementation would               while (termination criteria not met)         1) Compute gradient at a current location         2) Move in that direction until the objective function begins to increase         3) Set current location to the location before the "increase" in step 2         4) GOTO 1      Is this understanding correct?


title: Is deep belief network a belief network

text: I am new to deep learning. Please help me.. What is belief network. What is the use of it? How to learn belief networks? Is deep belief network a belief network with multiple layers?


title: What is the distribution of q-values under the null

text: It's well known that under the null hypothesis, the p-value has a uniform distribution. However, how can we determine the distribution for q-values (false Discovery Rate adjusted p-values)? I guess it is no longer uniform since the lowest is always increased. Thanks!


title: Which groups of results is the closest to a central point?

text: I'm building an application where a specific location is chosen, multiple services are polled to return results for that specific location and shown on a map. I have the results from the different sources shown in different colours on the map, so you can get a feel for which results are the "Best". In this instance "Best" is hard to precisely define, but means which are generally the closest to the starting point Because everything is lat/lon, I can then calculate the distance each one is from the starting point. I was thinking calculating the average of the sum of the squares of the distances, as this would favour those result sets who have the closest points. This might skew otherwise good results that have a just a couple of distant outliers, but it's worth a go Are there any good algorithms out there that would help in this sitution? I've looked at "clustering" algorithms, but they are good for organising data into groups. However, in this case, I already know what the groups are, I just want to know which group is closet to my target point thanks


title: Issues with SEM and CFA in AMOS

text: **Firstly:** I want to run a CFA with 7 latent variables and 39 observed variables. For a guide of how I constructed my CFA, I used this Youtube Video: http://www.youtube.com/watch?v=JkZGWUUjdLg It looks to me to be exactly correct according to the process explained in the above mentioned clip. However the CFA will not run, and the output reads this message: "The model is probably unidentified. In order to achieve identifiability, it will probably be necessary to impose 15 additional constraints". My first question is, why is the CFA not running, and why does it say I need to add 15 more constraints? To my understanding this CFA is correct. **Secondly:** I am running a SEM with 9 latents and 49 observed variables. This model will also not run, and when I check the output it reads: "Iteration limit reached, the results that follow are therefore incorrect". I do not understand what this means. Does anyone know what is going wrong with this SEM? A calculation of the sample size needed suggests I need 289 participants. I have 325 participants, although they are composed of a smaller data set multiplied by 3 or 4 times to reach sufficient sample size (could this be an issue?).


title: Clustering cloud fitting inverse functions

text: I have a cloud of points that can be clustered so that each set of points in the cluster can be fitted with an inverse function : $f(x) = cte / x$. What would be an approach for clustering my dataset when discrimining criteria is fit to an inverse curve. I am on the computational side, I am looking for a simple fitting algorithm for a beautiful display with curves that are fitting the data thus demarcating several cluters of points.


title: How to estimate a stochastic time delay

text: Say I have irregular samples of two random processes that are both based on an unobserved process, $Z(t)$. Their SDEs are $\mathrm{d}X(t) = \mathrm{d}Z(t-u(t)) + \epsilon(t)$ $\mathrm{d}Y(t) = \mathrm{d}Z(t-v(t)) + \eta(t)$ where $\epsilon(t)$ and $\eta(t)$ are white noise and $Z(t)$ is any Levy process. The delays, $u(t)$ and $v(t)$ are independent, continuous, and $\ge0$. I am interested in the relative delay, $u(t) - v(t)$, between the two observed processes. Any suggestions on how to estimate its expectation? What about its 2nd moment?


title: Question about asymptotics of steepest descent method in the context of adaptive filtering

text: The model which will be used is defined as $e(n) = d(n) - y(n)$ with $y(n) = x(n)^Tw(n)$. where $e(n)$ is the error term of the n-th observation, $x(n)$ the input vector of the n-th observation, $d(n)$ the target output of the n-th observation and last but not least $y(n)$ the observed dependend variable of the n-th observation. The weight difference in the steepest decent approach is defined as $\Delta w(n) = - \mu \frac{\partial \mathbb{E}[e^2(n)]}{\partial w(n)}$ where $w(n)$ is a $m\times 1$ vector of weights. Hence we can write down the following equation $w(n+1) = w(n) - \mu \frac{\partial \mathbb{E}[e^2(n)]}{\partial w(n)}$ If we assume that $\lim_{n \to \infty} x(n)x(n)^T = R_{xx}$ which is a finite matrix with full column rank and $\lim_{n \to \infty} x(n)d(n) = R_{dx}$ which likewise is a finite matrix with full coloumn rank we can compute the derivative $\frac{\partial \mathbb{E}[e^2(n)]}{\partial w(n)}$ as follows: Since $e^2(n) = [d(n) - x(n)^Tw(n)]^2$ we get $\frac{\partial \mathbb{E}[e^2(n)]}{\partial w(n)} = \mathbb{E}[-2x(n)[d(n)-x(n)^Tw(n)]]$ So if we split up the rightmost term we get $-2\mathbb{E}[x(n)d(n)] = -2R_{dx}$ and $-2\mathbb{E}[x(n)x(n)^Tw(n)] = -2R_{xx}w(n)$ (*) by the previous assumptions. The thing I dont get is why (*) is true. Since both $x(n)$ and $w(n)$ are dependend on random variables I dont see how we can get the equaltiy which is obviously achieved by taking $w(n)$ to be constant or to be independend of $x(n)$? For reference see http://www.ee.cityu.edu.hk/~hcso/it6303_4.pdf slide number 19. A hint in the right direction is very much appreciated :)


title: Given a historical disease incident rate of x per 100,000, what is the probability of y per 100,000?

text: Excuse the rather basic question, but I was reading this article on thyroid cancer in Fukushima and it was reported that 3 cases in 38,000 children were found in the previous fiscal year. Another web site tells me that for 15-19 year olds in the UK (I cannot find Japanese figures) the rate is about 1.6 per 100,000 per year, which would be 0.6 per 38,000. The rational person in me thinks that three cases rather than zero or one in a year is just a statistical hiccup, but the Think Of The Children person in me says it's a five times higher incidence. So, can someone statistically assure me/worry me about these figures? (PS: I realise that increased testing specifically looking for thyroid problems may skew the stats for last year by finding tumours earlier, but can we ignore that here, please?)


title: Methods of dealing with non-independent data in ordinations

text: Forgive if this question is too broad, but i'm wondering if there is a way to deal with spatial pseudoreplication in ordinations? I have a split-plot treatment design, and am wanting to run ordinations (maybe using NMS) for each year of my experiment to see how relationships between different plots receiving different treatments change. My concern is that due to the nature of the split plot design, these data are not truly independent, and this might affect the outcomes of the ordination...


title: Is PCA the same as CCA for the same random vector?

text: Reading Wikipedia about canonical correlation analysis (CCA) for two random vectors $X$ and $Y$, I was wondering if PCA is the same as CCA with $X$ and $Y$ are the same? Thanks and regards!


title: Classes of asymptotically normal statistics

text: Consider an iid sample $X_1, \dots, X_n$, from a "well behaved" distribution with pdf $f(X)$. My question is the following: what are the classes of statistics $S(X_1, \dots, X_n)$ that are asymptotically normal? I know about the asymptotic normality of specific statistics (such as the sample mean), but think it would be nice to have a list of all the families or classes of statistics that are asymptotically normal. I think M-estimators are an example of such a class, but are there others?


title: Finding correlations between aggregated and non aggregated data

text: I am trying to correlate server usage data with CPU utilization. E.g., X number of transactions on the server correlates with Y amount of CPU utilization. The transactions are reported hourly and the CPU utilization is reported on a per minute basis. The current approach is to average the CPU utilization over the hour and then plot that against the hourly transaction data. My concern is that there can be a spike in the load that lasts only a few minutes which would drive the CPU utilization up to a critical level during that time, but that spike (and the server performance degredation) would be lost when the CPU is averaged. Is there a better way to approach this? More granular load data may not be available so any help is appreciated.


title: Multinomial logit with different sample sizes in three groups

text: I have 3 groups in my dependent variable. The sample sizes on the 3 are very different. I have about 10000 in the first group, 35000 in the second and 100,000 in the third. My question is: Should we sample down group 2 and 3 before modeling? Paul Allison talks about this topic briefly in his book on logistic regression using SAS, where he says if you have two groups with different sample sizes, we could sample down one and finally adjust the intercept of the regression equation. I want to know how this will work in a multinomial situation. Any help is appreciated. If anyone can help me with what the equations would look like for a multinomial logit model, that would be great with intercepts being adjusted. I believe there are two equations. Wonder how the intercepts would get adjusted. Thank you


title: Are bootstrapped parameter standard errors unbiased in time series data?

text: I'm analyzing multivariate panel data which includes autocorrelated time series variables. The response variable is binary and I'm considering a random intercept to account for between subject variation. So a mixed effect logistic regression seems to be in order. Is it safe to say if I ignore the autocorrelation structure that my parameter estimates will still be unbiased, but their standard errors will be biased too low? To remedy this, can I use block bootstrapping to estimate unbiased standard errors for my model parameters?


title: Gibbs sampling product of normals as conditional

text: I am deriving a gibbs sampler for a joint distribution, where the conditionals of various parameters are product of two non-standard normal distributions. Usually, I have seen that in Gibbs sampling people sample their conditionals from standard distribution. But if I sample from the two normals and then multiply the samples should I be fine or should I use Metropolis Hastings step inside Gibbs?


title: Accuracy Measurement in Biased Dataset

text: Lets say we have data-set consists of 100 instance, each instance belong to one of 6 classes {1,2,3,4,5,6}. Furthermore, 70% of the instances belong to class 3. For classification problem, if the output of the classifier is always 3 we will get accuracy of 70%! Is there a way for measuring accuracy with heavily biased data like the one above?


title: Why does KL divergence show up in the proof of Hoeffding's inequality?

text: In some textbook the KL divergence shows up in the proof of Hoeffding's inequality (e.g., eq. (5) of this material). In contrary, most other textbooks seem not mention this fact. I know that KL divergence is essentially Bregman divergence, i.e., the difference between the true value of a convex function and tis locally linear approximation. But I have no idea why there is a place for KL divergence in Hoeffding's inequality? Any intuitively explanations?


title: Mathematical formulation of correspondence analysis?

text: From a website > you can think that Correspondence Analysis is a categorical data version of > PCA. But the main usage of Correspondence Analysis is different from that of > PCA, and it is more like clustering or Factor Analysis. With Correspondence > Analysis, we can analyze and visualize the relationships among your observed > data, and see which parts of the data are associated with another part of > the data. Wikipedia also says > It is conceptually similar to principal component analysis, but applies to > categorical rather than continuous data. In a similar manner to principal > component analysis, it provides a means of displaying or summarising a set > of data in two-dimensional graphical form. But I don't understand its details process, because I don't know what kind of problem Correspondence Analysis is trying to solve? Is there a clear mathematical formulation of the problem? Thanks and regards!


title: way to determine best number of clusters weka

text: When using weka library for clustering ,Is there any way to find best number of clusters. The EM methods are defined but how to use these methods in java code. Please suggest


title: Simulating values from an estimated value and confidence interval in R

text: I am using the estimated county-level poverty measure from the Small Area Income and Poverty Estimates (SAIPE) as the dependent variable in a regression analysis. This value is itself the result of a model and comes complete with upper and lower 90% confidence interval bounds. To be clear, each of my 3000+ observations has an estimated value and its own confidence interval based on the sample size for that county (targeted at 2.5% of population) I am wondering about the best way to incorporate the uncertainty in my dependent variable into my regression model. One way I have imagined is doing a random draw from a distribution using the estimated value and confidence interval for each observation. I would then re- run my regression using this simulated value for the dependent variable in my regression analysis and compare model outcomes to the model using the estimated value. By simulating new values and comparing many times I would gain an understanding of how sensitive my findings are to the estimates on the dependent variable. A key component of this is pulling a random number based on the estimated value and the upper and lower confidence intervals. The data is left and right censored at 0 and 100 and the estimated value is not centered within the confidence interval that is: abs(estimate-cl) != abs(estimate-cu) I was intrigued by the discussion here: Sampling random numbers from a distribution with asymmetric confidence intervals generated by a bootstrapped estimate but a modified version of the code just generates the estimated value Here is an example using the first 6 records from the 2008 SAIPE               sample.size<-c(1258.850,4405.300,745.900,539.725,1444.850,273.02)     POV08L90<-c(8.77,8.24,18.08,13.90,10.55,22.45)     POV08H90<-c(12.54,11.18,24.82,20.87,15.27,33.83)     POV08<-c(10.7,9.9,24.5,18.5,13.1,33.6)     test.data<-data.frame(sample.size,POV08L90,POV08H90,POV08)          gammaGenerate<-function(dat){       for(i in 1:length(dat$sample.size)){         n<-dat[i,"sample.size"]         cl<-dat[i,"POV08L90"]         cu<-dat[i,"POV08H90"]         barx<-dat[i,"POV08"]         talpha = qt(p=0.95,df=n-1)         s = (cu - cl)*sqrt(n)/(2*talpha)         kappa = 6*s*s*n*( cl - barx + talpha*s/sqrt(n) )         gamma.shape = 4/(kappa*kappa)         gamma.scale = s/sqrt(gamma.shape)         gamma.shift = barx - gamma.shape*gamma.scale         print(c(barx,(rgamma(n = 5, shape = gamma.shape) + gamma.shift)))           }     }     gammaGenerate(test.data)      Any help you can offer--either directing me to a better method of dealing with the uncertainty in my dependent variable, or an explanation for why my rgamma always lands at 0 would be very welcome.


title: Bootstrapping with replacement

text: If I want to nonparametrically bootstrap the following matrix (with replacement):               1   2   15  12  14  22     3   5   1   9   29  19     2   22  12  11  9   13     14  3   6   16  17  23      5   34  22  12  13  16     7   14  13  19  17  3      Would I "jumble" up only the rows, to obtain matrices like:               5   34  22  12  13  16     14  3   6   16  17  23      2   22  12  11  9   13     3   5   1   9   29  19     2   22  12  11  9   13     5   34  22  12  13  16      Or would I jumble up every element to obtain rows with any combination of all the values in the original matrix? The columns represent different species of animal and the rows different areas of woodland, the elements are the number of animals observed in each area. I'm looking to form confidence intervals for the diversity index of the whole woodland. I'm defining the diversity index (Shannon index) by summing the rows to get a species count for the whole woodland, then performing a diversity index on these rows. It returns a single value.


title: Recommendation for peer-reviewed open-source journal?

text: I have a manuscript on a bootstrap method for testing hypotheses of one mean, and I would like to send it for publication, but I have a moral dilemma. I have signed on to the protest against Elsevier for their unethical business practices, and reading up on the whole issue really made me question the ethics of other for-profit academic journals. I would therefore like to publish in a journal which is non-profit, preferably open-source, even though I understand that such journals aren't yet considered on par with the more established journals as far as prestige is concerned. Fortunately, since I already have tenure, that's not a big consideration for me. I would appreciate any recommendations.


title: Mean squared error for data with skewed distribution

text: I am doing regression task and the response variables in my dataset have a skewed distribution. Say, for the sake of simplicity, that I have a model Y~X and Y(response variable) is in [1,5] but there are many more values in the [4,5] range than the [1,2] range. As a result a majority predictor that predicts 4.5 for Y regardless of the value of X can complete with linear regression if I just look at MSE error. I am wondering if there is any methodical way to correct MSE to consider this case and maybe penalize the majority predictor more when it gets a Y=1 wrong than when it gets Y=5 wrong. Basically I am looking for a fair error measure for skewed data. **Update:** for the sake of simplicity let's say the majority predictor predicts 4.5 for everything regardless of the value of X. My predictor predicts 1 acuratly but always predicts 4.4 for 5's. the test set for Y is one 1 and the rest are 200 numbers and each equals 5. Based on MSE the majority predictor is better than my predictor but it does not make any sense. I want to modify MSE to favor accurate predictions for 1's compared to accurate predictions for 5's. Perhaps I can multiply each residual by the inverse of the frequency of the actual Y? also how can I use such an MSE for cases that Y is continuous? **update 2:** So some folks suggested that maybe Y should be sampled more properly and I need to find a subsample of my data that provides a uniform distribution of Y. This is not possible in my case. Say I am crawling amazon and most of the ratings that I see are 5 (since amazon removes under performing products) but there are some items with rating 1 as well. Now if I use a majority predictor that predicts 5 everywhere it is going to beat my SVM in terms of MSE value but the majority predictor brings no value to my system. Moreover I do not want to throw out my data just to make the distribution of ratings uniform. I believe it should be one by properly selecting the metric (error measure) Attached is an example distribution. ![enter image description here](http://i.stack.imgur.com/nTK14.png)


title: Plot with 2 conditions on the y axis

text: First, I want to clarify this is a very basic question. I'm new using R (and programming in general). I have a data.frame with 3 columns: block, average condition 1, average condition 2. The blocks run from 1 to 9 and the average are response times averaged across subjects. What I want to do is to plot both averages, one for each condition, while having in the x axis the progression of blocks. Therefore, I would expect to have one line for each condition running above the x axis. My question is what plot function should I use. I've tried with nearly everyone but what I thought would serve (simply plotting a vector with two columns as entry, against the block column) doesn't work, it gives the error:                'x' and 'y' lengths differ      So I thought I should first plot only one condition and then somehow add the other line with a second command. I successfully plotted one condition as I wanted to:               plot(meanrts2[,3],meanrts2[,2],ylab="RT",xlab="block",type="l")      So here I'm choosing the 3rd and 2nd columns to do the plotting, but what I want is to have two columns running parallel, each one with different colors. I tried abline, using the other columns I haven't plotted yet as an argument, but I get the "invalid a=, b= specification" Any suggestions are very welcome. Any other information I missed you can always ask.


title: Linear Discriminant Analysis for Training and Testing of data sets

text: I browsed through lots of literature on **Linear Discriminant Analysis** on the web. I can say that, I have understood most of it. This is what I understood. You take the data, calculate mean of the **data** , **mean per class** , **scatter between classes** , **scatter within a class**. Use this information to project multi-dimensional data onto a single line and hence reduce dimensionality. However, I don't understand how do I use it for **training and testing of the data sets**. Also, I know that **Y = X dot W** Y is projected data X is original data and W is the calculated vector used for projection. So, do I save **W** and use it to project new data? I hope I am on the right track!


title: Mathematica SmoothKernelDistribution equivalent in R

text: @vucko gave me excellent answer on my question, unfortunately using Mathematica code. I'm trying to rewrite it in R and I'm lost in R functions providing kernel density estimation. I have bivariate dataset with more than 46,000 rows (so I am also looking for a high performance solution--@vucko's solution is very time consuming). I would like to apply kernel density estimation and decide if some point lies in area with some density estimation level (confidence level respectively). ![kernel estimate](http://i.stack.imgur.com/lDgfo.jpg) @vucko in his answer selected two groups. I need only to know if some point lies in the green group or not. And that should be done with R. I experimented with `kde` and `bkde2D` functions but they don't provide me desired functionality as Mathematica `SmoothKernelDistribution`. Can you please me show the direction? For normal distribution I found the `ellipse` function which approximated data with some confidence level and used `inside.owin` function.


title: Which is the correct hypothesis test formulation for this problem?

text: The problem I am working on can be formulated as follows: $S$ is the set of events $X_i$'s where each $X_i$ has occured different number of times. Also, if you think this process in turns; in each occurance of $X_i$, it is followed with one of the events $Y_{i,j}$'s. What I want to investigate is whether one -or any numbers- of $Y_{i,j}$ is occuring significantly more often than others. I first thought of calculating the mean and variance of the number of times that $Y_{i,j}$'s occur, seperately for each $X_i$. After that, I could perform a t-test where $t={( {mean}_i -Y_{i,j} )}/ {stddev}_i $. And if any $Y_{i,j}$ was significantly more then this mean with an alpha of 0.05, then I would conclude that the proportion of occurance of this $Y_{i^* , j^*}$ is significantly more often than others. Then I noticed that it might be required to divide stddev with root of $n$, where $n$ is the number of occurance of any given $X_i$. After that, I again doubted this conclusion, thinking that a hypothesis tests based on proportion formulas might be better representative of the problem at hand. Since I am not a statistician, I thought of asking for advice on this issue. What approach do you think, would be appropriate for this problem?


title: Does it make sense to measure recall in recommender systems?

text: Assume I've built a recommender system that (given say movie rankings or whatever of many users) will produce a list of 10 recommended movies for each user to watch. Imagine that I also have some large pool of movie items, together with a log of user ratings together with movies that they actually decided to watch. So I want to use this data set to evaluate my system. I've seen in the literature that these "suggest some good items" tasks are usually evaluated using precision, recall and F1-scores (e.g. see [1]). I guess that I should be interested, in particular, on "precision at 10". However I'm not quite sure how one is supposed to compute these measures (or if they make any sense) in the scenario that I've described above. Apparently, the preferred thing to do is to randomly break the sample into a "training" and a "testing" part. And then feed the training data to my algorithm so that it can come up with a list of 10 predictions. Now precision sort of makes sense, I can check from the 10 predictions how many of these are actually found in the movies watched by the user in the test data. However for recall, if the user watched a lot of movies in the test data, say 50 or so; there is no way to obtain a "good" recall score, simply because my system was constrained to produce only 10 movies and I would get at most a 1/5 = 0.2 of recall. Alternatively, if I constrain the test only to guess the "next 10 watched" movies of the user (so that there is a chance to get a "perfect recall"), then precision and recall will always be exactly the same number (if the number recommended and the number relevant for the user is the same, precision and recall are also always the same). Am I doing something wrong? Or these metrics simply don't make much sense in the considered scenario?


title: Simple problem formulation for EM algorithm

text: I understand the EM algorithm, I understand for example how we get $Q(\theta, \theta^t)$, but I have trouble translating a real-world problem into the EM framework. For example, I'm given this problem where each document $x$ is classified into a class $y$, where there are $M$ documents and $K$ classes. $V$ is the vocabulary, and $C(w,x)$ is the number of times a word $w$ appears in document $x$. The joint probability distribution is: $$ P(x, y) = P(y)\prod P(w\mid y),\quad w\in V$$ I have to derive the update equations for the EM algorithm which maximizes the log likelihood of the incomplete data. The parameters are $P(y)$ and $P(w\mid y)$. But I don't understand where to start. For example, how is the $Q$ equation for this problem formulated?


title: Dimension reduction issues in self-organizing maps (SOM)

text: Self organizing maps are claimed to be able to visualize/cluster high- dimensional data in a smaller dimensional space. I have some difficulties in understanding this statement. Consider a six-dimensional data set; the codebook vector/reference vector is also six-dimensional. According to the SOM algorithm, updating these reference vectors is also conducted in the six-dimensional vector space. If we are considering a two dimensional map, how should I understand the map between the six-dimensional data space and two-dimensional map space?


title: Comparing two methods of sampling from a bivariate normal distribution

text: Is it possible to have two different methods of sampling from a bivariate normal distribution with a non-identity correlation matrix in such a way that one method would "consistently" result in a sample whose sample correlation matrix is closer to the correlation matrix of the distribution? I do not have a precise definition of what "consistently" means. I heard from a professional statistician that a method which gives a sample whose correlation is always very close to that of the underlying distribution may not be desirable, because there should be some randomness in the sample correlation matrices computed from different samples using the same method. I found this statement surprising as I would have expected that it is better to have a sample correlation as close as possible to the correlation of the original normal distribution. But I am not knowledgeable enough in this area to argue against that statement. EDIT: The "sampling" in my discussion involved using a Gaussian copula and a correlation matrix to generate the sample. In this case, would I be doing something "bad" if I computed many samples in a loop and returned the one with the closest match to the desired correlation as _the_ sample.


title: What variables do I keep in the model? multiple linear regression

text: What do I do when none of the predictor variables are significant? All the p-values are greater than the significance level of 0.05.


title: Given a mean and 95% confidence interval, do I need to know the sample size to calculate the standard deviation?

text: This question regards the basic statistics of a normal distribution, but I can't figure it out. I have been given the mean and 95% confidence intervals for a distribution, but would like to know the standard deviation. In my example: $$ \mu=53.4\quad 95\%\ c.i.=(52.3, 54.3) $$ I had thought that the solution for $\sigma$ would be something like: $$54.3=53.4+(SE*1.96)$$ $$SE=(54.3-53.4)/1.96=0.46$$ and then, $$SE=\frac{\sigma}{\sqrt{n}}$$ $$\sigma=0.46*\sqrt{n}$$ So, if I don't know the $n$, is this possible?


title: How to combine many Goodman-Kruskal gamma measures?

text: I have 5 samples (from the same population). I want to calculate the Goodman- Kruskal gamma ($\gamma$) for the population, based on the gamma of the 5 samples ($G_i$). Each data has 2 ordinal variables (X and Y), but the X's and Y's of different samples are not comparable, so I cannot put all the samples in the same bag and calculate the gamma for this new set. I can calculate the $G_i$ for each sample, and the standard error for the $G_i$ and I need to combine these $G_i$. Furthermore, I would like to test the hypothesis $ \gamma \ne 0$ and thus I need the standard error for the combined G. Just to add information, I am working in R, and I use the function `GKgamma` from the vcdExtra package to calculate the gamma itself and the standard error.


title: How to compute the maximum a posteriori probability (MAP) estimate with / without a prior

text: I am a newbie in this area so I hope someone could explain the following problem to me in plain English. Assume I want to use MAP to estimate some parameters on the basis of some observations. I know the method of computing MAP is: $$ \theta(x) = {\rm argmax} \ f(X|\theta) g(\theta) $$ where $g$ is the prior. However, I cannot find any answers online on how to compute this using a real world example. So here is my proposed question: Assume you asked 100 people of who they are going to vote for in an election (out of 2 candidates A and B), and assume the end result is 60% of them saying they will vote for A. How do you estimate the result of an election using MAP if:   1. candidate A is known to have a popularity of 40% and candidate B 60% (assume this to be the prior distribution)    2. the popularity is unknown. I also looked at this answer but I'm still confused: Example of maximum a posteriori estimation


title: Test for granger causality after fitting a GARCH(1,1)

text: I have two time series, where i wish to test for Granger causality of lagged values of $x$ on $y$, $y$ is changed to "rate-return" and $x$ is the positive or negative "rate-return", that is everywhere positive, and 0 where negative, where both variables are stationary. I fit the model by some autoregressive order p (using 2 for shorthand): \begin{equation} y_t = \beta_0 + \beta_1 y_{t-1} + \beta_2 y_{t-2} + \gamma_1 x_{t-1} + \gamma_2 x_{t-2} + \epsilon_t \end{equation} Where i wish to test $H_0:$ $\gamma_1 = \gamma_2 = 0$ Checking the $\epsilon_t$ term for normality, reject the hypothesis of normally distributed errors. Perform a heteroskedasticity test with IID. Fail to reject the null hypothesis of constant variance (p = 0.07). Test for autocorrelation using Breuch-Godfrey LM, and Durbin's alternative test, where both fail to reject the null hypothesis of no serial correlation. Perform a ARCH-LM test, this fails to accept the null hypothesis of no ARCH effects. To the question: After fitting a GARCH(1,1) (Which was found by theory to be most accurate using this type of time-series data, where i fit the variables to "rate-return"), of the form : \begin{equation} y_t = \delta_0 + \delta_1 y_{t-1} + \delta_2 y_{t-2} + \rho_1 x_{t-1} + \rho_2 x_{t-2} + u_t \end{equation} Finding the ARCH and GARCH errors, and interpreting them. Can i still do a granger causality test on this model,$H_0:$ $\rho_1 =\rho_2 = 0$? Or should i disregard the ARCH-LM test, and fit the model using newey west s.e? I have to point out that fitting the 3 different models, only GARCH find that X granger causes Y $(p= 0.0473)$, while newey west se and normal ols fails to reject the null of X does not granger cause Y $(p= 0.3)$.


title: How to interpret results from an experiment where covariate is influenced by the experimental setup?

text: I am doing study investigating the effect of company's heritage on consumer attitude towards the company/brand. For this purpose I have created two texts:   * one revealing the history of the company,    * and the other one - the control text - revealing some general information about the firm.  My results indicated that the two texts differ not only regarding heritage perception, but also regarding information credibility and believability. So the text with history details is also perceived as more credible. And this is also said to influence the attitude. My questions are:   1. Can I use the perceived text credibility as covariate in the ANCOVA, although it is measured after the treatment and caused by it?   2. How can I interpret the ANCOVA output in this case? Can I separate both effects - the one of the manipulation and the one of the covariate?    3. Would mediation analysis be one possible solution to check whether the effect of the treatment on my dependent variable is fully mediated by text credibility?


title: Is chi square the best approach for looking at number of cases by year

text: I have the number of people diagnosed with a condition in each quintile of deprivation over five different years. We are interested in whether numbers of diagnoses are going up faster in more deprived quintiles. The data looks like this: ![enter image description here](http://i.stack.imgur.com/QZwzV.png) It looks to me very much like it rises much faster in the more deprived quintiles, on the left (Q1 and Q2). I thought chi square was the best approach but:               test = structure(list(Q1 = c(98L, 109L, 263L, 323L, 312L),                           Q2 = c(90L, 113L, 199L, 237L, 247L),                           Q3 = c(70L, 83L, 133L, 166L, 182L),                     Q4 = c(20L, 39L, 60L, 87L, 90L),                           Q5 = c(38L, 50L, 75L, 101L, 115L)),                      .Names = c("Q1", "Q2", "Q3", "Q4", "Q5"),                      class = "data.frame", row.names =                        c("2008/09", "2009/10", "2010/11", "2011/12", "2012/13"))          chisq.test(test)              Pearson's Chi-squared test          data:  test     X-squared = 17.285, df = 16, p-value = 0.3674      The thing that I can see is missing is that the years are really ordered, we are looking at an increase, but of course the chi-square just treats them as nominal. I don't think I can use logistic regression because I'm just looking at a caseload- I don't have loads of people who don't have the disorder to compare with, just the increase over time. Or is it the right approach and it really just isn't significant? Thanks.


title: What happens when you initialize EM with a consistent estimate?

text: I have a certain family of models with latent random variables (not observed in the data), which I have a consistent estimator for. I now run EM on top of it, meaning, I get the consistent estimation, and initialize EM with it and run it for another few iterations. I know that EM has problem with local maxima, and also specifically for the model family I am handling as well. But, however, I think when you initialize EM with a consistent estimate, you can only improve your estimate (make it more efficient) and you won't have issues with local maxima. Is there anything in the literature to support that?


title: Trees and Cross Validation - # misclass

text: I am cross validating a classification tree and am able to plot the number of observations misclassed by different sizes of trees. My question is, what does it mean to return "the" number misclassed, for a given tree size, when there were _k_ different runs of that sized tree (where each of the _k_ runs has a presumably different number misclassed)? Several texts that I have read say that it is the average over the _k_ folds that should be returned for each size of tree, but I do not think this is what I am getting since the numbers I see plotted for the "number misclassed" are always perfect integers. Here is an example:               library('ISLR')     attach(Carseats)          High = cut(Sales, c(-Inf, 8, Inf), labels=c("Small", "Large"))     Carseats = data.frame(Carseats, High)          set.seed(2)     train = sample(1:nrow(Carseats), 200)          library('tree')     tr0 = tree(High ~ . -Sales, data=Carseats, subset=train)          set.seed(3)     tr0.cv = cv.tree(tr0, FUN=prune.misclass)     plot.tree.sequence(tr0.cv)          tr0.cv$dev      #[1] 55 55 53 52 50 56 69 65 80     # These don't look like averages across k folds      **Update** From my example, here are the different tree sizes and their corresponding dev value (here, meaning number misclassed)               tr0.cv$size     [1] 19 17 14 13  9  7  3  2  1     tr0.cv$dev     [1] 55 55 53 52 50 56 69 65 80      So we have for a tree of size 19, the number misclassed is 55. So is this saying that across the 10 runs, the sum total of all misclasses was 55? Thus about 5.5 misclasses on average, for each individual tree fit of size 19? This seems suspect, because when I do fit a tree of size 19, I see the number misclassed as 21.               tr0.Prune19 = prune.misclass(tr0, best=19)     summary(tr0.Prune19)          Classification tree:     tree(formula = High ~ . - Sales, data = Carseats, subset = train)     Number of terminal nodes:  19      Residual mean deviance:  0.4282 = 77.51 / 181      Misclassification error rate: 0.105 = 21 / 200       Doing this 10 times I would expect the total sum of number misclassed to be roughly 200, which is quite different from the reported 55. Granted, I know across 10 different runs there will be some variability in the number misclassed, but this sounds like too large of a discrepancy. Am I missing something?


title: Birthday paradox with a (huge) twist: Probability of sharing exact same date of birth with partner?

text: I share the same birthdate as my boyfriend, same date but also same year, our births are seperated by merely 5 hours or so. I know that the chances of meeting someone who was born on the same date than me is fairly high and I know a few people with whom I share my birthday although for the little I've read about the birthday paradox, it doesn't take same year into account. We've argued before about the probabilities and I am still not satisfied. My point was that the chances are tiny if you consider the probabilities of being in a relationship (+ being successful at it for X amount of time). I find the amount of factors to take into account quite vast (up to a point, gender and age, availability, probabilities of separation in our region, etc.) Is it even possible to calculate the probabilities on something like this? How would you go about it?


title: Principal Component Analysis for portfolio replication

text: I'm trying to implement Principal Component Analysis (PCA) in a portfolio replication procedure. (The replication procedure looks like regression: there is a vector representing payoffs of an asset under different ecocnomic scenarios, I need to find a linear combination of vectors of another assets that fits the most closely to an initial asset). Coefficients I get from regression should tell me how many assets I need to buy and sell to get approximately the same payoff as the asset I'm trying to replicate. There is a problem with regression as the candidate asset matrix is ill- conditioned: assets show high correlation. There's a hope that orthogonal principal components could resolve this problem. I need to do a PCA decomposition of candidate assets matrix, take only first $n$ components, do optimization and get components coefficients. Then I need to transform the coefficients back into an original basis. Now the problem: PCA usually works with mean-centered data, but if I subtract means from the original data, I don't know how to interpret resulting coefficients in my case and don't know how to reverse the operation. So far I'm doing eigen-decomposition of a covariance matrix, then using eigenvectors to make an orthogonal transformation of the data that is not mean-centered. Then I'm running a regression (actually L1 norm optimization) to get coefficients and transform them back into an original basis. The results are not bad, but I can't stop thinking about the problem with mean- centering, if I'm doing it completely wrong. I was hoping to find a detailed math reasoning for this problem, but unfortunately failed. I'm very much a noob in this and my math skills are far from being good, so I would really appreciate your help if you can share some insights on the problem of mean-centering in the PCA.


title: Looking for ways to compare between coxph models

text: I'm running Cox proportional hazards regression in R, and would like to test the option of categorizing one of my continuous variables to factor (I'm aware of the loss of data issue, just checking). Another thing that I'd like to check is putting the difference between 2 cont. variables inside the regression instead of putting them both. [It's actually testing if the pulse pressure, e.g, the difference between the systolic and diastolic blood pressure is more significant than each of them separately] My question is: What is the best way to compare between the different variations of the regressions (lets assume that I'll use step() in each of the attempts). No missing values in the dataframe whatsoever. I'm pretty confused between AIC, R2 (of coxph) and concordance of coxph. Can anyone clear things up for me? Is there any other option of comparing between different models on the same data? Thanks!


title: "Fully Bayesian" vs "Bayesian"

text: I have been learning about Bayesian statistics, and I often have read in articles > "we adopt a Bayesian approach" or something similar. I also noticed, less often: > "we adopt a **fully** Bayesian approach" (my emphasis). Is there any difference between these approaches in any practical or theoretical sense ? FWIW, I am using the package `MCMCglmm` in R in case that is relevant.


title: Coefficients change signs

text: I have a dataset where there is a high degree of multicollinearity, with all variables correlating positively with each other and the dependent variable. However, on some of the models I run I get a couple of significant negative coefficients. Basically there are two coefficients that depending on what variables I include in the model, I can manipulate their signs. My understanding is that if the variance-covariance matrix only contain positive values, then all coefficients should also be positive. Is this correct?


title: What is a good textbook on statistical genetics for someone without a lot of biology background?

text: I am looking for a good book to study from to learn about statistical genetics. I have no background in biology, some learning in chemistry, reading some of Griffith's book on genetic analysis but I am interested in just statistical genetics. I have ordered the book reviewed in the most current edition of AmStat by the fellow from TAM; but I would like to know if there are any others to use as a backup. Is Lange too difficult?


title: Which test to use for comparing changes in computer source code?

text: I am not sure which test I have to use to analyze my data. The situation is as follows: I have two computer programs which are doing exactly the same. Both share a large amount of source code, but some parts were developed with different technologies. I want to test the effect of the different technology. I performed the same changes on both programs and measured some source code properties after every change. The changes were performed in a chain. I performed change 2 on the result of performing change 1 on the initial version, and so on. Now, with the help of the properties measured after every change, I want to test if the program which was written with technology A were easier to change than the program written with technology B. Which statistical test do you think is the best choice? I thought of the Student's T-Test for paired samples or the Wilcoxon signed-rank test. But I am unsure if theses tests are applicable to this situation.


title: AR1 Modelling using dlmModReg

text: I'm trying to model AR1 using `dlmModReg()`. The main purpose is to keep phi a variable so that if phi >1, I know that mean reversion is not occurring. Below is my code:               buildFun <- function(x) {dlmModReg(e1, addInt = FALSE, dV = exp(x[1]), dW = exp(x[2]))}     fit      <- dlmMLE(e0, parm = c(0,0), buildFun)     dlmRes   <- buildFun(fit$par)       ResFilt  <- dlmFilter(e, dlmRes)     R        <- unlist(dlmSvd2var(ResFilt$U.R, ResFilt$D.R))     Q        <- R + as.vector(V(dlmRes))       "e is some residues from OLS and I'm trying to model e using AR1" My Question:   * Is my Q calculated wrongly as I try to built a confidence interval for `ResFilt$f` and the Q values are at least order 20 bigger than `ResFilt$f`?


title: Gamma-normal distribution as prior

text: When I study the Bayesian econometrics, the book firstly introduces Gamma- Normal distribution as (conjugate) prior, then the posterior will have the same distribution as the prior. But my question is, why can't we use other simple distribution as prior for the introduction example, such as normal distribution, why do we have to use a such as strange 'compound' distribution?


title: Difference time series before Arima or within Arima

text: Is it better to difference a series (assuming it needs it) before using an Arima OR better to use the d parameter within Arima? I was surprised how different the fitted values are depending on which route is taken with the same model and data. Or am I doing something incorrectly?               install.packages("forecast")     library(forecast)          wineindT<-window(wineind, start=c(1987,1), end=c(1994,8))     wineindT_diff <-diff(wineindT)          #coefficients and other measures are similar     modA<-Arima(wineindT,order=c(1,1,0))     summary(modA)     modB<-Arima(wineindT_diff,order=c(1,0,0))     summary(modB)          #fitted values from modA     A<-forecast.Arima(modA,1)$fitted          #fitted from modB, setting initial value to the first value in the original series     B<-diffinv(forecast.Arima(modB,1)$fitted,xi=wineindT[1])               plot(A, col="red")     lines(B, col="blue")      **ADD:** Please note I am differencing the series once and fitting arima (1,0,0) then I am fitting arima (1,1,0) to the original series. I am (I think) reversing the differencing on the fitted values for the arima(1,0,0) on the differenced file. I am comparing the fitted values - not the predictions. Here is the plot (red is arima(1,1,0) and blue is the arima (1,0,0) on the differenced series after changing back to the original scale) : ![enter image description here](http://i.stack.imgur.com/mQjAb.jpg) **Response to Dr. Hyndman's Answer:** 1) Can you illustrate in R code what I would need to do in order to get the two fitted values (and presumably forecasts) to match (allowing for small difference due to your first point in your answer) between Arima (1,1,0) and Arima(1,0,0) on the manually differenced series? I assume this has to do with the mean not being included in modA, but I am not entirely sure how to proceed. 2) Regarding your #3. I know I am missing the obvious, but are not $\hat{X}_t = X_{t-1} + \phi(X_{t-1}-X_{t-2}) $ and $\hat{Y}_t = \phi (X_{t-1}-X_{t-2})$ the same when $\hat{Y}_t$ is defined as $\hat{X}_t - X_{t-1}$? Are you saying I am "undifferencing" incorrectly?


title: Using STAN (related to BUGS/JAGS) to do linear regression with with ARMA(1,1) noise?

text: EDIT: I've modified my STAN code and it looks like I am getting numbers close to using R's `arima`. The original code, now moved to the end, was incorrect. I've been using STAN for simple things, but want to do a linear regression that has ARMA(1,1) noise. I've heard that this is basically equivalent to AR(1) + White Noise. I've made an attempt that gives reasonable numbers but I can't tell if I'm actually doing what I think I'm doing. The reason I'm attempting to do this is that the data is a climate temperature time series, which is thus serially correlated, and so I believe OLS regression -- used to calculate a trend -- will underestimate the uncertainty in the regression, which is reflected in the SE's of the slope and intercept coefficients. (And also in the residual error?) My code, using `rstan`:               temps <-         structure(c(-0.59, -0.17, 0.05, -0.7, -0.27, -0.94, -0.69, -0.96,          -0.58, -0.35, -0.58, -0.54, -0.48, -1.41, -0.82, -0.73, -0.48,          -0.37, -0.07, -0.16, -0.58, -0.43, -0.16, -0.19, -0.81, -0.37,          -0.52, -0.55, -0.51, -0.85, -0.43, -0.72, -0.43, -0.63, 0.16,          -0.26, -0.14, -0.48, -0.61, -0.36, -0.05, 0.22, -0.34, -0.23,          -0.2, -0.18, 0.51, -0.2, 0.28, -0.53, -0.07, 0.09, 0.45, -0.27,          -0.12, -0.35, -0.21, 0.11, 0.37, 0.09, -0.18, 0.14, 0.35, -0.16,          0.62, 0.04, 0.32, -0.12, 0.41, 0.28, -0.4, -0.34, 0.2, 0.26,          -0.27, 0.44, -0.11, -0.15, 0.68, 0.24, 0.12, 0.16, 0.33, 0.12,          -0.04, -0.01, -0.22, -0.14, -0.26, -0.42, -0.01, -0.1, -0.57,          0.16, -0.31, 0.12, 0.06, -0.19, 0.1, -0.01, 0.16, 0.81, -0.13,          0.53, 0.31, 0.06, 0.34, 0.23, 0.57, 0.07, 0.41, 0.51, 0.52, 0.39,          0.34, 0.73, 0.3, 0.38, 0.66, 0.6, 0.35, 0.55, 0.98, 0.89, 0.67,          0.86, 0.6, 1.31, 0.2, 0.75, 0.72, 0.53, 0.48), .Tsp = c(1880,          2012, 1), class = "ts")          stan.code1 <-     "data         {         int<lower=1> N ;         real x[N] ;         real y[N] ;         }          parameters         {         real alpha ;         real beta ;         real kappa ;              real<lower=0> sigma0 ;         }          model         {         real sigma[N] ;              alpha ~ cauchy (0, 5) ;         beta  ~ cauchy (0, 5) ;              kappa ~ gamma (1.2, 1) ;         sigma0 ~ gamma (3, 1) ;         sigma[1] <- sigma0 ;              y[1] ~ normal (alpha + beta * x[1], sigma[1]) ;              for (n in 2:N)             {             sigma[n] <- sigma0 + kappa * sigma[n-1] ;                  y[n] ~ normal (alpha + beta * x[n], sigma[n]) ;             }         }"              stan.list1 <- list (N=length (temps), x=1880:2012, y=temps)              stan.model1 <- stan (model_code=stan.code1, model_name="GISS NH Jan", data=stan.list1, iter=15000, chain=4)              print (stan.model1, digits_summary=8)      The alpha, beta, and sigma appear to be reasonable if I use them to throw lines onto the data using `abline`. QUESTIONS:   1. Have I really done a linear regression with AR(1) noise plus white noise? That is ARMA(1,1) noise? I believe $\kappa={{1-\theta}\over{1-\phi}}$.   2. Given that it is done properly -- either I lucked out, or someone will post a correct version of the STAN code -- how would I use kappa in calculating the larger uncertainty interval around the regression? EDIT: Original code I used, with `rstan`:               temps <-     structure(c(-0.59, -0.17, 0.05, -0.7, -0.27, -0.94, -0.69, -0.96,      -0.58, -0.35, -0.58, -0.54, -0.48, -1.41, -0.82, -0.73, -0.48,      -0.37, -0.07, -0.16, -0.58, -0.43, -0.16, -0.19, -0.81, -0.37,      -0.52, -0.55, -0.51, -0.85, -0.43, -0.72, -0.43, -0.63, 0.16,      -0.26, -0.14, -0.48, -0.61, -0.36, -0.05, 0.22, -0.34, -0.23,      -0.2, -0.18, 0.51, -0.2, 0.28, -0.53, -0.07, 0.09, 0.45, -0.27,      -0.12, -0.35, -0.21, 0.11, 0.37, 0.09, -0.18, 0.14, 0.35, -0.16,      0.62, 0.04, 0.32, -0.12, 0.41, 0.28, -0.4, -0.34, 0.2, 0.26,      -0.27, 0.44, -0.11, -0.15, 0.68, 0.24, 0.12, 0.16, 0.33, 0.12,      -0.04, -0.01, -0.22, -0.14, -0.26, -0.42, -0.01, -0.1, -0.57,      0.16, -0.31, 0.12, 0.06, -0.19, 0.1, -0.01, 0.16, 0.81, -0.13,      0.53, 0.31, 0.06, 0.34, 0.23, 0.57, 0.07, 0.41, 0.51, 0.52, 0.39,      0.34, 0.73, 0.3, 0.38, 0.66, 0.6, 0.35, 0.55, 0.98, 0.89, 0.67,      0.86, 0.6, 1.31, 0.2, 0.75, 0.72, 0.53, 0.48), .Tsp = c(1880,      2012, 1), class = "ts")          stan.code1 <-     "data         {         int<lower=1> N ;         real x[N] ;         real y[N] ;         }          parameters         {         real alpha ;         real beta ;         real kappa ;              real<lower=0> sigma ;         }          model         {         alpha ~ cauchy (0, 5) ;         beta  ~ cauchy (0, 5) ;         kappa ~ cauchy (0, 5) ;         sigma ~ gamma (3, 1) ;              y[1] ~ normal (alpha + beta * x[1], sigma) ;              for (n in 2:N)             {             y[n] ~ normal (alpha + beta * x[n] + kappa * y[n-1], sigma) ;             }         }"          stan.list1 <- list (N=length (temps), x=1880:2012, y=temps)          stan.model1 <- stan (model_code=stan.code1, model_name="GISS NH Jan", data=stan.list1, iter=15000, chain=4)          print (stan.model1, digits_summary=8)


title: How to interpret a QQ plot

text: I am working with a small dataset (21 observations) and have the following normal QQ plot in R: ![enter image description here](http://i.stack.imgur.com/kLcP7.jpg) Seeing that the plot does not support normality, what could I infer about the underlying distribution? It seems to me that a distribution more skewed to the right would be a better fit, is that right? Also, what other conclusions can we draw from the data?


title: Boosting in Weka with single PC

text: Given that I want to run Weka in a single PC, is it a good a idea to use boosting with a three layer feedforward neural network as a base learner? Why or why not? I think boosting requires a lot of computational power as it trains several base learners and averages them (am I right?), but how can I know it will be efficient or not taking into account the base learner (again for a single PC)?


title: Mann-Whitney test

text: I am carrying out a study in which I have $3$ groups $n= 2$, $n= 5$, $n=17$. Each group has a different inner ear abnormality, following the fitting of an ear implant I have measured their speech (ordinal data) and hearing scores (continuous data) at 1 and 2 years. The group with only $2$ patients in must be excluded from stats I believe. For the two remaining groups, most of the data are normally distributed and most variances are similar for speech (ordinal data), but not for the hearing scores (continuous data). I wonder if a Mann-Whitney $U$-test should be used for the following reasons:   1. Ordinal data DV: although the data are normally distributed and most variances are similar, ordinal means that only a M-W can be done.   2. Scale Data DV: the data are only normally distributed with homogeneity of variance in 2/3 of the measures, therefore again M-W? Or should a $t$-test be used (but doesn't a $t$-test assume homogeneity of variance)?


title: How to estimate a distribution?

text: Suppose I have estimated the discrete income distribution for Maryland, and the average distribution for the U.S. What is the algorithm or method to estimate the income distribution for Pennsylvania, for example, based on these two pieces of information? I know bringing in some additional information that helps explain the two different income distributions will help. How would I incorporate this information into the algorithm?


title: Meaning of marginal prior distribution

text: I have a variable $\mu$ with distribution, given $\phi$ of $ \mu|\phi \sim N(\phi, \sigma^2)$, and the distribution of $\phi$ is $\phi \sim N(\nu, \tau^2)$. (Here, $\sigma^2,\nu$ and $\tau^2$ are known). I need to find the marginal prior distribution of $\mu$, but am not sure what this means. Any help would be really appreciated. Thanks!


title: Curvature of a pencil's stroke

text: I'm trying to evaluate the curvature of an image of a pencil stroke based on the image's pixels and their shade of gray. I'm trying to get the curvature of the line at every point of the stroke. The part I'm not sure is the calculation of the curvature itself considering the noise. Here is how I'd process:   1. Image treatment. Higher the contrast and drop the pixels under a threshold of gray. I want to keep a scatterplot of the stroke and discard the white in the image.   2. Calculate the position of a pixel on the line. I'd use the diffusion space (R. R. Coifman and S. Lafon) to evaluate the distance between pixels, thus their position of the line.   3. Calculation of the curvature at each pixels. Here is the part I'm not so confident in. Most expression of the curvature of a line I've seen don't consider the noise (here the thickness of the stroke). I'd use an approximation based upon the curvature of a function according to it's derivative. I'd evaluate the derivative based on an average of the points far enough on the diffusion space.   4. Now I have for each pixel the distance on the diffusion space and their curvature. I could average that, but ideally there shouldn't be too much noise. So really, I'm not so sure how to go here. I know I can obtain something meaningful in the way I described, but the result will be highly dependent on the scale of the evaluation of the curvature. Having bias in the result is not so bad, because I'm more interested in comparing the curvature at different points.


title: Generate data based on existing data and their distribution

text: Consider the following dataset:               patient id, gender, age, systolic blood pressure, diastolic blood pressure     1, male, 66, 120, 80     1, male, 66, 119, 83     2, female, 45, 119, 90     3, female, 55, 120, 19      Using this dataset, I can calculate, for example, what the distribution of systolic blood pressure measurements is for males aged 60-70 etc. I am trying to figure out a way to create a separate synthetic dataset, with more patients than the original one, that still conforms to the distribution of data elements of the original dataset. I am guessing the problem has to be approached step-wise i.e. first generate synthetic data on all males, for age group 1, then age group 2 etc, then start with a new group and follow that through. I am thinking that this must be called something, but I am not sure what the proper term is so any hints would be appreciated. Are there any tools that could help with this?


title: Discrepancy between lagged lm() and lagsarlm()

text: If we have a spatial autoregressive process, we can estimate a model to control for the autoregression with a spatial lag, $$y=\rho W y + X\beta + \epsilon$$ Where $\rho$ is the strength of the spatial correlation, and $W$ is a matrix of spatial weights. The `spdep` package for R contains the `lagsarlm` command which is designed to estimate precisely this model. The package contains methods for creating the weights. But there seems to be some discrepancy between the model fit between `lagsarlm()` and `lm()` fitted to what should be a similar model. As an example, consider the example given with `?lagsarlm` in R.               library(spdep)     data(oldcol)     COL.lag <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,                     nb2listw(COL.nb, style="W"), method="eigen", quiet=TRUE)     summary(COL.lag)     Residuals:           Min        1Q    Median        3Q       Max      -37.68585  -5.35636   0.05421   6.02013  23.20555           Type: lag      Coefficients: (asymptotic standard errors)                   Estimate Std. Error z value  Pr(>|z|)     (Intercept) 45.079251   7.177347  6.2808 3.369e-10     INC         -1.031616   0.305143 -3.3808 0.0007229     HOVAL       -0.265926   0.088499 -3.0049 0.0026570          Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588     Asymptotic standard error: 0.11768         z-value: 3.6626, p-value: 0.00024962     Wald statistic: 13.415, p-value: 0.00024962      We can estimate what (I think) should be the same model by computing the actual spatial lag variable,               crime.lag <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME)     linearlag <- lm(CRIME ~ crime.lag + INC + HOVAL, data=COL.OLD)     Residuals:         Min      1Q  Median      3Q     Max      -38.644  -6.103   0.266   6.563  21.610           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 38.18099    9.21531   4.143 0.000149 ***     crime.lag    0.55733    0.15029   3.709 0.000570 ***     INC         -0.86584    0.35541  -2.436 0.018864 *       HOVAL       -0.26358    0.09136  -2.885 0.005986 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 10.12 on 45 degrees of freedom     Multiple R-squared: 0.6572, Adjusted R-squared: 0.6343      F-statistic: 28.75 on 3 and 45 DF,  p-value: 1.543e-10       The two models, which I think should be identical, are in fact significantly different from each other in every parameter and in model fit (with the `linearlag` model providing significantly lower AIC). Are there reasons why this should be? Why should I just not use the second model and abandon the special methods?


title: Simulate constrained normal on lower or upper bound in R

text: I'd like to generate random data from a constrained normal distribution using R. For example I might want to simulate a variable from a normal distribution with `mean=3, sd= 2` and any values larger than 5 are resampled from the same normal distribution. Thus, for the general function, I could do the following.               rnorm(n=100, mean=3, sd=2)      I then had a few thoughts:   * Iterate an `ifelse` function with a loop that repeats until all values are constrained to lie within the bounds.    * Simulate many more values than required and takes the first `n` that satisfy the constraint.   * Avoid vectorised normal variable simulators and instead use a for loop with a do while inside to simulate each observation one at a time and loop where required. All of the above seem a little bit clunky. ### Question   * What is a simple way to simulate a constrained random normal variable in R from normal with mean = 3,sd = 2 and max = 5?   * More generally what is a good general way of incorporating constraints into simulated variables in R?


title: P-value 2-way ANOVA vs 1-way ANOVA

text: I am comparing several treatments (8) at three time-points. For ease of analysis, I combined these data on a a two-way ANOVA table in Prism and ran the analysis. The first two time-points had stat sig differences between grouped variables, and the latter none. When I reran the exact same data, one time-point at a time, the second and third groups had stat sig diffs and none in the first. The strength of the significance was typically higher in one-way than two-way. I am not sure I understand why the sig found for the first group via 2-way disappeared in the 1-way test. Do I report both stating the test the values were sig in? Do I choose one or the other and present those results? Can anyone clarify? I'm using Prism v5.


title: Formulating a model statement using lmer for repeated measures in a blocked design in R

text: I am in the midst of analyzing data that I thought were rather simplistic, but I find myself in need of guidance. The function `lmer()` has been suggested, but I'm unclear about groups to be included when identifying random effects despite consulting every resource I could get my hands on. I think I know just enough about R to get myself into trouble, so please accept my apologies in advance if this should be obvious! The design: Three blocks [to which experimental units & treatments were randomly assigned], each containing every combination of 3 species (`spp`) and 4 treatments (`trt`). Response (`PN`) was measured weekly for 10 weeks (`Week`). I expect the species will differ, so - to simplify - I am running a separate model for each species. Within each of the `spp` x `trt` "tubs", I have three units which allowed a mean response (`meanPN`) and variance (`varPN`), so I'm using `meanPN` as my response and will weight the model by `1/varPN`. my model:               meanPN ~ Trt + Block + Time (+ interactions) + error      maybe the statement (excluding interactions for now as I try to wrap my mind around this) would be:               lmer(meanPN ~ Trt + (Trt|Week) + (1|Block) , data=SumExpt, weights=1/varPN)      Can you help me to understand why I would use `(1|Week)` versus `(Trt|Week)` \- or perhaps, when would each case be appropriate? Also, does one need to identify `Week` as the repeated variable, or is R interpreting that based on the pipe character? Finally, if I understand what I've read correctly, the interactions are handled in the denominators calculating the F-statistic. Is that proper understanding?


title: Calculating user preferences based on purchases: how to incorporate different variables for comparison

text: I am trying to come up with a metric to calculate user preference correlation for my final project (a web-shop for shoes) at school. Originally I intended to include user ratings and use Pearson's correlation coefficient but building a rating system might be part of the exam so I have to leave that to the side for now. I decided to compare purchases directly, taking binary values based on the number of purchases from the first user. This, of course, isn't terribly accurate. So, after looking at the various attributes of the shoes, I am trying to incorporate preferences in brands and style. For each shoe that is bought by both a 1 is awarded (0 if the shoe is bought only by user1), for brands and styles I was thinking 0.8 and 0.6 respectively. But then I lose the thread. I am not sure how to combine the variables in a relevant way. I would really appreciate any hints or tips. (Also, this is a part of the project that I decided to build in myself and is most definitely not required. I am not attempting to cheat by asking here; I just want to learn.)


title: Is maximin a machine learning algorithm?

text: Is maximin a machine learning algorithm, or simply another concept which is used for increasing the probability of winning?


title: Determining confidence intervals: using partial information on possible outcomes

text: Let's say we have a mathematical model that provides the probability of finding oil at a location in terms of a system of 10 bins with probabilities going from very low, say 2%, to 20% for the best bin. Let's also assume that _one_ of the measures we use to monitor the performance of this model is the Powerstat (Gini coefficient). [To be sure: This question is not about the Powerstat's suitability for this task]. Presumably, it would provide an indication how well the model is able to differentiate future successful oil strikes from failures. Let's say we have had a drilling campaign with 1000 trials. Before drilling, those 1000 locations were divided by the model over the 10 bins, for instance like {50,75,100,125,150,150,125,100,75,50}. Let's say we have found 100 oil wells among those 1000 trials distributed like {0, 2, 6, 14, 13, 18, 11, 13, 10, 8}. We can now calculate the Powerstat given this result. We want to determine whether or not the Powerstat of the model for this given data set of1000 drillings significantly deviates from what we would expect to find if the model works as described above. We could exactly calculate (at least in theory) the Powerstat for all possible outcomes and their probabilities or estimate that by means of Monte Carlo simulation and determine a confidence interval in this way. The main discussion we have is whether in this process we should use the information that we have in fact found 100 oil wells or whether we should not. So, in the MC case, should we generate random drilling outcomes based on a binomial model, the bin probabilities and the bin populations and throw away all results not totaling to 100 or should we use all samples regardless of their total number of oil strikes? In one view, the 100 oil well number is a more or less arbitrary outcome of the campaign and is just one realization of a large range of conceivable numbers (in fact anything between 0 and 1000, with varying probabilities). Using the given total of 100 would be just as arbitrary using a bit of information as would be, for instance, using information on bin 10 (containing 8 oil strikes), or the totals of bin 1-5 and 6-10 (each of which would lead to different intervals). The other view is that we should use all information we have; that having an outcome of 100 says something already about the model quality and that we should take that into account. So, the question is: which view is correct? PS: I had difficulties coming up with a descriptive title. Any improvement is welcome.


title: Comparison of two values in a statistically sound manner

text: For a university paper, I have to compare whether two program versions are statistically significantly different. The comparison in terms of runtime performance is straightforward -- I determine whether the confidence intervals for the two sets of measurements overlap. However, in the case of memory footprint, I have no sets of measurements, because the memory footprint does not vary from measurement to measurement (e.g., version 1 has a constant footprint of 920 kb and version 2 consumes 1040 kb). Is there any way to compare both values in a statistically sound manner?


title: When is Generalized Belief Propagation exact?

text: It is known that belief propagation gives exact result on trees, are there interesting examples when Generalized Belief Propagation is exact? ( **edit** junction tree is not interesting because it is exactly solvable without GBP) On the surface, Belief Propagation passes messages between cliques and separators, whereas GBP allows more general region hierarchy. That helps with convergence rate, but I wonder if this also extends the class of exactly solvable inference problems. **Edit** : as Thomas Minka points out, junction tree algorithm can be viewed as a version of generalized belief propagation. But it can also be viewed as a version of (cluster)belief propagation. What I'm wondering specifically is if GBP can give exact solution for any problem for which BP can't. The motivation is that with exact solution GBP gives result in a finite number of steps and you can view the result as a kind of algebraic factorization of the problem, in the spirit of Generalized Distributive Law paper


I would like to gather input from people in the field about the Yates continuity correction for 2 x 2 contingency tables. The Wikipedia article mentions it may adjust too far, and is thus only used in a limited sense. The related post here doesn't offer much further insight. So to the people who use these tests regularly, what are your thoughts? Is it better to use the correction or not? And a real world example which would yield different results at the 95% confidence level. Note this was a homework problem, but our class does not deal with the Yates continuity correction at all, so sleep easy knowing you aren't doing my homework for me.               samp <- matrix(c(13, 12, 15, 3), byrow = TRUE, ncol = 2)     colnames(samp) <- c("No", "Yes")     rownames(samp) <- c("Female", "Male")          chisq.test(samp, correct = TRUE)     chisq.test(samp, correct = FALSE)      Thanks!


title: Calculation of incidence rate for epidemiological study -- prevalence rate this time

text: Thanks for all the answer for the question Calculation of incidence rate for epidemiological study in hospital. And here come's the second part of the question: What about the prevalence rate then? I have read `The new public health` as suggested by Chi, the book says that prevalence is usually not available when using ordinary incidence rate, but I saw another formula here: `total case count in that period of time`/`total patient bed days during that period of time` It puzzled me again, what is it? I have never heard of prevalence calculated using denominator as `patient-bed days`. Thanks!


title: Formal definition of asymptotically uncorrelated

text: I have a big problem grasping the definition of "quickly" in the following sentence: "We need it to be weakly dependent, which occurs when the correlation between $x_{t}$ and $x_{t+h}$ tends towards zero sufficiently quickly as $h\rightarrow \infty$." So in the text it said that formally it means: "$corr(x_{t},x_{t+h}) \rightarrow 0$ as $h\rightarrow \infty$". Given that the correlation of a random walk process is $corr(y_{t}, y_{t+h} = \sqrt{t/(t+h)}$ and is violating this condition, I am wondering what does quickly means? Or am I making a mistake? Thanks in advance!


title: $P(X_1 < \min(X_i,\ldots, X_n))$ across different normal random variables

text: I have a set of mutually independent normal distributions $X_1$ to $X_5$ (with means and standard deviations) which represent finishing times for swimmers over a certain distance. The actual data is as follows: $$X_1(60, 3.0)$$ $$X_2(61, 1.0)$$ $$X_3(58, 2.3)$$ $$X_4(63, 2.4)$$ $$X_5(61, 1.7)$$ So swimmer 1 ($X_1$) has a mean finishing time of 60 seconds with a standard deviation of 3.0 seconds. Question 1: What is the probability of an event where $X_i$ finishes first. e.g. $$P(X_1 \lt X_i, i=2,\ldots,n)$$ Question 2: If I calculate this for all swimmers, can I simply order the results to determine the most probable finishing order? This is not homework. Based on the answers to this Cross Validated question, I have tried to solve this problem _based on_ the first answer. i.e. $$\Pr(X_1 \le X_i, i=2,\ldots,n) = \int_{-\infty}^{\infty} \phi_1(t) [1 - \Phi_2(t)]\cdots[1 - \Phi_n(t)]dt$$ Where $\phi_i$ is the PDF of $X_i$ and $\Phi_i$ is its CDF. Based on this formula, the results I obtained were: $$\Pr(X_1 \le X_i, i=2\ldots5) = 0.259653$$ $$\Pr(X_2 \le X_i, i=1,3\ldots5) = 0.214375$$ $$\Pr(X_3 \le X_i, i=1\ldots2, 4\ldots5) = 0.611999$$ $$\Pr(X_4 \le X_i, i=1\ldots3, 5) = 0.0263479$$ $$\Pr(X_5 \le X_i, i=1\ldots4) = 0.0697597$$ However, the probabilities add to 1.182135 when they should add to 1.0. I’m not sure if the formula is incorrect or my implementation of the integral (I used Excel and the trapezoidal method). I also attempted to solve the problem using Dillip’s method (from the above mentioned question) as follows: \begin{align*} P(X_1 < \max X_i) &= P\\{(X_1 < X_2) \cup \cdots \cup (X_1 < X_n)\\\ &\leq \sum_{i=2}^n P(X_1 < X_i)\\\ &= \sum_{i=2}^n Q\left(\frac{\mu_1 - \mu_i}{\sqrt{\sigma_1^2 + \sigma_i^2}}\right) \end{align*} However, the probability results were much greater than 1 in most cases so abandoned this approach. By the way, what exactly does $\max X_i$ mean? Any assistance in calculating the probability would be appreciated.


title: Analyzing repeated rank data.

text: I have a data set of N people, T items. Let's say N=100, and T=10. Each person goes through the following exercise.   1. She is shown 2 random items from the set of T=10, and ranks them as rank 1 and 2.    2. She is next shown 2 more random items from the remaining 8 out of 10 items, and ranks as rank 1 and 2.  At the end, the data set is of size 100x10, where each row has 4 numeric entries (two of which will be 1, and the other two will be 2) and 6 empty entries. My goal is to compare the 10 items against one another, and come up with an estimated rank value for a given item. What is the best way to analyze such data ? Thank you.


title: Understanding d-separation theory in causal Bayesian networks

text: I am trying to understand the d-Separation logic in Causal Bayesian Networks. I know how the algorithm works, but I don't exactly understand **_why_** the "flow of information" works as stated in the algorithm. ![enter image description here](http://i.stack.imgur.com/gPnim.png) For example in the graph above, lets think that we are only given X and no other variable has been observed. Then according to the rules of d-separation, the information flow from X to D:   1. X influences A, which is $P(A)\neq P(A|X)$. This is OK, since A causes X and if we know about the effect X, this affects our belief about the cause A. Information flows.   2. X influences B,which is $P(B)\neq P(B|X)$. This is OK, since A has been changed by our knowledge about X, the change at A can influence our beliefs about its cause, B, as well.   3. X influences C,which is $P(C)\neq P(C|X)$. This is OK because we know that B is biased by our knowledge about its indirect effect, X, and since B is biased by X, this will influence B's all direct and indirect effects. C is a direct effect of B and it is influenced by our knowledge about X. Well, up to this point, everything is OK for me since the flow of the information occurs according to intuitive cause-effect relationships. But I don't get the special behavior of so called "V-structures" or "Colliders" in this scheme. According to the d-Separation theory, B and D are the common causes of C in the graph above and it says that if we did not observe C or any of its descendants, the flow information from X is blocked at C. Well, OK, but my question is **_why?_** From the three steps above, started from X, we saw that C is influenced by our knowledge about X and the information flow occurred according to the cause- effect relationship. The d-Separation theory says that we cannot go from C to D since C is not observed. But I think that since we know that C is biased and D is a cause of C, D should be affected as well while the theory says the opposite. I am clearly missing something in my thinking pattern but can't see what it is. So I need an explanation of why the flow of information blocked at C, if C is not observed.


title: Nonlinear principal component analysis MATLAB code

text: I'm working on multivariate process monitoring using NLPCA method. I need MATLAB code to implement this method.


title: Coding an interaction between a nominal and a continuous predictor for logistic regression in MATLAB

text: So our data is structured as follows: We have $M$ participants, each participant can be categorized into 3 groups ( _G_ $\in {A,B,C}$), and for each participant we have $N$ samples of a continuous variable. And we are trying to predict values that are either 0 or 1. How would we use matlab to test for an interaction between the continuous variable and the categorical variable in predicting these values?


title: How to best visualize one-sample test?

text: We are currently writing a paper with several one-sample Wilcoxon tests. While visualizing two-sample tests is easy via _boxplots_ , I was wondering whether there is any **good way to visualize one-sample test results?**               # Example data     pd <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64,             0.73, 1.46, 1.15, 0.88, 0.90, 0.74, 1.21)          wilcox.test(pd, mu = 1.1)          #   Wilcoxon signed rank test     #     # data:  pd     # V = 72, p-value = 0.5245     # alternative hypothesis: true location is not equal to 1.1      ...and also: I would like to get Z-value instead of V-value. I know that if I use `coin` package instead of basic `stats` I will have z-values, but `coin` package seems not to be able perform one-sample Wilcoxon test.


title: Plot a linear discriminant analysis

text: Is there a package in `R` that can do this plot on LDA data graph plot: ![enter image description here](http://i.stack.imgur.com/x1wle.jpg) original question: How does linear discriminant analysis reduce the dimensions?