title: Does the first line of interaction output from lmer contain all levels of variables?

text: I am using LMM from lmer to interpret my data. As the system of its output is dummy, so it should mean that the number of levels is always the same in every line of output. For example, if I have three variables: color (red and white), sex (male or female), height (tall, low). The base line is always the second one in each variable. The output is like this:               Red     Male     Tall     Red:Male     Red:male:tall      Then does it mean that on the first line, the DV of red and female and low is higher than white and female and low. Or it simply means that DV of red is higher than white?


title: False discovery rate of a Bayesian classifier: scaling based on prior odds?

text: I am trying to assess the performance of my Bayesian classifier. One measure that I calculate is the false discovery rate (FDR): FP / (FP + TP), where FP = False Positive and TP = True Positive. The ratio (and sizes) of my positive and negative training sets do not accurately reflect the estimated 'real' ratios (and sizes). To clarify, I am working with sets of genes that have a certain property. It is estimated that there are approximately 18000 negative (not having that property) and 400 positive (having that property) genes (ratio 45:1). However, I only have 6000 negative and 200 positive examples (ratio 30:1) for which this property has been confirmed (my training set). I've read that if the sizes of the training sets do not reflect the prior odds, the FP an TP rate must be scaled in order to avoid over- or underestimating the false discovery rate. My questions are:   1. Do I have to scale/correct the FP and TP? If so, is this for the reason stated above (underestimation of the FDR) or some other reason?   2. If scaling is necessary, how can I achieve it?   3. Other performance measures, such as specificity and sensitivity, also use FP and TP (and TN&FN). Should I use corrected FP TP TN FN for these too? A more general question: I believe that what I refer to above as "estimated 'real' ratios (and sizes)" of the positive and negative sets are the prior odds. Is that correct? I am not a statistician and I have a hard time trying to find good resources on the subject. This is mainly because false discovery rate seems to be primarily associated with Benjamini-Hochberg multiple testing correction. I'm not sure if the two applications of false discovery rate are in fact the same.. Any help would be greatly appreciated.


title: Suggestion about morphometric analysis in R using lm or lme

text: Maybe this is an over-the-top question but I have many doubts regarding my recent analysis about deer skull measurements and how to proceed with the analysis. This is a sample of my dataset:                  Factor1 population manage foraging height biome abundance area  forest plough       -0.6033788 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000       0.3250981 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000       0.5577059 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      -0.1596194 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      -1.3089952 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      -2.1693392 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      -0.9669080 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      -1.8857842 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000       0.7242678 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000       1.6815373 ADA_BEC   best    fields   plain  agS    1500    73154  61154 12000      Factor1 are factor scores of all individuals (567) divided into 12 populations (population column). Others are either factor (manage (4lvl), foraging (3lvl), height (2lvl) and biome (4lvl)) or continuous, different for every population, (abundance, area (in ha), which totals forest+ploughland (also in ha)). Now I tried with all traditional uni-variate statistics such as anova, ancova, lm but, off course, my design is unballanced. My question is are there any general modeling solutions to incoporating all of these in a maximal model, simplifying it further and decide what is the most influential factor. Could mixed-effects model be used? Basic data are in fact, 50 measurements on every individual skull.


title: How to do binary logistic regression on people (couples) clustered within homes?

text: I am looking at the relationship between housing characteristics and a health outcome. To make the example simple, I have data for a continuous predictor (exposure) collected from 1000 homes and health outcomes S (a binary outcome) for 2000 people (1000 couples) living in each of those homes. I would like to look at the relationship between S and E using binary logistic regression. Apart from sharing the same exposure, there is no mechanistic reason to believe that status of partner 1 in the couple can affect the status of partner 2 e.g. its not a transmissible disease etc. Can I do an ordinary logistic regression? Or must I take into account the fact that people are clustered within homes? If so, why? What syntax would be appropriate in Stata, `xtlogit` with `i(house)`? or some kind of `xtmixed`? Many thanks


title: Correspondence analysis on a table of means

text: I have a table of the following kind:                    Cat1 cat2  cat3 ...     Var1 6.3  5.3   8.3     Var2 5.2  5.7   6.1     Var3 2.2  3.9   7.6     .     .     .      It is a table of means for different continuous variables across a categorical variable. My question is: I want a visualization of how the categories relate to the variables. For this purpose I have been using Correspondence analysis. I get a result that makes sense to me, but I am wondering if it is OK to use correspondence analysis on a table of non-categorical data?


title: Labeling boxplots in R

text: I need to build a boxplot without any axes and add it to the current plot (ROC curve), but I need to add more text information to the boxplot: the labels for min and max. Current line of code is below (current graph also). Thanks a lot for assistance.               boxplot(data, horizontal = TRUE, range = 0, axes=FALSE, col = "grey", add = TRUE)      The other solution is add the line from 0 to 1 (instead of x-axis), but I want it to go through the central line...for example like this graphic ![enter image description here](http://i.stack.imgur.com/goNiC.jpg)


title: Improving data analysis through a better visualization of data?

text: I ran four programs `a, b, c, d` _parallely_ on two different machines `X` and `Y` separately for 10 times. The below is a sample of the data. The running- times (milliseconds) in `10` runs of each program are given under their respective names.               Machine-X:          a b  c d     29 40 21 18     28 43 20 18     30 49 20 28     29 50 19 19     28 51 21 19     29 41 30 29     32 47 10 18     29 43 20 18     28 51 30 29     29 41 21 19               Machine-Y:          a b c d     16 24 19 18     16 24 19 18     16 23 19 18     16 24 19 18     16 24 19 18     16 22 19 18     16 24 19 18     16 24 19 18     16 24 19 18     16 24 19 18      I need to create graphs for visualizing the following:   1. Compare each program's performance (i.e. running-time) on both the machines X and Y.    2. Compare the variation in the running-times of each program on both the machines X and Y   3. Which machine is fair in providing computing resources to each program?   4. Compare the total running times (a+b+c+d) of the four programs in each run on both the machines X and Y.   5. Compare variation in the total running-times of the four programs in the 10 runs. For 1 and 2, I made Figure A, Figure B is for 3, and figure C is for 4 and 5. However, I am not satisfied because there are three graphs and it is difficult to fit all the three graphs in my paper. Moreover, I believe that we can produce better than these. I really appreciate if someone helps me to draw one or two nice graphs instead of three in R while satisfying my requirements. Please see below for the R code I have used to produce these graphs. Figure A: ![Figure A](http://i.stack.imgur.com/ZigVt.png) Figure B: X-axis shows the runs, Y-axis shows the running-times of the four programs in a particular run. ![Figure B](http://i.stack.imgur.com/GURrU.png) Figure C: ![Figure C](http://i.stack.imgur.com/ck18m.png) R Code               > pdf("Figure A.pdf")     > par(mfrow=c(1,2))     > boxplot(x,boxwex=0.4, ylim=c(15, 60))     > mtext("Time", side=2, line=2)     > mtext("Running times of each program in 10 runs", side=3, line=2, at=6,cex=1.4)     > mtext("Machine X", side=3, line=0.5, at=2,cex=1.1)     > boxplot(y,boxwex=0.4, ylim=c(15, 60))     > mtext("Machine Y", side=3, line=0.4, at=2,cex=1.1)     > dev.off()                    > pdf("Figure B.pdf")     > par(mfrow=c(1,2))     >  boxplot(t(x),boxwex=0.4, ylim=c(0,50))     > mtext("Run Number", side=1, line=2, at=12, cex=1.2)     > mtext("Fairness", side=3, line=2, at=12,cex=1.4)     >  mtext("Machine X", side=3, line=0.5, at=5,cex=1.1)     > boxplot(t(y),boxwex=0.4, ylim=c(0,50))     > mtext("Machine Y", side=3, line=0.4, at=5,cex=1.1)     > dev.off()               > pdf("Figure C.pdf")     > par(mfrow=c(1,2))     > barplot(t(x), ylim=c(0,150),names=1:10,col=mycolor)     > mtext("Run Number", side=1, line=2, at=14, cex=1.2)     > mtext("Total Running-Times in 10 Runs", side=3, line=2, at=14, cex=1.2)     > mtext("Machine X", side=3, line=0.5, at=5,cex=1.1)     > barplot(t(y), ylim=c(0,150), names=1:10,col=mycolor)     > mtext("Machine Y", side=3, line=0.5, at=5,cex=1.1)     > legend("topright",legend=c("a","b","c","d"),fill=mycolor,cex=1.1)     > dev.off()


title: Factor analysis: regression equation and residuals

text: I have used factor analysis with regression method. I was trying to develop an index using this regression equation by using the output of the component score coefficient matrix (SPSS) generated for each factor. Should I include the residuals found under the "Reproduced Correlations Table" (SPSS) in the regression equation?


title: I want to show $E(X)=\sum_{n=1}^{\infty}P(X\ge n)$

text: Let $X:\Omega \to \mathbb N$ be a random variable on probability space $(\Omega,\mathcal B,P)$ .show that $$E(X)=\sum_{n=1}^{\infty}P(X\ge n).$$ my definition from $E(X)$ is equal $$E(X)=\int_{\Omega}XdP.$$ Thanks.


title: Posterior distribution

text: I read the following in some document: > Let $Y$ be a random variable with distribution > $\mathcal{N}(\theta,\sigma^{2})$. The variance $\sigma^{2}$ is known. Let > $p(\theta) = 1$ a flat prior on $\theta$. Then, the posterior distribution > $p(\theta \vert Y)$ is $\mathcal{N}(Y,1)$. From Bayes' theorem, I know that: $$ \begin{eqnarray*} p(\theta \vert Y) & \propto & p(Y \vert \theta) p(\theta) \\ & \propto & \exp \Big( -\frac{1}{2\sigma^{2}} (\theta-Y)^{2} \Big) \\ \end{eqnarray*} $$ To me, the posterior distribution would be $\mathcal{N}(Y,\sigma^{2})$. I don't understand why the document says it is $\mathcal{N}(Y,1)$.


title: Rephrasing the risk assessment of a project

text: An energy saving project is being considered and its benefit is computed by subtracting the "post-project" energy cost from the "status quo" energy cost, e.g if the post-project energy cost is estimated to be £1m and the status quo energy cost is estimated to be £1,5m then the benefit is £0.5m. The net benefit is just the benefit minus the cost of the project. The decision makers are considering whether or not to fund the project. There is some uncertainty in the parameters and it turns out that **if the project is funded, there is a 90% chance of it leading to positive savings and 10% change of loosing money**. Without giving any insight into the actual numbers, does it make sense to rephrase the result in the way the **"if the project isn't funded, there is a 90% chance of loosing money"**?


title: Learning from unordered tuples?

text: EDIT: I reduced _my_ problem to a more specific question: http://math.stackexchange.com/questions/26573/ But I am still interested in other ideas. Let's say our data is generated by $$Y_i = f(X_i) + \epsilon_i$$ where $X_i$ are observed vectors, and $f$ is an unknown function. We know that $f$ is invariant with respect to _permutation_ of the elements of $X$. For example, if $X_i=[x_{i1},x_{i2},x_{i3}]$, then we have $$ f([x_{i1},x_{i2},x_{i3}]) = f([x_{i1},x_{i3},x_{i2}]) = f([x_{i2},x_{i1},x_{i3}])=\cdots $$ Are there modified versions of linear regression, support vector machines, forests, etc. which can be used to estimate $f$? I'm specifically interested in the case when $X_i$ are eigenvalues of matrices (so they have complex- valued entries). EDIT: A desperation move would be to make replicates of each data point with all permutations of each $X_i$ vectors and then apply standard methods, but this is clearly computationally impractical.


title: How to interpret results from a zero-inflated negative binomial model?

text: I have a zero-inflated negative binomial model (ZINB) for highly skewed, high zero-count data, $n=6800$. The null model does not reject the ZINB model, and Stata `count fit` also indicated ZINB over other count models. Yet, when I run the model with more than the minimum of covariates, the inflation intercept has a coefficient but no error, no $t$, no $p$-value. Does this mean the model is incorrect? Or just the zero-inflated assumptions? Any suggestions?


title: Basic R script/package with lots of automated calculations

text: Does that exist? A basic R script or package that just needs the data, and it spits lots of “basic” statistics (like the mean, the median, quartiles, make a plot, a boxplot, a histogram, a cumulative histogram, stem-leaf plot, $α_1$, $α_2$, $β_1$, $β_1$, the variance, etc.). I mean, you give the script lots of numbers (in a vector), and it gives you everything you need and more in a .txt, and a few pdfs with small plots. This could be useful to people who needs to do basic data treatment but doesn't want to click one thousand times to get everything with SPSS or Statgraphics. But (also) doesn't need to learn R language. In case it doesn't exist, I'm not asking you to make it, a plain “no” would be accepted as an answer. But in that case, it would be great idea if someone wants to write it :P


title: Main effect of the first independent variable in two-way ANOVa lost depending on the second independent variable

text: My problem is this: I have one dependent variable and 4 independent ones: one is age and the other three are temperament dimensions. I did 3 sets of two-way ANOVAs. The first independent variable is always the same (age) and the second is always different - one of the tempeament dimensions. In one case I get that age has significant effecet, the temeprament1 does not, and no interaction. in other case I get that age doesn't have the significant effecet, but still no influence of temperament2, and no significant effect of interaction. one-way ANOVA for age shows there is significant effect. My question is how do I interpret the data? My plan was to say if there is effect of age and temperament and interaction of age and temeprament dimensions, but now - the effect of age is sometimes there, and sometimes not?!


title: How to estimate vector autoregression & impulse response function with panel data

text: I am working on vector auto-regression (VARs) and impulse response function (IRFs) estimation based on panel data with 33 individuals over 77 quarters. How should this type of situation be analyzed? What algorithm's exist for this purpose? I would prefer to conduct these analyses in R, so if anyone is familiar with R code or a package designed for this purpose that they could suggest, that would be especially helpful.


title: What do you call $p(Y|X_1)$?

text: Given random vectors $Y= [Y_1, \dots, Y_m]$ and $X = [X_1, \dots, X_n]$, $p(y|x)$ is the conditional distribution of $Y$ given $X$. We can call $p(y_1|x)$ the _marginal_ conditional distribution of $Y_1$ given $X$. What would you call $p(y|x_1)$, the conditional distribution of $Y$ given $X_1$ then? Can it still be called " _marginal_ conditional distribution of $Y$ given $X_1$"? This is some terminology in Bayesian statistics, where you can think of the conditional distribution of $Y$ given $X$ as either the posterior distribution or the sampling distribution.


title: Probability distribution of income

text: I create an agent-based simulation of some economic phoenomenon, and I need to assign each citizen a random income level. I want to select the income level using a probability distribution that approximates the realistic distribution of income in a modern economy. What distribution should I use? Does it make sense to use exponential distribution? NOTE: I am not looking for actual census data of a certain economy, but for a more general approximation.


title: Testing for weak instrument: include intercept in regression of instrument?

text: When you want to use the IV (instrumental variable) estimator, you typically first test if you have a strong instrument. You do so by regressing the (endogenous) predictor against the instrument. With the regression coefficient, you can calculate the F statistic. My question is: do you have to include an intercept in this regression model? In all the formulas I find, there is no intercept included. But this would not be meaningful? Below is the regression of my endogenous variable on the instrument, with no intercept, which leads to a wrong conclusion about the regression coefficient? Thanks in advance. ![enter image description here](http://i.stack.imgur.com/3EC1X.png)


title: Imputing categorical variables before binarization

text: I wish to replace the missing values with mode of that categorical variable. In scikit-learn, we can something like               Imputer(strategy="most_frequent", axis=0)      but if I apply it on categorical column, an error message will be given since scikit-learn will try to convert the categories (expressed in strings) into float numbers. How can I do it?


title: Optimism bias - estimates of prediction error

text: The book Elements of Statistical Learning (available in PDF online) discusses the optimisim bias (7.21, page 229). It states that the optimism bias is the difference between the training error and the in-sample error (error observed if we sample new outcome values at each of original training points) (per below). ![enter image description here](http://i.stack.imgur.com/g1JCR.png) Next, it states this optimism bias (w) is equal to the covariance of our estimated y values and the actual y values (formula per below). i have trouble understanding why this formula indicates the optimism bias; naievely i would thought that a strong covariance between actual y and predicted y merely describes accuracy - not optimism. Let me know if somehow can help with the derivation of the formula or share the intuition. ![enter image description here](http://i.stack.imgur.com/eaBXO.png)


title: Modelling a random walk with varying distributions

text: I have some financial data i'm trying to fit a random walk too, but the daily change increments have different distributions when studied over the last month, last year, last 5 years etc, along with heteroskedascity. What is the best way to proceed with modelling this? Techniques and interesting areas to look into further?


title: Short term road traffic forecast modeling using neural networks toolbox in matlab

text: I have a hourly time series data of road traffic (i.e. count of the number of vehicles passing on a particular segment of road) collected over 7 days a week (Mon to Sun) for two weeks starting from 9 am in the morning to 10 pm in the night. I have (91 X 2 = 182) data points in my time series for two weeks duration. I wish to build a forecast model using neural networks which could forecast the count of no of vehicles for a given any one hour duration during any day of the second week(e.g. 9 am to 10 am on Tue) using the time series data of the first week. I am using neural networks toolbox inbuilt in MatLab. I wish to know the following things.   1. How to classify the data set that I have into training test, testing set and validation set (i.e the appropriate ratio which can be used initially)?    2. Is it the first week data set which i should classify into the 3 above mentioned sets?   3. Since I have only one array of time series in my problem, is it appropriate to use NAR model for my problem which is present in NN Toolbox or should I go with NARX model because I am splitting my data set into two parts (i.e. first week data set and second week data set)?    4. What should be the target set in my problem? Is it the traffic data set of the second week which is to be used? Since I am a beginner in this field, any help would be highly appreciated.


title: Is there an intuitive characterization of distance correlation?

text: I've been staring at the wikipedia page for distance correlation where it seems to be characterized by how it can be calculated. While I could do the calculations I struggle to _get_ what distance correlation measures and why the the calculations look as they do. **Is there a (or many) more intuitive characterisation of distance correlation that could help me understand what it measures?** I realize that asking for _intuition_ is a bit vague, but if I knew what kind of intuition I was asking for I would probably not have asked in the first place. I would also be happy for intuition regarding the case of the distance correlation between two random variables (even though distance correlation is defined between two random vectors).


title: Expressing a sum as a compound Poisson Distribution

text: I am given two Poisson distributions, $X_1$ and $X_2$ with parameters $\lambda_1=1$ and $\lambda_2=2$ respectively. Also given that $Y = -3X_1$ + $2X_2$, how do I represent $Y$ as a compound Poisson distribution? I know that I have to use the methods of convolution. However, I am entirely lost on this one.


title: Effect of one independent variable of several dependent variables – best strategy?

text: I have a question regarding which analysis strategy is best suited for our objective. In an exploratory study based on data from a survey we conducted ourselves in India, we are analyzing the consequences of a lack of suitable brides as a result of sex-selection. The determinants of sex-selection are well-researched and we ourselves also followed a separate objective where we analyzed these in separate analyses and regressions. What the consequences of a lack of suitable brides might be is still only limitedly researched as there is a lack of data. This means that we are essentially looking at an independent variable and are interested in exploring what other variables/aspects it has an effect on. We find bivariate correlations that make a lot of sense theoretically even though they are generally not very strong. Our variable of interest is a binary variable of whether there is a shortage of suitable brides as a result of sex- selective abortions. My main question is; is there any other suitable analysis strategy we could use where we are not simply limited to bivariate correlations? Furthermore, for some of the correlations we find there are issues of directionality and endogeneity. One example is crime, where crime could be a reason for parents to prefer and sex-select sons as sons could be thought to offer protection in an insecure environment where daughters could instead be perceived as liabilities. However, crime could also be a consequence of a shortage of brides and we do find such a correlation. From our separate analyses of the determinants of sex-selective behavior where we used logistic regressions with odds of births being sons as outcome, we know that crime is not a determinant of sex-selection. Our interpretation of this is that crime is to a certain degree a consequence of bride shortages but not a “cause” of sex-selective abortions. However, one of our variables that is negatively correlated with bride shortage is a measure of restrictions of women’s mobility, i.e. a measure of disempowerment of women. This suggests that a shortage of brides is related to a better situation for women. The difficulty is that this same variable also has, unlike our other variables, a significant effect on the odds of births being sons (once again from other regressions). Its effect on the odds of births being sons is positive, however, whereas it is negatively correlated with our variable for bride shortages. Our interpretation of this is that greater disempowerment of women is related to greater likelihood of sex-selection but that less disempowerment (and thus a better ‘situation’ for women) is a consequence of a bride shortage (and thus by extension a consequence of sex-selection)? My second question is; are we correct in drawing these conclusions? Is there a better way for us to make sure we are drawing the right conclusions? I would be very grateful for any help we could get. KML


title: Content-control (web filtering) using machine learning

text: I'm trying to build a content-control (web filtering) application using machine learning (just for training purposes). For example define gaming sites. I'm somewhat familiar with machine learning algorithms, but I have almost no idea how to get features from HTML pages. Can anyone suggest any good material? Maybe some resource with training data? I will be glad of any advice!


title: statistical analysis between two time series

text: I have two time series, which by looking at they appear to be correlated i.e. when one is high the other is high, and when one is low the other is low. Besides looking at the figure and saying these are correlated, I am looking for a statistical approach to assess whether these two time series are correlated or not. The image below shows the time series plotted against one another and a corresponding scatter plot. The data can be found here.![enter image description here](http://i.stack.imgur.com/fNgKa.png) I work primarily in MatLab but also sometimes use R. If anyone could provide some support in either programs that would be great. I basically want to analyse what is the relationship between the 2 series and is the series greater during some stages compared to others, but I cannot determine which method is most suitable. Thanks in advance.


title: Model position tracking error

text: I am building a whiteboard to track the position of its pen. I have `n` ultrasound receivers placed on the periphery of the whiteboard, and an ultrasound emitter in the pen. The pen emits pulses the receivers detect. A microprocessor gathers the arrival times of the pulses for each receiver, and makes an estimate of the position of the pen following an algorithm. The emitter/receivers are imperfect. They emit/receive the ultrasound with an error (described in the relevant datasheet), which causes an error in the measured arrival times. These timing errors can be modeled by a random variable. I want to optimise the placement of the receivers around the whiteboard, and optimise the algorithm, so that the computed position of the pen is close to the actual position. Of course, I could do live testing with a prototype, various receiver positions, and various algorithms, but that would be tedious. Is there statistical software where the timing error random variables can be specified, and where I can experiment with different receiver positions and algorithms to optimise the accuracy of the calculated pen position?


title: Compute probability of a grouping being correct

text: I have an exemplar grouping of objects (each with their own feature vector) into categories. I am then given a new grouping of compeltely different objects, and Iw would like to compute the probability of this grouping being 'correct'. I realise this is a tricky problem, so I have tried to break it down a little. My plan so far is to:   * For the training grouping, compute the mean variance of each feature within each cluster.   * Next, re-assign the training groups into random clusters, and again compute the mean variance of each feature.   * This should allow me to compute the probability of an unseen clustering being a true clustering, by comparing the variance of each feature with those in the 'true' and 'random' training examples Does this sound like a sensible method? I am a little wary of the concept of taking means of variances.


title: Correlation of bivariate grouped data?

text: Which test should I use, if I want to test the correlation between 2 bivariate grouped variables? The case is: I have asked several hotel owners about their feelings about the occupational rate of their hotel and the propagation of the region, where the hotel is. Both of the questions had ordinal scales, spreading from "low" to "high". Because of few participants of the survey I had to categorize the answers to both of the questions to two categories: "low" and "high". My hypothesis is, that there is a correlation between the occupational rate and the propagation of the region, so for example, the owner, who feels like his hotel has a high occupational rate, also feels like the propagation of the region is high. My question is, which test could I use to test this hypothesis. I would appreciate any kind of help! Thank you!


title: How to generate binary data using Cholesky decomposition in R?

text: m <- matrix(nrow=5, ncol=5)       m <- ifelse(row(m)==col(m), 1, 0.4)       ch <- chol(m)    # Choleski decomposition       u <- matrix(rnorm(2000*5), ncol=5)       uc <- u %*% ch       cr <- pnorm(uc)       cr <- qbinom(cr,1,0.5)       cor(uc)       cor(cr)         The cor(uc) is 0.4. However, I also expected the cor(cr) to be 0.4, but it is around 0.2.   Why the correlation is shrunken?


title: Problems with modeling a cumulative dependant variable

text: I am building a .NET program. One of its functions is to provide a predictive model for a vehicles life-to-date maintenance costs, basically what is the cumulative cost(Y) for a vehicle at specific year(X). I decided to use a 2nd degree polynomial least squares fit and for the most part it does a good job. Sometimes though the curve will peak and start trending downward which doesn't make sense for life-date-cost since its a cumulative value...(X,Y) > (X-1,Y). This negative trend happens when the difference in cost for say, year 2 to year 3 is less than year 1 to year 2. Some sample data that gives me a negative trend: (1,328.76) (2,1133.12) (3,1366.07) My solution for now is to check for a negative trend and if its found use a linear best fit instead but I feel like that's a messy fix. I've thought about implementing some sort of minimum value for the change from year to year...essentially turning the curve into a linear line at a certain X value but that seems complicated to implement. Does anyone see a better way of doing this or a better model to use? I'm not very knowledgeable with statistics so go easy on me :-p **Edit** Each vehicle has a varying amount of data depending on how long its been in service, with a soft max at 15 years. So the last data point for each vehicle is for the most recent year(2011 in this case) and we are really only interested in extrapolating 5 years beyond that point. As we use the model year to year, we will get more data for the vehicles which require the model to be altered. Thats why I choose the polynomial least squares fit because its easy to just run the new data back through that function and get a new equation.


title: GAM cross-validation to test prediction error

text: My questions deals with GAMs in the **mgcv** R package. Due to a small sample size I want to determine the prediction error using leave-one-out cross- validation. Is this reasonable? Is there a package or code how I can do this? The `errorest()` function in the **ipred** package does not work. A simple test dataset is:               library(mgcv)     set.seed(0)     dat <- gamSim(1,n=400,dist="normal",scale=2)     b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)     summary(b)     pred <- predict(b, type="response")      Thank you very much for your helping hand!


title: Appropriateness of Wilcoxon signed rank test

text: I've poked around a bit in the Cross Validated archives and haven't seemed to find an answer to my question. My question is the following: Wikipedia gives three assumptions that need to hold for the Wilcoxon signed rank test (slightly modified for my questions): Let Zi = Xi-Yi for i=1,...,n.   1. The differences Zi are assumed to be independent.   2. (a.) Each Zi comes from the same continuous population, and (b.) each Zi is symmetric about a common median;   3. The values which Xi and Yi represent are ordered...so the comparisons 'greater than', 'less than', and 'equal to' are useful. The documentation for ?wilcox.test in R, however, seem to indicate that (2.b) is actually something that is tested by the procedure: "...if both x and y are given and paired is TRUE, a Wilcoxon signed rank test of the null that the distribution ... of x - y (in the paired two sample case) is symmetric about mu is performed." This sounds to me as though the test is performed for the null hypothesis that "Z is distributed symetrically around median mu=SomeMu" -- such that rejection fo the null could be _either_ a rejection of the symmetry _or_ a rejection that the mu around which Z is symmetric is SomeMu. Is this a correct understanding of the R documentation for wilcox.test? The reason this is important, of course, is that I am conducting a number of paired-difference tests on some before-and-after data ("X" and "Y" above). The "before" and "after" data individually are highly skewed, but the differences are not skewed nearly as much (although still skewed somewhat). By that I mean that the "before" or "after" data considered alone has skewness ~7 to 21 (depending on the sample I am looking at), while the "differences" data has skewness ~= 0.5 to 5. Still skewed, but not nearly as much. If having skewness in my "differences" data will cause the Wilcoxon test to give me false/biased results (as the Wikipedia article seems to indicate), then skewness could be a big concern. If, however, the Wilcoxon tests are actually testing whether the differences distribution is "symmetric around mu=SomeMu" (as ?wilcox.test seems to indicate) then this is less of a concern. Thus my questions are:   1. Which interpretation above is correct? Is skewness in my "differences" distribution going to bias my Wilcoxon test?   2. If skewness _is_ a concern: "How much skewness is a concern?"   3. If the Wilcoxon signed rank tests seem grossly inappropriate here, any suggestions for what I should use? Thanks so much. If you have any further suggestions about how I might do this analysis I am more than happy to hear them (although I can also open another thread for that purpose). Also, this is my first question on Cross Validated; if you have suggestions/comments on how I asked this question, I am open to that as well! * * * A little background: I am analyzing a dataset that contains observations on what I'll call "errors in firm production." I have an observation on errors occuring in the production process before and after a surprise inspection, and one of the goals of the analysis is to answer the question, "does the inspection make a difference in the oberved number of errors?" The data set looks something like this:               ID, errorsBefore, errorsAfter, size_large, size_medium, typeA, typeB, typeC, typeD     0123,1,1,1,0,1,1,1,0      2345,1,0,0,0,0,1,1,0     6789,2,1,0,1,0,1,0,0     1234,8,8,0,0,1,0,0,0      There are roughly 4000 observations. The other variables are catagorical observations that descrie characteristics of the firms. Size can be small, medium, or large, and each firm is one and only one of those. Firms can be any or all of the "types." I was asked to run some simple tests to see if there were statistically significant differences in observed error rates before and after the inspections for all firms and various sub-groupings (based on size and type). T-tests were out because the data was severely skewed both before and after, for example, in R the before data looked something like this:               summary(errorsBefore)     # Min.  1st Qu.  Median   Mean  3rd Qu.    Max     # 0.000  0.000    4.000  12.00    13.00  470.0      (These are made up -- I'm afraid I can't post the actual data or any actual manipulations of it due to proprietary/privacy issues -- my apologies!) The paired differences were more centralized but still not very well fit by a normal distribution -- far too peaked. Differences data looked something like this:               summary(errorsBefore-errorsAfter)     # Min.   1st Qu.  Median   Mean  3rd Qu.    Max     # -110.0  -2.000   0.000  0.005   2.000   140.0      It was suggested that I use a Wilcoxon signed rank test, and after a brief persusal of ?wilcox.test and Wikipedia, and here, this seems like the test to use. Considering the assumptions above, I believe (1) is fine given the data generating process. Assumption (2.a) is not strictly true for my data, but the discussion here: Alternative to the Wilcoxon test when the distribution isn't continuous? seemed to indicate that this wasn't too much of a concern. Assumption (3) is fine. My only concern (I believe) is Assumption (2.b). **One additional note** , some years later: I eventually took an excellent non-parametric stats course and spent a lot of time on the rank-sum tests. Embedded in assumption (2.a), "Each Zi comes from the same continuous population", is the idea that both samples mush come from populations with _equal variance_ -- this turns out to be _extremely_ important, practically speaking. If you have concerns about differing variance in your populations (from which you draw the samples), you should be concerned about using WMW.


title: Logarithms and regression

text: I have a regression where I am trying to predict Y, using several variables. My variables are A, B, C. ![enter image description here](http://i.stack.imgur.com/9Qg6q.png) I'm using `R` for my analysis. My equation is as follows:               mod <-lm(  log(Y/B) ~ log(A/B) + log(A/C), data = dataframe)      and a summary of the model is as follows:               Coefficients:                           Estimate Std. Error t value Pr(>|t|)         (Intercept)          -1.326698   0.080683  -16.44   <2e-16 ***     log(A/B)              0.756779   0.007507  100.81   <2e-16 ***     log(A/C)             -0.292373   0.012762  -22.91   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.6925 on 17063 degrees of freedom     Multiple R-squared: 0.4043, Adjusted R-squared: 0.4042      F-statistic:  5790 on 2 and 17063 DF,  p-value: < 2.2e-16       I understand that `log(j/k) = log(j) - log(k)`; so `log(Y/B)` can be written similarly. Is there a way that I can include the `B` variable in the dependent part of my regression equation? Also, it is acceptable to use a variable for normalizing, on both dependent and independent variables? I have seen other examples where the same variable is used in this way, but I was not certain if this makes sense.


title: Sample size calculation for ROC/AUC analysis

text: As a background, I am not familiar with stats except on a basic level. I have been tasked with doing some analysis that is out of my comfort zone. I am trying to figure out how to compute necessary sample sizes for an ROC analysis based on desired statistical power. I have read Hanley and McNeil's paper (http://www.med.mcgill.ca/epidemiology/hanley/software/hanley_mcneil_radiology_82.pdf) and understand it fairly well. The relevant portion is the (over)estimate of standard error: $$ SE(W) = \sqrt{\frac{A(1-A) + (n_a - 1)(Q_1 - A^2) + (n_n-1)(Q_2-A^2)}{n_a n_n}}$$ where $A$ is the expected AUC, $n_a$ is the number of abnormal samples, $n_n$ is the number of normal samples, and $Q_1 = \frac{A}{2-A}$, $Q_2 = \frac{2A^2}{1+A}$. Using this, I can estimate confidence intervals fairly easily ($A\pm 1.96\cdot SE$ for $95\%$ confidence level). The part I am unsure about is how the extra parameter of statistical power comes into play. Let's say, for example, that I want a statistical power of $80\%$ when I have an expected $A$ value of 0.8. How would I go about finding suitable $n_a,n_n$ pairs?


title: How to measure the similarity among categories?

text: Suppose we have several documents. These docs are classified into several categories. But there could be issues like these:   * The categories may not be properly defined.   * Or The categories are properly defined, but the docs are not properly classified. Given the category definitions and some classified docs, I want to tell **how well the situation is**. I have tried to use the probability distribution of words to describe **the current status** of each category. And calculate the dissimilarity of different distributions through Kullback-Leibler Divergence. But I went into some issue here: Calculate the Kullback-Leibler Divergence in practice? I am wondering if the above way is proper. Is there a canonical way to do it? (Pardon me if the tags I assigned are not proper.)


title: What are good references for the different methods of combining p-values?

text: I need a good reference for the methods and/or the difficulties that arise when attempting such a combination. I've found Loughin, TM, A systematic comparison of methods for combining p-values from independent tests, _Computational Statistics & Data Analysis_ (2004) 47(3):467–485, for example. Is there a better/more precise one?


title: Distributed Median

text: What is the most optimum method to calculate the median in a distributed setting? I want to optimize with respect to data transmission i.e. an algorithm that finds the median across data spread across multiple nodes by sending the least amount of data


title: How to calculate a confidence interval for estimating a proportion when the population is very small?

text: I plan to execute a little study on whether dentists tend to grant privately insured patients (PRIV) a sooner appointment compared with publicly insured patients (PUB). I want to restrict this study to the dentists of a single town in Germany. In this town there are about 500 dental practices and I am willing to call maximum 100 practices (twice, once as a PUB and once as PRIV). (A sample without replacement of size up to 100) So the result of this examination would be a list of up to 100 practices with an associated value of 1 - if the appointment for the PRIV was more than two days sooner - and a 0 - if not. The proportion of 1s is what I want to estimate. Now the sources I have been reading always assumed a much larger population with a finite population correction close to 1. This is why I doubt that blindly using a conservative or approximate confidence interval would lead to a statistically meaningful result.   1. So my first question is, of course, whether it is reasonable at all trying to estimate the proportion using a confidence interval in this case?   2. And if yes, then what statistical technology/method (leading to a CI) would be a good choice for this use case? One idea that comes to mind is using a sample size of 25 f.x., in that case the FPC would be 0.95 - "close enough" to 1. But it seems stupid trying to improve the quality of an estimation by reducing the information used for the estimation.


title: Is there a difference between an index score and a composite score?

text: Following on from my previous question on forming scale scores, my sample has 100 people who have answered 10 likert scale questions each on two academic subjects: Maths and English. The likert scale is a 5 point scale, going from not important to absolute importance. I have aggregated the scores:   * The maximum score a person can get is 40 (which is 10 responses x 4 which is the maximum score on the likert scale e.g. 4 = absolute importance).    * The minimum is 0 (which is 10 responses x 0 which is the minimum score on the likert scale e.g. 0 = not important)  If a person gets 20 points (total of 10 likert scale responses for Maths) and 30 points (total of 10 likert scale responses for English), my questions are: **1. Is 20 or 30 a composite score or an index score?** (given that they are both based on subjective weighting aligned with the likert scale response) **2. Can I do a correlation analysis on the two scores?** (given they are composite/index scores, hence represent a latent variable) Many thanks


title: Kaiser's eigenvalue or MAP/parallel analysis?

text: i ran exploratory factor analysis and according to eigenvalue > 1, 41 factors were extracted out of 142 items. but when i ran MAP and parallel analysis, 16 factors were prescribed as a proper number of factors to retain. now what should I do? With 16 factors, my total variance is 51% and with 41 factor, total variance is 75% but with too many trivial items/factors. Which method should I adhere to? Eigenvalue or MAP and parallel analysis? My second question is that i read about another method named comparison data (Ruscio and Roche, 2012) which outperform other methods, but i cant find any information about how it works, which software should I use or any manual about it?


title: How to test whether variables are more frequently associated than expected by chance?

text: I have this presence/absence matrix (S2 is missing on purpose).               dat <- read.table(text="S1  S3  S4  S5  S6  S7  S8  S9  S10     0   0   1   1   0   0   1   0   0     0   0   1   1   0   0   1   0   0     0   0   1   0   0   0   0   0   0     0   0   0   1   0   0   0   0   1     0   0   1   1   0   0   1   0   0     0   0   1   1   0   0   0   0   0     0   0   1   1   0   1   1   0   0     0   0   0   0   0   0   0   1   0     0   0   1   1   0   0   1   0   0     0   0   1   0   0   0   0   0   1     1   0   0   0   1   0   0   0   0     0   0   1   0   0   0   0   0   0     0   1   0   0   0   1   0   0   0     0   0   0   1   0   0   1   1   0     0   0   0   0   0   0   0   1   0     0   0   1   1   0   0   1   0   0     1   1   0   0   0   1   0   1   1     1   1   1   0   0   0   0   1   1     0   0   0   1   0   0   1   0   0     0   0   0   1   0   0   1   0   0     0   0   0   1   0   0   0   0   0     0   0   1   1   0   0   1   0   0     0   0   0   0   0   0   1   0   0     1   0   0   0   0   0   0   0   0     0   0   1   0   0   0   0   0   0     1   0   0   0   0   0   0   0   0     0   0   0   1   0   0   1   0   0     1   0   0   1   0   0   0   0   0     0   0   0   0   0   0   1   0   0     1   0   0   0   0   0   0   0   0     0   0   0   0   0   0   1   0   0     0   0   0   0   0   0   1   0   0     0   0   1   1   0   0   1   0   0     0   0   1   1   0   0   0   0   0     0   0   0   1   0   0   0   0   0     0   0   0   0   0   0   1   0   0     0   0   0   0   0   0   1   0   0     0   0   1   1   0   0   0   0   0     0   0   0   1   0   0   0   0   0     0   0   0   1   0   0   1   0   0     1   0   1   1   0   0   0   0   0     0   0   0   1   0   0   0   0   0     0   0   1   1   0   0   0   0   0     0   0   1   0   0   0   0   0   0     0   0   0   1   0   0   0   0   0     0   0   1   1   0   0   0   0   0     0   0   1   0   0   0   0   0   0     0   0   0   0   0   0   1   1   0     0   0   1   0   0   0   0   0   0     1   1   0   0   0   0   0   1   1     0   0   1   1   0   0   0   0   0     0   0   1   0   0   0   1   0   0     0   0   1   1   0   0   0   0   0     0   0   0   0   0   0   1   0   0     0   0   1   1   0   1   1   1   0", header=T)      On looking at the number of times a particular variable is found with another, it seems to me that S4, S5 and S8 are frequently found together. Is there a way of testing this statistically either for pairs of variables or for triplets? Cheers


title: How to combine unbound variables with very different frequency distributions?

text: I want to combine three unbound variables. Each variable is the score provided by three different algorithms. Each algorithm predicts the likelihood (score) that a specific interaction between two molecules takes place. I have a gold standard with which I can rate the performance of each algorithm (that is, how many predicted interactions actually take place in real life). My goal is to come up with a combined score for each interaction, using the scores from the individual algorithms. Ex: interaction1 is predicted by algorithm1: $756$ score, algorithm2: $0.7$, algorithm3: $6.5$ Each has different samples sizes and frequency distributions. Example:               score1: from 1000 to 0, huge peak around 895     score2: from 1 to .5, looks like a standard distribution     score3: from 10 to 0, huge peak around .1     ...      My first idea was to normalize each score using $\frac{X-min(X)}{max(X)-min(X)}$, and then do a weighted average with each score: $combinedScore=weight1*score1 + weight2*score2 + weight3*score3$. I get the weights from another variable. Should I use z-scores, or the logistic function to normalize? Am I doing something wrong?


title: Find posterior distribution

text: Let $X_{1},..,X_{n}$ be a sample from a poisson$({\lambda})$ distribution. Let the prior be ${\pi}({\lambda})=1/{\sqrt{\lambda}}$. Find the posterior distribution. My work: We have $f(x|{\lambda})=\frac{{\lambda}^{x_{1}+..+x_{n}-0.5}e^{-{\lambda}n}}{x_{1}!....x_{n}!}$. To find the marginal density I believe I need to sum this expression through ${\lambda}=0,1,2,......$. How to do this? Thanks


title: How to determine the proper sets or subsets of outcomes when autocorrelating frequency of each from year to year

text: I am trying to determine whether a baseball player has the ability to influence the outcome of an "at bat" by calculating the correlation coefficient between the percentage of at bats which result in a given outcome from one year to the next, for each player-year in the sample. For example, suppose that an at bat can result in one of the following outcomes:   * The player strikes out (SO)   * The player puts the ball "in play" by hitting a ground ball, a line drive, or some other type of batted ball (BIP) Let's further suppose that the correlation coefficient between strikeouts per at bat (SO%) from one year to the next is found to be 0.85. Based on this, we can conclude that a player's ability to put the ball in play and avoid striking out is a "skill", of which different players have varying degrees. ### My question is, how can I take it a step further and determine whether players have the ability to influence which _type_ of batted ball they put in play? For example, suppose instead of starting with a set of two possible outcomes (strikeout or ball-in-play), we started with four:   * The player strikes out (SO)   * The player hits a ground ball (GB)   * The player hits a line drive (LD)   * The player hits a fly ball (FB) The difficulty arises when we try to calculate the correlation coefficient between line drives per at bat (LD%). Since a LD is a subset of BIP, and BIP% is well correlated from year to year (recall we found r = 0.85), LD% will also be well correlated from year to year (although possibly to a lesser degree). In fact, we would expect to see some correlation even if the number of line drives is a totally "static" percentage of the number of balls in play and has nothing to do with a player's skill or propensity to bat "frozen ropes" around the ball park. The only solution I can think of is to calculate r for LD per BIP instead of per at bat. However, the question now becomes: how do we know whether to treat line drives as a percentage of _all_ balls in play? What if line drives are a random occurrence of "air" balls, any number of which can end up a fly ball or a line drive? Or what if we find that the correlation for LD per at bat is higher than that of SO%? Do we then treat LD as the independent outcome and calculate SO as a percentage of SO+GB+FB? That doesn't even make sense. I guess what I'm trying to ask is, is there a systematic way to "tease out" which sets or subsets of outcomes are a result of player skill and, if so, how do we decide which of these sets to include in our correlation? I'm not even sure if what I'm saying makes sense so please let me know if anyone has any questions.


title: Model parameters, statistical significance and errors

text: I am fitting a basic generalized linear model (binomial) to some survey data that has a small sample size (N=376). The model my colleague has fit has 23 model parameters. Two or three have p-values less than 0.1.   There are theoretical reasons to expect that those particular independent variables are linked to the dependent variable, so that's good. But, with so many independent variables included in a model, particularly with a reasonably small sample size, is it not the case that the risk of a type I error increases? In other words, putting so many model parameters into the model, are we not increasing the chances that some parameters will present themselves as statistically significant, just out of random chance? To test this, would it be generally adviseable to begin removing non- statistically significant independent variables from the model to see if the remaining parameters remain significant? I think I'm sensitive to the danger of chasing the stars of statistical significance, but I want to be sure my logic and modelling strategy is correct. S.


title: Neural Network: What if there are multiple right answers for a given set of inputs?

text: For a given input into the input nodes, there are multiple correct values for the output nodes. In the training set, there are times when the inputs result in a certain output, and other times when the inputs result in a completely different (but equally valid) output. Will a neural network still be able to "figure out" a pattern? From what I know about backpropagation, it seems like the different right answers would prevent it from functioning properly. If that's so, are there any solutions? I don't need the neural network to predict all possible correct solutions, but I do need it to output _a_ correct solution.


title: Evaluation measures of overlapping clustering

text: I have a dataset of Facebook users and a set of different clustering algorithms. The project goal is to draw up a rank between these algorithms in order to understand which of them are the good ones. The only information I have is the users dataset and nothing else: no ground-truth, no class / truth labels or whatever. Moreover, each algorithm is designed to create an overlapping clustering, i.e., a given user can belong to more than one cluster. Since the ground truth is not available, I must use _internal indices_ to evaluate the clustering quality. According to this post: Evaluation measure of clustering (without having truth labels), we can find a lot of indices of this kind in the literature (see this paper for more details). The problem is: all these indices are not useful when the algorithms return an overlapping clustering. Do you know some indices to evaluate the quality of overlapping clustering?


title: What ANOVA non-superiority test corresponds to an omnibus set of TOST equivalence tests?

text: I am not sure if such a beast exists, but is there some non-superiority test for a oneway ANOVA of $k$ groups, probably based on a $\Delta^{2}$ equivalence threshold, for $m=\left[k(k-1)\right]/2$ _post hoc_ pairwise tests for equivalence (using the two one-sided tests procedure) based on an equivalence threshold of $\Delta$? (Where $\Delta$ is in units of the score being tested; my question could be reframed in terms of an equivalence threshold in units of the test statistic distribution, as in $\varepsilon = \Delta/s_{\bar{x}}$.) **Edit: turning this into self-study** On the one hand, I might argue that the _F_ statistic is a ratio two variances of the variable being tested, and that since variances are in squared units, that I could simply test something like H$_{0}\text{: }F \ge \frac{\bar{n}\Delta^{2}}{s^{2}_{\text{within}}}$, and reject H$_{0}$ if P$(F \le \frac{\bar{n}\Delta^{2}}{s^{2}_{\text{within}}}) < \alpha$. Is this naïve? On the other hand: Relationship between _t_ and _F_ tests $\begin{eqnarray} t_{\nu} & = \frac{\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}}\text{, where }\nu=N-1.\\ t^{2}_{\nu} & = \frac{\mathcal{N}^{2}\left(0,1\right)}{\chi^{2}_{\nu}/\nu}\\ & = \frac{\chi^{2}_{1}/1}{\chi^{2}_{\nu}/\nu}\\ & = F_{1,\nu}\\ \end{eqnarray}$ TOST for a _t_ test for equivalence, given symmetric equivalence threshold $\Delta$: $\begin{eqnarray} t_{1\nu} & = & \frac{\Delta - \mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}}\\ & = & \frac{\Delta}{\sqrt{\chi^{2}_{\nu}/\nu}} - \frac{\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}}\\ t^{2}_{1\nu} & = & \left(\frac{\Delta}{\sqrt{\chi^{2}_{\nu}/\nu}} - \frac{\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}}\right)^{2}\\ & = & \left(\frac{\Delta}{\sqrt{\chi^{2}_{\nu}/\nu}}\right)^{2} - \frac{2\Delta\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}} + \left(\frac{\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}}\right)^{2}\\ & = & \frac{\Delta^{2}}{\chi^{2}_{\nu}/\nu} - \frac{2\Delta\mathcal{N}\left(0,1\right)}{\sqrt{\chi^{2}_{\nu}/\nu}} + \frac{\mathcal{N}^{2}\left(0,1\right)}{\chi^{2}_{\nu}/\nu}\\ & = & \frac{\Delta}{\sqrt{\chi^{2}_{\nu}/\nu}}\left(\frac{\Delta}{\sqrt{\chi^{2}_{\nu}/\nu}} - 2\mathcal{N}\left(0,1\right)\right) + F_{1,\nu} \end{eqnarray}$ But I feel like I am kinda falling down a rabbit hole with this approach, and don't know how to translate that to $k$ number of groups >2.


title: Taking the expectation of Taylor series (especially the remainder)

text: My question concerns trying to justify a widely-used method, namely taking the expected value of Taylor Series. Assume we have a random variable $X$ with positive mean $\mu$ and variance $\sigma^2$. Additionally, we have a function, say, $\log(x)$. Doing Taylor Expansion of $\log X$ around the mean, we get $$ \log X = \log\mu + \frac{X - \mu}{\mu} - \frac12 \frac{(X-\mu)^2}{\mu^2} + \frac13 \frac{(X - \mu)^3}{\xi_X^3}, $$ where, as usual, $\xi_X$ is s.t. $|\xi_X - \mu| < |X - \mu|$. If we take an expectation, we will get an approximate equation which people usually refer to as something self-apparent (see the $\approx$ sign in the first equation here): $$ \mathbb{E}\log X \approx \log \mu - \frac12 \frac{\sigma^2}{\mu^2} $$ _QUESTION_ : I'm interested in how to prove that the expected value of the remainder term is actually negligible, i.e. $$ \mathbb{E}\left[\frac{(X - \mu)^3}{\xi_X^3}\right] = o(\sigma^2) $$ (or, in other words, $\mathbb{E}\bigl[o(X-\mu)^2\bigr] = o\bigl(\mathbb{E}\bigl[(X-\mu)^2\bigr]\bigr)$). _What I tried to do_ : assuming that $\sigma^2 \to 0$ (which, in turn, means $X \to \mu$ in $\mathbb{P}$), I tried to split the integral into two, surrounding $\mu$ with some $\varepsilon$-vicinity $N_\varepsilon$: $$ \int_\mathbb{R} p(x)\frac{(x-\mu)^3}{\xi_x^3} \,dx = \int_{x \in N_\varepsilon} \ldots dx + \int_{x \notin N_\varepsilon} \ldots dx $$ The first one can be bounded due to the fact that $0 \notin N_\varepsilon$ and thus $1/\xi^3$ doesn't bother. But with the second one we have two concurring facts: on the one hand $$ \mathbb{P}(|X - \mu| > \varepsilon) \to 0 $$ (as $\sigma^2 \to 0$). But on the other hand, we don't know what to do with $1/\xi^3$. Another possibility could be to try using the Fatou's lemma, but I can't figure out how. Will appreciate any help or hint. I realize that this is sort of a very technical question, but I need to go through it in order to trust this "Taylor-expectation" method. Thanks! P.S. I checked out here, but seems it's a bit of another stuff.


title: Looking for a good followup to the Stanford AI course

text: I took the original Stanford AI course and got an immense amount out of it. I got a poor score since I did not have a great deal of time to work out all the problems.. But what I really enjoyed about the course was getting a grasp of the new techiques, way beyond the conventional statitistcs.. I'm looking for a way to get to the next level. The course really did not provide a good reference for reviewing the subject. The big book by the course authors (Artificial Intelligence: A Modern Approach by Norvig and Russell) (which I have) is too complex for me. I try to do all my work in Python and I am making some headway in getting the tools, Sci-py, scikit learn, NLTK. Pandas (am currently reading Data Analysis with Python, by Mckinney, excellent). and more.. But I'm in need of references to help me leap from the conceptual level to actual implementation , including proper tool selection and problem selection and definition, that is, sort of a cook book approach. This question may be too open ended,, but it expresses my dilemma. The whole space is quite open ended, and I'm looking for references to help me navigate through it. What references (accessible) might you suggest?


title: 10-fold cross validation

text: (5 points) Assume that we are interested in generating a model (e.g., a decision tree) from a sample of examples of a specific size drawn from some distribution. Assume further that we would like to investigate how sensitive the resulting model is to the actual choice of training examples (i.e., how the performance varies over different sets of training examples of the specific size). Assume that we have access to 100 training examples drawn from the underlying distribution. If we are interested in investigating how the performance varies for models generated from 90 examples, would we obtain a reliable estimate of the variance of the model performance by performing a 10-fold cross-validation? Motivate your answer. Is the 10-fold cross validation reliable for estimating the variance of model performance?


title: Independence in a sum relationship

text: I have a model of total reaction time T, which is a composite of a selection time S and a discrimination time D. So a person first finds something, this takes the tS. Then he discriminates and reports an attribute of this thing, this takes the time D. We only measure the total time T, though. So, my model for T is: T = S + D In considering, I try to think about how these random variables behave in simplified settings; and these are my questions: 1.) If I would define S and D to be independent, could I infer that S and T, or D and T are independent? If not, could I infer that D and T are conditionally independent, given S? 2.) Defining S and D to be independent and normally distributed, I know I can infer that T is normally distributed. But what if I define T and D as normally distributed, can I infer that S is normally distributed? I think this boils down to the question 1.) whether I can say that T and D are independent.


title: Is it possible to have a variable that acts as both an effect modifier and a confounder?

text: Is it possible to have a variable that acts as both an effect (measurement) modifier and a confounder for a given pair of risk-outcome associations? I'm still a little unsure of the distinction. I've looked at graphical notation to help me understand the difference but the differences in notation are bewildering. A graphical/visual explanation of the two and when they may overlap would be useful.


title: Efficiency as a dependent variable

text: I have a dataset, consisting of two independent variables each with two levels. The independent variables are individuals (1 vs. 2) and light (day vs. night) The dependent variable is the efficiency to detect an object (in other words the proportion of objects detected in a known population of objects). I want to find out whether there is a difference between individuals and between amount of light. I wonder which statistical test is should use for the analysis, is it possible to use a two-way ANOVA?


title: Hyperplane problem in linear congruent generator

text: From Wikipedia > if an LCG is used to choose points in an n-dimensional space, the points > will lie on, at most, $m^{1/n}$ hyperplanes (Marsaglia's Theorem, developed > by George Marsaglia). This is due to serial correlation between successive > values of the sequence $X_n$. ![enter image description here](http://i.stack.imgur.com/BPxy4.gif)   1. Mathematically, what is the serial correlation between successive values of the sequence $X_n$?   2. I was wondering how the points within one period are distributed to different hyperplanes, as they are generated one after another in a sequence? I guess it is unlikely that the points firstly fill out one hyperplane, and then fill out the next one, and will never visit the hyperplanes visited before in the single period?  Is the time interval between every two consecutive visits to each fixed hyperplane fixed, and same for all hyperplanes? Is it connected to the period of the LCG? Does taking a small front portion of a sequence of a full period somehow overcome this drawback of LCG? That is the reason why I ask the above questions.   3. Note I think there is a typo. "at most, $m^{1/n}$ hyperplanes" should be "at most, $(n!m)^{1/n}$ hyperplanes". Am I correct? Thanks!


title: How to get odds ratio in a table in SPSS?

text: I have survey data from students from two Universities. One of the questions in the survey had 4 response categories and students were allowed to check only one (categorical variable). I want to make a table that shows the odds ratio and 95% CI for OR with University arranged in columns and all 4 possible answers in rows. I also want to see if there is a statistical difference between ORs in rows. Is it possible? If so, how should I do it?


title: How to perform Fisher exact test in SPSS?

text: I've collected around 10 numbers of ants from a colony and introduced the same numbers of ants to the same colony (say Colony-1) at their entrance to see their behaviour (whether they accepts or rejects). Again I introduced 10 ants from different colony (say Colony-2, i.e. non-nestmate) to Colony-1 and checked their rate of acceptance and rejectance. Likewise I introduced another 10 ants from another colony (say colony-3) to colony-1 to see the same.   The data I got looks like that:                                *In Colony-1*         Accepted         Rejected      Colony          8                  2         Colony-1          3                  7         Colony-2          2                  8         Colony-3      Now my question is that how I should arrange these data in SPSS 11.5 and run Fischer exact test in order to get the P-value and percentage of acceptance and rejection in Colony-1?


title: How to explain the connection between SVD and clustering?

text: Is there an intuitive explanation for how SVD is related to co-clustering when performing SVD on a covariance matrix? (i.e. the SVD is performed on the matrix $E[X Y^{\top}]$ where $X \in \mathbb{R}^n$ and $Y \in \mathbb{R}^m$.) How does the projections $U$ and $V$ in the decomposition $U \Sigma V^{\top}$ relate to clustering? If I take some instance of $X$ and project it through $U$, what do I get? (note I am talking about thin svd, so that $\Sigma$ is of dimensions lower than $n$ and $m$, and $U^{\top}x$ for some instance $x$ of $X$ is projected into that lower dimension.)


title: Nonparametric Pearsons's r-value and p-value

text: Pearson's R Value shows the correlation of two arrays of real numbers. I would like to calculate an correlation value for text which is associated with numbers then a P value to judge the statistical significance of this relationship. Is this possible with non-parametric statistics? Thank you


title: Why are there so many random generators in R?

text: A quick look at the help page of the random number generator of R shows that you can choose among 7 predefined generators (`Wichmann-Hill`, `Marsaglia- Multicarry`, `Super-Duper`, `Mersenne-Twister`, `Knuth-TAOCP-2002`, `Knuth- TAOCP`, `L'Ecuyer-CMRG`).               ?Random      The default is Mersenne-Twister, which seems to be very good. So why would you ever need to use another one?


title: What books to read for MCMC theory?

text: Suppose one is interested in stochastic processes for the purpose getting a theoretical understanding of MCMC. They already have a decent understanding of probability theory (let's say at the level Billingsley) and want to develop enough tools to potentially prove theorems about the Markov Chains they construct at a research level. What book or sequence of books would be reasonable for getting up to speed on the theory?


title: Checking assumptions for random effects in nested mixed-effects models in R / S-Plus

text: I am modelling a nested mixed-effects model with just the intercept in the random part, of the form:               fit4<-lme(fixed = Stdquad~factor(LayoutN)+factor(nCHIPS.fixed), random = ~1|Class.Ordered/student)      When trying to check the assumption on the independent, identically distributed (iid) random effects in SPlus using `ranef()`, I keep getting the error:               Problem in sort.list: ordering not defined for mode "list": sort.list(x, partial)      I suspect the nesting is the problem since if I remove _student_ , the plots come up fine. I would like to know how I can check that the random effects are iid and are independent for different groups using S-Plus (or R since its very similar).


title: Model selection using multivariate normal as input and multivariate negative binomial as response?

text: I'm trying to determine a good model to use to predict multivariate count data given a row of multivariate normal as inputs.   * The training set is $N\times D$ and the response set is $N\times P$, where $N$ is the number of observations and {$D,P$} are the number of features in the training and response respectively.   * Each number in the training set are continuous, and Box-Cox was used for each column to transform to normality with reasonably good fit.   * The columns in the response are count variables for different labels (# apples, # oranges, etc.).    * From training data it seems the counts for each column follow a negative binomial distribution very well, but with different parameters for each column.    * They are correlated, so that if you have a lot of apples, for example, there's a good chance you also have lots of oranges, etc. Any good choices for model selection in this case? I tried predicting each response column individually, but it missed some obvious relationships on a whole; for example predicting 1,000 Apples and 0 everything else.


title: Large Sample Theory

text: In large sample theory, I'm told that as $n$ grows larger and larger ($n$ being the number of samples in a dataset) that $\sqrt n(\hat \beta_1-\beta_1)$ gets closer and closer to normal distribution. What exactly is $\sqrt n(\hat \beta_1-\beta_1)$?


title: Gaussian "messages" for Belief Propagation

text: how can I find the marginalized "messages" $m(x_0)$... $m(x_n)$ ($x_1$ as example) $$m(x_1) = \int{ P(x_0|x_1,x_2,...) P(x_0)P(x_2)... dx_0 dx_1}$$ where $P(x_0|x_1,x_2,...) \sim N(m_0+x1+x2;v_0)$ and $P(x_1) \sim N(m_1;v_1)$, $P(x_2) \sim N(m_2;v_2)$, ... ? So I have the product of a normal distribution for $x_0$ whose mean is the sum of several other normal distributed variables $x_1$, ..., $x_n$ plus $m_0$. What I need to find is all the marginal distributions for the different variables ($x_0$ was easy...). I want to use this for the belief propagation algorithm with gaussians instead of discrete random variables, so I would guess that somebody has seen a solution before ;-). (all still univariant so far) Thanks for your help!


title: what if the residual series from arima() not looks normal?

text: My time series is here and my code is here. I fit a model to a series with `arima()` in R with ARMA(5,5) and regression on some covariates.               fit5 = arima(x, order=c(5,0,5), xreg=covaraites, include.mean=F)      I am now checking if the fitted model is adequate. The residual series looks like: ![enter image description here](http://i.stack.imgur.com/FcNpS.png) The residual series pass the Ljung Box test with a p value 0.3859. The ACF and PACF of the residual series are as follows (seems it is uncorrelated, right?): ![enter image description here](http://i.stack.imgur.com/ANqWq.png) ![enter image description here](http://i.stack.imgur.com/qcTPR.png) The qqplot for the residual series is as follows:               qqnorm(fit5$residuals, asp = 1)         qqline(fit5$residuals, asp=1)      ![enter image description here](http://i.stack.imgur.com/zYC3k.png) It looks okay within (-2, 2). I wonder if it is not okay to be Gaussian? IIC, arma model requires its residual series to be a white noise but not necessarily Gaussian. But `arima()` fits the model using MLE (assuming Gaussian residual series?). So if my residual series can't be taken as Gaussian, how shall I revise my model and what function in R I can use to fit to my time series?


title: How to perform regression calibration with only one variable in Stata?

text: I am attempting to perform regression calibration on some test data to familiarize myself with the idea behind it. I don't have much of a stats background. The method rcal in Stata seems to be what I am looking for, however I am having trouble when trying to regress on only one variable. In the following example, there are two independent variables (x1,x2) with an error in x2 in the rcal section. Ordinary regression with no measurement error ` gen x1 = uniform() gen x2 = uniform() gen y = 1+2*x1+3*x2 regress y x1 x2 ` This outputs the coefficients of x1 and x2, 2 and 3 respectively, with _cons=1 Now if I was to add some error on x2, and perform a regression calibration ` gen a1 = x2 + .3*invnorm(uniform()) matrix suu = 0.1 rcal(y=x1) (w1:a1), suuinit(u) ` This outputs an estimate of the coefficients of x1, and w1. Thats good. However, if I only have one variable ` gen x1 = uniform() gen y = 1+2*x1 regress y x1 ` This outputs the coefficient of x1. I can however not use the rcal method now if I was to add noise on the x1 variable. The rcal method seems to require at least one variable without error for its calculations. e.g., the following is not syntactically correct ` rcal(y) (w1:a1), suuinit(u) ` I can't seem to figure out how to only use only one independent variable. I'm thinking there might be a mistake in my understanding of regression calibration.


title: How to compare multiple proportions between multiple groups

text: I measured prey capturing success (=proportion of successful captures divided by the total number of prey captured) in four fish species (species 1: n=19, species 2: n=18, species 3: n=4, species 4: n=3) and would like to compare the prey capturing success between the four species. However, as there might be an individual effect within species (some fish might be better due to slight differences in morphology) I do not think one can simple sum everything up and compare four values (proportions). Therefore I want to know if it is possible to compare multiple proportions between the four species. Note that the proportions do not show a normal distribution within the species. Many thanks in advance!


title: How to update weights in RBM (Restricted Boltzmann Machines)?

text: Related Question: Learning weights in a Boltzmann Machine I'm trying to understand RBMs and how they are applied in training of Deep Architecture. Being new to the field of statistics, I stumbled upon various notations. Question 1: What does this symbol represent: <$s_is_j$> It is said while updating weights of RBMs that, the binary states of each node in hidden unit is set to 1 with a probability **p** equal to a function (a sigmoid function with the input: the sum of product of visible vector and corresponding weights vector connecting a hidden node). Question 2: Why is every hidden unit being set to 1 (with probability **p** ) and what is the function of the probability **p**? Question 3: As far as I understand weights are updated once the system goes to equilibrium, so when is the system said to be in equilibrium? Also is there any step-by-step procedures for training a network using RBMs


title: Prior predictive with discrete prior

text: I'm working with a uniform distribution as a prior, defined as: $\pi(\theta) = \begin{cases} \frac{1}{7} & \text{if } \theta\in\{0,\frac{1}{6},\frac{2}{6},\ldots,1\} \\ 0 & \text{otherwise } \end{cases}$ I've to compute the prior predictive, but I don't know how to deal with this discrete distribution. With a continuous one, I'd compute as, $P(\tilde{y})=\int_{\Omega}P(y|\theta)\pi(\theta)d\theta$ I know how to deal with this integral with a continuous distribution, but not with a uniform. As a model I'm using a $P(y|\theta)=Bin(20,\theta)$ I guess I've to do something with a sum $\sum$, but I'm not proficient in maths. Thank you in advance! Kenneth.


title: Difference between iid data and non-iid data for a simple regression problem

text: I am trying to understand the difference between iid and non-iid data. Let's consider a given time series, and say it's reasonable to assume that at each time point the random variable $X_t$ depends on $X_{t-1}$. Now say someone gives me the dataset $ D = \{ ( t_i, x_i ) \} $ for $i=1 \dots n$. I might be tempted to apply some regression method, which usually assumes that the corresponding pairs of random variables $(T_i, X_i)$ are iid. Would that be a bad idea if we really believe that $X_{t}$ depends on $X_{t-1}$?


title: Distance between a transition matrix and an instance

text: I am trying to put a number to the distance of a sequence and how close it is to the original training corpus. From the original training data, I got a markov transition matrix (TM). So from the sequence I am trying to evaluate, I have all the transition probabilities. I could calculate a new TM from the generated sequence (though much sparser), and for each element calculate the euclidian distance. I can calulate this in -log probs perhaps instead of probabilities, to be able to add them. Would there be another approach to evalute how much the new sequences looks like the training data/model?


title: Two stage model with a randomly transformed independent variable

text: Hi all and thank you in advance for helping! I have three variables: $y$ the dependent count variable, $x$ the independent positive continuous variable, $U$ is the observation unit (and I expect a random effect with it). I would like to test the hypothesis that the slope for $x$ is positive. I would normally use the following model: $y = \beta_0 + \beta_1 x + \phi + \epsilon$, where   * $\phi$ is the observation unit-level random error and   * $\epsilon$ is the observation-level random error. But I have a complication: at each observation unit, the original $X$ is transformed in some unknown order-preserving way, but not necessarily the same way at each observation unit. For example:   * At $U=$"A", $x=X/2$,   * At $U=$"B", $x=X^2$,   * At $U=$"C", $x=3X+3$ How would I test my hypothesis that the slope for $x$ is positive? Is it feasible?


title: Document Similarity Gensim

text: I am trying to get related documents for a list of 10,000 documents from the same set of 10,000 docs. I am using two algorithms for testing: gensim lsi and gensim similarity. Both give terrible results. How can I improve it?               from gensim import corpora, models, similarities     from nltk.corpus import stopwords     import re          def cleanword(word):         return re.sub(r'\W+', '', word).strip()          def create_corpus(documents):              # remove common words and tokenize         stoplist = stopwords.words('english')         stoplist.append('')         texts = [[cleanword(word) for word in document.lower().split() if cleanword(word) not in stoplist]                  for document in documents]              # remove words that appear only once         all_tokens = sum(texts, [])         tokens_once = set(word for word in set(all_tokens) if all_tokens.count(word) == 1)              texts = [[word for word in text if word not in tokens_once] for text in texts]              dictionary = corpora.Dictionary(texts)         corp = [dictionary.doc2bow(text) for text in texts]          def create_lsi(documents):              corp = create_corpus(documents)         # extract 400 LSI topics; use the default one-pass algorithm         lsi = models.lsimodel.LsiModel(corpus=corp, id2word=dictionary, num_topics=400)         # print the most contributing words (both positively and negatively) for each of the first ten topics         lsi.print_topics(10)          def create_sim_index(documents):         corp = create_corpus(documents)         index = similarities.Similarity('/tmp/tst', corp, num_features=12)         return index


title: diagnosing bias, variance from learning curve

text: I was doing machine learning course at coursera and there was a lecture about diagnosing high bias, or high variance from learning curves. If anyone interested here is the lecture - https://class.coursera.org/ml-003/lecture/64. I just can't find it enywhere else, I mean in some serious literature. So my question is, if this is really good, or commonly used approach to diagnose bias - variance in **classification** task? Or are there any commonly used, probably preferable way how to determine bias, variance error part in **classification** task? Any comment, point to useful resource would be very appreciated.


title: product of normal and lognormal variates

text: If x and y are uncorrelated normal variates, x*exp(y) will have a symmetric unimodal distribution with positive excess kurtosis. Has this distribution been named and studied?


title: lognormal distribution, standard-deviation and (physical) units

text: I have a distribution of microparticles that follows a lognormal distribution. The cumulative distribution function thus is given by: $$ F_X(x;\mu,\sigma) = \frac12 \operatorname{erfc}\!\left(-\frac{\ln x - \mu}{\sigma\sqrt{2}}\right) $$ $$ \mu = ln(M) + \sigma^2 $$ Now, the plot of the distribution function should be exactly the same no matter if the particle diameter $x$ is given in micrometers or meters (as long as I adapt the x-axis accordingly of course). However, this only works if I only convert $x$ and $M$, while not touching the numerical value of $\sigma$, and I don't understand why. $F_X$ has to be unitless, so $x$, $\mu$ and $\sigma$ should all have the same unit, right?


title: I have a source of numbers. Is there a way to estimate the probability for the given number to originate from it?

text: The source in question is some blackbox which output looks completely random to me, yey I haven't tested if for randomness quality. I don't exactly know insides of blackbox but can easily caputure it output 160bit integers, millions if necessary.


title: Goodness of fit for a Poisson random effect regression model

text: I have run a classical Poisson regression model and a Poisson random effect model in R. I want to know if there is any goodness of fit for the Poisson random effect model?


title: Several items for one construct

text: I am a dummy in statistics but need to accomplish some simple survey analyses. My question is: if I have several questions (items) for one construct (say, questions JS1-JS3 for job satisfaction) with the same 5-point scale, how do I do descriptive stats for the construct? Is it ok to calculate the mean & standard deviation across all questions for the items? (imagine I have 20 respondents, then the table of answers is 20 x 3 columns of questions JS1-JS3 - do I select the whole area for excel functions?)   Do I use Std Deviation or Std Deviation for population function in excel? How do I get cronbach alpha for this set? Also, if I have several constructs within the survey that are measured with multiple items, how do I perform correlations analysis? Please, help! Thanks a lot, experts!


title: How to calculate the probability that the arithmetic mean of a multinomial process exceeds some value

text: I've been mulling over a problem that has something like the following form. I don't have a math or stats background so advice and answers at various levels, from terminological to strategic, would be helpful. I'll start off with the example that I've been thinking about and then move on to the specific case I'm working towards. Suppose you have a jar full of coins, say pennies, nickels, dimes, and quarters. They are worth, respectively, $1$, $5$, $10$, and $25$ cents. You're holding a particular coin, in this case let's say a dime. You draw a coin from the jar, write down its value, and put it back in the jar. You do this some number of times and take the average (arithmetic mean) of the values. What is the probability that the average of the values of the coins you've drawn is strictly greater than the coin you're holding (a dime)? What about for the case when you're holding any other coin? For the simplest case of drawing a single coin, it's easy enough. If you're holding a penny, it's the probability that any single coin is greater in value than a penny. So, it's the probability of drawing a nickel, dime, or quarter. If you're holding a nickel, it's the probability of drawing a dime or a quarter. If you're holding a dime, it's the probability of drawing a quarter. If you're holding a quarter, there's no coin that you can draw that is greater in value than the one you're holding. For the case of two coins, it's also easy enough. If you're holding a penny, it's the probability of drawing at least one coin that isn't a penny. If you're holding a nickel, it's the probability of drawing at least one dime or quarter. If you're holding a dime, it's the probability of drawing at least one quarter. If you're holding a quarter, then no coins will yield a higher value. As we increase the number of coins drawn, however, the number of different combinations increases exponentially. Let's call the number of coins in the jar $C$, the number of coins you draw from the jar $n$, and the probability of drawing a penny $p_P$, nickel $p_N$, dime $p_D$, or quarter $p_Q$. Where $p_i = \frac{\text{number of i coins}}{C}$. Is there a simple, but general way to express the probability that the average of your draws will be greater than the coin you're holding? Now, instead of monetary values, let's say the coins are worth $-b$, $0$, $a - b$, and $a$. Where $a > b$. Again, is there a nice way of expressing this?


title: Piecewise linear regression with knots as parameters

text: I would like to fit a piecewise linear regression with knots as parameters. I would like to know what's the best solution.   * Should I run a set of regressions with all the possible knots and choosing the knots which could minimize an information criterion such as AIC (Akaike Information Criterion)?    * If that's the best solution, how can I compute standard errors for my estimates?


title: Nested cross validation for model selection

text: How can one use nested cross validation for **model selection**? From what I read online, nested CV works as follows:   * There is the inner CV loop, where we may conduct a grid search (e.g. running K-fold for every available model, e.g. combination of hyperparameters/features)   * There is the outer CV loop, where we measure the performance of the model that won in the inner fold, on a separate external fold.  At the end of this process we end up with $K$ models ($K$ being the number of folds in the outer loop). These models are the ones that **won** in the grid search within the inner CV, and they are likely different (e.g. SVMs with different kernels, trained with possibly different features, depending on the grid search). How do I choose a model from this output? It looks to me that selecting the best model out of those $K$ winning models would not be a fair comparison since each model was trained and tested on different parts of the dataset. So how can I use nested CV for model selection? Also I have read threads discussing how nested model selection is useful for analyzing the learning procedure. What types of analysis /checks can I do with the scores that I get from the outer K folds?


title: Using correlate in Stata to compare within a variable

text: I am trying to compare correlations between subgroups of a variable. For instance, I am trying to compare the income of group 1 v. income of group 2 v. income of group 3 v. income of group 4, where the groups are specified by another categorical variable. What I did was create new variables separating out the income variable into four new variables (i.e. income_1, income_2, income_3, income_4). Naturally, there are missing values that look something like this (each bracketed portion represents a new line where the '.' represents missing values):               [income_1   income_2    income_3   income_4]       [1 . . . ]       [1 . . .]       [. 2 . .]       [. 2 . .]       [. . 1 .]       [. . 1 .]       [. . . 2]       [. . . 2]        So when I try to perform `correlate` on these four values, I get an error that says 'no observations' probably because `correlate` performs listwise or pairwise deletion. How can I fix this problem? I don't know how to get rid of the missing values.


title: Is the linear probability model (LPM) less inappropriate when the independent variable of interest is binary? (and controls are continuous)

text: I need to use an LPM for a paper I'm doing (for reasons orthogonal to this question). I have about 6000 observations in a time series cross-sectional structure (countries and years) and my dependent variable is binary. Now, I know all the standard arguments against the LPM, but I am wondering whether the case for using it is better when the independent variable of interest is binary? As far as I can tell, if all independent variables are binary the LPM is completely general, but I have some controls that are continuous. So, is the LPM equally bad if the independent variable of interest is binary?


title: SVM with sparse Gram matrix

text: I have a corpus of many short entries. Because of this it is very rare that for two entries $x_i, x_j: \ <x_i,x_j> \neq 0$. Therefore for almost all the entries, $x\in TeS$, in the test set I get $\forall x_i\in TrS, \ <x,x_i>=0 $ where $TrS$ is the training set. Therefore SVM is unable to classify most of the entries in the test set. Is SVM the wrong approach to take with short entries? What other algorithms do not face this problem?


title: Use of log-linear analysis with percentage data

text: Although my question sounds like I have data suitable for log-linear analysis, I may not. I have a large number of students (about 1000) of two levels who solved a test of 8 different types of problems (which have specific features and categorized at two different levels). Thus I have three types of independent variables: `type of problem`, some other `features of problems`, and `student level`. I want to see if there are main and interaction effects along these dimensions. Since the dependent variable is dichotomous, I need a nonparametric test. Also it should be a repeated design, since same students answered 8 different questions. So the first question is: **Is there such a test?** What I did was calculate the percentage of correct solutions for each case and ran a log-linear analysis, which I am not quite sure about (as it violates the independence assumption).


title: Why is ordering values necessary for nonparametric regression

text: Suppose we are want to consider the asset price $P_t$ (business daily) of some stock. The log return is defined as $$X_t:=\log{\frac{P_t}{P_{t-1}}}$$ suppose we are considering the prices between some interval, eg. June 1986 and March 1990). Log returns can be modelled by $X_t=\sigma_t\epsilon_t$, with $E[\epsilon_t]=0,Var(\epsilon_t)=1$. We assume Markov-property, thus $\sigma^2_t=v(X_{t-1})$. The model can be fitted by nonparametric regression of the function $v$ in $$Z_t:=X^2_t=v(X_{t-1})+\phi_t$$ where $\phi_t=\sigma^2_t(\epsilon^2_t-1)$ can be seen as error. We want to use `ksmooth`, `smooth.spline` and `loess`. In my notes they say: "First we sort the values, in order not to get problems with `ksmooth`, which orders the values internally and gives back results corresponding to the ordered values." I do not understand why we have to do that. Basically we have a function $v$, with values on the $y$-axis and the $x$-axis is our time interval between June 1986 and March 1990. So why do we have to order the data?


title: RESET Test and Heteroscedasticity Test for Hierarchical Multilevel Model

text: Are there implemented tests in R to check heteroscedasticity and linearity for multilevel regressions? I am using a two Level Regression Framework (lmer function in R). Is it necessary to do those tests if i use a MCMC Approach (MCMCglmm in R)? How do i perform those test when using MCMC to estimate multilevel regressions (gaussian distribution)? Best regards, Linus holtermann@hwwi.org


title: Hausman test: Why requirement of positive definite covariance matrix?

text: Can someone please explain to me why the covariance matrix needs to be positive definite for the Hausman test? And why we can skip this restriction by applying the "test of overidentifying restrictions"? Isn't the issue easily explainable when simply comparing the formula of the Hausman test statistic and the formula of the Sargan test statistic?


title: Guttman's smallest space analysis

text: I am trying to replicate an analysis as conducted by Raven, et al. (1998). A scale is analyzed with Guttman's smallest space analysis. First of all I wasn't able to find a suitable procedure in R and it appeared to me as if SPSS might be the solution. However even though I use the same scale it is impossible to reproduce the graph as printed in the original publication: ![SSA graph as published by Raven, et al. \(1998\)](http://i.stack.imgur.com/dgjNj.png) The graph I manage to get out of SPSS looks like the one below: ![SPSS SSA graph](http://i.stack.imgur.com/Z2BJx.png) I suppose there must be some settings to adjust the dimension or scales. Any recommendation to an R package or some guidance to the SPSS settings would be appreciated. * * * ## References Raven, B. H.; Schwarzwald, J. & Koslowsky, M., Conceptualizing and Measuring a Power/Interaction Model of Interpersonal Influence, _Journal of Applied Social Psychology, Blackwell Publishing Ltd,_ **1998** , _28_ , 307-332 [Link]


title: Geographic regression

text: I'm working on a project to estimate real estate and started with some classique techniques, such as linear regression etc. The obtained results are already going in the good direction, but to get some more precise results, I've started to lookup some more advanced models. First I've looked into a spatial autoregressive model, but these models are a bit hard to use if we want to predict new values, or am I wrong? The next one is the kriging technique. In general, these models use only two coordinates and one prediction value (like temperature or height). In my project I'm using way more variables to predict the outcome. My question is which technique/model I would best look into to obtain some better results.


title: Wald test and F distribution

text: Assume that the sample $X_i, i=1,\dots, n$, is iid normally distributed with unknown mean and variance $\theta = (\mu, \sigma^2)$. The Wald test statistic of testing if $\mu = \mu_0$ is $$W=\frac{n(\hat{\mu} - \mu_0)^2}{\hat{\sigma}^2}$$ where $\hat{\sigma}^2 = \frac{\sum_i (X_i - \hat{\mu})^2}{n}$ is the maximum-likelihood estimator of $\sigma^2$.   1. I was wondering if it is correct that $\frac{(n-1)W}{n^2}$ has a F distribution with $1$ and $n-1$ d.f. under null?   2. From a note, it seems that $W$ is an F distribution with $1$ and $n-1$ d.f. under null. Isn't this incorrect? Let $p_2$ be the dimension of $\mu$, and $p$ be dimension of $X_i$. > Under the assumption of normality we have a stronger result. The > distribution of $W$ is exactly chi-squared with $p_2$ degrees of freedom if > $σ^2$ is known. In the more general case where $σ^2$ is estimated using a > residual sum of squares based on $n-p$ d.f., the distribution of $W/p_2$ is > an $F$ with $p_2$ and $n-p$ d.f. Thanks and regards!


title: Problem regarding the scale and range of factor variables generated in CFA (Stata)

text: I ran a confirmatory factor analysis in Stata using the `sem` command. The idea was to extract two factors: TRUST, extracted from four variables measuring the trust towards different groups of people (friends, strangers etc.) on a scale from 1 to 7; SELF-RELIANCE, extracted from four variables measuring the consent of the test-subject regarding statements declaring socio-economic self-reliance as a virtue in our society on a scale from 1 to 4. I intend to use SELF-RELIANCE as a dependent variable in a multiple linear regression. TRUST will be one of the predictors in said regression. I get pretty decent factor scores, have a decent model fit, everything seems fine. But then I try to generate the factor variables TRUST and SELF-RELIANCE from the PCA using               predict <varname>, latent(factorname)      I expected the factors to have the same scale as their endogenous variables (1-7; 1-4), but they have a range from -2.944997 to 1.157046 (in the case of TRUST). What does that even mean? Did I use the wrong predict option? My professor said the variables with the constrained factor scores in `sem` serve as the scale template for the factors, so I don't understand why my generated variables would have such a weird scale. Does anybody know what this means? Here's my syntax for the CFA (names are in German, as am I, so I'm sorry if I formulated something wrong):               sem (vertrverw vertrfreund vertrbek vertrfremd <- VERTRAUEN) (einkuntmot ranguntakz sozuntger <- EIGENV)     estat gof, stats(all)     estat mindices      predict VERTRAUEN, latent(VERTRAUEN)     predict EIGENV, latent(EIGENV)


title: "grouped error distribution"?

text: I can't figure out the meaning of emphasized phrase in the following: > Representing a _grouped error distribution_ by a multinomial distribution, > Simpson (1756), Langrange (1776), and Laplace (1776) used the classical > methods for deriving the distribution of the mean rectangularly and > triangularly distributed errors. I think I know what the author means by an "error distribution", but I can't figure out what the "grouped" refers to.