title: Combining $\chi^{2}$ tests

text: I have a problem with statistics. I'm looking at whether variation in the concentration of a certain substance affects how strongly insects are drawn towards it. I did this by baiting insect traps with the substance in different concentrations. I spread the traps out as pairs with the traps in the pair having different concentrations. Now the data I have looks a bit like this:               More : Less        6 : 0          1 : 0           3 : 1          1 : 0         15 : 3        What I need to do is to show that insects prefer higher concentrations of bait. Pooling the data and using $\chi^{2}$ is an option, but I don't know if it's the best option. Trapping was done in various locations during one month so there might be perturbation factors involved. Thanks in advance! EDIT: That is only a small portion of my full data. In total I have around 40 pairs of values to compare and around 300 insects caught. That is just a small, but representative portion.


title: Why there is usually no correction for multiple comparisons in Bayesian statistics

text: How would you explain to an only moderately statistically literate scientific reviewer that adjustment for multiple comparisons is not necessary when sampling from a Bayesian posterior?


title: How to train SVM correctly on a 1D dataset

text: I am trying to use svmtrain (Statistic Toolbox) to train a linear (2 class) SVM on a 1D feature vectors. The features are not fully separable and the classes intersect. The naive approach would be to search for the optimal threshold for a simple classification, but one can expect to get similar (or even the same) results using a linear SVM. I am using all the default parameters and the SVM gives a different threshold and performs much worse than searching for a threshold (0.73 accuracy with threshold search vs. 0.66 accuracy with SVM). The motivation for the process is as follows: While 1D vectors can be separated into 2 classes using a simple threshold, when the dimension is > 1, the threshold approach is not feasible and requires an SVM (or similar methods) to find a separation. To make sure I configure and use the SVM correctly I first expect to get very similar results using threshold search and SVM on a 1D vectors. Would be happy to know what am I doing wrong to get this kind of difference in performance.


title: Regression model where output is a probability

text: I am trying to fit a model, and have the suspicion that what I am doing is not quite right. The data tracks what proportion of people made a decision, and what factors were active when they made their decision, i.e. something like this:   [1,0,1,0, 23%]   [1,1,0,1, 41%]   etc...   I also know how big each group is. The goal is to predict the percentage, based on the binary input. My initial thought was, the model cannot be a straight linear combinations, if only because the output is bounded. It is not exactly a logistic regression either, because the output is not a label, but the average for each group. My first take has been to transform the output in a fashion similar to logistic regression, into log (p/(1-p)), and fit a linear regression. This has given me some decent results, but I have a nagging feeling this is not quite right. Besides the output transformation, I am also concerned that treating the input as numbers, when they really represent binary values, is probably not the best way to go.   So my question is, if this is not the right way to go, are there models that address this specific type of situation? What should I be looking for? [Edit for clarification] From the comments / responses it seems my description of the data was a bit lacking, so here is a bit more, as well as why I am uncertain about using logistic regression. I'll illustrate on similar data. Suppose products had a set of binary characteristics, and various products were presented to customers, and the result (buy/no buy) was recorded. Then the raw dataset would look like:   F1,F2,F3,... Fn, Buy/No Buy   1, 1, 0, .. 0, 1   1, 1, 1, .. 1, 0   where each row is a specific product, and what the customer did. Now I can aggregate these by identical products, which will have the same characteristics, and simply record the proportion of buys, as well as the number of customers that were presented with that choice. This is essentially what I have.   I could disaggregate back into the original dataset, and run a logistic regression on it, but the groups themselves are very large, and of very different sizes. On top of that, I have 2 problems. First, I could reconstruct synthetic groups that have the same proportion as the original (i.e. if 4% purchased, construct 4 Buy, 96 No Buy rows), but the the Buy/No Buy ratios are pretty small, which would mean reconstructing large groups to properly approximate. Second, the groups have very different size, and I believe the group composition in the complete sample should have similar composition to the original one, which would mean creating potentially very large groups. Which is why I was essentially wondering if there was a way to work off the much smaller dataset directly, without having to reconstruct an artificial gigantic dataset.   My current approach has been to use a gradient descent approach, weighting observations by the group sizes, but I was wondering if there was a smarter way to handle this!


title: How to test if a value is over-represented in one sample vs another

text: I have two multinomial data samples that both have `N` discrete categories. I know that a Kolmogorov-Smirnov test will let me know if the distributions of the two samples are different. But how can I statistically test if a particular category is significantly over-represented in one sample compared to the other?


title: SD larger than mean, non-negative scale

text: I was given an article reporting a study very similar to one my lab wishes to run. But, I noticed that for the variable of interest, Duration, the SDs are larger than the mean...since this is duration measured in minutes it can never be negative and this seems very strange to me. This happened in 2 studies reported, below is one. Beyond that, this is a mixed-design. Control v Treatment (between groups), and Time1,Time2,Time3 (repeat measures). Here are the means (SDs), N >200                                      Time1                Time2                  Time3      Control               15.1 (14.6)          14.4 (14.8)            13.3 (15.7)     Treatment             14.8 (13.2)          10.0 (12.2)            8.2 (9.9)      ...they ran an ANOVA and reported a p<.001. I was asked to use this as a basis for a power analysis to determine sample size for our study. I am pretty sure this indicates that the data is non- normal or has outliers and I don't feel comfortable determining sample size based on this. Am I just way off base?


title: Where to find data on the use of statistical methods in Psychology?

text: I'm trying to find out, how the use of statistical methods in Psychology have been changing through history. For example, are regression models more popular now than ten years ago? When did simple non-parametric statistics, such as Mann-Whitney U, become less used (they seem to be used much less now than in 50s)? **Do you know any papers or resources that provide such data?**


title: What is the meaning of using a compound variable in PCA

text: I am using PCA as a dimension reduction method before feeding the data into a classifier. I am typically selecting only the first two or three PCs. Typical data includes factors such as length, width, height, weight, energy output, and other key characteristics. Ultimately I am using only the first several PCA- scores instead of the original factors. **My question is** : If I were to include compound factors such as `area = width*length`, `energy per unit volume` or `surface area = (length*width + length*height + width*length)` in addition to the simple factors that go into the PCA algorithm would this be expected to add, subtract or be neutral in terms of the information expected to be contained in the first several principal components? What about including multiple variations of the same factor such as `length^2`, `length^3`?


title: Applying ARIMA in R

text: Small question as a follow up from this topic: Relationship between two time series: ARIMA I am trying to reproduce the above through R. In the first steps, IrishStat applies AR(1) on doubly differenced x and single differenced y. However, after applying this filter in R, we obtain the following:               > ar1_12 <- arima(df1$x, order=c(1,2,0))     > ar2_12 <- arima(df1$y, order=c(1,1,0))     > ar1_12          Call:     arima(x = df1$x, order = c(1, 2, 0))          Coefficients:     ar1     -0.3761     s.e. 0.1306          sigma^2 estimated as 0.3199: log likelihood = -42.53, aic = 89.06     > ar2_12          Call:     arima(x = df1$y, order = c(1, 1, 0))          Coefficients:     ar1     -0.4815     s.e. 0.1213          sigma^2 estimated as 0.0328: log likelihood = 14.64, aic = -25.29      What can we conclude from this? Could you provide maybe some elaboration on the thought process regarding the cross-correlative structure of these two arima filters? Why is it that both autoregressive factors are -.383 in the example by IrishStat?


title: Forecasting a time series based on a collection of other time series

text: I am new to Time Series analysis, but I read a lot of questions here which deal with how you can forecast values of a time series based on its history alone. Now, the data I have is of this kind : TimeSeries1:               timestamp value     ..        ..     ..        ..      Now, I have a huge number of other timeseries, which follow the same schema. Also the data is not seasonal. If plotted, the time series will look like a rectangular wave. What I want to do is forecast values of TimeSeries1, based on the aggregate behavior of all the remaining timeseries data that I have. I had a look at this question, but my problem involves incorporating a large number of other time-series rather than just one like in the above question. What are the possible approaches to this problem?


title: how to check for robustness for categorical variables in STATA

text: I performed logit regression on my variables. I have 10 variables, all are categorical varibales. After performing logit, I want to check for robustness. How do I got about it in STATA


title: Best way to represent data

text: What is the best way to represent the following data graphically? Can I use a histogram               Year    Output per Person       Capital Employed     2010    16.3 units/annum        £40,000 p.p     2011    15.1 units/annum        £38,000 p.p     2012    14.4 units /annum       £35,000 p.p     2013    11.7 units per annum    £33,000 p.p     2014    10.8 units per annum    £30,000 p.p


title: Should grades be assigned to students based on a normal distribution?

text: Assume you have a class of approximately 800 students and following a set of assessments each student has a raw grade.   * How should these raw grades be converted into a final grade?   * Is it a good idea to scale the raw grades to a normal distribution?


title: Identifiability in generalized linear random effect model?

text: Suppose I observe binary $Y_{ij}$ for $i = 1, ..., N$ and $j = 1, ..., J$ and I want to model $$\Pr(Y_{ij} = 1 \mid \lambda_{i}) = \Phi(\lambda_{ij}), \qquad [Y_{ij} \perp Y_{ij'} \mid \lambda_i]$$ where the vector $\lambda_{i} = (\lambda_{i1}, \ldots, \lambda_{iJ})$ is a multivariate-normal random effect, $\lambda_i \sim \mathcal N(\mu, \Sigma)$ and $\Phi(\cdot)$ is the probit link function (in general this could be any link function and we would run into approximately the same issues, but probit is easier to analyze with normal random effects). $\Sigma$ is not necessarily full-rank so this specification also covers random effects models of the form $\lambda_{ij} = x_j^T b_i$ where $\dim(b_i) \ll J$. My question is: under what conditions on $(\mu, \Sigma)$ are they identified? My concern is due to the fact that it is well-known that the related multivariate probit model is unidentified in the absence of further restrictions (such as requiring $\Sigma$ to be a correlation matrix). The model as written is equivalent to a multivariate probit with mean $\mu$ and variance $\mathbf I + \Sigma$ so similar concerns should apply here. For example, I feel confident heuristically that setting $$ \Sigma(\theta)_{jj'} = \theta_1 e^{-\theta_2|t_j - t_{j'}|} $$ for known $t$ leaves $\theta_1$ unidentified. On the other hand $\Sigma = X \Sigma_b X^T$ is a standard random effects covariance matrix and so is apparently identified as long as $\Sigma_b$ has small dimension. What can one say about (say) $$ \Sigma = \Sigma(\theta) + X\Sigma_b X^T $$ where $\Sigma(\theta)$ is defined as above?


title: Measuring quantization error for clustering - squared or not?

text: * When measuring the quantization error of a clustering, should the distance between samples/centroid be squared or not? I found both variants in the literature.   * Furthermore, is (squared) quantization error not the same as squared sum of errors?


title: Computing $(X^TX)^{-1}X^Ty$ in OLS

text: Let $A\in\mathbb{R}^{n \times n}$ be a dense symmetric positive-definite matrix (the $X^TX$ from here) and $b$ a vector in $\mathbb{R}^n$. I need to compute $A^{-1}b$. Two questions:   1. Could you recommend an efficient and numerically stable algorithm for computing $A^{-1}b$ for $n \approx 1000$?   2. Let $\tilde{A_i}$ denote the matrix obtained from $A$ by removing its $i$-th row and $i$-th column. Is there an algorithm that, having been allowed to pre-process $A$ in some way, would enable me to quickly compute $\tilde{A_i}^{-1}\tilde b$ for any $i \in \\{1,2,\ldots,n\\}$ and any $\tilde b \in \mathbb{R}^{n-1}$?


title: Is it logical to stand on the chance-line (50-50 %) when we don't know a-priori probabilities?

text: Let we have two hypotheses $H_0$ and $H_1$ and we don't know their _a-priori_ probabilities. If we wish to calculate the average probability of error, does it makes sense to assume 50-50 % chance of occurrence of any of two, in the context of Neyman-Pearson criteria?


title: Is PCA unstable under multicollinearity?

text: I know that in a regression situation, if you have a set of highly correlated variables this is usually "bad" because of the instability in the estimated coefficients (variance goes toward infinity as determinant goes towards zero). My question is whether this "badness" persists in a PCA situation. Do the coefficients/loadings/weights/eigenvectors for any particular PC become unstable/arbitrary/non-unique as the covariance matrix becomes singular? I am particularly interested in the case where only the first principal component is retained, and all others are dismissed as "noise" or "something else" or "unimportant". I don't think that it does, because you will just be left with a few principal components which have zero, or close to zero variance. Easy to see this isn't the case in the simple extreme case with 2 variables - suppose they are perfectly correlated. Then the first PC will be the exact linear relationship, and the second PC will be perpindicular to the first PC, with all PC values equal to zero for all observations (i.e. zero variance). Wondering if its more general.


title: Relationships between two variables

text: Comparing two variables, I came up with the following chart. the x, y pairs represent independent observations of data on the field. I've doen Pearson correlation on it and have found one of 0.6. My end goal is to establish a relationship between y and x such that y = f(x). What analaysis would you recommend to obtain some form ofa relationship between the two variables? ![Graph](http://koopics.com/ask_math_chart.jpg)


title: How can I obtain null and residual deviance/degrees of freedom for assessing model significance?

text: I am trying to determine model significance using:               1-pchisq(null deviance-residual deviance, null df- residual df)      I have 5 models:   * Four models were estimated with GLMs, which gave me null and residual DFs in the summary.    * The fifth model was estimated with the `lmer` function because of the nested structure of my data. This does not give me degrees of freedom for null or residuals.  How do I get my null and residual deviances and degrees of freedoms? Thanks!


title: What is Affine Transformation?

text: What is affine transformation? Which distribution families are closed under affine transformation?


title: Bias-variance tradeoff in the paired t-test

text: Suppose we have $K$ subjects and a treatment with two levels, "Before" and "After". A paired t-test is equivalent to fitting a fixed effects ANOVA: $Y = Subject + Treatment +\epsilon$ It is also equivalent Repeated Measures ANOVA or Mixed ANOVA, where Treatment is fixed and Subject is random. It tried all four methods for this simple dataset and the Treatment p-value is exactly the same. For the sake of this question, let us consider the fixed effects version only. The total number of observations is $2K$, while the number of parameters is: $ p = intercept + (K-1)$ subject effects + (2 - 1) treatment effects = $K + 1$ That is, there are about 2 observations per parameter for any $K$. To me it suggests that the model is likely to be very overfitted, unless there is a very substantial subject effect. In practice, how often have you seen the subject effect so large that it justifies pairing the observations by subject?


title: Mean squared error / Max error squared

text: I'm use this metric for variation that is similar to CV, but it is guaranteed to be on a 0 to 2 scale for positive numbers, where 2 is "all over the place" and 0 is "perfectly even". Since its dimensionless and bounded, it's useful to me in ways that CV really isn't. Does this have a standard name? It's the mean squared error divided by the maximum error squared. I use it all the time, but don't know if it's called anything:               function(v) {(sum((v-mean(v))**2))/((sum(v)-mean(v))**2)}      And a skewed version that deals with low-counts:               function(v) {(min(2,(max(0,sum(floor(abs(v-mean(v)))**2)-2*sum(v)))/max(1,((sum(v)-mean(v))**2-2*sum(v)))))}       It only reports things as "out of whack" when they are out of sync/not balanced enough. It works perfectly. When the data are fairly uniform, you get zero. When there's a spike of some kind (size of the spike relative to length of vector), you get a larger value. I need to explain _why_ I'm using it. I know that it works, empirically, in my situation, because we get better results than using things like maximum likelihood, or a cutoff with CV alone. Essentially it signals outliers, which are nicely removed. But is there a better way of describing these functions for dimensionless coverage summary statistics? IE: Is there a better function to use to detect a spike proportional to the size of the vector?


title: rebuild model based on what columns a record has

text: I've built a model in R using glm, and in the new dataset that I need to predict, there exist some new levels for the columns that are non-numeric. I know there are so many approaches to deal with missing values, but seems like most of them are trying to impute the missing values, in another way, it will try to fill the missing values in an "intelligent" way (KNN..etc.), however, here is my concern: Say my model is based on 10 columns that are all discrete (factors like manufacturer, product family name, grade, origin country.). In my new dataset, some product has certain columns and the other products have that few columns and the missing columns are inconsistent. And in an extreme scenario, some record only has 2 columns. My goal is to predict all of them. Instead of trying to impute the missing value, filling in tons of values by guessing, I actually prefer an approach which will rebuild a model based on what columns that certain record has and predict only using the columns they have. Here I wrote some code using the `mtcars` dataset, where I manually insert a new level in column `vs` and then rebuild a model without `vs` column to predict the record.   1. Is there any general mathematical method to help me solve my problem.   2. Is there any implementation in R / Python handling this properly, thanks a lot!


title: Quality of PLS Regression at different interaction levels

text: I am fairly new to multivariate statistics and have run into the following situation: I have a data set of 12 response sets based on a Likert scale (1-5), data which is commonly (in social research) treated as more-or-less interval scale / quantitative. Only about half of the responses show a normal distribution, many show a right-skewed distribution. The aim is to identify correlations within this response data set. Using XLSTAT for Excel, I started out with multiple regression and identified two main influential IVs (R2 over all variables = 0,32) for a particular DV (selected from among the 12 variables). I then re-ran the calculations for 2 and 3 interaction levels. While the R2 of the results understandably increased considerably (to 0,85 and 0,94), the variables that now were most influential (according to standardized coefficient value for significant results) changed with each interaction level. While for no interactions, multicollinearity was low (tolerance 0.6-0.8, VIF always under 2), with 2 or 3 interaction levels multicollinearity soared (2 levels: VIFs over 400, 3 levels: VIFs in the millions). Assuming that multicollinearity and probably suppressor variables were the problem, I then changed to PLS regression, again 11 IVs for one DV. With no interactions accounted for, I got a Q2 cumulated of 0.12, R2Y cum=0.26 and R2X cum=0.25. The most influential IV according to VIP value was the same as at the non-interactions multiple regression (with VIP value of 1.4). The same with 2 interaction levels again showed that same IV (VIP now 2.0) at the top (and the previously second-most significant IV still at second place), but Q2cum slightly down at 0.09, R2Y cum improved at 0.41, R2X cum lower at 0.09. Trying now for 3 interaction levels, the rank of most influential IVs shifted completely, with a combination of three IVs showing a VIP of 2.4, the previously most important IV still being in third place (VIP of 2.1). However, model quality is now terrible, with a Q2cum of -0.03 (R2Y/X cum about the same as for 2 interaction levels). What I would like to understand from this:   1. I assume that higher VIP values for increased interaction levels come from more complex modeling due to more interactions being accounted for, but the Q2cum for 3 interaction levels shows that model predictive quality has dropped. Why would this be the case? Which of the three PLS models (0/2/3 interaction levels) will serve me best, or does the low Q2 mean that none of them is significant at all? Assuming that Q2 and R2 were the same for all three, would the higher interaction model have the greatest explanatory power?   2. How should I interpret the model quality values (Q2/R2x/y cum)? In particular, is Q2cum (predictive power of model) the primary quality factor to consider, or are reasonably R2x/y cum levels (actual explanative power of model) sufficient? Since R2y reflects the variation of the DV explained by the IVs (if I understood correctly), then I assume that the R2y value is more significant for the quality of my model than R2x, correct?   3. How bad is it that my IVs do not all show a good normal distribution for all this? How robust is PLS regression in this case? Thanks for your patience with a newbie, I appreciate simple and understandable explanations, thanks a lot!!


title: Treating multiple observations per object

text: I am working on a project whose aim is to analyze the relationship between machine elements and their price. My data consists of thousands of machine elements, their price, as well as technical and non-technical attributes. Most of the data is categorical The technical data contains information such as weight, dimensions, material. Each element (identified by unique code) has distinct technical attributes. The non-technical data contain information such as supplier name, contract type, INCOTERMS, first day of validity of the contract, last day of validity, minimum order quantity, etc. In some cases the prices depend on the order quantity. One of the problems I am facing is that often each element has a different price depending on the financial data, i.e. element X will cost 100 dollars if Incoterms is A, 200 dollars if Incoterms is B, etc. In other words, there are rows that contain price information on the same element, but one of the columns has a different value and so the price is different. In other cases the price is 100 if 50-100 elements are ordered, 80 if 100-200 elements are ordered, 60 if 200-500, etc. I am planning to do some correlations as well as regression. I will probably also try data mining (using Rattle and R). I need advice on how to treat the (rather) similar observations that each object has in case of correlation, regression and in general for data mining. Should I try to select a single observation per element or do the analysis for “all in”. I guess the last option won’t work for correlation at all. Perhaps including a sample screenshot of the data with clarify my problem. So what is visible in the picture happens quite often: one object, different attributes, different price. I am wondering whether I should (try to) have only one row per object with only one set of attributes (we don’t need to discuss the criteria) or whether I should keep the multiple rows for the same object when running a correlation analysis or a regression. Thank you in advance. ![enter image description here](http://i.stack.imgur.com/DEOqg.png)


title: Scale parameter MLE scheme known but how to find according distribution PDF?

text: For known location, we can find the scale parameter of a normal distribution by calculating the sum of squared differences to the location, then dividing by n-1 and taking the square root. This is the maximum-likelihood and most- efficient scheme for parameter estimation for the normal distribution. However, for the double-exponential distribution (Laplace) using the sum of absolute values is the MLE scheme. So the power p is not 2 but 1! And for the uniform distribution the maximum difference to the center gives the most- efficient parameter estimation scheme. Also that can be described as sum of power, by let the power p go to infinite (Chebyshev norm). So for these 3 distributions with $p=1,2,\ldots,\infty$ we know both the MLE and the PDF well. Question: If we generalize the estimation scheme by using any positive power (like p=4, 6, or whatever), what is the according distribution PDF(p,x)? This would lead to quite a general PDF and can form a generalized normal distribution (like the Student-t extends the normal for more long-tailed distributions).


title: Testing constraints/restrictions in multiple regression using R

text: I have a Cobb-Douglas equation I am modelling in `R`: $$ \ln(Y) = β + β_a \ln(A) + β_l \ln(L) + β_f \ln(F) + β_o \ln(O) + u $$ Is there an easy way of testing the hypothesis about whether the betas sum to one; less than one; more than one?


title: Randomisation/permutation test for paired vectors in R

text: I'm not an expert, so forgive me if some of the terminology is a little clumsy. Happy to provide more information where required. I have two vectors of 50 paired numeric values in R. I want to perform a two- tailed randomisation or permutation test to determine whether their differences are due to chance or not. > A permutation test (also called a randomization test, re-randomization test, > or an exact test) is a type of statistical significance test in which the > distribution of the test statistic under the null hypothesis is obtained by > calculating all possible values of the test statistic under rearrangements > of the labels on the observed data points. I want to do this type of test because I believe the distributions of the values in the vectors violate the assumptions of other tests such as the t-test (for example, many of the numeric values in the vector are 0). The `permtest` function in the BHH2 library, almost does what I want, but it operates on all $2^{50}$ permutations, which will take too long. Instead, I want to estimate the p-value, by sampling a large number of the possible permutations. I had a look in the coin package, but nothing in there seems to do a permutation test with sampling from paired numeric vectors. Some googling lead me to this email, which suggests that the reason I can't find a package to do it is that it's a one-liner in R. Unfortunately, I'm not experienced enough with R to be able to produce that one-liner. **Is there a package or method that will perform a two-tailed paired permutation test using only a sample of the permutation space?** **If not, would somebody be able to share a short bit of R code to do it?**


title: CART Decision Tree

text: I am running a CART Decision Tree, with 5 independent variables and I wonder how can I interpret the results when the same variable splits the terminal node and also splits the terminal node closest parent. For instance.. The independent variable is Income and one terminal node says x <= 11.000 "Income" and its parent says x > 5.0000, so it seems that they overlap (between 5.000 and 11.000) which can be a way to interpret or on the other hand should I only consider the terminal node figure x <= 11.000? ![enter image description here](http://i.stack.imgur.com/QGO5o.jpg) The question is: Node 18 should be “read” like the customers with income between 10.001and 15.000 OR just the customers with income above $15.000?


title: Ranking different variables by standardization

text: I have seven numeric variables and 12 cases. All of them are measures for how well 5 different casinos are performing. The higher the variable, the better the performance. Four of the variables are: `attendance figures`, `no. of new members`, `visit frequency count`, and `spend per head`. I need to combine these four variables into single performance score which will allow me to successfully rank in order of performance. My idea is very simple: Standardize all four variables, sum them up and whatever you get is the score which should be rank-order. Do you see any problem with applying this approach? Is there any other (better) approach that you would recommend?


title: Wilcoxon-Mann-Whitney as a loss function

text: I'm reading a paper where the authors are using Wilcoxon-Mann-Whitney loss function while minimizing an objective function. As the authors say in the paper, the role of the loss function is to give a positive value (penalty) when the input to the function is positive and it gives a zero otherwise. Here's the equation of the Wilcoxon-Mann-Whitney loss function: $\frac{1}{(1+e^{-x/b})}$ However when I used the function, it gave a positive value whether I gave it a positive value or a negative value! So what's wrong?! Could it be that the authors mis-wrote the equation? I wonder if someone knows this loss function and can clarify. By the way, it's also mentioned in the paper that when someone aims to maximize AUC then Wilcoxon-Mann-Whitney is used.


title: Confidence interval for the difference of two means using boot package in R

text: I have two samples, one of size 52, and one of size 31, that are obtained at different times. I'd like to get a 95% bootstrap confidence interval for the difference between the means of the populations that these samples represent. I've been trying to use the "boot" package in R, and I'm getting an error that I can't figure out. I was hoping someone here could help me out. This is how my data look like (in a dataframe named "totalData"):                     X              samplingTime     1  -0.29              initial     2   0.3               initial                 ....     52  -1.2              initial     53   0.7              final     54  -1.2              final                ....     83   1.52             final       This is what I did to get my bootstrap CI:               meanDiff = function(dataFrame, indexVector) {          m1 = mean(subset(dataFrame[, 1], dataFrame$samplingTime == "initial"))         m2 = mean(subset(dataFrame[, 1], dataFrame$samplingTime == "final"))         m = m1 - m2         return(m)     }          totalBoot = boot(totalData, meanDiff, R = 10000, strata = totalData[,2])     totalBootCI = boot.ci(totalBoot)      and in the last line I get the error:               Error in bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t.o, : estimated      adjustment 'w' is infinite.       I'd very much appreciate any comments. Thanks!


title: Maximum likelihood estimation for a sequence of observations

text: Assuming $A_1, A_2, \ldots, A_n$ are independent exponential random variables (each having the same parameter and, for the sake of simplicity, let's assume the value of the parameter is 1). Define $B_i = A_i + k$ (where $k$ is a constant of unknown value). What would the maximum likelihood estimate (MLE) of $k$ be, if we're also provided with a sequence of observations $(b_1, b_2,\ldots, b_n)$?


title: Logistic regression with time series predictor data

text: I like to know if we can model binary outcome with time series predictors. For example lets say Y is binary. $X_1, X_2, X_3,...,X_n$ is the same predictor variable but is a historical snapshot over time period $1,...,n$. I am interested in predicting $Y$ but include the auto correlation among $X_1,...,X_n$ and seasonality if any among $X_1,...,X_n$.


title: What is an appropriate technique for measuring correlation between classes?

text: I am a school teacher doing a small project that involves collecting data on 5 different teaching methods and the resulting test score of children who are taught was each teaching method. I have 50 students equally split among the 5 teaching methods (10 students/method). The data I have is as follows. Each line below is one student where the first column of data is the teaching method and the second column is the student's score:               Method A,88     Method B,72     Method A,90     Method D,55     Method D,60     Method A,84     ...      I would like to somehow determine if test scores correlate with the different methods. I've been reading about Intraclass Correlation and I want to know if Cohen's kappa or a similar technique would be an appropriate way to look at the correlation of test scores and teaching methods. Any help or direction is much appreciated.


title: How to specify the Wishart distribution scale matrix

text: I'm running the below Bayesian mixing model in R using the rjags package, but I am having difficultly in specifying the scale matrix for the Wishart distribution. Essentially, I want Sigma.inv to be a covariance matrix that accounts for instrumental measurement error as determined from repeat analysis of a control sample. I have chosen the Wishart distribution as the Sigma.inv prior, with the scale matrix taking the form of the covariance matrix of the control samples.However, this results in the model failing to converge, even after 1 million runs with 500,000 burn-in and 125 thinning length, and tiny errors are generated for the estimates of p. So my question is should the scale matrix be in a different form? Or is there perhaps an alternative to the Wishart distribution that would be more appropriate in this case? I should add that the model works quite well if a use the identity matrix to scale the Wishart distribution, but of course this doesn't represent the measurement error that I want to incorporate. Any help will be most appreciated.               model     {      # Likelihood     for(i in 1:N) {     Y[i,1:J] ~ dmnorm(mu[i,],Sigma.inv)     mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]     }          # Priors on s     for(i in 1:N) {     for(k in 1:K) {       s[k,1:J,i] ~ dmnorm(mu.s[,k],Sigma.s.inv[,,k])       }     }          # Prior on p     for(i in 1:N) {     p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])     phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV     for(k in 1:K) {       expphiV[i,k] <- exp(phiV[i,k])     }     }          # Prior on phi     for(i in 1:N) {       for(k in 1:(K-1)) {       phi[i,k] ~ dnorm(0,kappa.inv[k])     }     philast[i] <- -sum(phi[i,1:(K-1)])     }          for(k in 1:(K-1)) {     kappa.inv[k] ~ dgamma(2,1)     kappa[k] <- 1/kappa.inv[k]     }          # Prior on Sigma     Sigma.inv ~ dwish(R.Sigma,k.Sigma)     }      ### # # R Code to run model ### #               # Y values     Y<-matrix(c(5.8565, 18.5979,    0.0045, 5.7579, 1.1401, 0.4945, 0.2122, 0.2613, 0.3984,               14.0696,  10.6925,    0.0084, 8.5126, 2.1779, 1.1406, 0.4051, 0.3472, 0.5971,             13.7417,    10.0652,    0.0084, 8.2977, 2.2172, 1.1204, 0.3929, 0.3583, 0.6143,             11.7867,    15.1036,    0.0072, 8.3222, 1.9592, 0.9697, 0.3454, 0.3067, 0.5695),           byrow=TRUE, ncol=9)     colnames(Y)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")          # mu.s values     mu.s<-matrix(c(6.97463333,  6.657444444,    10.560500000,   13.944108333,               35.47273333,  13.221666667,   6.794500000,    4.149395833,               0.00358000,   0.004533333,    0.008404167,    0.008891667,               5.04523333,   7.128111111,    6.149708333,    6.932004167,               1.19620000,   1.121222222,    2.086083333,    2.446995833,               0.61910000,   0.485333333,    1.010833333,    0.880383333,               0.19653333,   0.199222222,    0.491125000,    0.406904167,               0.07113333,   0.312111111,    0.339916667,    0.298254167,               0.45840000,   0.307222222,    0.616708333,    0.658416667),              byrow=TRUE, ncol=4)     rownames(mu.s)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")     colnames(mu.s)<-c("CB","FD","RV","TS")          # Covariance matrix accounting for measurement error of control samples      Sigma.inv<-structure(c(7.516889e-03,  4.719111e-04, 5.117778e-06,   4.896822e-03,   2.683444e-03,   5.368000e-04,   4.266667e-05,   -8.364444e-05,  4.630667e-04,                       4.719111e-04, 8.097005e-05,   1.132174e-06,   8.006300e-04,   3.065903e-04,   4.421449e-05,   -4.869565e-06,  -5.677295e-06,  7.178551e-05,                       5.117778e-06, 1.132174e-06,   9.927053e-08,   1.090671e-05,   4.145217e-06,   2.788406e-07,   -4.459903e-07,  -2.015942e-07,  9.478261e-07,                       4.896822e-03, 8.006300e-04,   1.090671e-05,   9.350603e-03,   3.260499e-03,   4.947188e-04,   -3.441932e-05,  -4.135604e-05,  7.779478e-04,                       2.683444e-03, 3.065903e-04,   4.145217e-06,   3.260499e-03,   1.390332e-03,   2.255014e-04,   -9.909179e-06,  -3.122995e-05,  2.906319e-04,                       5.368000e-04, 4.421449e-05,   2.788406e-07,   4.947188e-04,   2.255014e-04,   7.714783e-05,   1.326377e-05,   -6.255072e-06,  4.525217e-05,                       4.266667e-05, -4.869565e-06,  -4.459903e-07,  -3.441932e-05,  -9.909179e-06,  1.326377e-05,   1.484058e-05,   2.859903e-06,   -3.663768e-06,                       -8.364444e-05,    -5.677295e-06,  -2.015942e-07,  -4.135604e-05,  -3.122995e-05,  -6.255072e-06,  2.859903e-06,   4.060386e-06,   -5.611594e-06,                       4.630667e-04, 7.178551e-05,   9.478261e-07,   7.779478e-04,   2.906319e-04,   4.525217e-05,   -3.663768e-06,  -5.611594e-06,  7.074783e-05), dim=c(9,9))     Sigma.inv<-as.positive.definite(Sigma.inv)               # Sigma.s.inv values - a precision matrix     Sigma.s.inv<-structure(c(53.2149530065154, 7.78555807206446, -23455.5492262134,                           3.98226035437711, 181.147593576025, -554.144707627592, 810.959913329385,                           -222.049892052889, -676.160229178196, 7.78555807206446, 1.66555095380872,                           -2875.24335449319, 0.59057208979672, 35.2056764704436, -90.5513686454713,                           122.419754836375, -24.1464188840069, -110.985869289511, -23455.5492262134,                           -2875.24335449319, 26623099.2220065, -551.965187580738, -109728.586986558,                           276663.105767497, -386431.378216831, -82170.2904382779, 259038.414073531,                           3.98226035437711, 0.59057208979672, -551.965187580738, 0.918694682125983,                           9.31656679803296, -36.0723479510435, 38.501457515572, -39.2951967395694,                           -43.3516097843844, 181.147593576025, 35.2056764704436, -109728.586986558,                           9.31656679803296, 1046.08776695798, -2268.37491730957, 2936.52705816089,                           60.3905421599093, -2791.52606704426, -554.144707627592, -90.5513686454713,                           276663.105767497, -36.0723479510435, -2268.37491730957, 6336.11679874626,                           -8941.85134505843, 1263.96712724819, 7151.07272932156, 810.959913329385,                           122.419754836375, -386431.378216831, 38.501457515572, 2936.52705816089,                           -8941.85134505843, 14092.5591797243, -2587.63879760624, -10498.261488095,                           -222.049892052889, -24.1464188840069, -82170.2904382779, -39.2951967395694,                           60.3905421599093, 1263.96712724819, -2587.63879760624, 5072.16198995712,                           2677.23974494659, -676.160229178196, -110.985869289511, 259038.414073531,                           -43.3516097843844, -2791.52606704426, 7151.07272932156, -10498.261488095,                           2677.23974494659, 10377.5013426694, 39071331365.1695, -857206689.09299,                           13250725351142.8, 4293968882.52526, -138526471018.157, -346806050176.474,                           -185731096148.676, -40154997287.1756, -13832086879.5978, -857206689.09299,                           45035996.2737049, -291049796239.663, -94316253.2767046, 3042708991.03323,                           7617532441.71197, 4079550081.74548, 881997284.581011, 303819297.524648,                           13250725351142.8, -291049796239.663, 4499056063320832, 1457943322445.14,                           -47034282019304.7, -117752032879430, -63061801053651.4, -13633939080436.7,                           -4696447331898.41, 4293968882.52526, -94316253.2767046, 1457943322445.14,                           513861374.526114, -15241712287.5009, -38158180147.8331, -20435516111.352,                           -4418151354.45332, -1521909040.23616, -138526471018.157, 3042708991.03323,                           -47034282019304.7, -15241712287.5009, 491753442721.598, 1231009860799.41,                           659264193054.762, 142532685331.566, 49097861284.5838, -346806050176.474,                           7617532441.71197, -117752032879430, -38158180147.8331, 1231009860799.41,                           3081922920585.39, 1650491845606.43, 356835753178.233, 122918278499.906,                           -185731096148.676, 4079550081.74548, -63061801053651.4, -20435516111.352,                           659264193054.762, 1650491845606.43, 883961734443.138, 191102478025.129,                           65828570726.718, -40154997287.1756, 881997284.581011, -13633939080436.7,                           -4418151354.45332, 142532685331.566, 356835753178.233, 191102478025.129,                           41361275161.6651, 14232113705.0418, -13832086879.5978, 303819297.524648,                           -4696447331898.41, -1521909040.23616, 49097861284.5838, 122918278499.906,                           65828570726.718, 14232113705.0418, 4947128692.00388, 123.244158717727,                           15.1472733328015, -29469.8213424378, -162.933686694468, -830.16396214664,                           94.7853096969282, -68.8560246643179, 741.612274164014, 2264.45755928194,                           15.1472733328015, 5.35479274130578, -399.822164505876, -26.2908447667995,                           -48.3456870303433, -20.9894871349876, -7.0699744151633, 67.0416328725067,                           116.600590769204, -29469.8213424378, -399.822164505876, 25454078.4009846,                           34253.206150171, 238183.903396863, -33771.2583884048, 38523.5333162132,                           -288045.203177399, -1016390.47602203, -162.933686694468, -26.2908447667995,                           34253.206150171, 237.878595765458, 986.561801969833, -68.5166440497585,                           209.1509212556, -915.217874817308, -2932.38505457965, -830.16396214664,                           -48.3456870303433, 238183.903396863, 986.561801969833, 6669.29110276054,                           -1078.91990285319, 352.917095152379, -5609.73421176863, -17946.4520240348,                           94.7853096969282, -20.9894871349876, -33771.2583884048, -68.5166440497585,                           -1078.91990285319, 1052.57768971324, -913.546632402807, 484.410845136349,                           2431.79560931122, -68.8560246643179, -7.0699744151633, 38523.5333162132,                           209.1509212556, 352.917095152379, -913.546632402807, 2990.61062774375,                           525.038417043516, -3685.88997836577, 741.612274164014, 67.0416328725067,                           -288045.203177399, -915.217874817308, -5609.73421176863, 484.410845136349,                           525.038417043516, 7256.45767695138, 17207.7716622753, 2264.45755928194,                           116.600590769204, -1016390.47602203, -2932.38505457965, -17946.4520240348,                           2431.79560931122, -3685.88997836577, 17207.7716622753, 62661.2158837379,                           7.0546793436162, -0.0754213178285816, -2357.02875362113, 2.50749046482452,                           12.5166502527345, -86.6941792229212, 1.66923366612312, 7.93904545083449,                           -126.082128066627, -0.0754213178285816, 0.618360785191319, 1538.6589391227,                           -0.31799222497181, -1.26341570245756, 2.43556703274056, 18.6459686649981,                           -1.52020047556644, -1.33047266147014, -2357.02875362113, 1538.6589391227,                           7450231.40845542, -1467.47518565875, -8484.65408805661, 38482.6779430954,                           36677.1639056923, 5037.57654299884, -16285.3259202736, 2.50749046482452,                           -0.31799222497181, -1467.47518565875, 1.63035346006115, 1.13929112479727,                           -17.7642468677628, 2.98912100156106, 7.74881374361886, -54.902708360705,                           12.5166502527345, -1.26341570245756, -8484.65408805661, 1.13929112479727,                           106.706180072993, -244.571429668905, 24.7937851688084, -114.000942411739,                           -690.553652091766, -86.6941792229212, 2.43556703274056, 38482.6779430954,                           -17.7642468677628, -244.571429668905, 1534.27886445762, 549.368288309068,                           64.0901199201054, 980.291164452268, 1.66923366612312, 18.6459686649981,                           36677.1639056923, 2.98912100156106, 24.7937851688084, 549.368288309068,                           1781.88219182038, -107.322815612669, -1733.09607995999, 7.93904545083449,                           -1.52020047556644, 5037.57654299884, 7.74881374361886, -114.000942411739,                           64.0901199201054, -107.322815612669, 507.111649502059, 661.404359682356,                           -126.082128066627, -1.33047266147014, -16285.3259202736, -54.902708360705,                           -690.553652091766, 980.291164452268, -1733.09607995999, 661.404359682356,                           9442.32984843676), .Dim = c(9L, 9L, 4L))          #N     N<-nrow(Y)     #K     K<-4     #J     J<-9          myjagsfile<-"JAGSmodel.txt" # the above jags model     tV = t(ilrBase(D=K))     mydata <-list(N=N, J=J, K=K, Y=Y, mu.s=mu.s, Sigma.s.inv=Sigma.s.inv, R.Sigma=Sigma.inv, k.Sigma=J, tV=tV)     mypars <-c("Sigma.inv","p","kappa","phi", "s")     myiterations<-25000     myburnin<-12500     mythin<-175     set.seed(5000)     run.mcmc <- jags(data=mydata, inits=NULL, parameters.to.save=mypars, model.file=myjagsfile, n.chains=2,                   n.iter=myiterations, n.burnin=myburnin, n.thin=mythin, DIC=TRUE)


title: Expected value of q given y is weighted average of mean q and and y

text: It is assumed that: 1) $y=q+u$ Where $q$ is productivity and $y$ a testscore that measures true productivity. $u$ is a normally distributed error term, independent of $q$, with zero mean and constant variance; $q$ is also assumed to be normally distributed with a mean equal to $\alpha$ and with a constant variance. The outcome of this is: 2) $E(q | y) = (1-\gamma)\alpha + \gamma y$   where $\gamma=Var(q)/(Var(q)+Var(u))$ How do you get equation 2? The equation can be expressed as a group effect and a individual effect. [This is a model of statistical discrimination, see: Dennis J. Aigner and Glen G. Cain. Statistical theories of discrimination in labor markets. Industrial and Labor Relations Review, 30(2):175{187, January 1977. URL http:// ideas.repec.org/a/ilr/articl/v30y1977i2p175-187.html.]


title: How to calculate confidence bands for a regression line, given intercept, slope, and standard errors

text: I have the following set of model-averaged fixed effects from a set of binomial GLMMs: ![model parameters image](http://i.stack.imgur.com/yN4wR.png) I would like to plot the predicted effect of "NBT", along with confidence bands, while holding all the other variables at their baseline levels. My attempt to do this in ggplot:               Xvars <- seq(from=0, to=100, by=0.1)  #NBT range is 0-100       binomIntercept <- 1.317       binomSlope <- -0.0076            binomSE <- 0.009         Means <- logistic(binomIntercept + binomSlope*Xvars)                   loCI <- logistic(binomIntercept + (binomSlope - 1.96*binomSE)*Xvars)     upCI <- logistic(binomIntercept + (binomSlope + 1.96*binomSE)*Xvars)     df <- data.frame(Xvars,Means,loCI,upCI)     p <- ggplot(data=df, aes(x = Xvars, y = Means)) +      geom_line() +               geom_line(data=df, aes(x = Xvars, y = upCI),col='grey') +     geom_line(data=df, aes(x = Xvars, y = loCI), col='grey')     p                                                  ![graph image](http://i.imgur.com/eMJBxQQ.png) I'm assuming that the confidence bands are cone shaped because I'm not accounting for uncertainty in the estimate for the intercept. Maybe this is okay (?), but it does look different from every regression line I've ever seen with confidence intervals plotted. Can someone please tell me how I should be writing my equations to get the correct confidence intervals, given the intercept, slope, and standard errors from my model output? (I know I can use the predict function to do this in R, but would like to know how to do it by hand.)


title: k nearest neighbors with repeated points

text: An example illustrates my question. Here I use the `FNN` package in `R` to find the indices of the 400 nearest points for each point in the set 1 to 10000:               require(FNN)     data = 1:10000     x = knn.index(data, k = 400, algo = "kd_tree")      Now suppose that I use a different dataset that still has 10000 points, but only 100 are unique:               data = sample(1:100, 10000, replace = TRUE)     x = knn.index(data, k = 400, algo = "kd_tree")      Both of these examples take about the same amount of time to run. But isn't it be possible to perform the second task much more quickly than the first, since there are only 100 unique points in the second task? If so, where can I look to find such an algorithm? Table 1 in this paper is potentially relevant, although I'm not sufficiently familiar with the various methods to be able to recognize any of these as the solution.


title: What is $\nu$ in $\nu$-SVM controlling?

text: I have read the $\nu$-SVM formulation, but still cannot understand the intuition behind introducing $\nu$. Can anyone give a interpretation of the parameter? I read about it from this paper. This in turn was referred in the book Learning with Kernels.


title: Regression analysis (intercept help needed)

text: I am having a situation out here. I got data from system of ODE's and want to estimate a quantity. The problem is this For linear estimation, I used this relationship               Y = C + C_a A + C_b B + C_d D      where `C` is the intercept and the rest of the `C's` are coefficients. I want this estimation to work for a void range but the problem is with the intercept. I can find intercept before hand. When, I approximate A, B, D in a small range say [0.01, 0.05], the intercept matches perfectly but when I increase the range say [0, 0.5] there is a big different between the intercept I get through the estimation and the one I can find before hand. I believe that there should not be that much difference between the estimated intercept and the one we can get before estimation. What should I do about this? ![enter image description here](http://i.stack.imgur.com/BHf3Q.png) For 625 sample in the range [0, 0.5], I get this estimation `YY = 1.50651907424 - 0.295152822848*A - 0.931922695296*B - 7.59056025338e-09*DD - 0.656723555386*EE` For 108 sample in the range [0, 0.5], I get this estimation `YY = 1.62459346269034 - 0.397600596769767*A - 1.302971457454*B - 7.39245088646173e-09*DD - 0.827832435224737*EE` For 256 sample in the range [0.01, 0.05], I get this estimation `YY = 1.75844218786707 - 0.729590804005662*A - 2.56393407417859*B - 7.39245088646173e-09*DD - 1.29876329140818*EE` The exact intercept before estimation is `1.77094768051425`


title: How can I show $(X_{1} / r_{1}) / (X_{2} / r_{2})$ and $(X_{3} / r_{3}) / ( (X_{1} + X_{2}) / (r_{1} + r_{2}) ) )$ are independent?

text: In the book by Hogg and Craig 6th edition, there is a problem numbered 3.6.14. saying that If $X_{1}$, $X_{2}$, and $X_{3}$ are three independent chi-squre variables with degrees of freedom $r_{1}$, $r_{2}$, $r_{3}$ respectively, show that $X_{1} / r_{1}) / (X_{2} / r_{2})$ and $(X_{3} / r_{3}) / ( (X_{1} + X_{2}) / (r_{1} + r_{2}) ) )$ are independent. Actually, it is a part (b) of 3.6.14., but from (a), we know that $X_{1}/X_{2}$ and $X_{1}+X_{2}$ are independent. Of course, we can solve (b) by tedious calculations, but I want to know the neat solution. But even if I tried, I cannot know the solution. I guess that if there are three random variables $X$, $Y$, and $Z$ s.t. $X$ and $Y$ are independent, and $X$ and $Z$ are independent, then for any functions $f$ and $g$, $f(X)$ and $g(Y,Z)$ are independent. If this is true, then the question above would be obvious fact. But I don't know whether the statement is true or not. Can anybody help me?


title: I want to show $Y_{n} \xrightarrow{ p} c$ ,then show that $f(Y_{n}) \xrightarrow{ p} f(c) $

text: Let $c \in \mathbb R$ . Let be a sequence of random variables . let $f:\mathbb R \to\mathbb R$ be continous function . if $Y_{n} \xrightarrow{ p} c$ ,then show that $f(Y_{n}) \xrightarrow{ p} f(c) $


title: Why such a poor result from sparse PCA R package?

text: I'm preparing to use R to perform sparse analysis on my data. I tried to get started with an ad hoc example, but the reconstruction result turned out really poor. I'm wondering if I was making any mistake in the code, or was I encountering any intrinsic drawbacks of the method? The data is ad hoc:               > prod          [,1] [,2] [,3] [,4] [,5]     [1,]    1    1    1    1    2     [2,]    1    2    1    2    3     [3,]    0    1    0    1    1     [4,]    0    0    1    1    1     [5,]    0    1    0    1    1      and I use this               rst = nsprcomp(prod, retx=T, nneg=T)      to perform non-negtive sparse PCA. I reconstruct the data like this:               recon = rst$x %*% t(rst$rotation) + matrix(1,5,1) %*% rst$center      The final result turned out to be really poor:               > abs(prod - recon)                 [,1]      [,2]      [,3]      [,4]      [,5]     [1,] 0.386523767 0.2432819 0.1186152 0.3535836 0.3770144     [2,] 0.184944985 1.1214659 0.6230812 1.2173770 1.8819417     [3,] 0.096185913 0.4301642 0.2249522 0.4825192 0.6942987     [4,] 0.009206957 0.5044193 0.2917920 0.6059221 0.8703587     [5,] 0.096185913 0.4301642 0.2249522 0.4825192 0.6942987      Any idea or comment is appreciated. ## Edit The package is called nsprcomp


title: Which statistical test to use? Paired-sample or independent-sample?

text: I have 1 subject rating 40 participants on the big five personality traits and other variables. these 40 participants are in 2 conditions (20 each). I want to compare whether there is a difference between how the subject rates the two conditions. I know it sounds like paired samples, but doesn't paired samples compare the difference of 2 answers for the same participant? I have one subject giving 40 answers. Also practically, how should I construct my data set? Paired samples cannot test the difference between conditions. I would have to make new variables for condition 2 and then compare these. In SPSS, however, won't it look as if participant 1 and participant 20 are the same person, since their columns would be next to each other on the data set.   Any tips on what analysis to use? I have to hand in the study in a few days, so fast advice would be appreciated!


title: Contradictive results of correlation and intergroup comparisons

text: I got strange results in spss and I don't know how to explain them statistically. I made correlations between: 1) personality and emotional inteligence (EI) 2) personality and emotional competence (EC) Correlations (r-Pearson) showed that personality has many connections with EC, but almost none with EI. Everything was OK until I made intergroup comparisons (ANOVA). I made 4 groups of people:   * with high EC and high EI   * with high EC and low EI   * with low EC and high EI   * with low EC and low EI I checked if there are significant differences between these gropus, in personality profiles. And results showed that:   * If two groups have the same level of EI they don't differ with personality (even if they have different level of EC).   * And when groups differ with EI, they also differ with personality (even if EC are the same). These results are unexplainable for me because correlation showed that personality is connected with EC , NOT with EI! What do you think... What does it mean? I've checked if I did all analysis in proper way.


title: How to calculate prevalence ratios?

text: In a reviewer request from a submitted manuscript, they requested that I use **prevalence ratios** instead of odds ratios. The reviewer says: > Using logistic regression when prevalence of outcome is greater than 5-8% > produces OR that are an exaggeration of the true relative risk. In this > study, the prevalence of the entity was 50%. It may be best to consider > prevalence ratios. This clinical study was of measures of a specific area of the body (K, continuous variable) and its association with disease X (present/not-present). The observation of the study is that the larger K is, the higher the OR of having disease X. The number of subjects is 200. I'm stumped in how to calculate prevalence ratios as requested. I used the logistic fit model in JMP to calculate OR however am yet to find instructions on how to obtain prevalence ratios for this particular sotware. I'm aware that SAS has a routine for this but wonder if anyone can give me a direction on how to calculate prevalence ratios by hand or other method. I'll be glad to provide additional info/data.


title: Can we reconstruct the hidden (latent) variables after executing EM?

text: The question is in the title. I know that EM algorithm could do maximum likelihood estimation for models that have latent variables. I would like to know can we get the (estimated) value of these latent variables once the algorithm is finished ? I know the question is vague and it depends on problems so I try to describe mine as follow: I have a random point process (failure process) with ${T_i}_{ \lbrace i \geq 0 \rbrace}$ denote the failure times. At each failure we have a flag $Z_i \sim \mathcal{B}(p)$, where $\mathcal{B}(p)$ is a Bernoulli distribution, to indicate the realized action, though in the end we can not observe $Z_i$ but only $T_i$. I might note that $T_i$ is not iid because it depends on the realized action at $T_i$ to jump from $T_i$ to $T_{i+1}$. With the help of EM algorithm, I can estimate $p$ and others parameters BUT I couldn't know the value (or estimated) of $Z_i$ since this algorithm don't mention. So I would like to know if there are any articles discuss this one? Any idea is welcome and thanks in advance !


title: Gaussian mixture regression in higher dimensions

text: **Problem:** I have a discrete representation of a surface/height-map $z = f(x,y)$ that i want to model as a mixture of gaussians (please take probability distributions out of your mind for a moment). More precisely, $ z = f(x) = \sum\limits_{i=1}^{n_g} a_i \; g(x; c_i, \Sigma_i)$ where   * $x \in R^d$   * $g(x; c_i, \Sigma_i) = exp\\{-\frac{(x-c_i)^T\Sigma_i^{-1}(x-c_i)}{2} \\}$,    * $a_i \geq 0$ is the max-height of the ith gaussian component,    * $c_i \in R^d$ is the center of the ith gaussian component, and    * $\Sigma_i$ is the covariance of the i'th gaussian component that determines its support in space and its orientation. Given the data of the surface, i want to fit this model i.e estimate the parameters $[...,a_i, c_i, \Sigma_i, ...]$ given some reasonable initial guesses by minimizing the least squares objective function given below: $E(a,c,\Sigma) = \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right)^2$ _Constraints:_   * $a_i \geq 0$,    * $\Sigma_i$ should be symmetric ($\Sigma_i = \Sigma_i^T$), and positive-definite ($u^T\Sigma_i u > 0 \; \forall{u} \in R^d$). I would like to emphasize that i want to optimize for the full covariance matrices $\Sigma_i$ (even diagonal wont suffice). **Approach so far** I will describe below the approach i have been venturing on so far but i _would really appreciate any other ideas in fitting this model_ to data or _atleast point me to the right literature_ to read (googling takes me to so much literature on fitting a mixture model to a probability distribution using EM). I might have a huge number of components in the mixture. So it will be nice if you could _supplement your answers with comments on computational complexity, the quality of the solution (how far from global minimum) and sensitivity to initialization (how close should the initial guesses be)_. My initial idea is to use iterative minimization algorithms such as the Gauss- Newton algorithm or the Levenberg–Marquardt algorithm. So i started off by deriving the derivatives of the error/objective function $E(a,c,\Sigma)$ with respect to each of its parameters: $\begin{eqnarray} \frac{\partial E}{\partial a_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial a_i} \\\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) g(x_j; c_i, \Sigma_i) \end{eqnarray}$ $\begin{eqnarray} \frac{\partial E}{\partial c_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial c_i} \\\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \frac{\partial g(x_j; c_i, \Sigma_i)}{\partial c_i} \\\ & = & -2 \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \; g(x_j; c_i, \Sigma_i) \; \Sigma_i^{-1} (x_j - c_i) \end{eqnarray}$ $\begin{eqnarray} \frac{\partial E}{\partial \Sigma_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial \Sigma_i} \\\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \frac{\partial g(x_j; c_i, \Sigma_i)}{\partial \Sigma_i} \\\ & = & -2 \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \; g(x_j; c_i, \Sigma_i) \; \frac{\Sigma_i^{-T} (x_j - c_i)(x_j - c_i)^T\Sigma_i^{-T}}{2} \end{eqnarray}$ Now while i minimize the error function $E(a,c,\Sigma)$ using one of the iterative minimization algorithms, _i want to enforce that the covariance matrices $\Sigma_i$ are symmetric and positive definite. And i dont know how to enforce that constraint?_ (should i try to optimize for elements on one side of the diagonal which is sort of projecting it to the domain of symmetric matrices but not necessarily positive definite?) I am guessing this is an already solved problem atleast reasonably. It would be great if you could point me to the literature for this problem so i can move forward. Alternatively, coming from a computer-vision + machine-learning background, i can't resist but to think about applying the Expectation-Maximization algorithm. Intuitively, if i knew which guassian component each (x,y,z) location belongs to (in a soft sense) then i can compute the parameters of the gaussian components. And, if i know the parameters of the guassian mixture, I can compute the degree to which each (x,y,z) location belongs each component in the mixture. And, Expectation-Maximization algorithms in the end boils down to minimization of an objective function using coordinate descent. So it would be nice if you could _point me to any literature you can point me to that takes the EM approach to this problem_.


title: Comparison of t-value and F-value

text: I have conducted a t-test via SPSS, followed by an ANCOVA (Analysis of Covariance). Now I am wondering if it is correct to report the t-value and then the F-value of the ANCOVA. But are those two values really comparable? Or would I have to run an ANOVA first, instead of a t-test? I know that technically there shouldn't be a difference between the results of a t-test and an one-way ANOVA if I only have two groups, but this really confuses me. Additionally the criteria of homogeneity of variances isn't met - is there a way to correct that when conducting an ANCOVA?


title: Exploratory data analysis

text: I have two data sets: one has a continuous response variable whereas the other discrete. I will use a linear model for the continuous data and a logit and probit model for the discrete data. What I want to do is to conduct some exploratory data analysis for both data sets. I am wondering what plots can be used to assess the model assumptions. I am also open to new thoughts. Any help would be greatly appreciated.


title: Analyzing split-plot design in lme4 in R

text: I have the data from a split-plot design where A is my whole plot fixed factor with two levels and B is my subplot fixed factor with 2 levels and C is my random block factor. How do I analyse these data to get appropriate the DF (=1 for each fixed factor and interaction) using GLMM in the lme4 package in R? I tried `response~A*B+(A/B)+(1/C)` which gives the same result as when leaving the (A/B) out of equation completely; thus I am not sure how to incorporate into the equation that the factor B is nested within A. I also tried `response~A*B+(1|C/A/B)` although this is probably wrong, but all the codes are giving me different results and I cannot tell which one seems correct as there are no DF in the summary and the visual output of my data is not so straight forward. Can anyone help with nesting? Also if I should decide my fixed factors are in fact crossed rather than nested (still have to re-think that), would the code look different?


title: Indicator variables and backward elimination with GLM

text: I am running logistic regression using `glm` in R on data, that has some indicator variables to it. Two of those have multiple levels and have been rewritten as (#levels-1) new dummy predictors respectively. As for variable selection, I use step (backward variable selection based on AIC), which is where I struggle: Even though all dummies of the initial indicator were kept or removed, I wonder what one would do, if a single dummy (i.e. one of a set of associated dummies) was eliminated during the variable selection? After all, this would mean that I lose the information encoded for the base-level (since it is "indicated" by 0s in all the other dummies). Could someone tell me what to do if such a case occurred? Thanks & Regards


title: Hotelling's T (or something similar) on results of relative importance analysis

text: I am using the R package relaimpo (Grömping, 2006) to analyze the relative importance of 16 regressors on one dependent variable (N=7000). The results show that regressors vary drastically from each other in explained variance, varying between 18.2% and 0.8%. That is, after controlling for each regressor some regressors explain a lot of variance, others only very little. For the paper, I need to report one overall p value, e.g.: "Overall, regressors varied significantly from each other (p<.001)." The package only provides me with information about significant differences between regressors A and B, regressors A and C, et cetera, not an overall test. Since I only have aggregate data here for 16 regressors (18.2%, 16.5%, ..., 0.8%), I don't know what test to use.


title: Binary classification of DNA motif sequences (bioinformatics)

text: I've been working on on a method for binary classification of DNA sequences. In more detail, here is what the method does. Given a family of DNA sequences, for example DNA sequence motifs, I try to predict whether other sequences belong to this family, by measuring their similarity to the sequences in the family. My method fits a distribution to the sequences, and then assigns a pvalue to sequences not in the family. In the data I'm using, the sequences are all of length 28. The sequences I am analyzing are here, specifically human 12 RSS and mouse 12 RSS. Since a comparison to existing methods is always a good thing, I am wondering what are the standard methods to beat, if any? I'm not very familiar with the available methods/algorithms. I am in the process of trying MEME. This does not seem to do exactly what I want, and I don't know if I will be able to persuade it to. Specifically, I'm not sure if I can tell it that the sequences of 28 length _are_ the motif. I got the impression from the documentation that it decides what the motifs are by itself, or something. This would probably be more on-topic on a bioinformatics site, but SE does not currently have one, so I'm trying here. I can give further details if necessary. NOTE: I think it would be nice if a bioinformatics tag could be created and used here. This does not seem to exist currently. UPDATE: I think my original question was poorly written and lacked sufficient detail, so I'm adding some more (hopefully useful) details. I'm analyzed two RSS data sets, each of which is a collection of sequences, as I state above. The main purpose of the analysis is to do prediction. So, I used a cross-validation method. I divided each data set into 5 parts, and use 4 of the five parts as a training set in turn. (The number 5 here is a bit arbitrary, but since I want to include the results per training set, I don't want the number to be too large.) After fitting a model to the training set, I then used this model for prediction as follows. The RSS data set is contained in gene segments, typically one or two RSS per gene segment. The gene segments are often much larger than the RSS. These are 12 RSS, so each RSS is of length 28. I took all the gene segments I could find that contained an RSS, and selected from them all contigous sequences of length 28. The current total number of these sequences is 449905 for one, and 624400 for the other. The corresponding number of RSS was 118 and 201. Note that these sets did not necessarily contain all distinct values. So, I used the model derived from the training set to calculate pvalues for all these approx 500,000 sequences. (I'm leaving out some details here, but I don't think it is important how exactly I calculated the values.) Then I ranked the sequences by order of decreasing pvalues. The idea was that the RSS sequences would rank highly in this ranking, and in the event they did. Now, I'd like to find an algorithm which can perform a similar procedure. So far I haven't had much luck - I've been wading through a morass of confusing papers and results, mostly written by biologists, and therefore hard for me to understand. Much of the related material seems to be about de novo motif discovery, which is about magically finding motifs apparently without a training set. Also, much of the software does not seem designed to deal with sequences as long or as numerous as in my example. Someone suggested tomtom from MEME, but so far I've not got it to work.


title: Confidence intervals for log normal responses

text: I got this assignment from Generalized Linear Model class. At first glace it looked like it is an easy task, but there are a lot of subtle (at least in my opinion) things, which I would like to clarify with somebody. So here is the problem: Assume that we are interested in random variable $Y=e^Z$. Variable $Z$ is obtained from simple linear regression i.e. $Z=a+bX+\varepsilon$, where $\varepsilon \sim N(0,4)$. Given 100 observations of $(Z,X)$ we have following empirical characteristics $\bar x=2$, $\sigma^2_x=1$, $\bar z = -1$ and $\operatorname {cov} (z,y)=-1 $.   1. How do we construct _approximate_ Wald's 95% confidence interval for parameter $\theta = E(Y|X=1/2)$, using the fact $E e^{\lambda U} = e^{a \lambda+\sigma^2\lambda^2/2}$, when $U \sim N(a, \sigma^2)$?   2. How do we construct exact Wald's 95% confidence interval for individual prediction of $Y$ when $X=1/2$?   3. Now assume that we do not know variance of $\varepsilon$ and from those 100 observations we know that the empirical variance of $Z$, $\sigma^2_z=4$. I started by plugging in the empirical statistics to specify simple regression part and I got $Z=1-X+\varepsilon$, so the $$\operatorname{\log} Y =1-X+\varepsilon. $$ As the $\operatorname {log} Y$ is normaly distributed then the distribution of the mean of $Y$ should be lognormal? If $Y$ is log normal then I can't write write it in canonical exponential family form $$\operatorname{exp}\left\\{ \frac{y \theta-b(\theta)}{a(\phi)}-c(y,\phi) \right\\},$$ because I have $\operatorname{log}y$ instead of $y$ multiplied by $\theta$. That means that I do not know what $b(\theta)$ and other parameters are, which I need for standard error estimation. So any tips on getting at least first question answered would be great!


title: How to put variables related by an IF condition in a model

text: The first variable is 'SMOKE' that indicate if the person is a smoker or not. The second variable 'Daily' is related to the first by an If condition: If the person is a smoker then does he smoke daily or not. How can I put these two variables in my model? Thanks in advance, EDIT: (1)I prefer to exploit the evident ordering from non-smoking to smoker to daily smoker. (If this will be significant in my model) (2) Smoking is an independant variable in my model (In fact, It's a risk factor for a disease)


title: SVM: non-linear versus linear models

text: In the context of classification on somewhat large datasets (say at least 50Kx50K), I am wondering in which cases non-linear models are superior to linear ones to warrant the added complexity. I often see in my own research that for these larger datasets, non-linear datasets cannot outperform linear ones (say for a linear kernel SVM and an RBF kernel SVM). But this might be biased due to my 'repository selection' of datasets which are all sparse and drawn from transactional data. My intuition says that specifically for an RBF kernel, a linear kernel should be a lower-bound for the performance that you can achieve with an RBF kernel, but my hope of attaining higher performance than this lower bound is not fullfilled because in the end they all achieve more or less the same performance. Specifically my question is this: have you encountered situations in which non-linear models were worth the effort? Or do you perhaps know about research confirming/rejecting my own observations?


title: Set hypothesis instead of point hypothesis

text: Statisticians are often interested on testing point null hypotesis, such as: $$H_0: \mu =0 \,\,\, vs. \,\,\, H_1:\mu\neq0.$$ As Jeffreys himself said, this actually corresponds to some presumption that $\mu$ is fairly small. My question is: has the alternative formulation: $$H_0: \vert\mu\vert \leq \epsilon \,\,\, vs. \,\,\, H_1:\vert\mu\vert>\epsilon,$$ for some $\epsilon>0$ been considered? If yes, could you please provide some references?


title: An up to date keyword set on global news

text: My question is not strictly binds to the topic of text mining, but maybe you can help. I am hunting for a keyword set, which has the following criterions: \- contains only english words/n-gramms or named entities \- manual (tagged by human) tags of global news \- has some main topics (5-10), e.g. tech, business, sport, ... \- has to be relatively big (10000+ tags) \- has be up to date \- it would be nice that the keywords have frequency weights. I thought that couple of the biggest news portals has tag cloud which (partly) fits on the criterions. But I didn't find anything on BBC News, CNN News, Reuters, ... Interestingly I found some portals in my mother language (hungarian), I cant believe that, there isnt anything on global level. I dont need API, I can parse the HTML if necessary. Maybe corpora can be useful. Thanks.


title: Model fitting to data by using machine learning algorithms?

text: I am trying to fit an equation to data. I know the form of the equation but I need to know constants parameters in the equation. I used non-linear fitting and optimization techniques but I could not get promising results. This can be a dump question but I wonder if it is possible to fit an equation to data by using machine learning algorithms such as guassion process, neural network and so. For example, I have data and I wanna fit this equation: a*sin(x) +b*cos(x). Is it possible to find the constants of "a" and "b"?


title: Hidden Markov model & statistical significance

text: I am using HMM to explore the language development of one individual with six variables as input, which are trained into three sequences. In the first state, the covariance of variable A and B is -0.25. In the second state, the covariance of variable A and B is -0.01. In the third state, the covariance of variable A and B is 0.10. My question is: Can I make a claim that these two variables (A and B) show a competing relation to a supporting relation without testing the "significance" as we do for the correlation test?


title: EM algorithms - confidence interval estimation

text: Does anybody know how to find the confidence intervals for estimated parameters of a mixture of Gaussians by using EM algorithm?


title: Is a simple time series one- or two-dimensional?

text: Say we have a simple monthly time series:               > ldeaths           Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec     1974 3035 2552 2704 2554 2014 1655 1721 1524 1596 2074 2199 2512     1975 2933 2889 2938 2497 1870 1726 1607 1545 1396 1787 2076 2837     1976 2787 3891 3179 2011 1636 1580 1489 1300 1356 1653 2013 2823     1977 3102 2294 2385 2444 1748 1554 1498 1361 1346 1564 1640 2293     1978 2815 3137 2679 1969 1870 1633 1529 1366 1357 1570 1535 2491     1979 3084 2605 2573 2143 1693 1504 1461 1354 1333 1492 1781 1915      Is it one- or two-dimensional? The values are of course one-dimensional, but since we have information about the time, do we consider a simple time-series as this one as **two-dimensional data**?


title: How to choose an appropriate covariance structure for a simple one-way within-subjects RM ANOVA design?

text: It is postulated that one of the main issue is to find an appropriate covariance structure for repeated measures designs [Ref1]. SAS' PROC MIXED contains a number of covariance structures. Despite _"many choices among models to fit to a given data set in the mixed model setting... [and] we must always remember that all models are wrong (because they are idealized simplifications of Nature), but some are useful_ [citation]." there are different recommendations for choosing among the covariance models are known [Ref2], [Ref3]. I do experiments with simple one-way within-subjects RM ANOVA (balanced and unbalanced) design without a between-subject factor described here. My question is   *how is it important if I have a " **condition** " ( _not time_ ) as a within-subject factor and how to choose an appropriate covariance structure in this case*? Perhaps this is a strange case as I have a native substance and its chemical derivatives with _similar_ molecular structures. So probably I have correlation caused by the substance itself?


title: Undefined real result in a zero-inflated negative binomial

text: I have run a zero-inflated Poisson model in WinBUGS without problems, and now I am trying to run its equivalent negative binomial. However, I get an "undefined real result" trap message over and over. Here is the data set http://www.megafileupload.com/en/file/474809/data1.html               #Data handling     library(R2WinBUGS)     setwd('') # working directory     dget('data1')     data1 <- data1[order(data1$C), ]         n = nrow(data1)         n0 = length(which(data1$C == 0))          sink("ZINB.txt")     cat("     model {          # Pr(Y=0|x,z) # probability of zeros     K <- 10000     for(i in 1:n0){     zeros0[i] <- 0     zeros0[i] ~ dpois(mu0[i])     mu0[i] <- -log(p00[i] + (1-p00[i])*pow(overd[i], r)) + K          overd[i] <- r/(r + lambda0[i])     p00[i] <- min( 0.999999, max(.000001,p0[i]) ) # avoid stackoverflow     logit(p0[i]) <- a0 + b0.A*A[i] + b0.M*M[i] + b0.cat*cat[i] + b0.d*d[i]     log(lambda0[i])<- a + logh[i] + b.A*A[i] + b.M*M[i] + b.cebo*cebo[i] + b.descart*descart[i] + b.d*d[i]     }          # Pr(Y>0|x,z) # probability of positive counts     for(j in (n0+1):n){     zeros[j] <- 0     zeros[j] ~ dpois(mu[j])     mu[j] <- -log((1-p1[j])*fnb[j]) + K          lfnb[j] <-  r*log(overd[j]) + C[j]*log(1 - overd[j]) + loggam(C[j] + r) - loggam(r) - loggam(C[j] + 1) # - log(1 - pow(overd[j], r))     fnb[j] <- exp( lfnb[j] )     overd[j] <- r/(r + exp(lambda[j]))     p1[j] <- min( 0.999999, max(.000001,p[j]) )     logit(p[j]) <- a0 + b0.A*A[j] + b0.M*M[j] + b0.cat*cat[j] + b0.d*d[j] + b0.descart*descart[j]     log(lambda[j])<- a + logh[j] + b.A*A[j] + b.M*M[j] + b.cebo*cebo[j] + b.descart*descart[j] + b.d*d[j] + b.dcos*dcos[j]     pred[j] <- (1 - p1[j])*lambda[j] # predicted values (see Zuur et al 2009)     var[j] <- (1 - p1[j])*(lambda[j] + lambda[j]*lambda[j]/r) + lambda[j]*lambda[j]*(p1[j]*p1[j] + p1[j]) # variance (see Zuur et al 2009)     res[j] <- (C[j] - lambda[j]*(1 - p[j]))/sqrt( var[j] ) # Pearson residuals     disp[j] <- pred[j]/var[j] # dispersion     }          # Vague priors for model coefficients     r ~ dunif(0,10)     a0 ~ dnorm(0,0.001)     a ~ dnorm(0,0.001)     b.A ~ dnorm(0,0.001)     b0.A ~ dnorm(0,0.001)     b.M ~ dnorm(0,0.001)     b0.M ~ dnorm(0,0.001)     b.d ~ dnorm(0,0.001)     b0.d ~ dnorm(0,0.001)     b.descart ~ dnorm(.0,0.001)     b0.descart ~ dnorm(0,0.001)     b0.cat ~ dnorm(0,0.001)     b.dcos ~ dnorm(0,0.001)        b.cebo ~ dnorm(0,0.001)     }     ",fill=TRUE)     sink()          # Data passed to WinBUGS     win.data <- list( C = data1$C, A = as.numeric(data1$A), M = as.numeric(data1$M), d = as.numeric(scale(data1$dcol)), dcos = as.numeric(scale(data1$dcos)), cebo = as.numeric(data1$cebo), descart = as.numeric(data1$descart), cat = as.numeric(data1$cat), logh = log(data1$h), n = n, n0 = n0 )          # Initial values     inits <- function(){ list( r = runif(1, 1, 5), a0 = rnorm(1), a = rnorm(1), b.A = rnorm(1), b0.A = rnorm(1), b.M = rnorm(1), b0.M = rnorm(1), b.d = rnorm(1), b0.d = rnorm(1), b.descart = rnorm(1), b0.descart = rnorm(1), b0.cat = rnorm(1), b.dcos = rnorm(1), b.cebo = rnorm(1) ) }          # Parameters     params <- c( "r", "mu", "lambda", "a", "a0", "b.A", "b0.A", "b.M", "b0.M", "b.d", "b0.d", "b.descart", "b0.descart", "b0.cat", "b.dcos", "b.cebo", "res", "pred", "var" )          # MCMC settings     nc <- 3     ni <- 30000     nb <- 3000     nt <- 3          # WinBUGS execution     out <- bugs( data = win.data, inits = inits, parameters.to.save = params, model.file = "ZINB.txt", n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = TRUE )          print(out, dig = 2)      The trap arises from the Pr(Y > 0) part of the model since the Pr(Y = 0) runs without problems. I've tried different inits and priors to no avail.


title: Fitting two sets of values by changing x-axis

text: I have two datasets. The first one consist of measurements with a known date. The second dataset consists of measurements along a distance axis. The second dataset is sequentially deposited (for example an ice-core) with an unknown growth rate. Values in the two datasets are expected to match within a certain random error. An example:               curve <- data.frame(date = as.Date(c('01/09/11', '01/10/11', '01/11/11', '01/12/11', '01/01/12', '01/02/12', '01/03/12', '01/04/12', '01/05/12', '01/06/12', '01/07/12', '01/08/12', '01/09/12'), format = "%d/%m/%y"), value = c(3.33, 3.80, 4.19, 4.43, 4.56, 4.76, 4.71, 4.96, 4.72, 4.66, 4.35, 4.12, 3.88))          points <- data.frame(distance = 1:10, value = c(3.5, 4.2, 4.5, 4.4, 4.7, 4.8, 5.1, 4.9, 4.1, 3.7))          library(ggplot2)     library(scales)          a <- ggplot(curve, aes(x = date, y = value)) + geom_line() + scale_x_date(name = "Month", breaks = date_breaks("months"), labels = date_format("%b"))     b <- ggplot(points, aes(x = distance, y = value)) + geom_point()          grid.arrange(a,b,ncol=1)      ![enter image description here](http://i.stack.imgur.com/fQNxg.png) I call the dataset with known date "curve" because the number of measurements is much larger than for the dataset with a distance axis ("point" dataset). In reality both of them consist of measurement points, but I suppose I could approximate the curve dataset to an actual curve by using loess or splines. I would like to fit the "point" dataset to the "curve" dataset by changing the x-axis (or growth rate) for the "point" dataset with a constraint that the growth rate > 0 (distance between "point"s). I assume that I would be making a growth model for the distance axis. What would be the best method to do this kind of fitting? Is there an R package developed for this kind of problem?


title: Contrast matrix definition

text: What _exactly_ is contrast matrix and how _exactly_ is contrast matrix specified? I.e. what are columns, what are rows, what are the constraints on that matrix and what does number in column `j` and row `i` mean? I tried to look into the docs and web but it seems that everyone uses it yet there's no defition anywhere. I could backward-engineer the available pre-defined contrasts, but I think the definition should be available without that.               > contr.treatment(4)       2 3 4     1 0 0 0     2 1 0 0     3 0 1 0     4 0 0 1     > contr.sum(4)       [,1] [,2] [,3]     1    1    0    0     2    0    1    0     3    0    0    1     4   -1   -1   -1     > contr.helmert(4)       [,1] [,2] [,3]     1   -1   -1   -1     2    1   -1   -1     3    0    2   -1     4    0    0    3     > contr.SAS(4)       1 2 3     1 1 0 0     2 0 1 0     3 0 0 1     4 0 0 0


title: Estimating the distribution from partially aggregated data

text: I have a dataset that consists of a list of tuples:               (1, .5), (3, .7), (1,.1), (23,.8),...      In each tuple, the second element is the average price, and the first element is the number of items. So for example, (3, .7) means that three items were bought for an average price of .7. From this list of data, I'd like to fit a power law tail (although really, I'm interested in the general version of fitting a distribution to partially aggregated data, in a nonparametric way). What are the common ways to split up, say (3, .7) into three separate observations, or more generally, how should I treat this partially aggregated data?


title: Missing input value during prediction of a generalized linear model

text: If I'm doing prediction with a generalized linear model and a new batch of inputs comes with some missing values, what strategies can I use to minimize the loss of information from the missing inputs? Let's assume that some of the predictors may be correlated and that the missing inputs are random. My first guess would be to use mean imputation from mean of the training data. I read some bad things about mean imputation, but they were all about when doing estimation, not prediction. Is it a better choice in this case?


title: Generating null distributions by a residual permutation procedure

text: I am trying to understand the method described in this paper which describes an hypothesis-testing framework for stable isotope ratios. The data are in a bivariate isotopic space and the metrics that are of use to me are the group centroid locations, mean distance of observations to the centroid and group eccentricity. These values are fairly straightforward to determine (I'm using Mathematica for data analysis); the problem arises when the authors compare the test statistics to _null distributions generated by a residual permutation procedure (RPP)_. The authors describe this procedure in what is probably a satisfactory manner to those who are comfortable with linear algebra and R; however I cannot include myself in either of those categories. I believe the general procedure proposed is to compare the vectors from each group centroid to that of an overall centroid using a multivariate test statistic, in this case, Hotelling's $T^2$. (Again, I'm delving into foreign territory here; however I'm going to assume for the moment that I can figure out the multivariate analogue of the t-test.) What confuses me is the use of the term _residual_ : > All test statistics were compared to null distributions generated by a > residual permutation procedure that works by shuffling residual vectors, > where each observation is described as a residual vector from the overall > centroid and also as a residual vector from each group centroid. Can someone assist me in understanding how to generate the null distributions? Also, can someone confirm that this type of analysis is, in essence a null hypothesis test analogous to the univariate t-test? Below is a subset of the data, if it is of help in answering the question.               sample1 = {{6.9, 11.3}, {10.5, 8.6}, {9.2, 12.2}, {10.5, 0.4}, {13.9, -1.7}}     sample2 = {{8.4, 13.7}, {10.2, 14.8}, {14.4, 14.3}, {11.6, 13.1}, {9.5, 7.9}}


title: Random parameters with endogeneous regressors by Stata SEM or GLLAMM?

text: I'm interested in estimating a linear model with random parameters (RP) and some endogenous regressors in stata. The main issue is to obtain estimates of **parameters variances** and not the parameters by themself (a two stage method provides unbiased parameter estimates). Is there any specific setup within SEM or GLLAMM modules to do an straight estimation of this model ?


title: An elementary question on binomial test: why should I take a sum?

text: I conducted an experiment to compare two computer reversi programs, namely $A$ and $B$, in which I had them play with each other 50 times and $A$ won $B$ 31 times. Now I would like to test if $A$ is stronger than $B$, i.e., the probability $p$ of $A$ winning $B$ (neither drawing or losing) is greater than 1/2, with significance level $\alpha$. According to Wikipedia's article on binomial test, I can compute $\frac{\binom{50}{31}+\dots+\binom{50}{50}}{2^{50}}$ and compare it with $\alpha$. What I do not see is why I should take a sum, and consider the probability of winning 31, 32, ... or 50 times. If I were dealing with a continuous random variable, it would surely make sense to take a sum (that is, to integrate) to obtain probability from probability density, because there is no such thing as the probability of the variable coinciding exactly with the observed value. But I am dealing with discrete variable and I can calculate the probability of winning exactly 31 times. So why should I take a sum (or do I really have to take a sum) ? EDIT: I corrected my mistakes and clarified my question.


title: wilson's adjustment for sample proportion

text: Could somebody provide a derivation for Wilson's adjustment for the CI of a sample proportion? The formula is: $$\tilde{p} \approx \frac{x+2}{n+4}$$ where $x$ is the number of successes, and $n$ is the sample size. Specifically, I want to know: **where do the 2 and the 4 come from?** Those numbers seem a bit arbitrary to me, and I've had no luck with Google. Just trying to scratch my curiosity here...


title: Naive Bayes classifier calculation

text: I'm trying to use naive Bayes classifier to classify my dataset.My questions are: 1- Usually when we try to calculate the likehood we use the formula: P(c|x)= P(c|x1) * P(c|x2)*...P(c|xn)*P(c) . But in some examples it says in order to avoid getting very small results we use P(c|x)= exp(log(c|x1) + log(c|x2)+...log(c|xn) + logP(c)). can anyone explain more to me the difference between these two formula and are they both used to calculate the "likehood" or the sec one is used to calculate something called "information gain". 2- In some cases when we try to classify our datasets some joints are null. Some ppl use "LAPLACE smoothing" technique in order to avoid null joints. Doesnt this technique influence on the accurancy of our classification?. Thanks in advance for all your time. I'm just new to this algorithm and trying to learn more about it. So is there any recommended papers i should read? Thanks alot.


title: 95% CI that contradict a paired samples t-test

text: I conducted a paired samples t-test to analyse the effect of an exercise intervention on blood pressure measurements. The test demonstrated a significant p-value. However, when I have calculated 95% confidence intervals for the mean difference between pre and post measures the result suggested there to be no significant difference between pre and post measures, i.e lower level minus and upper level positive. Is this possible?


title: Covariance estimation

text: Suppose I want to estimate the covariance of $n$ $p$-dimensional iid random vectors $X_i$, where $n>p$. I've read in several places that if $n-p$ is small then the MLE covariance matrix estimate will be poor (especially the small eigenvalues). However the MLE estimate of the mean will be good if $n$ is large (not necessarily bigger than $p$). Is there an intuitive reason why covariance estimation is harder than mean estimation?


title: How would one calculate a geographic estimate from a dataset of people?

text: Given you had a large dataset of people with a name and country. How would you calculate an estimate of the origin country from one of the names? If i looked up a name, Kim and just averaged the rows from each country, it seems to me that my numbers would be skewed since the dataset might contain a lot more names from a given country in the first place. How can i make up for this uncertainty?


title: Support Vector Machine(SVM) and log transformation

text: Why may log(natural logarithm) transformation improve results of SVM prediction( **regression** , eps-svm)? Is SVM based on the assumption of normal distribution or something else? update1. I use Radial basis function kernel.


title: nonparametric MARS regression

text: Most statistical methods assume homogeneous (outlier-free) data in which all data points satisfy the same model. However, real data are (not) NEVER homogeneous; and accurate identification of outliers plays an important role in statistical analysis. Moreover, the OLS classic regression needs to satisfy many assumptions (normality, homoscedasticity of the residuals, linearity of the regression function...). And with the OLS regression, we have as well to take care of the (multi)colinearity problem between the explanatory variables (measured by the Variance inflation factor (VIF)). My 2 questions are the following : 1) When we realize the nonparametric MARS (Multivariate adaptive regression splines) regression, do we have to care about the outliers problem, about the (multi)colinearity problem ? And do we have to care about the different assumptions (normality, homoscedasticity of the residuals and linearity of the regression function...) like for OLS regression ? Or when we practice the MARS regression, we don't have to care about all these ? 2) If there are some, what are the conditions of application of MARS ? I have read that there must be lots of observations ? For example, if I have a sample size of n=33 countries and 1 dependent variable and 3 explanatory variables, is it enough for MARS regression ? Best Regards, looking forward to reading You.


title: Ways of choosing test statistics: likelihood-based or not?

text: I just learned that there are three ways of constructing test statistics based on likelihood: likelihood ratio test, Wald test, and score test. I think it nicely categorizes various tests, and I liked it very much   1. I would like to know if the three can cover all or most of the specific tests. What are some commonly-seen examples of specific test tasks and test rules that belong to any of the three? Are there some test task and test rule which doesn't belong to any of them? Especially, do the nonparametric tests belong to any of them?   2. What are some other ideas or heuristics to construct test statistics, besides the three likelihood-based ones? By "other", they can be either likelihood-based or not. Thanks and regards!


title: Why Normalization (Standardization) values should be smaller than $1$?

text: The books gives some examples about content based recommendation. An example of what I understood is at below. > A movie's attributes are values between $1$ and $10$. The duration attribute > gets values between $1$ and $100$. If we use these raw values to calculate > distance, duration will dominate purely because of wider range, so we should > to normalize that value. Standardization formulas usually cause values smaller than $1$. It may between $1$ and $-1$. But if I scale values between $1$ and $10$ so how this normalization can be right? I expect that $1$ duration value should be represent $1$ and $100$ duration value should be represent $10$. But as you have know standardization formula cause smaller than $1$. Why it is so? Do I have to re-scale result for $1-10$ range again? For example if result is $0.43$ so that should represent $4$?