title: What is the relation between the estimated standard deviation of a normal distribution and the scale of a t distribution when applied to normal data?

text: If I generate normally distributed data and then use this data to estimate the parameters of a normal distribution using, say, maximum likelihood I would expect the estimated standard deviation to be close to the true standard deviation given enough data. If I instead use this normally distributed data to estimate the parameters of a t distribution with a fixed DF the scale does not approach the true sd. There should, however, exist a functional relationship between the SD of the data and the scale of the t-distribution. I haven’t manage to figure this relationship out so my question is: **What is the relation between the estimated standard deviation of a normal distribution and the scale of a t distribution when applied to the same (truly) normally distributed data?** I've done some numerical simulations, where I generated normally distributed data, fitted a normal and a t with a fixed DF and looked at the ratio between the estimated SD of the normal and the scale of the t. It definitely looks like there is a relation, but it would be nice to have a formula for this :) Here is the simulation in R code:               library(MASS)          scale_ratio <- sapply(seq(1, 10, 0.1), function(df) {       x <- rnorm(9999, sd=100)       fitdistr(x, "normal")$estimate[2] / fitdistr(x, "t", df=df)$estimate[2]     })          plot(seq(1, 10, 0.1), scale_ratio, xlab="DF", ylab="normal SD / t scale")      ![scale vs sd](http://i.stack.imgur.com/iMj53.png) _Edit: 1_ A possible use for the transform I'm looking for would be when using a t-distribution as a robust alternative to the normal distribution when modelling data. The t distribution is a consistent estimator of the mean of normally distributed data with the added benefit that it will be more robust against contamination/outliers in the data. As shown above, the t-distribution will not be a consistent estimator of the SD but it would be neat if it could be made consistent by transforming the scale somehow. _Edit 2:_ One way of transforming the scale of the t distribution so that it becomes more similar to the SD of the normal data is by taking the Interquartile range (IQR) of the estimated t distribution and then dividing the IQR by 1.349 which, _if_ the IRQ would have come from a normal distribution, would have resulted in the SD of that normal. This transform feels a bit like a hack to me and does not make the discrepancy go away completely as seen below. ![enter image description here](http://i.stack.imgur.com/0U8p0.png) Another possibility would be to use the SD of the t-distribution, instead of the scale parameter, as an estimumate of the SD of the normal data. This does not work well, however, but instead overestimates the SD of the normal data as seen below: ![enter image description here](http://i.stack.imgur.com/QNQm2.png)


title: Canonical correlation analysis

text: I have to do canonical correlation analysis between two multivariate datasets X and Y. One dataset contain numerical data and the other binary data.I would like to know what features are highly correlated with the features in second data set. Would the normal CCA available in Matlab be sufficient if one data set has numerical(integer ans floating point values) and the second dataset has binary vlaues(0 or 1) for all the valuese for the various features?


title: MCMCglmm for dummies

text: In R, there is the function MCMCglmm (in the package "MCMCglmm"). It allows one to perform a Markov chain Monte Carlo Sampler for Multivariate Generalised Linear Mixed Models. My question is: **How does`MCMCglmm` work (explained to dummies)?** There is no need to go into the details (neither talking about programming) but just the very basic understanding about the statistical process behind this function. * * * Note: I had a look to `vignette("Overview", "MCMCglmm")` but the article is long and quite hard for me. I have got few knowledge in Bayesian statistics and some experiences in Approximate Bayesian Computation (ABC).


What does "chi" mean and come from in "chi-squared distribution"? Does "chi" mean a random variable with a standard normal distribution? Thanks!


title: Transforming an estimator to overestimate

text: I have an estimator $\mu^*$ of a the mean $\mu$ of a certain distribution that I obtained using a variational technique (basically just establishing a bound on $\mu$ and finding a trial function that minimizes that bound). My estimator $\mu^*$ has therefore the property that: $\mu \leq \mu^*$ I want to compare the performance of my estimator to estimators obtained by other means (e.g. through MCMC) if possible by using the same variational bound. To do this, I want a transformation that converts any estimator into an overestimator. Specifically, I want some sort of transformation $f$ such that for any estimator $\mu'$ of $\mu$, we have that $\mu \leq f(\mu')$ and that for a loss function $L$ with respect to $\mu$ (e.g. least squares) $f$ has the property that $E[L(\mu', \mu)] = E[L(f(\mu'), \mu)]$. My question is, does this transformation $f$ even exist or if not is there an approximation such that $E[L(f(\mu'), \mu)]$ lies relatively close to $E[L(\mu', \mu)]$? Rao-Blackwellization may help here, but finding a sufficient statistic may be hard in my case and I don't know a priori if the estimator $\mu'$ is unbiased.


title: How to convert percentile rank into z score in R

text: I have a vector of percentile ranks. I want to convert them into z-scores, so it will be interval scale. I have to do it in R, but I could not find a function or package that can do this. Does anyone have any idea?


title: Is Latin hypercube sampling effective in multiple dimensions?

text: I am currently using a Latin Hypercube Sampling (LHS) to generate well-spaced uniform random numbers for Monte Carlo procedures. Although the variance reduction that I obtain from LHS is excellent for 1 dimension, it does not seem to be effective in 2 or more dimensions. Seeing how LHS is a well-known variance reduction technique, I am wondering whether I may be misinterpreting the algorithm or misusing it in some way. In particular, the LHS algorithm that I use to generate $N$ spaced out uniform random variables in $D$ dimensions is:   * For each dimension $D$, generate a set of $N$ uniformly distributed random numbers $\\{u^1_D,u^2_D...u^N_D\\}$ such that $u^1_D \in [0,\frac{1}{N+1}]$, $u^2_D \in [\frac{1}{N+1}, \frac{2}{N+1}]$ ... $u^N_D \in [\frac{N}{N+1}, 1]$   * For each dimension $D \geq 2$, randomly reorder the elements from each set. The first $U(0,1)^D$ produced by LHS is the a $D$ dimensional vector containing the first element from each reordered set, the second $U(0,1)^D$ produced by LHS is the a $D$ dimensional vector containing the second element from each reordered set, and so on... I have included some plots below to illustrate the variance reduction I get in $D = 1 $ and $D = 2$ for a Monte Carlo procedure. In this case, the problem involves estimating the expected value of a cost function $E[c(x)]$ where $c(x) = \phi(x)$, and $x$ is a $D$-dimensional random variable distributed between $[-5,5]$. In particular, the plots show the mean and the standard deviation of 100 sample mean estimates of $E[c(x)]$ for sample sizes of 1000 to 10000. ![LHS for $D=1$](http://i.stack.imgur.com/Ti9Qd.jpg) ![LHS for $D=2$](http://i.stack.imgur.com/aO7Cj.jpg) I get the same type of variance reduction results regardless of whether I use my own implementation or the `lhsdesign` function in MATLAB. Also, the variance reduction does not not change if I permute all sets of random numbers instead of just the ones corresponding to $D \geq 2$. The results make sense since stratified sampling in $D = 2$ means that we should sample from $N^2$ squares instead of $N$ squares that are guaranteed to be well spread.


Consider two binary observations, $x_1$ and $x_2$, which are independent given $\theta$. Now suppose there is uncertainty with respect to $\theta$. Find the correlation between two binary observations when the parameter is one of two possible values with priors $p(\theta_1)$ and $p(\theta_2)$.


title: Visualizing an interaction effect of a Logistic Regression Analysis (in SPSS)

text: What would be the best way to go about visualizing an interaction effect of a Logistic Regression Analysis (preferably using SPSS)? There is a significant interaction between one interval (HDI) and one binary (before/after Fukushima) variable as well as between that same interval variable (HDI) and another interval variable (share of non-renewable electricity consumption). The dependent variable is binary (ecological benefits of solar power mentioned in news article about topic yes/no). Failing any kind of more specific way of doing this, what would be more general steps to arrive at some visualization? Is a visualization even a good idea? (For context: Until now I have been looking at interactions using ANOVA where plotting interaction effects is straightforward. Those visualizations have been tremendously helpful and I would like something similar for this.)


title: How to analyze data from a 2x3 experiment where one cell is missing?

text: I have an experiment with two factors. The first factor, picture, has three levels: nice, ugly, absent (i.e., the picture is not presented). The second factor, information, has two levels: present, absent. The dependent variable is a continuous variable: liking. Now, I did not run the picture=absent + information=absent condition because it would mean not presenting anything and asking people how much they liked it, which doesn't make a lot of sense. So I'm left with 5 cells instead of 6. Any ideas on how to analyse these data?


title: Should my test set be balanced or imbalanced?

text: I have an imbalanced dataset (90% class 0 10% class 1), should I first split it in training and test set, then balance my training set (my test set would still be imbalanced), or could I randomly downsample the majority class in my dataset then split in training and test set (my test set will not contain any observations used in the training set because I down sampled)? I am getting very different results for each approach.


title: In cluster analysis, can you use Gower's coefficient of similarity with a k-means clustering method?

text: I am researching cluster analysis, and I am interested in variables that are both categorical and continuous, for which I have read that a Gower's similarity coefficient is a good proximity measure. I have read that Gower's similarity coefficient is generally not compatible with Ward's method, so I was planning to initially cluster using average linkage, but I was also seeking to compare the cluster structure (for content validity purposes) with another clustering method, specifically the k-means method, using the number of clusters and initial centers obtained in the average linkage method. Is Gower's coefficient of similarity a compatible proximity measure for k-means method?


title: Pearson correlation test for non-gaussian data

text: I want to perform significance test for correlation coefficient. However my data is not normally distributed. In my books I have found tests only for the normally distributed data. Is there such test I could apply to data having any distribution?


title: Looking for examples or alternatives to R RuleFit ensemble package

text: Does anyone know of any good example code illustrations for the rulefit Rule Based Learning Ensembles package? The documentation is incredibly lacking. I was guided to the package by this paper. If anyone who is familiar with the model could give a brief interpretation of how the Rule based method works (and why it works well or not), it would also be appreciated. The results in the paper look very promising. Also, are there any other alternative ensemble packages that perform as well or have similar functionality? I'm familiar with Caret, although, I don't know if any of the Caret methods are necessarily comparable to the above. For instance, Caret uses linear correlation for the variable importance selection, whereas the Rule Based Method seems to use interaction statistics and partial dependence functions and plots to determine importance.


title: Correlation of non-uniform time-series data?

text: I have two datasets of time-series data, each with arbitrary, non-linear sampling intervals. I'm planning to analyze them using either Pearson or Spearman correlation (determined by the detection of outliers using the 3-sigma rule) but I believe that the correlation between times at which each sample was taken may be relevant. Does anyone know of a good algorithm to correlate this type of time-series data, taking into account sample time if it's relevant? Pseudo-code is particularly welcome... A million thanks, in advance!


title: Suggestion for MATLAB statistical pattern recognition toolbox?

text: I am looking for a pattern recognition toolbox in MATLAB capable of doing well and have enough number of possible routines. What can be good candidate?


title: How do I find fields related to a class in a naive bayesian classifier?

text: I've built a naive bayesian classifier which classifies a binary category (either in or out). I have a number of variables which the classifier gets trained on v1, v2, v3, v4, v5, etc. How do I identify which variables contribute more or less to the overall classification? EDIT: I'm still working on this problem, but I think I've found an appropriate answer. Information gain ratio. Wikipedia says: "Information gain ratio is between 0 and 1, cf. ordinary information gain". [1] Which seems like exactly what I want. I can calculate IGR, and then use that to weight each feature. I also get a very intuitive list of which features are related to my target class, which I can use to make decisions in the real world. [1] https://en.wikipedia.org/wiki/Information_gain_ratio#Advantages


title: Comparing two mean values to check for statistical significance

text: This is a very simple question but it has been years since I have done statistics. I am also new to this website (was recommended by a friend) so apologies if this question has been done 'wrongly'. I am comparing several criteria between two different treatment methods of a certain disease. I have the following table as an example of results.                    Treatment A                 Treatment B               14                          23               8.2                         8.1               17.4                        17.1               0                           2      The first row is the amount of patients I have on each treatment. The next 3 rows are the average results of different criteria. Basically, how do I check if the difference in results of each criteria are statistically significant? Currently, from searching around on the interwebs I believe I am supposed to do a T-test. I have been using this website: http://www.quantitativeskills.com/sisa/statistics/t-test.htm I input the mean result of treatment A and B, the amount of each data, the standard deviation (calculated using excel with the =stdev() function) and with a width of C.I. set to 95%. The website then gives me a lot of jargon which I do not understand. However, for some of my results it has given the conclusion that 'Difference not significant at 5%' and some of my results do not show that line, which signifies to me that the results are indeed significant. Have I gone about this the right way?


title: Statistical comparisons of multiple classifiers performance?

text: If the accuracy of $classifier1$ is statistically significantly better than $classifier2$ as per some hypothesis test, and likewise the accuracy of $classifier2$ is statistically significantly better than $classifier3$, does it follow that $classifier1$ is statistically significantly better than $classifier3$. This question relates to my question asked here.


title: How to compute a probability threshold for a linear combination of two variables ~ N(0,1)?

text: I have a variable which is a linear combination of two other variables, each one following an N(0,1) distribution. I need to compute the threshold of the distribution of this combination variable (to compute the Z for a specific probability). Does anyone know which function should I use, i.e. which package is available to perform this task?


title: Marginalisation on conditional probability

text: I'm afraid this is an extremely simple question, but I didn't understand completely why if both $O_1$ and $O_2$ can be marginalised over $R$ the following holds: $ P(O_2 = o_i | O_1 = o_j) = \sum_r P(O_2 = o_i | R = r) P(R = r | O_1 = o_j)$ Thanks


title: Statistical method for modelling of disease incidence for NGO

text: I am trying to provide help to an NGO in the non-profit sector that is running a disease screening program: The program visits thousands of villages a year. A village has a population (on average around 800), with an unknown number of infected individuals. Each village is visited by the NGO for mass screening. For a given population in a given year:   * An infected person will attend the NGO screening and get cured with a probability of S (e.g. 50%)   * An infected person might independently self-diagnose himself, seek treatment and get cured in a health clinic outside the screening program with probability P (e.g. 10%)    * An infected person will infect another person per year with probability R (e.g. 50%)    * An infected person can live for 4 years.    * There might a fixed probability that a new infection will jump over from another village (J) These probabilities (S, P, R) will differ per village (within a estimated prior distribution), but are constant per village over time. The NGO has data going back 7 years for each village. ![Example data for 2 villages](http://i.stack.imgur.com/22AJD.png) The NGO would like to predict the number of infected people in the population for each village in year 8. I am trying to help the NGO, but have struggled to apply statistical learning to this problem, due to the stochastic nature and path-dependency. Most of my approaches try to fit a max. likelihood estimate to a forward projection of a starting population variable in year 0 with sampling of parameters from the prior distributions. I understand we are limited by the many hidden parameters and very limited data per village, but I feel there must be at least some predictive value in the data, but struggle to find the right approach. Hence – I was wondering whether someone can give their perspectives on what is the best approach / are the best tools to perform predictive modelling on such a problem, and help me understand which might be possible in terms of predictive accuracy?   Any tips, links or advice much appreciated;


title: How to find the inflection point in logistic regression using R?

text: I have conducted a binary logistic regression (with 7 predictor variables) in R using "year" as a random effect, i.e. a generalized linear mixed model (GLMM). For a few of the predictor variables, I would now like to find the inflection point, i.e. the value of the predictor variable for which the probability of the response variable being "1" equals 0.5. I would like both to produce a curve of this and to calculate the exact value, so that I won't have to read it from the curve (with the inaccuracy this would cause). A post from November 2013 made by Francis and edited by Scortchi touches on this topic (with a different purpose, though) and shares a specific code for producing such a curve (Logistic Regression and Inflection Point). However, I don't understand the 4th line of the code (`newdata`) and can't make the code work for my data. Is the 4th line something that has to be done in order to produce the curve? Might someone help me by uploading a more general code, or gently translating the referred one, for me to understand exactly how to fit my own data into the code? Also, if someone knows how to make R calculate the exact value, I would very much like to know. As you might be able to tell, I'm still learning R. I hope that someone can help.


title: What can cause a "Error() model is singular error" in aov when fitting a repeated measures ANOVA?

text: I have the following data set:                                |       Scenario 1       |     Scenario 2         |                      |Trial 1|Trial 2| Trial 3|Trial 1|Trial 2| Trial 3|      -------------------------------------------------------------------                   S1 | ...      Condition 1  S2 | ...                   S3 |      -------------------------------------------------------------------                   S5 |      Condition 2  S6 |                   S7 |      Thus the `Trials` are nested in the `Scenarios` and all of them are within subject. I am trying to run an ANOVA on this data set. Here is the model without defining that `Scenarios` (and `Trials`) are within subject.                my_data.aov <- aov(value~Condition*Trial%in%Scenario,data=my_data) #works fine      But when I specify that these are within subject:               my_data.aov <- aov(value~Condition*Trial%in%Scenario+Error(Player/(Trial%in%Scenario)),data=my_data)       I get the following error               In aov(value ~ Condition * Trial % in % Scenario + Error(Player/(Trial %in%  :     Error() model is singular      The closest set-up I could find was Split plot in R but there the subjects are nested inside each `Trial` not in each `Condition`. **EXAMPLE FILE** Here is an example file in long format. **What about this approach?** If I treat each `Trial` as a sample, then I can collapse across `Scenarios` by averaging them, so I will have a simpler model, where each `Subject`'s behavior is described per `Scenario`. And since I need to analyze the relationship of `value~Condition*Scenario` I can do so by defining the `Error` like `Error(Subject/Scenario)`. Will this approach invalidate my analysis?


title: R regression with categorical response variable

text: I have four variables, two are categorical and two are numeric:               a<-c("yes","yes","no","no","no",NA,"yes","no")     b<-c("high","low","medium","medium","medium","low",NA,NA)     c<-c(12,23,23,12,23,34,12,NA)     d<-c(1.2,1.3,4.5,3.4,NA,5.4,9.4,7.4)          df<-data.frame(a,b,c,d)      I want to replace the `NA` values in each variable by using a model fit. For example:               a~b+c+d; b~a+c+d       and so on. I know that if all the variables are numeric, a simple multiple linear regression model can be built to estimate the value of the `NA` data. However, I am not sure how I can achieve this if the response variable is categorical.


title: Von Mises distribution to detect outliers

text: I am working out the difference between two angles from a circle, and I work out the mean difference across 96 trials in 10 separate samples. In order to detect outliers for statistical analysis, Barret & Lewis ( _Outliers in Statistical Data_ , 1984) suggest the use of a von Mises basic model.   1. Is a von Mises distribution is appropriate in my case? I'm not interested in the raw angle values _per se_ , but the difference between them.   2. I understand how outliers are calculated from _z_ -scores, standard deviations from the mean, etc., but I don't understand how how they are calculated from von Mises – can anybody offer simple clarification? Thanks.


title: Categorizing personality outcomes

text: I measured the five factor model (extraversion, neuroticism, openness to experience, and agreeableness) and the results gave me a continuous value between 1-5. Using these values as an IV is giving me problems as the frequency of the values (can be anything between 1-5) is way to low to produce any results. So, I believe the only option for me is to cluster them. Would clustering them using the quantiles be a good approach here? I can imagine that this will differ across the 5 factors of the personality and thereby also influencing the interpretation. What would here be the best approach to deal with this data.


title: Customer discovery as a mathematical problem

text: The last couple of days I tried to find systematic approaches to customer discovery, i. e. how to answer the question "what kind of people are interested in product X?". I didn't find anything useful (lots of tactics, but it's no strategy), so I started to think about it myself. Imagine, there is a population of people, each of which has following attributes   1. age (several age groups),   2. gender (male, female),   3. region,   4. income level,   5. education level. Each of these variables can take a fixed set of values. In order to find out whether people with a certain combination of attributes (e. g. 35-40 years old, male, US, X dollars annual income, bachelor degree) are interested in the product, I can run a test, which is associated with certain costs (think of showing an ad to 100 people from that target audience and counting how many of them reacted to that ad). I know for sure that the product I want to sell is bought by at least one type of people in that population (I know this because similar products are sold there already). I want to search for the perfect customer in as effective as possible. Effective means that the number of experiments is minimal. If I use brute-force (show the product to all possible types of people) approach, the number of experiments to run is equal to `number of age groups * 2 * number of regions * number of income levels * number of education levels` (where 2 is the number of genders). This is the worst-case option. So I need to find a way to find out the type of people (characterised by those 5 attributes), such that the number of experiments is minimal. By "find" I mean that the response to the ad is above some threshold. Probably I'm not the first person, who needs to solve a problem like this. AFAIR the backpropagation algorithm in neural networks searches for best possible values of weights, for example (like I'm searching for the best combination of attributes). Long time ago I read about robots, who learn how to find their way out of a maze using a series of "experiments" (the robot gets feedback from the outside world whether it has managed to go from the maze, in my case I get feedback from the response to the ads). I'd appreciate, if you told me what existing models or algorithms I can look at for inspiration (which solve similar problems).


title: How to specify and calculate conditional probabilities with trinomials?

text: Could someone give an example of a simple BN using trinomial nodes? I am a bit confused, particularly on how to specify and calculate conditional probabilities with trinomials.


title: How to test the randomness of residual plot?

text: Let's say in linear regression, I got a fit and I can plot residuals to see whether there is any systematic trend in such a plot. How to quantitatively determine whether the residues are really random? Is Durbin-Watson test used for this purpose? How to interpret such test if so? Please provide an example, preferably in R.


title: How do I calculate if the degree of overlap between two lists is significant?

text: If I have two lists A and B, both of which are subsets of a much larger list C, how can I determine if the degree of overlap of A and B is greater than I would expect by chance? Should I just randomly select elements from C of the same lengths as lists A and B and determine that random overlap, and do this many times to determine some kind or empirical p-value? Is there a better way to test this?


title: Fisher's exact test and large tables with pairwise

text: I found that it is possible to look at adjusted residuals and treat them as a z-score and then use a two-tailed test to find a p-value. The problem is how do you calculate adjusted residuals on what I realise is a 1xr table (as it is a goodness of fit test)? In addition if anyone is interested in seeing what I was doing in regards to my pairwise and more details on my issues then please look at a forum question here. I also have tables of habitat use vs. availability. Some of the values are small but only two groups have a mean value below 5. I use a $\chi^2$ as I did not know how to do a Fisher's on such a large table (9 or 5 habitats). The issue is, I then did pairwise comparisons for each habitat combo to add the p-values to come up with which habitats were contributing to the significant result. Should I somehow have used Fisher's for these 2x2 tests? I have not found a way to do these pairwise for chi-squared so I set up an Excel file to do it and then compared each p-value (smallest first) with Holm's Bonferroni p-values i.e. `0.05/9` then `0.05/8` etc., for each in turn. Does this sound feasible rather than delving into Fisher's within the test? The results found did look valid and ironically the habitat type with `O=4` and `E=2` was not significant. In addition, if I am doing the Holm's correction by hand (i.e., comparing each with the respective p-value) should the real p-value be shown in the results, or should you show a p-value that has been multiplied by the respective correction so that they can all look as if they are being compared to `alpha=0.05`? For example, my fifth smallest of 9 p values is `0.02`. According to Holm's I should compare this to `0.05/5`, i.e., `0.01` so it is not significant. If I need to report the p-values, should this one, for example, be shown as `0.10`, i.e., it has been multiplied by the 5? Or should I show the original `0.02` and then have a reference such as "significant using Holm's correction?


What is the difference between a population and a sample? What common variables and statistics are used for each one, and how do those relate to each other?


title: What is the distribution of the sum of non i.i.d. gaussian variates?

text: If $X$ is distributed $N(\mu_X, \sigma^2_X)$, $Y$ is distributed $N(\mu_Y, \sigma^2_Y)$ and $Z = X + Y$, I know that $Z$ is distributed $N(\mu_X + \mu_Y, \sigma^2_X + \sigma^2_Y)$ if X and Y are independent. But what would happen if X and Y were not independent, i.e. $(X, Y) \approx N\big( (\begin{smallmatrix} \mu_X\\\\\mu_Y \end{smallmatrix}) , (\begin{smallmatrix} \sigma^2_X && \sigma_{X,Y}\\\ \sigma_{X,Y} && \sigma^2_Y \end{smallmatrix}) \big) $ Would this affect how the sum $Z$ is distributed?


title: How do you do a power calculation for a 2-phase prevalence survey?

text: I'm trying to work out what sample size I need to calculate a prevalence rate within certain confidence intervals, using a two-phase design. There is a screener test. Everyone that screens positive is diagnosed. A proportion of those that screen negative are also diagnosed to check how well the screener worked. This will allow us to weight back the number of cases back to the full sample. Given:   * How well the screener is likely to perform (the positive and negative predictive values)   * An estimate of the overall prevalence   * How many of the screen-negatives will be diagnosed   * The desired standard error How do we calculate the total sample size? Thank you!


title: k folds cross validation on a multi-class dataset

text: **Cross validation** is one of the most important tools because it gives us an honest assessment of the true accuracy of our system. In other words, the cross-validation process provides a much more accurate picture of our system's true accuracy. **If we have for example a dataset that contains one class (lets consider a face dataset):** In this case, we divide our dataset into `k folds` (or `k portions`). A common value of k is 10, so in that case, we would divide our dataset into 10 parts. We will run `k rounds` of cross validation. In each round, we use one the folds for validation, and the remaining folds for training. After training our classifier, we measure its accuracy on the validation data. Average the accuracy over `k rounds` to get a final cross-validation accuracy.               Prepare our dataset     Divide it into 10 folds.     for i=1:10 % ten times     fold(i) for testing     the remaining for training      end     Final accuracy = Average(Round1, Round2, ...., Round10).       **Else if we have a dataset that contains multiple class (lets consider 3 classes : faces, airplanes and strawberry categories):** I don't know if my opinion is correct or not : each of the three categories is split into 10 folds. Does we measure the final accuracy as above? for example:               For the first round: take the first folds from the first, second and third categories and use them as testing and all the remaining folds (from all the categories) as training.            For the second round: take the second folds from the first, second and third categories and use them as testing and all the remaining folds (from all the categories) as training.          etc. until the round 10...      Does that is correct ? Does my opinion is correct? please I need your help and explanation. Any help will be very appreciated.


title: Categorical variable with a very large number of categories as a predictor

text: I am trying to use a categorical variable as a predictor in a supervised learning setting, but there are too many categories for the classification algorithm to handle, something like over a 1000 categories. What are some ways to get a manageable number of categories, is there a standard way of binning these categories? I suppose this binning should be performed on a training set disparate from a test set, to get a truer measure of out-of-sample error? If cross-validation is employed, I suppose the procedure should be run on on each fold.


title: Probabilistic latent semantic analysis vs bag of words

text: I am going through probabilistic latent semantic analysis (pLSA) for image classification. I did some image classification experiments with BOW (bag of words) and pLSA but based on my experiments BOW and pLSA performs similar. I could not find any improvement using pLSA. If pLSA does not perform better than BOW what is the purpose of it? (other than doing some image segmentation) Could you please suggest any computer vision references where pLSA significantly performs better than BOW.


title: Binary logistic regression with 3 similar outcomes

text: I was given three binary dependent variables, which are the following:               1) Paid in Full (Subject paid their full balance)     2) Settled in Full (Subject paid 80% or more of their balance in 1, 3, or 5 payments)     3) Rehab (Subject made 10 or more payments)      With these three binary dependent variables, the goal of this project is to build a predictive model to predict the odds of Paid in Full (PIF), Settled in Full (SIF), or Rehab. I have 250k records and ~500 potential predictor variables to work with. My question is how to best capture these 3 types of payers. I was told to simply combine all three dependent variables into one dependent variable, where: Paid = 1 (PIF, SIF, or Rehab) and 0 = Not Paid (All other records), and use the pool of predictive variables to predict this single outcome using binary logistic regression. However, I do not believe this is the best approach as the factors that result in PIF likely differ from factors that result in SIF and/or Rehab. Is there a better way to model this data than using only the single outcome variable/binary logistic regression? Any references or explanations are greatly appreciated!


title: How to develop intuition about parameter estimation in AR-GARCH model?

text: consider a typical ar-garch model:                y = sum( bi*xi ) + epsilon, epsilon~garch(p,q)      In such case a typical textbook says that we should first estimate all bi's and then estimate the garch parameters, because regression parameters and garch parameters are independent to each other (asymptotically.) Although it should be right, I do not understand the logic behind that. In general, if 2 sets of estimates (of model parameters) are independent, then can we estimate first set first then the second set? Can somebody please teach me the intuition over that?


title: Looking for correct distribution for modelling foosball scoring

text: I'm trying to adapt something like this association football prediction approach to 2v2 foosball games at our office. For football, they basically have an offense score and a defense score for each time. The game is a fixed length of time. They then use Poisson regression based on the scoring history to try to find out each team's offensive/defensive values. The problem is different for foosball, and I think Poisson may not be the correct distribution to use. For the following reasons:   1. We play first to 10 points, not for a fixed length of time   2. Players switch positions after their team scores 5 points Can anyone recommend a better distribution or a way of approaching the problem? We don't need a 100% optimal solution, but we are aiming for something that is reasonably ok. Bonus points for ideas on how to split the offense/defense ratings among the players.


title: How are the p-values of the GLM in R calculated?

text: I have been running some binomial logistic regressions in R on a data set and I realised that the p-values of the estimated coefficients are not computed based upon a Normal distribution. For e.g. I have the following result from the glm() function:                   Coefficients:                                             Estimate Std. Error z value Pr(>|z|)          (Intercept)                    -1.6127     0.1124 -14.347  < 2e-16 ***         relevel(Sex, "M")F             -0.3126     0.1410  -2.216   0.0267 *           relevel(Feed, "Bottle")Both    -0.1725     0.2056  -0.839   0.4013             relevel(Feed, "Bottle")Breast  -0.6693     0.1530  -4.374 1.22e-05 ***      I previously thought that the beta hats (estimated coefficients) are asymptotically normal and the p-value should be calculated using the normal distribution. But it seems from the p-values here that a t-distribution is used to compute them (I think). I know that the asymptotic condition does not hold when sample size is small but if so, what is the process that causes the estimator to be distributed with a Students distribution? And how do I find out the degrees of freedom for such a distribution?


title: Was Hitler a "black swan" event?

text: This question was inspired by another question on another SE site: http://history.stackexchange.com/questions/7931/has-any-other-book-in-history- been-as-prescriptive-as-mein-kampf If someone in 1909 had pointed to a certain 20 year old AUSTRIAN and said, "Within 25 years, this homeless high school dropout will make himself absolute dictator, not of his home country (Austria), but a more powerful neighboring one (Germany)," the statement would have been dismissed as "statistically improbable." Yet, unfortunately, it actually happened. Suppose the birth of someone like Hitler were a "six sigma" event: I've taken the liberty of "rounding" this to "one in a million." Then, in a world of seven billion people, wouldn't there be 7,000 potential Hitlers around? Then wouldn't the question be not, "do these people exist? but rather, "Will one of these people ever get into power?" Put another way, does a large enough population offset a "statistical improbability" to a great enough degree so that it becomes at least "somewhat probable"?


title: Finding the distribution and/or fitting a model (to a biological problem)

text: I am a bioinformatician and work on _RNA-Seq data_. The data contains a lot of **reads** (of length 80 bp in my case). These reads are fragments of those _genes_ that were _expressed_. I map them back to my _reference genome_ which helps answer many questions (gene expression, snp variants, alternative splicing etc...). I work on **alternative splicing** , where I am interested in characterizing events where same genes _express_ in more than 1 ways leading to different _transcripts_. Let's say I my data consists of 4 libraries (or biological replicates) for each of the 2 plant subspecies (closely related ones belonging to same species). And, I would like to look at _exon skipping events_ (where certain exons from coding genes sometimes gets excluded). _My question is about the statistical analysis towards characterizing these exon skipping events between the two subspecies over those replicates_. The basis of my question (or line of thinking and therefore the complexity), however, comes from statistical analysis on RNA-Seq data for gene expression. So, I'll frame the problem from there. Please bear with me. For every gene, $g_{i}$, its expression is = number of reads mapped to that gene (or _count data_ ). For one gene, $g_{1}$, say, the data would look like this.                    Lib1    Lib2    Lib3    Lib4     sA   400     420      600     250     sB   180     229       60     125      Since its _count data_ , given the concentration $q_{1}$ of gene $g_{1}$, biostatisticians usually model the number of reads, $r_{1}$ that map to a $g_{1}$ as a poisson model. It can be shown that $r_{1} \propto q_{i}$, or $r_{1} = sq_{1}$ (poisson parameter). That is, $p(Y = r_{1} | sq_{1}) \sim poisson$. **One question here** : This means, that if I were to measure for the same gene, the read count and note down the concentration, then, I would end up in a poisson model, right? From here, since we can not know the concentration $q_{1}$, they make use of the biological replicates, to estimate the variation. If $Q_{1}$ represents the concentrations over (all) the biological replicates, then under the assumption that this follows a **gamma distribution** (for the sake of mathematical convenience), the read counts $R_{1}$ for gene $g_{1}$ would then follow a **negative binomial distribution**. This is an overdispersed poisson model, it seems. Then a statistical test is performed to estimate if the difference in expression is statistically significant under a given $\alpha$. There are R-packages available that implement this. Now, **coming back to my problem** , as opposed to looking at every gene, for every _exon_ , $E_{i}$, I find out the number of times this exon is _spliced out_ (= not included in the transcript) and the number of times it is included, in all libraries, in both subspecies. So, for an exon, $E_{1}$, the data would look like this, for example: The first number in each entry is the number of _exon skipping events_ and the second is the _normal events_.                    Lib1    Lib2    Lib3    Lib4     sA   2, 80   1, 65   0, 40   2, 66     sB   10, 120 0, 22   8, 90   4, 90      The difference is that I have two values for each entry. I would like to find out, for each exon, if there is a difference between the two subspecies. My immediate thought was to calculate the ratio, for each entry, and then use existing model, however, they are on based on count data. So, my questions are these: 1) In general, is it known/ possible to estimate what distribution the ratio of two counts would follow?   2) Are there any other approaches ( for example, **generalized linear models** ) which could from these data (including dispersion from replicates) help me calculate if the occurrence of exon skipping events between the two subspecies are statistically significant? PS: In case somethings are not clear enough, I'd love to clarify them (I am by no means a statistician). I would appreciate any thoughts you guys have on this problem. Thanks again!


title: Vertical line graphs in R

text: I want to draw vertical line graph (like the first graph on this page): ![enter image description here](http://i.stack.imgur.com/UXcb2.gif) The deal is I have about ~15k data points whose frequencies (along Y axis) I want to present. I also want to be able to to present the frequencies in different colors : for example for the 1st 5k points, I want yellow colored line graphs for the next 5k I want red etc. How can I do this in R? Pointers to packages (and to helpful source code where possible) will be very useful. Beginner in R here. Thank you.


title: Survey questionnaire to find out the important factors?

text: For my masters dissertation I am currently preparing a Survey Questionnaire. The dissertation main aim is to find the factors which have an impact on failure of IT related projects in SMEs. I have done some extensive reading in the respective area and shortlisted around 30 factors. My current plan is to use each of these factors as a Likert item, with 5 options. Then the idea was to find the mean and the standard deviation for each of these item. Then from the survey result shortlist the 10 most selected factors with the highest mean. But googling around especially on this site I have read that it is advisable not to calculate the mean of the likert item, mostly because they are ordinal. Is likert approach the right one? if yes, are there methods (other then mean calculation) to find the 10 most selected factors from the survey? Help on this would really be appreciated. Thank You


title: Markov chain from Poisson: null or positive recurrent?

text: Let $K_t$ be a Poisson process with rate $1$ and $X_n=K_n-n$ $, \ \ \ n\in \mathbb{N}$ Now I want to determine whether it is null or positive recurrent, we already know it is recurrent. I tried to compute the mean return time to zero, but, using the law of total probability, I got huge sums... If $T=\min(n \geq 1 : X_n=0 |X_0 = 0) $ then I get: $$ \mathbb{P}(T=k) = \sum_{t_1\not = 1}\sum _{t2\not = 2} \cdots \sum_{t_{k-1}\not = k-1} \mathbb{P}(K_1 = t_1)\mathbb{P}(K_2 = t_2)\cdots \mathbb{P}(K_{k-1}=t_{k-1})\mathbb{P}(N_k=k) $$ It does not look right.... How can I go on? thank you very much


title: Longitudinal relationship between chocolate consumption and happiness: repeated measures ANOVA?

text: I am using SPSS and having some trouble with a research question which is analogous to the hypothetical question: _Is there a longitudinal Relationship between Happiness and Chocolate Consumption?_ Let’s say I take a sample of people and contact them when they are aged 14 and aged 18 and ask them: a) What is your chocolate consumption in grams per day. b) Are you happy? I have my fictitious data in the following wide format: ID HAPPY.14 HAPPY.18 CHOC.14 CHOC.18 1 YES YES 100 5 2 YES NO 50 30 3 NO YES 30 50 etc. _I would like to know if the mean chocolate consumption per day is higher among happy people than those who are not happy while accounting for the fact that I have taken repeated measures of both chocolate and happiness at the two time points._ **Approach 1** I suppose one way of doing this would be to do an ANCOVA, using time (before/after) as a grouping variable and controlling for happiness status. However, I think this may be inadvisable as correlation between the two time points would be neglected. **Approach 2** I understand that one valid approach for this should be a repeated measures ANOVA. I’m just not sure how to do this correctly in SPSS. I have specified my within-subjects factor as chocolate consumption with age 14 data as level one and age 18 data as level 2. What I’m uncertain about is the next step – specifying covariates and between individual factors – I have the option of adding HAPPY.14 and OR HAPPY.18 as a between individual factor. If I add both the output tells me about the effect of HAPPY.14 and HAPPY.18 as you’d expect, not about the “effect” of happy (YES/NO)per se. I realise it’s a basic question. Any feedback would on either of the two approached would be greatly appreciated.


title: Expected number of shipments and its standard deviation

text: Recent history suggests that one supplier fails to meet this new specification 20% of the time. Assume that the next 15 batches of this alloy are a random sample. How can I find the expected number of shipments that do meet the new specifications, and the standard deviation?


title: What's the minimum number of individuals in a group for repeated-measures ANOVA?

text: I'm in the process of running a repeated-measures ANOVA on two groups of participants from an experiment. In one group, there are 10 participants, and in the other group, there is only 1 participant. This is because the 1 participant is a patient who is being compared to the 10 other participants. Now, my question is: is it fair/allowable to run an ANOVA comparing the patient with the other participants? I wasn't sure, so ran it in SPSS, and it seems to work, with a significant result. Is this naughty?


title: A smart way of clustering a collection of sets based on an inherent hierarchy

text: Given a collections of sets, which have an inherent but unknown (at runtime) hierarchy, I would like to cluster them based on the sub/super-relationships with respect to their elements. Let me try and illustrate this with a overly simplified example:               Set 1 = {a, b, c, d, e, f}     Set 2 = {a, b}     Set 3 = {a, b, c ,d}     Set 4 = {a, c, d, f, g, h}     Set 5 = {d, f}      In this example, there would be two main clusters with the following relations; cluster 1: Set 1 $\supset$ Set 3 $\supset$ Set 2;   Set 1 $\supset$ Set 5 ... and cluster 2: Set 4 $\supset$ Set 5 The way I see it, the complications here from a standard clustering approach are; 1) I can not come up with a good measure of correlation between sets that are to be clustered. I was initially thinking of using the number of common elements but then the following scenario (which is essentially rather likely) complicates things: $_s(Set1 \cap Set2) = 10$ $_s(Set1 \cap Set3) = 10$ $_s(Set3 \cap Set2) = 0$ 2) In theory there is no reason for a small set to not be sub-set under more than one superset. This effectively makes any Tree-based data structure unusable, or am I mistaken on this point? I did some googling, checked on both StackOverflow and here briefly but havent really found something that is useful. Before I start implementing something in Java from scratch I was wondering if anyone had ideas or previous experiences on something like this. If there are libraries/functions one can use for this purpose it would be pretty cool, although I doubt there is something like written in Java. I know that most of you use R, but as I said, the rest of the software is written in Java so I'd prefer to keep things there, if at all possible. Thanks, **EDIT:** Following @whuber's comments I'll try and clarify the question further. I believe a significant portion of reasoning behind the question got lost when I tried to generalize and abstract the concept. So here it goes: The sets mentioned above are gene/protein sets, and the elements are then genes/proteins. As these entities work in connection with one another, one speaks of functional groups/sets. However the databases that hold this data usually have a high degree of redundancy, in the sense that Set A usually has all the elements of Set B, C .. etc. My whole project is based on analyzing these sets, and when I am done with the analysis and present my results I have a long set of these sets with associated scores. However highly scoring sets sometimes cluster, they may or may not be in the same super-sets. Thus the need/desire to cluster these in a structure like a dendogram. Thus one can overlay the scoring data, with the hierarchy data. On a side-note: I was recommended by a colleague of mine to consider spectral clustering, on which I will read more in the coming days to see whether or not the method can be used here or not. I hope these notes make things more clear now, I'd do my best to further develop the ideas if necessary. Thanks again!


title: Why is R giving a bad result of addition?

text: I can't seem to do a comparison with a sum of numbers and get the right answer (lines alternate between my input and R output): > (.6 + .3 + .1) == 1 > > [1] FALSE > > (.6 + .3 + .1) > > [1] 1 > > 1 > > [1] 1 I've tried with and without parenthesis, tried comparing to 1.0, tried using as.numeric, but can't get it to work.


title: Error Correction Model & Trend Stationarity

text: I'm familiar with a traditional Error Correction Model (ECM) of the form $$Y_{t}-Y_{t-1}=\Pi Y_{t-1}+\sum_{i=1}^{p-1}\Gamma_{i}\left(Y_{t-i}-Y_{t-i-1}\right)+\varepsilon_{t} $$ where $\Pi=\alpha\beta'$, $\alpha$ is the long-run adjustment, $\beta$ is the matrix of cointegrating vectors, and $\beta'Y_{t-1}$ is stationary. For the data set I am looking at (includes economic variables, interest rates, and stock prices), I was testing cointegrating vectors and found that some are stationary and some are trend stationary. I think I should include the trend stationary ones, but I'm not entirely sure what's the best way to proceed. For instance, one approach is to estimate the model $$y_{t,i}=a_{i}+b_{i}t+x_{t,i} $$ to de-trend each variable and then apply the above ECM to $X_{t}$. However, I wonder if it is better to keep everything in terms of $Y_{t}$. Hence, I was thinking about estimating $$Y_{t}-Y_{t-1}=A+Bt+\alpha F_{t-1}+\sum_{i=1}^{p-1}\Gamma_{i}\left(Y_{t-i}-Y_{t-i-1}\right)+\varepsilon_{t} $$ where $F_{t}=\beta' Y_{t}$ and $\beta$ includes all significant cointegrating vectors (including trend stationary ones). Whereas the traditional ECM is stationary since you are adding stationary series, this is a linear combination of a deterministic time trend, trend stationary series, and stationary series. Is this what should be done? Should I de-trend $F_{t}$ if I include the deterministic time trend, and vice-versa if left out?


title: P-value for multiple regression changes as more points are added

text: I am working on a multiple regression model that will forecast the value of loans being granted within the current month. The data points are broken down per day ( Jan will have 31 data points etc) and I will be updating the model every week. I have 4 independent variables driving the model. I am currently testing the model and trying to understand what is happening in particular to the $p$-values. On Day 8, 2 of my Variables (loans and declines) have $p$-values of               loans   0.014030324     declines    0.980464984      On day 15 when I run the regression I get these               loans   0.003114471     declines    0.023498327      I am just wondering why the "declines" $p$-value is now starting to show as a significant value. Is it because the model is working over more data points or is there something in the data that is suggesting that this variable is becoming more significant?


title: Several intraclass correlations (ICC) when modeling heterogeneous variances?

text: In a nested design (e.g., students in schools, ants in colonies, individuals in a species, leaves in a tree), ICC is defined as the proportion of variance among observations from the same subject (over total variation, i.e., variance within + residual variation). See, e.g., answer here When in the design there are one or more categorical factors grouping several subjects (between-subjects factors, e.g., type of school, species of ants or trees), possibilities for heterogeneous variances appear. Residual variances can be modeled as heterogeneous in a mixed, random-intercepts model, for example as in (R code, nlme package):               lme(fixed = Y ~ A, random = ~1 | Subject, weights = varIdent(form = ~1 | A))      where there is then a single estimated variation among subjects, but as many residual variations estimated as levels in A. (1) Is it valid/appropriate to estimate and report several ICCs, one for each level of the fixed factor? This would represent the correlation among observations from the same subject for subjects within each of the groups defined by the factor A. (2) Moreover, if variances among subjects are modeled as heterogeneous among groups (this is, random intercepts for each level in A come from different normal distributions, see here), the answer to (1) is the same?


title: Puzzled by derivation of time series prediction based on its log

text: From _Introductory Time Series with R_ : > If the random variation is modelled by a multiplicative factor and the > variable is positive, an additive decomposition model for $\log(x_t)$ [where > $x_t$ is the observation at time $t$] can be used: > >> $\log(x_t) = m_t + s_t + z_t$ >> >> [where $m_t$, $s_t$ and $z_t$ are, respectively, terms for the trend, seasonal effect and an error term at time $t$] > > Some care is required when the exponential function is applied to the > predicted mean of $\log(x_t)$ to obtain a prediction for the mean value > $x_t$, as the effect is usually to bias the predictions. If the random > series $z_t$ are normally distributed with mean 0 and variance $\sigma^2$ , > then the predicted mean value at time $t$ based on Equation (1.4) is given > by > >> $\hat{x}_t = e^{m_t + s_t} e^{\frac{1}{2}σ^2}$ The discussion about the caution that is necessary with this sort of model that follows makes sense and, unless it is necessary, I will omit it. What is confusing me is why $e^{\frac{1}{2}\sigma^2}$ is used to represent the error term. There is apparently some point about mathematical statistics that I'm not aware of, but would like to know. Thank you.


title: MCMC and terrible convergence

text: The code and data I am borrowing come from http://www.perossi.org/home/bsm-1 under CS 5 from the book _Bayesian Statitics and Marketing_. I tried applying their model to another dataset and am getting failed convergence/terrible results. I am simply looking for conceptual reasons why I might be getting what I am seeing. I am able to duplicate the results from the book using the code. I then took the code, modified it to run on my own set of data and get a failed convergence (almost perfect reject rate, log likelihood still decreasing, parameters stuck or widly moving across their support space). Almost immediately, my beta_ij parameters (see model description below) get stuck at zero and never budge. When I only run the code on a subsample of records, I get reasonable results. There are a lot of zero quantities in the data where x_i's would be zero for many of the products. The model is of the form: Utility of product i and household j = beta_ij + delta_i log( x_i + 1), where beta_ij are drawn from a normal distribution, delta_i from a uniform (-1, 0) and x_i are the quantities purchased of product i. The random effects model parameters are given priors, the mean, beta_bar is modeled as a multivariate normal and the variance-covariance matrix is modeled as an inverse-wishart. The code draws first the prior parameters given the household beta_ij's, then draws for the delta_i's and finally updates the beta_ij's given these new parameters. For the book reference, these are formulas CS5.4 - CS5.7 on page 271 of _Bayesian Statistics and Marketing_ by Rossi, Allenby and McCulloch. Thank you and any thoughts would be appreciated.


title: What is the correct way of comparing two very small groups (N = 3)?

text: I have 3 samples of water contamination (bacteria/L) from a lake which I want to compare to another, more "pristine" lake. So I have 6 values in total, 3 belonging to each lake:               contamination,  lake type     7576.48,        pristine     7472.81,        pristine     7511.55,        pristine     6803.65,        test lake     6692.63,        test lake     7128.27,        test lake      What is the best test to compare them? In the boxplot they look different. Initially I thought about a Student-t test, but with _N_ = 3 it looks like it will be rather poor, right? Should I go then directly for the Wilcoxon rank sum test instead? Or what do you do when you have just 3 cases for each group?


title: Metrics for covariance matrices: drawbacks and strengths

text: What are the "best" metrics for covariance matrices, and why? It is clear to me that Frobenius&c are not appropriate, and angle parametrizations have their issues too. Intuitively one might want a compromise between these two, but I would also like to know if there are other aspects to keep in mind and maybe well-established standards. Common metrics have various drawbacks since they're not natural for covariance matrices, e.g. they often don't especially penalize non PSD matrices or don't behave well w.r.t rank (consider two rotated low-rank covariance ellipsoids: I'd like the same-rank intermediate rotation to have lower distances than the componentwise average, which is not the case with $L_1$ and maybe Frobenius, please correct me here). Also convexity is not always guaranteed. It'd be good seeing these and other issues addressed by a "good" metric. Here is a good discussion of some issues, one example from network optimization and one from computer vision. And here's a similar question getting some other metrics but without discussion.


title: How can I create a Graph of a Probability Density Function from Percentiles?

text: I play a game online (Heroes of Newerth) which has a large ladder of players, each player having a couple of different ratings. I've manually gathered the rating data for all percentiles of players but am blanking at how to turn this into a histogram and compute the Mean and Standard Deviation. I'd also like to have a graph to show the curve. I could easily do so using my TI-83+ and I consider myself very proficient at Excel, but I don't see an obvious way to do this. Below is a link to the data, or if you'd prefer to give me instructions on how to create the distribution I'd be happy to do so myself. https://spreadsheets.google.com/spreadsheet/ccc?key=0Atwn_lcLizk9dDZ4WW1YV3BZQnFSaGNoQVdRa0JVZ3c&hl=en_US&authkey=CKyI38QJ#gid=0


title: Why is R plotting standardized residuals against theoretical quantiles in a Q-Q plot?

text: In R, why do the default settings of `qqplot(linear model)` use the standardized residuals on the y-axis? Why doesn't R use the "regular" residuals?


title: How to compare two datasets with Q-Q plot using ggplot2?

text: As both a stats and R novice, I have been having a really difficult time trying to generate qqplots with an aspect ratio of 1:1. ggplot2 seems to offer far more control over plotting than the default R plotting packages, but I can't see how to do a qqplot in ggplot2 to compare two datasets. So my question, what is the ggplot2 equivalent of something like:               qqplot(datset1,dataset2)


title: What does the gene set score tell in GSA analysis?

text: I carrying out a gene set analysis to investigate if a list of genes are differentially expressed before and after a medical treatment (paired analysis). The list of genes was made (a priori) by a promoter analysis, where genes that are likely to bind to a specific DNA promoter region was computationally identified. This gene set consists of approximately 4000 genes, and is likely to contain many genes that is not a target (false positives). I have used GSA (Brad Efron and R. Tibshirani) on microarray data (24000 genes, maxmean statistics, paired, permutation 5000), but the result puzzles me:               Gene set score: -0.0325     p-value: 0.0238       Note that FDR is not important here, since only one gene set is being studied. Why is this gene set score so low? What does the gene set score actually tell? Is the gene set score an averaged T score value? If so, why is p significant (<0.05)?


title: Estimating the total number of customers based on a subset of orders

text: Assume there are a number of customers, each of which has made one more orders. You now see a sample/subset of orders (lets say 1,000,000 out of 10,000,000) but do not know the total number of customers. How to estimate the total number of customers in the full dataset, based on the sample? I can assume that the number of orders per customer follows a power law, but it would be better if I could do the estimate without this assumption, but base it purely on the sample.


title: Computing Jensen-Shannon Divergence between discrete and continuous distribution

text: Is it possible to compute the Jensen-Shannon divergence between a discrete and a continuous probability distribution, e.g. between a standard normal distribution and a distribution taking the values 1,2,3 each with probability 1/3? Or are there any other divergences which can be used as a similarity measure of discrete and continuous distributions.


title: improving fitting results of neural network

text: I've trained fitnet network for prediction steel's yield stress with MATLAB `ann` toolbox. The neural network should predict yield stress. I have about 250 vector inputs. Every input vector has 12 parameters as inputs plus one parameter as target (yield stress). `ann` should predict yield stress. The inputs are weight percent of chemical elements that are used in steel. Here is my code snippet:               inputs=transpose(inputs);     targets=transpose(targets);          net =  feedforwardnet(20);     %net.trainParam.max_fail=1000          % Setup Division of Data for Training, Validation, Testing     net.divideParam.trainRatio = 70/100;     net.divideParam.valRatio = 15/100;     net.divideParam.testRatio = 15/100;          % Train the Network     [net, tr] = train(net,inputs,targets);      Here are the `ann` training results: ![enter image description here](http://i.stack.imgur.com/I6zv9.png) Here is the regression plot: ![enter image description here](http://i.stack.imgur.com/wusUA.png) Here is the actual output (green) and neural network output (red) charts in one diagram: ![enter image description here](http://i.stack.imgur.com/hVSW9.png) As it can be seen the charts are not overlapping. This means that network has not worked well. How can I improve results? Should I use other network or topology. Is my data very noisy or incorrect or they are less to train `ann`? I've taken data from here: http://www.msm.cam.ac.uk/map/data/materials/welddb-b.html


title: Frequency weights in R (multilevel using lme4)

text: I'm building a multilevel model with lme4 (two levels: a class and a student level with educational data). I'm using the ICCS 2009 data: http://www.iea.nl/iccs_2009.html). ICCS 2009 advise the use weights. I made a “within-school student weight” as the product of the class and student level weight factors. I also made the sums of the weights equal to the sample size. My question: how can I specify in R that I want to use this weight as a frequency weight. Do I have to use the "survey" package?


title: Explosion in stochastic gradient descent

text: I am using stochastic gradient descent for minimizing an objective function of the following form ![enter image description here](http://i.stack.imgur.com/b2km3.jpg) with such an R script               require(rootSolve)          zz <- c(86.02,86.85,87.48,88.69,88.82,87.46,87.4,87.6,             88.38,86.91,88.02,89.4,89.96,90.52,91.07,91.03,             90.14,90.1,91.11,91.24,90.7,91.76,92.52,92.8,             92.19,91.79,91.53,92.38,92.08,92.09,90.47,91.46,             90.79,89.65,90.23,91.45,91.12,114.67,114.97,115.64,             116.54,117.33,117.56,119.08,120.35,121.03,119.31,             121.06,122.5,122.26,122.06,122.37,122.65,122.24,             121.87,124,123.01,122.76,123.42,124.78,125.94,124.93,             125.4,125.45,126.43,125.49,125.11,122.75,124.09,122.96,             121.57,120.59,121.7,122.69)     zz <- matrix(zz,ncol=2)     nu <- 0.000002076897     lambda <- 2     m <- 5     xprev <- c(0.5471319,0.4390035)     yy <- NULL     g <- function(x,y1,y2) sum(y1*x[1]-y2*x[2])^2-lambda*sum((y1*x[1]-y2*x[2])^2)     for(i in m:nrow(zz)){         y <- zz[(i-m):i,]         xprev <- xprev - nu*gradient(g,xprev,y1=y[,1],y2=y[,2])         yy <- rbind(yy,xprev)     }      Usually it works fine and found parameters are in the reasonable bounds, but sometimes they start to explode for some reason. For example here I used about 1300 data points during the in-sample period to learn the model (nu,lambda,m) and got the following x's values dynamics (which is normal) ![enter image description here](http://i.stack.imgur.com/Q81TF.jpg) but after about 600 data points during out-of-sample period where data didn't undergo any significant changes (as I can visually detect it) x's values started to explode ![enter image description here](http://i.stack.imgur.com/tPNJI.jpg) for this particular data set ![enter image description here](http://i.stack.imgur.com/0t9fw.jpg) Could you please help me to understand what can be the reason of such an explosion and maybe how to struggle with it. Thanks!


title: Statistical models that incorporate factors stepwise?

text: I am working on a dataset in R that essentially follows 20 groups of individual organisms through time. I am interested in assessing what factors impact their mortality. I have 5 factors that I have measured that I believe effect the final mortality (%) of the individuals in each group. Obviously as time goes on, some groups will have been completely eliminated and therefore some factors will have an NA (missing) value for that particular group while others will have a zero designating that nothing was measured for that factor. Any ideas on models that can take into account how predictors in succession can impact some final value (mortality in this case)?


title: Stochastic Diagonal Levenberg Marquardt in Convolutional Neural Networks

text: Can someone explain this method of optimized convergence used in CNNs? I understand it involves generating the hessian matrix for every epoch, but can someone outline the steps? Thanks!


title: Measuring real advertising effect: control group issue

text: **Background** : The goal is to measure impact of advertising technique - retargeting. Retargeting is about showing related ads to prospects who add something to shopping cart but abandon it. We would be paying a fixed % of all sales that are made by visitors that are successfully retargeted. The issue is that such program would cannibalize other marketing channels to some extent, i.e. some of those sales would have happened anyway. We would like to **place a cap on the payouts** to advertising company by measuring the real impact of their program. **Our Approach** : We are thinking or constantly splitting all website visitors into two groups:   * **Control Group** \- will not be retargeted (some of them will come back through other channels)   * **Target Group** \- will be a target of the retargeting program (some of the sales that are attributed to the program would have happened anyway - this is what it is all about) I came up with this formula: $$cap = (R_t - (R_c * \frac{C_t}{C_c})) * x$$ Where:   * $R_t$ - Target Group Revenue (last 30 days, e.g. \$100000)   * $R_c$ - Control Group Revenue (last 30 days, e.g. 20000)   * $C_t$ - Target Group Size ("seen" in last 30 days, e.g. 8000)   * $C_c$ - Control Group Size ("seen" in last 30 days, e.g. 2000)   * $x$ - Commission Rate (e.g. 10%) We are simply calculating how many times the target group is larger than control group, multiplying the result with revenue from the control group. This way we get an approximate estimate of how much revenue we could have made from the **target** group if they were **not** retargeted. I guess this logic is clear enough. **The Question** : I have been teaching myself statistics for the last few hours but still do not clearly understand all the confidence interval, confidence level and other things like that. I need to find out what % of all visitors to assign to the control group to have trustworthy data about which I could say something like this - "By paying out \$1000 we can be 99% sure that we are not overpaying and not underpaying by more than \$50". I need to prove it by calculations which should be described in plain English as well.


title: Predicting dichotomous outcome of temporal data set with covariates

text: I have a set of data, with outcome and time-varying variables, for patients during the course of their respective stays in the hospital. There is a dichotomous outcome on the last day. The length of stay may be different for each patient. I am most interested in making a predictive model, in the same spirit as one might do with standard regression. I believe that the patient may develop some "momentum" which would tend lead to one outcome or the other. My previous attempts to model this data have been limited to logistic regression. I am beginning to think about lmer and friends, but I'm not sure if it would be suitable for the directionality in time. Also, I have considered some stochastic approaches, but I haven't made much progress with hitching this into a binary outcome or even considering the covariates. Any suggestions would be most welcomed. Fictional Data:               Patient No. | Day of Stay | Outcome | Age | Gender | Blood Value | ...              1  |           1 |       - |  70 |      M |       123.1 | ...              1  |           2 |       - |  70 |      M |       134.2 | ...              1  |         ... |     ... | ... |    ... |         ... | ...              1  |          10 |       1 |  70 |      M |       148.3 | ...                ...           ...       ...   ...      ...           ...            100  |           1 |       - |  54 |      F |        98.3 | ...            100  |           2 |       - |  54 |      F |        95.2 | ...            100  |         ... |     ... | ... |    ... |         ... | ...            100  |           6 |       0 |  54 |      F |        54.1 | ...


title: How do I measure the statistical difference between two measures?

text: I am measuring rendering time on Stack Overflow, specifically pre/post changes. I have data as such               pre_post | RouteName      | count | avg | st_dev | min |  max     ---------+----------------+-------+-----+--------+-----+-----     1. pre   | Questions/Show | 24736 |  11 |     15 |   0 | 1004     2. post  | Questions/Show |   842 |  16 |     26 |   1 |  453      From this I would like to understand:   1. Are the requests statistically slower or faster after the change, or there is no discernible change?   2. How many measurements do I need before having reasonable certainty of point 1.? * * * Here is what the peaks look like (after cleaning, the data is not the same as above!) ![nice peaks](http://i.stack.imgur.com/pSFbi.png) Here is a link to the raw data. * * * Here's what the raw data looked like before cleaning and normalizing ![twin peaks](http://i.stack.imgur.com/y03N1.png)


title: Weighted least squares to correct for heteroscedasticity

text: I would like to use a weighted least squares (WLS) regression to perform tests on heteroscedastic spatial data. Each data point represents the mean of some variable over an area, and the sample sizes between the areas vary, so intuitively the things I'm measuring are more error prone in areas with a small sample size. The variance of a mean is inversely proportional to the sample size, so presumably I should weight the regression by the inverse of this i.e. weight each point by the sample size that point was derived from. But do I use, in this, the sample size from the dependent variable, the sample size from the independent variable - or both?


title: Interaction term in a linear log model

text: I am using a linear-log model to test whether overseas development assistance and remittances positively affect FDI in cases of good governance and financial market development. Let's say I want to interact the log of net official development assistance received (ODA) with governance. Governance is estimated as a score from -2.5 to 2.5. So I cannot apply the log transformation to governance, but can I still interact the log of ODA with governance? My hypothesis is that ODA will positively affect FDI in cases of good governance.


title: Pre-truncation moments for truncated multivariate normal

text: Suppose the random variable $Y$ has a multivariate normal (MVN) distribution, and consider truncating $Y$ in some way to create $T$. Given $T$'s mean and covariance matrix, I'd like to obtain $Y$'s mean and covariance matrix. Below I state the above problem more precisely, ask three specific questions, give an example, and close with several remarks. Suppose $Y$ is a random variable whose $d$-variate normal ($d$VN) distribution has mean $\mu = \mathrm{E}(Y)$ and covariance matrix $\Sigma = \mathrm{Cov}(Y)$, so $Y \sim \mathcal{N}_d(\mu, \Sigma)$. Let $T = Y | (Y \in \mathcal{A})$ be a truncated version of $Y$ that takes only values of $Y$ in the set $\mathcal{A}$ and has mean $\eta = \mathrm{E}(T)$ and covariance matrix $\Omega = \mathrm{Cov}(T)$. For instance, in the example below $\mathcal{A}$ consists of vectors that form valid correlation matrices. Collect $\mu$'s and $\Sigma$'s $d(d+3) / 2$ distinct parameters in the vector $\theta_Y$, and collect their counterparts from $\eta$ and $\Omega$ in $\theta_T$. Here are three questions about obtaining $\theta_Y$ for a specified $\theta_T$:   1. Under what conditions does a unique $\theta_Y$ exist for a specified $\theta_T$? This may be a question of invertibility: If $f$ is a function such that $\theta_T = f(\theta_Y)$, under what conditions is $f$ invertible?   2. In a condition from #1 above, how can we obtain $\theta_Y$ for a specified $\theta_T$, either exactly or by numerical approximation or simulation? For instance, can this be framed as a problem of optimization or solving a nonlinear system of equations?   3. If we approximate or simulate $\theta_Y$ in #2 above, how can we control the approximation's or estimate's accuracy? For instance, we might want an estimate of $\theta_Y$ to be near $\theta_Y$ with high probability. **EXAMPLE:** Consider a fairly simple version of the above problem where $d = 3$, so $Y \sim \mathcal{N}_3(\mu, \Sigma)$, and $\mathcal{A}$ used for truncation is the set of all 3-vectors containing Fisher- _z_ transformations of a valid Pearson- _r_ correlation matrix's distinct elements (e.g., below the diagonal in row-major order). For instance, $\mathcal{A}$ contains [1.20 1.20 0.50] but not [1.20 1.20 0.40]. Let's specify $T$'s desired mean and covariance matrix -- with commas separating rows -- as follows, along with $\theta_T$ (e.g., $\eta$ followed by $\Omega$'s lower triangle in row-major order):   * $\eta$ = [0.40 0.60 0.80]   * $\Omega$ = [0.25 0.12 0.06, 0.12 0.16 0.06, 0.06 0.06 0.09]   * $\theta_T$ = [0.40 0.60 0.80 0.25 0.12 0.16 0.06 0.06 0.09] An ad hoc method -- applying an iterative truncate-and-transform technique each of several large 3VN samples -- yields the following estimates of $Y$'s mean and covariance matrix and, hence, $\theta_Y$:   * $\hat \mu$ = [0.329 0.631 0.832]   * $\hat \Sigma$ = [0.292 0.101 0.043, 0.101 0.163 0.064, 0.043 0.064 0.095]   * $\hat \theta_Y$ = [0.329 0.631 0.832 0.292 0.101 0.163 0.043 0.064 0.095] These estimates are meant to be within about 0.1% of their estimands (in $\theta_Y$) with probability at least .95. **REMARKS:** Answering the above questions would probably be easier if $f$ in #1 above were available in closed form, but that's difficult for some truncation rules of interest (i.e., some sets $\mathcal{A}$). For now I'd mainly like an _accurate_ solution to these questions, but I'd eventually like to create a _fast_ implementation. I'm also interested in extensions of this problem, such as working with a mixture of $d$VNs for $Y$ or finding $\mu$ and $\Sigma$ for specified $\mathrm{E}(G)$ and $\mathrm{Cov}(G)$, where $G = g(T)$ is a nonlinear function of $T$. These questions arose while developing meta- analytic methods for functions of multivariate effect sizes, where between- studies heterogeneity in an effect-size parameter (e.g., correlation matrix) induces heterogeneity in a function of that parameter (e.g., path-model coefficients).


title: Why is Hedonic Regression used instead of Linear Regression

text: Why is Hedonic Regression used (especially in housing prices) instead of Linear Regression? There do not seem to be any libraries in Python (and R) for Hedonic regression, is it too niched a regression type? Are there regression methods similar to Hedonic regression that I can use instead?


title: Multivariate values to acf() in R

text: I am trying to understand the concept of auto correlation and I am looking for some help to clear some doubts regarding my data. I have a time series data and it has repeated experiments. Each sample has 4 time points and corresponding values for the genes under study. I have used `ind` to represent each sample in my data. So `ind` 1 means rat which is studied over 4 time points and the samples are in rows. I have 400 genes and corresponding values for each sample and each time point. The genes are in columns. my dataset: `M1`               No      ind  tme      gene_1      gene_2      gene_3      gene_4      gene_5  gene_7          A1T1:2  1   -64 0.0307  0.0022  0.0010  0.0001  0.0007  0.0035          A1T2:2  1   8   0.0365  0.0031  0.0003  0.0002  0.0009  0.0043          A1T3:1  1   48  0.0182  0.0014  0.0001  0.0001  0.0005  0.0018          A1T4:1  1   96  0.0134  0.0010  0.0001  0.0001  0.0003  0.0015          A2T1:1  2   -64 0.0387  0.0032  0.0003  0.0002  0.0010  0.0051          A2T2:1  2   8   0.0264  0.0022  0.0010  0.0001  0.0007  0.0032          A2T3:1  2   48  0.0205  0.0017  0.0002  0.0001  0.0005  0.0022          A2T4:1  2   96  0.0161  0.0012  0.0001  0.0001  0.0004  0.0018      What I want to do is identify is simple: if there is autocorrelation in my data? my code:               s1<-read.table("M1.txt", sep=" ",header=T)     s2<-s1[1:4,1]     s2.v<-as.vector(s2)     acf(s2.v)      I get a plot with 4 sub plots. I do not understand how to interpret the result. Also Can I give all the genes as input (for sample 1) to identify autocorrelation? I am new to these concepts hence my understanding is basic. Thanking you for all the help.


title: A good primer on how people deal with cluster-wise autocorrelation in mixed models?

text: The courses that I've taken have been primarily in the economics department, but I've been learning more of stats more broadly recently. One thing that I don't understand well is how exactly people account for cluster-wise autocorrelation within mixed-modeling frameworks. When I talk about cluster-wise autocorrelation, I'm talking about panel data, or repeated observations of an individual over time, or more generally repeated realizations of one semi-stochastic, logically-distinct unit. The standard approach in econometrics is to control for unit-wise dummy variables ("fixed effects" in econometrics-speak) and then estimate the VCV as $$ X(X'X)^{-1}\left(\displaystyle\sum_j X'_j \epsilon_j \epsilon_j'X\right)(X'X)^{-1}X' $$ What this ends up accomplishing is inflating the standard errors to the degree that one observation in a unit is predicted by other observations within that unit. So say I want to start estimating random coefficient models instead of "fixed effects" regressions. How do people deal with the problem of this sort of autocorrelation within non-economics stats? Recommendations for things to read are welcome. edit: or is it even needed when everybody gets their own slope? But then I guess that if some variables are not specified as random, then you'd still need to account for autocorrelation somehow, otherwise the SE's will be biased downward. Right?


title: Model notation correct? ARMA(p,q)-GARCH(r,s)

text: For my variable $l_t$ I want to use an ARMA(p,q)-GARCH(r,s). So the mean equation follows an ARMA(p,q) and the conditional volatility is modelled by a GARCH(r,s). Is my notation correct? ![ag](http://i.stack.imgur.com/io3dX.png)


title: How can I calculate a critical t value using R?

text: Sorry if this is a newb question; I'm trying to teach myself statistics for the first time. I think I have the basic procedure down, but I'm struggling to execute it with R. So, I'm trying to evaluate the significance of regression coefficients in a multiple linear regression of form $$ \hat y = X \hat \beta$$ I think the t-statistic for testing $H_0: \hat \beta_j = 0, H_a: \hat \beta_j \neq 0$ is given by $$t_0 = \frac{\hat \beta_j - 0}{\text{se}(\hat \beta_j)} = \frac{\hat \beta_j}{\sqrt{\hat \sigma^2 C_{jj}}} = \frac{\hat \beta_j}{\sqrt{C_{jj} SS_{Res}/(n-p)}}$$ where $C_{jj}$ is the $j^{th}$ diagonal entry in $(X'X)^{-1}$. So far, so good. I know how to calculate all of these values using matrix operations in R. But in order to reject the null, the book says I need $$|t_0| > t_{\alpha/2,n-p}$$ **How can I compute this critical value $t_{\alpha/2,n-p}$ using R?** Right now the only way I know how to find these values is by looking in the table in the back of the book. There must be a better way.


title: Normalization of three different variables for linear addition

text: I am a research student and stuck at a point in my work and want your help. I have three different variables for each node of my network (Energy, Traffic Load and Link quality). These three variables have different units. I want to add them linearly to get a grade or rating in range of (0,1) and assign it to a node. So what i did is i normalized these variables in (0-1) range by dividing them with their maximum values and and then added them together. I want to know whether its correct?               Energy = Energy(i)/Energy(max) [0,1]     Traffic Load = Traffic Load(i)/Traffic Load(max) [0,1]     Link Quality = Link Quality (i)/Link Quality(max) [0,1]     Grade = Energy + Traffic Load + Link Quality


title: What is a representative sample?

text: I understand the rationale that underpins representative sampling (i.e. to avoid bias so that the sample 'broadly' represents the target population). Suppose the size of the target population is 600 and the population comprises of four types of units as follows:   * 100 white units   * 250 black units   * 50 green units   * 200 yellow units **Questions**   1. What would be a representative sample (overall and subgroup) in this case?   2. Is the sample in (1) appropriate for only descriptive analysis or inferential analysis or both?   3. If you have a representative sample, is _a priori_ power analysis needed?


title: Regarding precision and recall for the highly unbalanced validation data set

text: > **Possible Duplicate:**   >  Optimising for Precision-Recall curves under class imbalance I built a classification model and tested it against a validation data set. The positive set is composed of 86 cases and the negative set is composed of 1256 cases. The confusion matrix is as follows                                    True positive  True negative   precision        Predict positive    55               338          13.99%        Predict negative    31               918          96.73%                Recall      63.95%       73.09%        The precision and recall for this classifier is not good, especially for the positive precision. However, the negative cases are much more than the positive case. I am not quite sure that, for this kind of unbalanced data, can we still use the precision and recall as the performance evaluation as usual?


title: What statistical test should I use to assess rate of photosynthesis in leaf discs?

text: My biology test is to measure the rate of photosynthesis in leaf discs. It's a simple experiment where I will infiltrate the leaf disc space with sodium- bicarbonate and cause a pressure difference to sink the leaf disc. Then I will shine a light on it causing photosynthesis and time it in minutes for it to rise to the top. EDIT** I will have three leaves of which I cut them all in half. So I will have 6 halves of leaves. Each half of a leaf would have been irradiated so in turn I will have 3 halves which are irradiated, and the other three are left as control. I will then use a cork borer to extract three sample circles from EACH half of a leaf. So I, in turn will have 18 samples. This is the results of ONE sample:               Time/minutes  Has it floated?     0                  No     0.5                No     ""     ""     7.5                Yes      I have worked out that it will be two tailed since my hypothesis is noting whether there is a change in photosynthesis if irradiated with UV, but other than that, I'm stumped.


title: How wise is the use of ANCOVA when groups differ on the covariate?

text: In this case I presume loss of ANCOVA power, so I don´t know what type of analysis should I run next. There was significant difference in covariate between groups (p=0,008). Is there some solution? Could you help me please? Can I run ANCOVA if my groups differ on covariate (I know, that I can, but is it right solution)?


title: Is it possible to plot factor scores in a graph?

text: I have 5 factor scores which were saved out from Mplus6 to SPSS using commands below.               Savedata;      SAVE=fscores;      I have the variables in SPSS and used them in some analysis but my professor has suggested that I take the 5 factor scores and plot a graph like a half- normal distribution plot showing all 5 factors. I'm wondering if anyone knows how I can achieve this?


title: Sampling from a Deep Belief Network: Treatment of biases in directed part of the model

text: When generating samples from a DBN, how do you handle the biases that have been learned for the layers below? I know that you normally perform a number of Block Gibbs sampling steps in the undirected top 2 layers of the DBN, then propagate the sampled values of the penultimate layer down to the visible layer. So when propagating downwards, do you also include the visible biases learned during greedy layer-wise training, when sampling a layer below given the states of the layer above? Hope my question makes sense...


title: Is it bad to leave in poor predictors in regularized multiple regression?

text: There are some variables that I measured but strongly suspect are useless because (for example) almost all my data points scored the same on that (binary) variable. It's been put to me that I may well leave them in, because in a regularized multiple regression it won't hurt to leave them in, even though it probably won't help either. Is that true? If so/not, why?


title: Binomial Events with unequal probabilities

text: Imagine I have a database of 1000 students. I test these students every week with questions that can be true/false. How do I rank the students according to the historical track records? The problems are:   1. Each student has answered a different number of questions in the past. Many have only answered a few, others have answered over 50.    2. The questions are of varying difficulty. Sometimes the answer is obvious and other times it really is a 50:50 guess. Imagine I have a difficulty score for each question ranging from 1 - 100 (100 meaning very difficult). How do I find the students who are truly geniuses, and separate them from the students who only perform at chance level (or only slightly above)? Am I right to think that these are binomial events with unequal probabilities?


title: Chi-squared test vs gamma correlation coefficent

text: I have two 5-point Likert questions x and y. I want to test whether they are associated using the chi-squared test of independence. I summarized the responses into the following frequency table:                     x     y      1     3    39     2     7    76     3    38    65      4    97    30      5    75     8       I'm not really sure if the table is appropriate for verifying whether there is an association between x and y. Anyway, calculating chi-squared leads to $X^2 =$ 184.7218, with $p < $2.2e-16. Now, when I calculate the gamma correlation coefficient for x and y, where both variables are lists of all responses (containing a Likert item from {1,2,3,4,5}), I get for the gamma coefficient -0.02189, and a $p$ value of approx. 0.781. According to the high $p$-value of the gamma coefficient, there is probably no association between these two variables ($p>$0.05). Chi-squared states the opposite. What is the right interpretation here? I'm using R for calculations. Since I'm not experienced with statistics, I would like to know whether this is plausible, or I am doing/understanding something wrong?


title: Integration step in Bayesian linear regression

text: I want to obtain the integral: $$\int_{{\mathbb R}^p} \frac{1}{(2\pi)^{\frac{n}{2}}\vert\Sigma\vert^{\frac{1}{2}}}\exp\left[-\frac{1}{2}({\bf y} - {\bf X}{\beta})^{\top} \Sigma^{-1}({\bf y} - {\bf X}{\beta})\right] d\beta,$$ where ${\bf y}$ is an $n\times 1$ vector, $\beta$ is a $p\times 1$ vector, $\Sigma$ is a symmetric positive definite matrix and ${\bf X}$ is an $n\times p$ matrix. This expression is very common in Bayesian linear regression but I do not know the trick used to separate ${\bf X}$ from $\beta$ in order to obtain the integral by using its resemblance to the multivariate normal distribution. I would appreciate any guidelines on this. Thanks.


title: Determine if two proportions from the same sample are different

text: Is there a test I can use to compare two proportions from the same sample that are not mutually exclusive? For example, I would like see if the proportion of people that answered "yes" to a particular question is significantly different from the proportion of people that answered "yes" to a different question.


title: Time series and stationnarity tests

text: I perform some time series fitting with the help of the forecast and urca packages. I have a question regarding the correspondance between results coming from statistical test such as KPSS, ADF or ERS, and the results arising from automatical selection procedures such as auto.arima from the forecast package. Suppose we are in a situation where the outcome of the KPSS test is that the series does not need to be differentiated. I then apply the ADF test, which shows that the series is stationnary up to a trend (at+b). In that situation, I could force the auto.arima function to consider only ARIMA models that include a linear trend, but I can also explore all the models (including or not the trend) and let auto.arima choose the best one. And there are cases where the outcome is in contradiction with the test result, e.g. the chosen model does not include a trend while the test advises to include one. What is the best practice in such a case ? Should the test result be prevalent to the selection of the ARIMA model, or is it best to consider all possible models and choose the best one ? Thanks.


title: Using Rpart to find which factor influence the outcome the most

text: I have 3 factors x1, x2, x3, and one outcome y (True, False). x1 has 3 levels, x2 has 40 levels, x3 has 2 levels. I would like to find which parameter (x1, x2, or x3) and associated levels influence the most the outcome y. Is this possible with rpart? I tried:               fit <- rpart(y ~ x1, model="class")     Error in if (model) { : argument is not interpretable as logical          > fit <- rpart(y ~ x1)     > fit     n= 181365           node), split, n, loss, yval, (yprob)           * denotes terminal node          1) root 181365 16370 FALSE (0.90974003 0.09025997) *      I am not sure which model to take: "anova", "poisson", "class" or "exp"? The example in the manual (Kyphosis) uses only numeric input, but I have factors. * * * Here is a reproduction of the problem in a small data set:               > a = data.frame(c("A", "A", "B"), c(FALSE, FALSE, TRUE))     > a       c..A....A....B.. c.FALSE..FALSE..TRUE.     1                A                 FALSE     2                A                 FALSE     3                B                  TRUE     > colnames(a)=c("x1", "y")     > a       x1     y     1  A FALSE     2  A FALSE     3  B  TRUE     > rpart(y ~ x1, data=a, model="class")     Error in if (model) { : argument is not interpretable as logical     > a$y = factor(a$y)     > rpart(y ~ x1, data=a, model="class")     Error in if (model) { : argument is not interpretable as logical     > a$y = as.factor(a$y)     > rpart(y ~ x1, data=a, model="class")     Error in if (model) { : argument is not interpretable as logical


title: When Restricted Boltzmann Machine "dreams" how do you sample the results?

text: I've tried to sample by analyzing single step of the "dreaming" but the output seems nonsensical. There are some patterns but they don't really look like the data the RBM was shown ( for example, simple lines on a black and white 10x10 matrix of pixels ). This is what it dreamt when I've trained it on pictures of circles: http://i.imgur.com/MjFm8jc.jpg \- doesn't really look like a circle to me. Do I need to perform sampling over bigger timeframe? Can't you just look at the results on one step? From what I've read, people use Gibbs sampling - is it the way to see the RBM dream?


title: What is a proper scale to do the social network analysis?

text: # Problem My problem is to profile the individual user (i.e. mine individual user's interest, location, and many more). What we have as input are the **network structure** (e.g. linked-in network) including user set $V$ and link set $E$ and the **attributes** for each user (e.g. affiliations, employment, skills). There is NO explicit weight on the link, of course. The attributes are only available for a subset of labeled user $V_l$. A user $u_i \in V$ has have partial attributes or no attribute available. Now we want to infer those attributes from his friends. The intuition is that if all friends of $u_i$ are PhD students in University B, then $u_i$ is very likely to be a student in University B. To solve this problem, the first step is to select a proper set of user $V_r$ to propagate their attributes to $u_i$. The criteria of $V_r$ is   1. They are relevant to $u_i$.   2. At least one of them has available attributes. ## Existing Solution for $V_r$ selection In social networks analysis, I only know the following two kinds of networks:   1. **Whole social networks** \-- all the nodes and links among them.   2. **Ego network** \-- a focal node and all the directly connected nodes (1-hop neighbors) and all the links among those nodes. In my opinion, using the whole social network is not feasible. From social network service, we can collect millions of users, most of whom may be completely irrelevant to our target. Besides, the full network structure is unavailable to public in general. Therefore, I want a proper sub-network. However, the ego network seems to be too small. Because, some nodes may have only a few _directly_ connected neighbors. I want to know what other $V_r$ can I use and why. ## Extra Problem It is clear that the ego network is one connected component. Is the whole social network necessary one connected component? If not, then there should be several connected components in this "whole" social network. Then given one focal node, we can find one connected component centered on this node. Is there a terminology (such as ego network) on this kind of network?