title: Whether to transform non-normal independent variables in logistic regression?

text: I recently received the following email, which I paraphrase below: > I want to do binomial logistic regression with the data and I have non- > normally distributed IVs. I tried doing square root transformation on the > non-normal IVs. This was successful in normalising the distribution; > however, when I did the logistic regression on the data I got ridiculous > results i.e. OR: 14, CI 6 - 180. So I transformed the variables back by > squaring them, and then ran the analysis again. ### Questions   * Why do the odds ratios look ridiculous after applying square root transformation?   * Should one apply the square root transformation to non-normal predictors when doing logistic regression?


title: Bayesian vs ML, overfitting problem

text: In Bishop's PRML book, he says that, overfitting is a problem with ML, and Bayesian can avoid it. But I think, overfitting is a problem more about model selection, not about the method used to do parameter estimation. That is, suppose I have a data set $D$, which is generated via $$f(x)=sin(x),\;x\in[0,1]$$, now I might choose different models $H_i$ to fit the data and find out which one is the best. And the models under consideration are polynomial ones with different orders, $H_1$ is order 1, $H_2$ is order 2, $H_3$ is order 9. Now I try to fit the data $D$ with each of the 3 models, each model has its paramters, denoted as $w_i$ for $H_i$. Using ML, I will have a point estimate of the model parameters $w$, and $H_1$ is too simple and will always underfit the data, whereas $H_3$ is too complex and will overfit the data, only $H_2$ will fit the data well. My questions are, 1) Model $H_3$ will overfit the data, but I don't think it's the problem of ML, but the problem of the model per se. Because, using ML for $H_1,H_2$ doesn't result into overfitting. Am I right? 2) Compared to Bayesian, ML does have some disadvantages, since it just gives the point estimate of the model parameters $w$, and it's overconfident. Whereas Bayesian doesn't rely on just the most probable value of the parameter, but all the possible values of the parameters given the observed data $D$, right? 3) Why can Bayesian avoid or decrease overfitting? As I understand it, we can use Bayesian for model comparison, that is, given data $D$, we could find out the _marginal likelihood (or model evidence)_ for each model under consideration, and then pick the one with the highest marginal likelihood, right? If so, why is that?


title: fast approximate maximum likelihood parameter estimation

text: Given a univariate probabilty density function depending on a set of parameters, is there a generic method to estimate these parameters approximately without having to run the full (numerical) maximization of the likelihood given a dataset ? In the simple case of estimating the mean and standard deviation of a Gaussian distribution, one can calculate the estimators for mean and standard deviation from the 'moments' of the dataset. Is there a generalization to arbitrary PDFs ? Are there methods based on basis function expansion (e.g. Fourier transform, Wavelet transform) to achieve this ? (The original problem is to perform maximum likelihood parameter estimates with a large number of datasets with the same PDF. Computing coefficients related to the PDF only can therefore be expensive because they only have to be done once but computation of coefficients for each dataset should be fast)


title: Interpreting coefficients from a VECM (Vector Error Correction Model)

text: I would like to ask a question about error correction terms from VECM if I may. I am currently working on a lot of time-series data and one of the questions I would like to address is whether there is some relationship between different time-series. Reading around the subject it seemed to me that testing for cointegration would be a good way of showing that two or more time-series were in a long run equilibrium. However, before embarking on such procedures I wanted to be clear on a number of things. Firstly, (assuming there is a cointegrating vector) I have been trying to work out how to interpret the error correction terms from a VECM. I have been using the vars package in R, which provides the error correction terms in the summary table for a vecm model. Now, as I understand it the error correction terms describe how the time-series adjust to disequilibrium. I was wondering what the best way to interpret negative and positive error correction terms was? Given two time-series Xt and Yt, from what I have read so far it seems that negative error correction terms would mean that when Yt-1 is above its long- run level then ΔY*t* will be negative, pulling Y back towards its long-run relationship with X. On the other hand, I am less sure on the interpretation of a positive error correction term and I haven't been able to find a simple description. Finally, I have been trying to discover how one could interpret ECTs when two or more cointegrating vectors exist. So far, the only thing I have found in the papers I have searched is that its very difficult. Thank you for help, any advice, clarification etc. would be much appreciated.


title: What is the most robust outlier detection and removal technique, particularly for classification datasets?

text: I'm trying to prepare my datasets for a classification algorithm that I'm testing. It's a tumour classification dataset to be precise, but I'd like to know of techniques that are dataset as well as learning-model agnostic. What would be the most robust technique to deal with outliers, but without having to lose a lot of information?


title: Performance metric for categorical outcome prediction

text: I have 3 independent subjects each of whom have one attempt at guessing the number of balls in a hidden tray with 20 compartments. Each subject therefore has to produce 20 guesses, one for each compartment. The number of balls in each compartment is an integer beginning at 0 with no upper limit (although the majority of the compartments have fewer than 10 balls). What are several ways of measuring the performance of each subject in predicting the actual number of balls. Note that I have no specific criteria in mind as to what is the most important feature in their predictive ability. I am open to suggestions. Some dummy data in R:               actual_numballs <- c(4,2,3,3,3,4,3,7,5,3,3,3,2,3,2,2,3,3,4,8) # number in each of the 20 compartments     s1_predict_numballs <- c(4,2,2,3,2,4,3,7,5,3,3,3,2,3,2,2,3,3,2,9)     s2_predict_numballs <- c(3,2,2,3,2,4,3,5,5,3,3,3,2,3,1,2,3,3,2,3)     s3_predict_numballs <- c(4,2,3,3,3,4,3,7,5,3,3,3,2,3,2,2,3,3,2,3)      Now, I want to extend this so that the same subjects are now trying to guess the number of balls in the 20 compartments of several different trays, all of which are unrelated (independent) of eachother. How can I summarize the performance of the subjects across trays? By the way, any solutions provided in R would be much appreciated.


title: Distribution of $R^2$ for pairs of random variables

text: if I have two uniform random variables $X$ and $Y$, and I sample $N$ values for each, what's the probability of getting an $r$ Pearson's correlation coefficient (or Spearman correlation) between them of at least $x$? In other words what is the distribution of the $r$ values in this case? Similarly, what if $X$ and $Y$ were normal RVs, each with its own mean and standard deviation ($X_\mu$, $X_\sigma$ and $Y_\mu$, $Y_\sigma$, respectively), then what would the distribution of $r$ be as a function of these parameters?


title: Correlation and coefficients in ols

text: I am trying to fit a simple linear model with OLS. I have say 4 terms and the number of data points is >>4, so no explicit overfitting. 1) My terms or my independent varibales do have a strong correlation amongs each other 2) They are all positively correlated with my y My question is that on building the OLS model, I actually see the coefficient on some independent variables to be negative. I was expecting all positive signs actually, given the positive correlation. What could cause this? What is a remedy, Non Neg Least Squares.


title: Confused about simple hypothesis testing question

text: I am using R to do some hypothesis testing on a large sample (n=311) and I am getting confused with the results. R returns very different results than I am getting by hand, using the formulas in my text book. My data is the difference between (what I think) are two paired sets of observations and is summarized as follows: $$ \bar y_d = 0.0006, s_d=0.02, n_d=311 $$ My hypotheses are: $$H_o:\mu_d = 0\\H_a: \mu_d < 0$$ Which, according to my text book, gives me this formula for my test statistic: $$t= {{\bar y_d - D_0}\over {s_d/\sqrt {n_d}}} = {{0.0006 - 0}\over{0.02/\sqrt{311}}} = 0.591$$ And this is the start of where I am confused. R gives me the following print out when I run a t.test on it:               data:  dVals      t = 0.6452, df = 310, p-value = 0.7404     alternative hypothesis: true mean is less than 0      95 percent confidence interval:             -Inf 0.002134155      sample estimates:     mean of x          6e-04       How is R getting this t-stat? Am I computing t using the wrong formula? Also, I can't seem to get the correct p-value no matter what formula/table I try to look in... And as a last bit, I'm not understanding how a value of t=0.6452 leads to a true Ha since to accept Ha: $$ t < -t_\alpha \\ 0.6452 < -1.645,\space where\space\alpha=0.05 $$ Which is clearly not true... So shouldn't that mean the null hypothesis is not rejected?


title: Confusion related to prediction

text: I was referring to this video lecture http://videolectures.net/mlss09uk_murray_mcmc/ where the speaker had something like this $$ p(x|D) = \int P(x|\theta,D)P(\theta,D)d\theta $$ equivalent to $$ \frac{1}{S}\sum_{s=1}^S P(x|\theta^{(s)},D), \theta^{(s)} \sim P(\theta|D) $$ I didn't get how the average was derived and why is $\theta^{(s)} \sim P(\theta|D) $ Can anyone explain?


title: Finding the number of trials

text: If the probability of success is 0.09, how many trials are required to have a probability of at least one success as 1/3 or more?


title: Auto-Regressional & Moving Average Model Formula Properties

text: I seeking help in understanding specific values underlying the formula's for the MA(p) model & the AR(q) model. I am attempting to implement the models (building up to the combined ARIMA model) in the programming language Java. I do not come from an overly mathematical (I'm fairly new to statistics at least) background so be gentle. Here is the formula I am using for the AR(p) model: $$X_t - μ = β_1(X_{t-1} - μ) + ... + β_p(X_{t-p} - μ) + Z_t$$ Where $X$ is the time series, $μ$ is the mean of the time series, $β$ is the auto-correlation coefficient at a specific lag, $p$ is the order of the model and $Z$ is white noise of mean $0$ and variance $σ^2$. I'm fairly certain I have the above figured out, however the term "$Z_t$" confuses me. How would I implement this in code? I understand it is "random" however what are its ranges? Surely there must be a maximum and minimum of the term $Z_t$. Is it somehow based on the variance of the overall dataset? How is the "$Z$" value calculated on implementation exactly? Here is the formula I am using for the MA(q) model: $$X_t - μ = Z_t - θ_1(Z_{t-1}) - ... - θ_q(Z_{t-q})$$ Where $X$ is again the time series dataset, $μ$ is the mean of the dataset, $Z$ is white noise with mean 0 and variance $σ^2$, $θ$ is the correlation coefficient at a specific lag and $q$ is the order of the model. Again I the same issue as the above model in regards to the "$Z$" term. Also $θ$ is also the correlation coefficient of the dataset at various different lags, correct? Any help on this matter would be extremely welcomed and if you have any questions I would more than happy to answer them. Any use of examples alongside a full dataset i.e., X = (1,2,3,4,5,6,7) would be also extremely welcomed as it helps me understand the concept much more easily. Also please try to keep the explaination as idiot proof and contained as possible.


title: Deriving the distribution of the sum of censored variables

text: I want to be able to calculate the distribution of $$Y = \sum_{i=1}^n\max\{0,X_i\}$$ where the random variable $X_i\sim N(\mu_i,\sigma_i)$. Is the calculation of $f_Y(y)$ possible and if so what is it?


title: Significance test across multiple simulated experiments

text: First time question on this site, so please bear with me, thank you: I have 6 coin-flip-type experiments for which I can calculate 6 binomial p-values. I would now like to calculate the significance of observing at least 4 p-values < 0.05 in six experiments total. In one approach, I used Fisher's method (http://en.wikipedia.org/wiki/Fishers_method) to compound the p-values, but I wanted to add an additional test based on simulating the original coin-flip data. To this end, I performed random coin-flips (P=0.5) for each of the 6 experiments; the total number of flips differs between these 6 experiments but is irrelevant. I then count how many times (out of 100 simulations) the binomial p-value < 0.05. Simulating the original data 100 times, I get the following number of significant binomial p-values ("false positive hits") from these 6 experiments: 12, 13, 9, 10, 7, 11 Or divided by 100 (= frequency of false positive hits in simulated experiments): 0.12, 0.13, 0.09, 0.1, 0.07, 0.11 How do I calculate the probability that 4 or more out of these 6 would be positive given these frequencies? I realize that for calculating the probability that 6/6 are positive, I would simply multiply 0.12 x 0.13 x 0.09 x 0.1 x 0.07 x 0.11. But for 1-5/6 it's more complicated. I'm leaning towards a hypergeometric test, since I have to draw 6 times and I _think_ there's no replacement, but I want to double-check with you experts. Thank you!


title: Binary classifications problem where number of negative instances is much larger than number of positive instances

text: I am trying to write a logistic classifier that will classify text tokens on a web page. One problem I am running into is the fact that their are many more negative instances of the tokens I am trying to classify then positive ones. My question is: what is the ratio of positive training data to negative training I should be using? It seems that the more neg I tried googling the problem to look for relevant papers, but had no success.


title: How to interpret basic output from a regression analysis?

text: I have been trying to interpret the results below, but I am finding it difficult. I wonder if someone could help me. All answers highly appreciated. > Number of obs = 30   >  F( 2, 27) = 19.73   >  Prob > F = 0.0000   >  R-squared = 0.5937   >  Adj R-squared = 0.5637


title: Hot deck imputation: validity of double imputation and selection of deck variables for a regression

text: **Background:** I had a data set containing 212 observations with a lots of missing values. Most of the IVs and DVs are categorical (DVs are ordinal) in nature. There are 3 DVs and about 30 IVs. My intention was to run an ordinal logistic regression. A list-wise deletion keeps only 42 observations, so I decided to use hot deck imputation to fill in the missing values. I chose similar variables as the deck variables during the hot deck imputation (the deck variables should always be categorical and as far I know there should be a maximum of 5 deck variables). **Here are my queries:** **1)** When I imputed via hot deck once, 169 of the observations were filled in completely. If I use these imputed values for another hot deck imputation, then all 212 observations will fill in completely. But I am not sure if it is valid to use the imputed values for a further imputation. Can anyone suggest? **2)** Someone suggested me (from his experience) to use the 3 DVs as the background or deck variables for imputing all the DVs and all the IVs, because that will probably facilitate my regression results. May I know your comment about it? **3)** If I see almost all of the values (except from a very few) of a continuous IV are 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90 etc. then isn't it better to impute them via hot deck rather than via EM (as the hot deck will impute a variable with it's existing values only)?


title: Collapse dimension when comparing model coefficients

text: I have a dataset that includes 2 dimensions (A, B), each of which has the same two categories (1,2). I fit a model to the data in each of the 4 conditions and get 4 coefficients for a parameter of interest as a result including various stats about them (SE, CI, etc). I want to determine if dimensions A or B produce significant differences in the coefficient in conditions 1 and 2, similar to what a two-way anova would do. How can I group model coefficients together in this way so that I can compare them across each dimension?


title: How much RAM is needed for simulation studies using R?

text: I run several simulation studies with R using a MacBook Pro, which has 8 GB of RAM. Unfortunately, some of the simulations cannot be done due to limited memory. My question is how much RAM is needed for large simulations? (For example, each data set has 10000 subjects and needs to generate 1000 datasets). Is 64 GB of RAM big enough for large R simulations?


title: Three-level hierarchical regression using lmer

text: I would like to fit a 3-level hierarchical regression in lmer, however, I don't know how to specify the grouping factor above the second level. the model would be:               lmer(depedent ~ independent 1 + independent2 + (1|group1)....      And I would like to specify another group nested within `group1`. I've tried `(1|group1/group2)` but this gives an error message and group1:group2 is an interaction. I've also tried separately `(1|group1) + (1|group2)` but i'm not sure if this is correct.   thanks


title: When do you know when modeling something with a binomial distribution is better than a Poisson?

text: I was trying to understand when its better to model some random variable as being distributed by a Poisson Distribution or when being modeled as a Binomial Distribution. I was reading the following article: http://www.umass.edu/wsp/resources/poisson/ but was not sure if I appreciated it. What kind of structure does your data need, for it to be better approximated by a Poisson rather than a Binomial? What assumptions are necessary to separate when to model by Poisson or by Binomial?


title: Robust multivariate Gaussian fit in R

text: I need to fit a generalized Gaussian distribution to a 7-dim cloud of points containing quite a significant number of outliers with high leverage. Do you know any good R package for this job?


title: Measure of spread of a multivariate normal distribution

text: What is a good measure of spread for a multivariate normal distribution? I was thinking about using an average of the component standard deviations; perhaps the trace of the covariance matrix divided by the number of dimensions, or a version of that. Is that any good? Thanks


title: Deep learning and curse of dimensionalty

text: From "] Learning Deep Architectures for AI". see: http://www.iro.umontreal.ca/~lisa/pointeurs/TR1312.pdf (see section 3.1) "The local generalization issue is directly connected to the literature on the curse of dimensionality, but the results we cite show that what matters for generalization is not dimensionality, but instead the number of “variations” of the function we wish to obtain after learning. For example, if the function represented by the model is piecewise-constant (e.g. decision trees), then the question that matters is the number of pieces required to approximate properly the target function. There are connections between the number of variations and the input dimension: one can readily design families of target functions for which the number of variations is exponential in the input dimension, such as the parity function with d inputs." Anyone care to explicate? Or try to reword this paragraph. Or give examples. (And is he saying that decision trees are piecewise-constant?)


title: Graphing and Analyzing Probit Regression

text: I am currently analyzing a data set having to do with oak tree mortality. I am trying to understand the correlation between dead crowns (dead = `0`, live = `1`) and the distance to **(1)** nearest dead crown and **(2)** distance to nearest bay tree. (See Below). I have run a probit regression model:               Living<-sample(0:1,10,replace=T)         Dead_Dist<-c(158.68,99.62,64.42,64.42,86.92,117.77,41.81,41.81,54.73,64.35)         Bay_Dist<-c(92.47,179.92,317.73,365.23,58.70,193.23,330.36,123.14,88.65,72.34)         mydata<-cbind(Living,Dead_Dist,Bay_Dist)         mydata1<-data.frame(mydata)         myprobit <- glm(Living ~ Bay_Dist + Dead_Dist, family = binomial(link = "probit"),data = mydata1)         myprobit      But I am at a loss for how to plot and interpret the results, as I am fairly new to R and regressions.   If anyone has suggestions with how to proceed with this analysis, I would appreciate it!


title: How do I statistically compare two time series of different length?

text: If I have two completely different time series at equally sampled times, but the length of each of these series is different. How do I statistically compare them? Dynamic Time Warping? Also, as a follow up question. How do I go about this if each time series is of say 12 dimensions?


title: Paired, repeated-measures ANOVA or a mixed model?

text: I have been asked to analyse some data from a clinical trial looking a two methods of measuring blood pressure. I have data from 50 subjects, each with between 2 and 57 measures using each method. I'm wondering how best to proceed. Obviously I need a solution that will account for the fact the measure of blood pressure is paired (two methods measured contemporaneously) and also a time varying covariate (with a varying number of observations per patient) as well as account for intra- and inter-patient variablility. I was thinking of somehow shoe-horning this into repeated measures ANOVA, but I'm thinking it might need to be a mixed model approch. I'd appreciate any helpful advice you could offer. I'm a complete R newbie but very excited to develop skills and I have a moderate experience in Stata so could always fall back on that.


title: Leaf ordering for hierarchical clustering dendrogram

text: Assuming merging process was completed and we have the history of `n-1` merged clusters (merge two clusters `p` and `q` to cluster `c = min{p,q}`), we are to draw the dendrogram with `k` leaf nodes (or number of records, whichever is smaller). If `k` is smaller then number of records then leaves represent cluster objects I've seen some papers with optimal leaf ordering algorithms that run in `O(n^4)` time and other fancy algorithms. I've also heard about "left sliding" ordering. I'm looking for a simple ordering algorithm to find an admissible ordering of `k` leaf objects for dendrogram, given history of `n-1` merged pairs. Pseudo- code or link would be very nice.


title: How to estimate the centroid of clustered sequences?

text: I have run a sequence analaysis using the Optimal Matching algorithm. Afterwards, I have clustered the resulting distance matrice using the Ward algorithm and calculated silhouettes as measures of cluster quality and to identify representative sequences. Now, I am curious whether it is possible to estimate the sequences of the cluster centroids which, to my knowledege, must not be an original data point. How can I estimate the sequence of a centroid? To get an idea of the different steps of the analysis, consider this manual example[1]:               library(TraMineR)      library(WeightedCluster)      data(mvad)      mvad.alphabet <- c("employment", "FE", "HE", "joblessness", "school", "training")      mvad.labels <- c("Employment", "Further Education", "Higher Education", "Joblessness", "School", "Training")      mvad.scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")           ## Define sequence objects     mvad.seq <- seqdef(mvad[, 17:86], alphabet = mvad.alphabet, states = mvad.scodes, labels = mvad.labels, weights = mvad$weight, xtstep = 6)          ## Computing OM dissimilarities     mvad.dist <- seqdist(mvad.seq, method="HAM", sm="CONSTANT")          ## Clustering     wardCluster <- hclust(as.dist(mvad.dist), method = "ward", members = mvad$weight)     clust4 <- cutree(wardCluster, k = 4)          ## Silhouettes     sil <- wcSilhouetteObs(mvad.dist, clust4, weights = mvad$weight, measure = "ASWw")          ## Sequence index plots ordered by representativeness     seqIplot(mvad.seq, group = clust4, sortv = sil)      In this example, it would be for example interesting to see whether the sequence of third cluster's centroid differes from the most representative, original sequences in the cluster which are printed at the very top of the sequence index plot. In other cases, the centroid sequence may even have a more idealtype character which does not exist in the original dataset but reflects certain typical structures. [1] See for the example Studer, Matthias (2013). WeightedCluster Library Manual: A practical guide to creating typologies of trajectories in the social sciences with R. LIVES Working Papers, 24.


title: Predictive website analytics with irregular data points using R

text: I am trying to use R to build a time-series model to predict a variable which is only reported on irregular intervals. I have data for independent variables for all time periods (days), however the key dependent variable, reported_UV, is only reported monthly (every 30 or so days). I am looking for help building this model so that reported_UV can be predicted based upon previous cases of reported_UV as well as the independent variables. The dataset is at the following link in .csv format: https://docs.google.com/file/d/0B42_3jCtgxwYRVp5dkJ2NGVnaWs/edit?usp=sharing A description of the data fields is as follows:   * Month_Year The Month & Year (e.g. Jun-11)   * Date The Date (e.g. 6/1/2011)   * DayInMonth The day in month (e.g. for June 1st, this would be "1")   * Weekday Shows as "1" if the day is a weekday   * Weekend Shows as "1" if the day is a weekend   * Holiday Shows as "1" if the day is a holiday   * Month Shows the month in the year (e.g. for June, this would be "6")   * daily_UV The total unique visitors for the day   * p28_UV A lagging indicator showing the total unique visitors for past 28 days. This is a deduped version of the previous 28 days of daily_UV.   * reported_UV This is the key dependent variable to predict, reported by an independent auditor. This is reported on the last day of the month and shows the total unique visitors for the month. For example, the value on June 30th shows the entire month of June. All other values are missing.


title: Is Fleiss's $\kappa$ as conservative as Cohen's $\kappa$?

text: I know Cohen's $\kappa$ for inter-rater reliability is often regarded as a conservative estimation--e.g., Strijbos & al. below, who say > The higher the number of categories, the smaller the likelihood of chance > agreement as opposed to possible deviations. In other words, kappa tends to > be more strict (a comparatively larger chance agreement correction) in the > case of fewer categories. Does this same problem apply to Fleiss's $\kappa$? Strijbos, J.; Martens, R.; Prins, F.; Jochems, W. (2006). "Content analysis: What are they talking about?". Computers & Education 46: 29–48.


title: Intraclass correlation coefficient interpretation

text: I'm having a look at the intraclass correlation coefficient in SPSS. Data: 17 participants rated two lists of 9 & 7 items from 0 to 5 (0 being unimportant and 5 being very important). All participants rated all the items and the participants are a sample of a large population. The following output has been produced in SPSS. ![7 items case](http://i.stack.imgur.com/c132G.jpg) ![9 items case](http://i.stack.imgur.com/s6FcW.jpg) I am struggling to find anything online which deals with interpreting this, nor does any book interpret this in the level of detail I need. There is surprisingly little information/examples of the interpretation online on this, the literature is about choosing an intraclass correlation coefficient and not interpreting it. One problem I foresee here is the F test. The data is only 17 responses, which do not follow the normality assumption.


title: Is it necessary to transform an econometric model in order to have only diagonal elements in the error covariance matrix?

text: ## Model and its assumptions I'm working on the methodology of a two-way error component model. Here is the model: $$y_{jis} = x_{jis} \beta + \upsilon_{jis}$$ $j$ refers to school, $i$ refers to individual and $s$ to tested area / topic (mathematics or English, ...). $J$ is the number of schools and $S$ is the number of tested fields. $N$ is the number of students among all schools. $\upsilon_{jis}$ is the error term and can be decomposed as follows: $$\upsilon_{jis} = \theta_{j} + \phi_{js} + \varepsilon_{jis}$$ $\theta_{j}$ is the random effect for students attending school $j$, $\phi_{js}$ is the random effect for students attending school $j$ for a topic $s$ and $\varepsilon_{jis}$ is the traditional idiosyncratic error. In multivariate form it is: $$Y = X \beta + \upsilon$$ where $\upsilon = R \theta + F \phi + \varepsilon$ $R$ and $F$ are matrices that make it possible to correctly distribute their respective random effects. $R$ is size $= NS \times J$, and $F$ is size $= NS \times JS$. $\theta$ follows a multivariate normal with mean $0$ and covariance matrix $\tau I_{J}$, where $\tau$ is a scalar and $I_{J}$ is the identity matrix (size $J \times J$). $\phi$ follows a multivariate normal with mean $0$ and covariance matrix $\gamma I_{JS}$, where $\gamma$ is a scalar and $I_{JS}$ is the identity matrix (size: $J S \times J S$). $\varepsilon$ follows a multivariate normal with mean $0$ and covariance matrix $\sigma I_{NS}$, where $\sigma$ is a scalar and $I_{JS}$ is the identity matrix (size $N S \times N S$). ## Estimation of unknown parameters I have to estimate the different variances: $\tau$, $\gamma$ and $\sigma$ before running FGLS method. I work step by step to delete the different random effects. First I made a within regression to get rid of $\phi$ effect. Since $\phi$ is an effect for a given school, that within operator ($W$) will automatically make the $\theta$ effect vanishing. $\varepsilon$ is the only variance term left and I can use the methodology in Hayashi to calculate the expectation of $((We)' We | X)$ in order to correct with the right degrees of freedom. It is easy since the variance of that element is only diagonal; we only need to calculate the trace. Now, I only want to make the term $\theta$ disappear with another within operator ($W_{2}$). By construction, it deletes $\theta$, so there are $\phi$ and $\epsilon$ left. Once again I try to perform the expectation of $( (W_{2}e)' W_{2}e | X)$ where $e$ is now composed of $\phi$ and $\epsilon$. The variance of that vector is not diagonal, thus we don't have a beautiful formula for the estimated variance of $\phi$. I am carefully following the appropriate methodology, but I don't have a good-looking formula. I can also provide more details about the development of Hayashi that I reuse. **Is it necessary to transform the model in order to have only diagonal elements on variance-covariance matrix of errors?**


title: Assessing similarities in "call repertoire" for animal vocalizations

text: I have count data for number of call types over an interval of time for multiple animals. What I want is (1) a summary of association between call types (i.e., if X1 is large, X2 is small) and (2) a measure of similarity between the animals "repertoire" or rather pattern in which they call (comparison of correlation structure of animals' calls?). My data looks like the following:               Animal    Time    X1    X2    X3    X4    ...    Xn      A            1    0      5     3     1    ...     0     A            2    5      1     1     8    ...     6     ...................................................     ...................................................     ...................................................     C            1    3      3     5     1    ...     3     C            2    0      0     2     1    ...     5      Where `X1, X2, ... Xn` represent call type (i.e., 'Scream', "Howl', 'Laugh', etc). I'm not really sure where to start with this sort of data.


title: Machine learning for dynamic pricing

text: I haven't been able to find a good answer for this (though I know that there are a few related questions): How does Amazon determine how to price a product? How does it figure out when to downprice, and by how much? I'm sure that there's some machine learning involved here, but what does it involve? Doesn't have to be Amazon. I have metadata and data on past demand, but a regression model and random forest neither showed that anything other than overall sales are significant. I also have a natural experiment, but I don't want to base my whole case on one limited time period.


title: Propensity score matching using R

text: I have 50 control subjects and I would like to get 150 treatment subjects (i.e., a 1:3 ratio between the control and treatment groups) using propensity score matching and a few covariates. After performing the matching I would like to get indices of matched subjects (i.e., indices for both control and treatment groups) - this is my main goal. I would like to do this in R. I see that R has the Matchit and Matching packages. However, I am not sure what the difference is between the two and which one is more appropriate for my purpose. I am also not sure about the code and what output I need to be looking at to get indices of the matched control and treatment subjects. Also, could I use the `match()` command without the outcome (i.e., `Y`) variable? The experiment is not yet complete and thus I do not have the outcome variable. Would it be appropriate to use `match(Y=NULL, Tr=.....)`?


title: What is causing degrees of freedom and negative variance issues when estimating CFA with uncorrelated factors?

text: Ultimately I want to estimate a second-order factor structure in AMOS. As such I am trying to compare the fit of a first-order model in which the first-order factors are (a) constrained to be uncorrelated, and (b) free to correlate. I am running into two problems when I try to run these models.   1. the chi-square and df's for the fully constrained and unconstrained models are coming out as the same (I am trying to constrain the model by naming the covariance parameters under the plugins tab which I think constrains them to be the same?) It is my understanding that they should be different and if the unconstrained model (factors free to correlate) has a smaller chi-square value then a second-order factor structure can be estimated.    2. If I constrain individual correlations (setting the covariance of the curved arrows between the latent variables in the measurement model to be 0) in several instances I get an error message that one of the variances is negative and the chi-square is not estimated.  **Additional Information:** I have 3 factors with 2 observed variables each. My sample size is 202. The second order factor is not included in these models. I am looking to compare a constrained model (constrain factors to be uncorrelated) to the same model with unconstrained correlations (allowing factors to correlate) to determine whether I can estimate a second-order factor (2 models), then I am also constraining correlations individually (there are 3 as there are 3 factors). Can I constrain the correlations to be uncorrelated by setting the covariances between each of the 3 factors to be 1?


title: Using LSA for dimension reduction on test data

text: I am looking to explore using LSA as a dimension reduction technique for some textual data. To be specific I would like to take a matrix of M documents and N n-grams (ie: variables) and then create a smaller set of D variables that accounts for the most of the variance (thus resulting in a MxD matrix). The problem I am having is that as far as I can see lsa's `fold_in()` function only seems to return a full MxN matrix (albeit one that was estimated using D number of singular values). Here is an example of what I mean:               > library(lsa)     >      > #make training set     > td = tempfile()     > dir.create(td)     > write( c("dog", "cat", "mouse"), file=paste(td, "D1", sep="/") )     > write( c("ham", "mouse", "sushi"), file=paste(td, "D2", sep="/") )     > write( c("dog", "pet", "pet"), file=paste(td, "D3", sep="/") )     > write( c("dog", "pet", "bird","cat"), file=paste(td, "D4", sep="/") )     > train = textmatrix(td)     >       > #make test set     > td = tempfile()     > dir.create(td)     > write( c("dog", "mouse","bird"), file=paste(td, "D1", sep="/") )     > write( c("ham", "mouse", "sushi"), file=paste(td, "D2", sep="/") )     > write( c("dog", "pet", "pet","cat"), file=paste(td, "D3", sep="/") )     > test = textmatrix(td)     >      > #run LSA withing training data, reducing dimension to 2     > train.lsa = lsa(train, dims=2)     >      > #apply LSA to test data     > (test.lsa <- t(fold_in(test, train.lsa)))          docs       cat       dog      mouse         ham       sushi        pet       bird       D1 0.6477117 0.7120160 0.85915463  0.43561856  0.43561856  0.3527842 0.22417560       D2 0.4613586 0.2370914 1.03194173  0.59632318  0.59632318 -0.4227943 0.02574004       D3 0.5074166 0.9633718 0.04236298 -0.08829291 -0.08829291  1.2886710 0.37676074      As you can see the resulting test.lda object is 3x7 whereas I am looking to reduce the number of number of variables to 2. I could apply `prcomp()` to `as.textmatrix(train.lsa)` and then `predict.prcomp()` to test.lsa to reduce dimension after the fact, but that seems redundant given that SVD is already being used to run the LSA and that only 2 singular values are being used to create test.lsa to being with. I assume that I could also, after reviewing a bit of linear algebra, use the SVD components stored within the train.lsa object to do some manual matrix multiplication to get what I'm looking for but the lack of an easy to use function in the lsa package to do this makes me weary that I'm missing something here. So my question is what am I missing? Is there a simple function that I just can't find or am I misunderstanding how LSA is supposed to be used for "dimension reduction"? Also I apologize if this question is better suited for cross-validated - I'm kind of hoping this just a technical issue of me not seeing some function or parameter in the lsa package so I'm going with stackoverflow to start with.


title: Explanatory variables with mixture of effects

text: I have a continuous normally distributed outcome variable ($Y$) and $k$ explanatory variables ($X_i$, $i$ is $1 ... k$). The observations in $X$ are rare, but they behave in an additive fashion, so for all the $X$'s in each observation, I sum the observations up, and then I do linear regression of `Y ~ the summation of X`. The number of $Xi$ could be anywhere from 500 - 10,000. The rarity of non-zero observations in any $Xi$ is about 1 - 10. The problem is, that while $X$ should behave additively (if it influences $Y$ at all that is), it can behave in three different ways. Any one $X$ variable can decrease, increase, or have no effect on $Y$. The summation procedure won't take account of this, so is less powerful. What I would like to do is be able to use a weighted sum of the $X$ variables, with the weights determined by examining the relationship of each individual $X$ on $Y$ in an independent sample, but I can't do that. I don't have a lot of data ($n=430$, each individual effect is going to be pretty mild, and given that each observation for each $X$ is rare, power is at an all time low anyway). So my idea is that if I can somehow figure out which of the $X$'s are having a positive effect on $Y$, which are negative, and which are neutral, I can then basically weight them as $-1, 0, 1$, or something like that. Unfortunately, I am not a statistician, and while I would have loved for one to be consulted at the outset of the study, this is the data I have been asked to work with. Any ideas or advice or a jumping point for me to start Googling something, as I totally don't even know what to search for.


title: Autocorrelation in panel data

text: I have a random effects panel data model with 4 variables, 10 years of data and 946 observations describing different financial values for different companies. The panel is balanced. The model looks something like (written without subscripts t for time and j for company): $ln(A)=b0+b1(ln B) + b2 (ln C)...+u$ However, when I run my model in Eviews, the given Durbin-Watson statistic shows that I have positive autocorrelation in my model (D-W statistic value is lower than d-lower (dL). I have found out that there are numerous ways to correct your model for pure autocorrelation, but I have run into trouble when implementing any of them. What I have tried so far:   1. First, I tried Generalized Differencing (also known as Generalized Least Squares). I didn't find any information as to whether this was a suitable method for panel data and had some considerable trouble retaining the natural-logarithmic functions of my model (since some values changed into negative values and taking logs from negatives is impossible), but even after implementing this to my best of knowledge the autocorrelation was not eliminated. In fact, instead of positive autocorrelation I now had negative autocorrelation.   2. Secondly, I found that Newey-West standard errors can eliminate autocorrelation. However, Eviews doesn't support this function for panel data, which made me believe that, perhaps, this method does not work for panel data. Would it work?   3. The third solution I found from some more complicated papers. Basically, the authors of the theory (Arellano and Bond 1991) include a lagged dependant variable (model turns into a dynamic model) and then estimate a final model with Generalized Method of Moments. However, I couldn't find much information on how to create such a lagged variable (and which one of my variables should be lagged), nor how big the lag should be be and overall how to do this in practice.  My questions:   * How do I eliminate autocorrelation in panel data?   * Which methods really work for panel data? All of the ones mentioned above or none? Maybe there is some additional way to do this?   * How do I bring one of those methods into practice? Not everyone uses Eviews, but it is easy to see a formula, but much harder to explain to the computer how you want this formula solved. * * * Reference: > Arellano, M. and Bond, S. (1991) Some Tests of Specification for Panel Data: > Monte Carlo Evidence and an Application to Employment Equations. _Review of > Economic Studies, 58_ (2): 277-297 [DOI]


title: Generate random correlated variable from known $X$

text: I want to use Excel to generate a random correlated $Y$ from a known $X$. From another thread, I found the equation $Y = r\cdot X + E$, where $X$ is standardized and $E$ is a random variable from normal distribution having mean $0$ and $\sigma = \sqrt{(1-r^2)}$. I assume $r$ is the correlation coefficient found using Excel's `CORREL` function. I also assume I can calculate $E$ by using Excel's `NORMDIST` function. Are my assumptions correct? If I have a known $X$, how do I "standardize" $X$? Thanks for any help.


title: Survival analysis with categorical variable

text: I have event time data for subjects with different categories (A, B, C etc.) yearly observed. To my understanding my data is both right and interval censored (?). Subjects' category can change from year to year, e.g. A for years 0 to 2, than B from year 3 onwards for subject 1 etc. Each subject can be exposed to the "event" independent of its category. For time being I am not interested in transition effects. Which estimators can be suggested for the 1 year failure probability per category?


title: Dividing and forecasting a normal distribution

text: I have a collection of real-world samples that I am trying to model and then forecast. There are 10 datasets – each has an increasing number of data points (from 1 to 10). The points are equi-distant and each has some percentage of the total allocation for that set. For example, the first set has 1 data point with 100% allocation; second has 2 points with 50% each, etc. The sets appear to have a roughly normal distribution – but with a small skew to the right of the curve. As the number of data points increases with each set, the x-axis expands and the distribution is spread thinner over the set (ie. the mean decreases and the curve becomes more flat (lower kurtosis?)). I am using MathNet Numerics for C#. I have tried importing one of the datasets into a normal distribution and having the library calculate the mean, SD, skew, kurtosis, etc. Firstly, given those parameters for a normal distribution, is it possible to reproduce the values for a given number of data points? e.g. Can a normal distribution be divided equally on the x-axis and the percentage allocation for each point be calculated? Secondly, given the sample datasets, is there any way to model or find a function for the change between the datasets so that I may forecast the distribution for any number of data points? (e.g. given the shape of the distribution for 5 and then 10 datapoints, what would the distribution be for 15 or 20?). I am not a statistician, so please forgive me if I have asked something obvious or if anything doesn't make sense.


title: Standard error of a ratio

text: I have a linear regression and two estimates, say A and B and their standard errors. I need to find the standard error of a ratio A/B [or A/(1-B)]. I guess the main problem is that I don't know the correlation between A and B... I guess that makes this unsolvable, since the correlation plays a major role... correct? If I knew it, what would be the way to calculate the standard error? Thanks!


title: When to use the standard error on the mean

text: I have a number of measurements of a data point which I want to compare to my model. I have the sample mean, sample standard deviation and therefore the SEM. It usually seems fairly obvious to me whether to use the std. dev. or SEM. I think in this case I want to use the std. dev. as I'm comparing my model to my data, not estimating a mean value. However it made me wonder, if I had no model and I was simply reporting an experimental error are there any hard a fast rules as to which one would use. I have probably been rather inconsistent in the past. I tend to prefer the std. dev, as I generally use systems which can collect thousands of discrete data points a minute, so SEM tends to be rather implausibly small compared to the uncertainty present in the measurement system. Often this is to the point where the screen/printer resolution isn't sufficient to distinguish the upper and lower bounds using the SEM. If I work on an experiment where repeated measurements are hard to come by my mind seems to go to SEM. Is it simply a case of honesty? If stats are your dominating uncertainty quote SEM and if its measurment error quote std. dev.? I feel like I should not be running on auto-pilot on what seems like an important question.


title: HOG Feature Implementation with SVM in MATLAB

text: I would like to do classification based on HOG Features using SVM. I understand that HOG features is the combination of all the histograms in every cell (i.e. it becomes one aggregate histogram). I extract HOG features using MATLAB code in this page for the **Dalal-Triggs variant**. For example, I have grayscale image with size of 384 x 512 x 1, then I extracted the HOG features at 9 orientations and a cell-size of 8. By doing this, I get 48 x 64 x 36 features and when I check the size of these HOG features, it becomes 48 x 2304. How can I make this a histogram and use it toward a SVM classifier?


title: rpart and the printcp function

text: I don't really understand how the columns "xerror" and "rel error" are calculated. I found out that the printcp() function "gives cross-validation estimates of misclassication error (xerror), standard errors (xstd) of those estimates and the training (resubstitution) estimates (error)". Here is an example:               car <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/car/car.data", sep=",")     set.seed(2)     tree <- rpart(V7~.,data=car)     printcp(tree)     # OUTPUT #####################################     Classification tree:     rpart(formula = V7 ~ ., data = car)          Variables actually used in tree construction:     [1] V1 V2 V4 V5 V6          Root node error: 518/1728 = 0.29977          n= 1728                   CP nsplit rel error  xerror     xstd     1 0.129344      0   1.00000 1.00000 0.036767     2 0.115830      2   0.74131 0.85135 0.034987     3 0.040541      4   0.50965 0.50965 0.028872     4 0.030888      7   0.38803 0.38803 0.025729     5 0.027027      9   0.32625 0.35328 0.024694     6 0.023166     10   0.29923 0.30888 0.023261     7 0.017375     12   0.25290 0.26641 0.021754     8 0.015444     14   0.21815 0.23552 0.020557     9 0.010000     16   0.18726 0.18726 0.018472     ############################################     # Missclassification Error     sum(predict(tree, type="class")!=car$V7)/nrow(car)     [1] 0.05613426      Our Missclassification Error seems to be 0.05613426. In the printcp-output the estimate is 0.18726. How could this happen? I also tried 10-fold cross- validation to estimate the expected Missclassification Error and got 0.05497715 as Error. So how is this 0.18726 from the xerror and rel error column of the printcp output calculated?


title: Convert back standardized linear predictor

text: I guess this is a quite basic question, but I have been struggling with this for quite some time, so I hope someone can help me with this. I have a model (type is not relevant for now) which includes a linear predictor (LP) that would normally be calculated as:               LP = sum(coefficient1*predictor1, coefficient2*predictor2, etc.)      However, the predictors were standardized by centering them at the mean values and rescaling them using the standard deviation (SD). Hence, the LP now looks like this:               LP = sum(coefficient1*(predictor1-mean1)/SD1, coefficient2*(predictor2-mean2)/SD2, etc.)      So far it is clear to me, but now I would like to present my models in an efficient way, so instead of writing out the mean and SD for each individual predictor I want to provide only the coefficients (see first LP) and then standardize the LP as a whole. If the LP was only centered at the mean values writing it out would look like this:                LP = sum(coefficient1*predictor1, coefficient2*predictor2, etc.)- sum(coefficient1*mean1, coefficient2*mean, etc.)      Since the latter part of this formula has fixed values I can present that as one number, which makes the whole look somewhat less complicated. However, I don’t know how to also account for the SD rescaling. Does anyone know whether this (algebraically) possible and if yes how could I do that? Many thanks!


title: How to sample a country

text: Let be Statistan a country of approximately 20 million farmers. I want to draw a representative sample of the farmer population in order to survey them, in particular with the goal of capturing their yield (as defined by kilograms per cultivated hectare) and understanding what makes yield varying across the country. Statistan is a very diverse country.   1. Its cultivated land is 50% rainfed and 50% irrigated;    2. Environmental conditions change a lot through the country with 33% of the cultivated land in mountainous regions (very steep fields), 33% in tropical regions and 33% in arid regions;   3. Social conditions also differ with 50% of the farmers living in regions with of high per capita income and 50% of the farmers living in region of low per capita income. Should I take into consideration these categories when drawing the sample? Of course these are all important factors in determining each farmer's yield. It seems that if I sub-group my total population before sampling it I dramatically reduce the overall representativeness of the sample. But also if I don't I won't have many observations from less populated areas (like arid regions) where the potential of relevant findings is higher. How can I better understand the trade-off between sub-grouping and representativeness of the sample population given the number of observation I will collect with the survey?


title: Standardized scores correlation positive instead of negative

text: I had to standardize the scores of a variable because the items composing it were measured on different Likert-type scales (5, 6, 7 points). After this, I ran a correlation between this variable and other variables, with which theoretically it is supposed to correlate negatively. This yielded though a positive correlation. Correct me if I'm wrong, but is this result due to the score standardization? Thank you!


title: Estimating parameters using Kullback-Leibler or Kolmogorov-Smirnoff via Nelder-Mead

text: I want to find the parameters of a model which specifies a set of classification probabilities, for say M classes. (I'll use the parameters in another model later.) Given a set of parameters $\theta$, I can calculate its implied classification probabilities, $\pi(\theta)$, using Monte Carlo integration - $\pi$ is a smooth vector function of $\theta$. I would keep changing the parameters until I get the $\theta(\mathbf{p})$ which generates the classification probabilities I'm looking for - fixed $\mathbf{p}$. (I'll set the seed to control the randomness.) So, I want to use Nelder-Mead to minimise the "distance" between the $\mathbf{\pi(\theta)}$ and $\mathbf{p}$. I would like the distance to be Kullback-Leibler or Kolmogorov-Smirnoff as opposed to least squares for my fitness function. Which fitness function will lead to the solution faster - KS or KL? Also, given that I can evaluate the function to be optimized, but not its derivatives, is Nelder-Mead my best option.


title: Repeating random forest models to get highest prediction accuracy

text: I am using a random forest code to run one random forest model and distinguish which variables are important for classification and then to run a second random forest model using only these variables to reduce noise. The problem is that different variables are added to the model on different iterations of this process. I decided to repeat this process 10 times and use the model that produces the highest prediction accuracy (repeated measures sub-sampling from withheld data) as the model I report in my results. Is this a valid method?


title: Simple Bayes network

text: Given the following Bayes network: ![BayesNet](http://i.stack.imgur.com/O7H7b.png) with   * $p(k=t)=.2$   * $p(o=t)=.1$   * $p(s=t|k=f,o=f)=.0$   * $p(s=t|k=f,o=t)=.2$   * $p(s=t|k=t,o=f)=.5$   * $p(s=t|k=t,o=t)=.95$ how would I calculate $p(s=t|o=t)$ and $p(o=t|s=t)$? I tried the following way: $p(s=t|o=t) = \dfrac{p(s=t \land o=t)}{p(o=t)} = \dfrac{p(s=t)p( o=t | s=t)}{p(o=t)}$, which obviously doesn't make sense, since I have to caluclate $ p( o=t | s=t)$ with Bayes rule.


title: Meta-analysis with known confounder

text: I’m performing a meta-analysis in which the main outcome of interest is a correlation coefficient between two variables, $X$ (a psychological measure) and $Y$ (a biomarker). $Y$ is known to be strongly correlated with covariate $Z$ (age). I originally intended to do an individual patient data meta-analysis (IPD), which would allow approaches like stratification by age or propensity scores. However, due to difficulties in obtaining raw data from authors and differences in the scale used to measure $X$ across studies, the IPD approach is unlikely to work. Instead, I am using traditional methods to meta-analyze effect sizes. As a first step, I’ve meta-analyzed the raw $r_{XY}$ and then the $r_{XY \cdot Z}$ (partialling out $Z$). But for reasons described in an earlier question, I don't want to rely heavily on the latter approach to “controlling for” the age covariate. I am also considering taking the 6 studies for which we have raw data (of 10 total) and performing a stratified meta-analysis on the raw $r_{XY}$. I would define, say, 4 age groups, compute effect sizes within each age group for each study, and conduct a separate meta-analysis within each stratum. Another similar approach, which could incorporate some of the studies without raw data, would be to meta-regress on mean age in the study. I strongly suspect both of these approaches might be badly underpowered. Given the limitations of the data we have, is there a better approach to assessing the extent to which the correlation between $X$ and $Y$ is confounded by covariate $Z$?


title: Odds of winning one of $p$ prizes with $t$ tickets entered

text: How do you calculate the odds of winning one of 10 prizes if you have $t$ tickets entered in a drawing with $T$ total tickets? Each person can only get one prize. A person with 20 tickets entered obviously has higher odds than someone with 1 ticket. There are 50 people with tickets in the jar. So if I had 2 tickets out of 500. I would have a 2/500 chance of getting the first ticket drawn. The odds of getting the second ticket would be 2/(500-however many tickets the first person drawn had). So I was thinking that the average person has 500/50 or 10 tickets each. So 2/500 chance for the first ticket drawn. Then 2/490 the second, 2/480 for the third and so on for 10 tickets drawn since you can't know who actually gets drawn first. Someone may get drawn first that has 50 tickets and then you have better odds for the next 9 because more tickets were removed. This is just a random example to try to figure out the idea of solving a similar problem.


title: What are some classification approaches for real-valued features that produce actual probabilities?

text: I am looking for a classification algorithm that emits probabilities for each label and supports real-valued features. 1) From what I gather, logistic regression may output a variable in $(0, 1)$, however, it does not describe the true classification probability. 2) Naive Bayes supports probabilities, however, it uses discrete features. I need this to compute expected utilities. There is no way to couple utilities with probabilities in the algorithm, because the features are not related to the utilities (they are external).


title: Bayesian analysis- Normal distribution with unknown mean and variance

text: Suppose we have some data points that we believe they follow $N(\mu, \sigma^2)$, and both parameters are unknown . I want to assign conjugate prior distributions on both $\mu$ and $\sigma^2$, and then obtain the posterior distribution. How this can be done in R?


title: How to test relative validity statistically?

text: I am trying to understand a) what relative validity is, and b) how to test it statistically. So far I have gathered, from sources such as this that relative validity is a larger reliability in causation associated with certain stimuli, given a set of stimuli. It seems that some statisticians test relative validity using regression models, such as in this article. However, in none of these sources it is explained _how_ regression (or any other statistical tool) is used to measure relative validity. Hence my question - how do you test relative validity statistically?


title: Coursera Neural Nets

text: If I go to https://www.coursera.org/course/neuralnets and click preview, there are many video lectures. Does anyone know if these are the same lectures you see if you enroll in a course?


title: Parameter estimation with GMM

text: I have estimated the parameters of normal distribution with GMM and got the follwing results : mean = -0.01168 , p-value = 0.83519 i'm bit confused in interpreting the result. can i say that the mean is zero as it is not statistically significant. thanks


title: Why does poisson regression need to assume observations are poisson distributed?

text: Zuur (2013) 'Beginners Guide to GLM and GLMM' states that if the Pearson residuals, when plotted against fitted values from a poisson regression, show the pattern below then the assumption of poisson distributed observations is violated. Zuur urges to try a negative binomial regression if the residuals show a pattern like below. If a poisson regression is used, then Zuur states the standard errors will be wrong. But why does the model need to make assumptions about the distribution of observations around the fitted values? Why can't it use information contained in the distribution of the observed residuals to calculate the standard errors? My best guess is that for every fitted values, there might only be a few observed values, and there is not enough information contained in these few observed values to calculate the standard errors. ![enter image description here](http://i.stack.imgur.com/V9ppj.jpg)


title: Multinomial logistic regression interpretation in R

text: I have run a multinomial logistic regression test for the interaction between species of deer, days a camera trap was in the field and type of reaction. The model with the best AIC value was:               Coefficients:        (Intercept) speciesmuntjac   speciesroe speciessika        days     r   -0.7471023      0.6263753 -0.005967869 -0.74253017 -0.05189515     sr   0.6909319      0.5552278 -0.355611180 -0.01622306 -0.03178001          Residual Deviance: 971.6464      AIC: 991.6464       But I also got the following                   (Intercept) speciesmuntjac speciesroe speciessika       days     r  0.0071741793     0.17599897  0.9865350   0.2095687 0.08669276     sr 0.0001402185     0.08536257  0.1331651   0.9574258 0.08829020      I'm just wondering if anyone would know what the interaction between species and days is in relation to r (reaction) and sr (strong reaction)? I know I can't reject the null hypothesis that there is no effect of species or days on r or sr but beyond that I'm lost!


What is the difference between statistical significance (based on probability) and an effect size? I am fairly new to stats and finding a bit of an up hill marathon at the moment. I have been asked to research the question posted above, yet I have exhausted my resources and I am still non the wiser on how to accurately answer the question.. I really do hope you guys/girls can help me out of this slowly deepening hole!


title: Parameter estimation of a power spectrum equal to a power law + white noise

text: Given $X_t$ a multivariate random gaussian variable of covariance matrix $N_{tt'}$ diagonal in Fourier space (sampling is equally spaced), I would like to parametrise its power spectrum as: $S_X(f) = \frac{a}{f^\alpha} + b$, for $f\in[f_{min},f_{max}]$ and estimate $a$, $b$ and $\alpha$. Any idea how to get these estimators (unbiased if possible), the goal being to use these parameters (and extrapolate to all f) to describe the covariance matrix? There are plenty of references on how to estimate the slope of a power-law distribution (see review by Clauset et al. 2007, SIAM Proceedings), but in the case I am interested in, the starting point is $X_t$, not $S_X(f)$


title: ANCOVA to check effect of covariate

text: I want to check whether a group of patients are significantly different from their control group. However, I also want to check if the p-value is still significant in case a covariate is taken into account. For that I created a simple data frame in R:               data <- data.frame(group = c(rep("CTRL", 10), rep("P", 10)),          response = c(10,11,14,16,17,17,19,20,21,22, 10,11,11,11,12,13,14,14,15,16),         age = c(40,41,45,43,50,51,55,57,60,62, 40,42,43,43,45,46,46,50,52,54))      First of all, I performed a normal ANOVA using the `lm` function. I tried 3 different ways to see if the results are the same. And they are:               anova(lm(data$response ~ data$group))          summary(lm(data$response ~ data$group))          summary(aov(data$response ~ data$group)      To check if the covariate `age` is correlated with the `response` variable I performed a correlation test:               cor.test(data$response, data$age)      Seeing a high and significant correlation I concluded that the significance between the control group and the patient group might be due to the effect of age. I am now unsure how to perform the ANCOVA analysis to check if the effect between the two groups is really there or if it is just because of the covariate age. To check this I did the following:               m1 <- lm(data$response ~ data$group + data$age)          m2 <- lm(data$response ~ data$group)          anova(m1, m2)      Comparing the two models resulted in a significant p-value. I would have thought of a p-value of much less significant then the one obtained with the normal anova. Is it actually right to compare the two models or can I achieve this by using only one model? I am really stuck here and hoping to find some helpful answers here.


title: Getting started with projection analysis

text: For one of my grad school classes I am working on a projection analysis. The data includes projected year-end balances. These balances are projected every month of the fiscal year (starting in July and ending in June) to see how much approximately is the Y-T-D amount going to be. I'd like to analyze the projected balances vs. the actual Y-T-D balance.   I'd like to find a way to estimate at which month we are comfortable at saying that this is approximately what the Y-T-D amount is going to be. I thought of using confidence interval and normal distribution, but I can't figure out how to approach this. Thank you.


title: Experimentially proving that random variables are independent

text: I just had a very simple (possibly very dumb) question. How would I, experimentally prove that two random variables drawn either from the same population or two different populations are independent (or for that matter of fact dependent). Thanks.


title: Fuzzy regression discontinuity design and exclusion restriction

text: In a fuzzy regression discontinuity design, what does the exclusion restriction look like in terms of a conditional expectation between the instrument in the first stage and the error term in the structural equation? In IV, we normally say that the first-stage must produce a non-zero estimate (and preferably one that is sufficiently statistically significant to avoid weak instruments). However, in RD, we modify that to say that the continuity of the first-stage must be non-zero. That is, we evaluate it in the limit. The same for the errors. In the limit, the difference in the conditional expectation coming from the left relative to coming from the right around the discontinuity must equal zero? How about the exclusion restriction? Is that also something that must be a continuous condition in that its difference is evaluated from both the left and right in the limit, and the resulting difference should equal zero?


title: ID3 and C4.5: how does “gain ratio” normalize “gain”?

text: The ID3 algorithm uses "Information Gain" measure. The C4.5 uses "Gain Ratio" measure which is Information Gain divided by SplitInfo, whereas SplitInfo is high for a split where records split evenly between different outcomes and low otherwise. My question is: How does this help to solve the problem that Information Gain is biased towards splits with many outcomes? I can't see the reason. SplitInfo doesn't even take into account the number of outcomes, just the distribution of records in the split. It may very well be that there is a low number of outcomes (say 2), and the records are split evenly between those 2 outcomes. In that case, SplitInfo is high, Gain Ratio is low, and a split with few outcomes is less likely to be chosen by C4.5. On the other hand, it may be that there is a low number of outcomes, but the distribution is far from even. In that case, SplitInfo is low, Gain Ratio is high, and a split with many outcomes is more likely to be chosen. What am I missing?


title: PCA weights vs. standardized

text: I have a number of timeseries on which I want to apply a PCA (using matlab). These time series have very different variances. My objective is: 1) all time series have the same weight for me, no one should get overweighted 2) I don't want the nature of my data (statistical properties) to get modified through transformations I have found two solutions on the internet: A) standardizing the time series B) using the pca command in matlab, while specifying that the weights should be the inverse of the variances Both approaches look ok to me, what are the advantages or disadvantages of A) vs. B)? Is there anything I should be aware of or are both approaches perfect substitutes? Thanks!


title: Properties of logistic regressions

text: We're working with some logistic regressions and we have realized that the average estimated probability always equals the proportion of ones in the sample; that is, the average of fitted values equals the average of the sample. Can anybody explain me the reason or give me a reference where I can find this demonstration?


title: Discrepancy between number of expected and observed events in logistic regression

text: > **Possible Duplicate:**   >  Properties of logistic regressions I have built a logistic regression model with two covariates. Let $y_i$ be the binary responses and $\hat{\pi}_i$ the ML estimated probabilities for $i=1,\ldots,n$. Now it happens that $\sum_i{y_i}\neq\sum_i{\hat{\pi}_i}$. So, computing a chi-square goodness of fit test, sum of expected frequencies is different from sum of observed frequencies. Is this normal, or there is something wrong?


title: Nelson-Aalen estimator and interval censoring/tied events

text: I just noticed that the Nelson-Aalen estimate of the cumulative hazard changes depends on the existence of tied events. As a toy example, consider a study with 3 patients. The patients die on day 1, 5, and 14 of the study. If we check the patients status daily, the Nelson-Aalen estimate for the cumulative hazard at the end of the study is $\frac{1}{3} + \frac{1}{2} + 1 = \frac{11}{6}$. If we check weekly (on day 7 and 14), the final estimate is $\frac{2}{3} + 1=\frac{10}{6}$. So the final estimate of the cumulative hazard differs, and I don't see why this is reasonable.


title: 2D object recognition using MATLAB

text: Have you any idea about implementing 2D object recognition with MATLAB? Which characteristics of objects can feed a neural network?   It's my training data-set (provided by ETH University of Switzerland)   What is the start point?   There are 5 view of each object, and whole objects are 66. ![enter image description here](http://i.stack.imgur.com/Yt6hN.jpg)


title: Variance based on given frequencies using SPSS

text: I'd be too happy, if someone could post a code snippet, which explains how to compute mean and variance of a set of records, which contain frequencies? Suppose we have records like **(FORMAT F)**               - GroupA, 6 x Grade 1, 5 x Grade 2, 10 x Grade 3     - GroupB, 2 x Grade 1, 7 x Grade 2, 18 x Grade 3     - GroupA, 23 x Grade 1, 5 x Grade 2, 1 x Grade 3     - GroupA, 10 x Grade 1, 15 x Grade 2, 12 x Grade 3      Goal: Compute for each Group's mean and variance of grades using SPSS syntax. Clarification: Using SPS, processing cases **(FORMAT C)** would be simple (as colleagues told me):               # Record number 1 (line 1 of the above data)     GroupA, Grade 1        ... 6 times     GroupA, Grade 1     GroupA, Grade 2        ... 5 times     GroupA, Grade 2     GroupA, Grade 3        ... 10 times     GroupA, Grade 3     # Record number 2 (line 2 of the above data)      GroupB, Grade 1     GroupB, Grade 1     GroupB, Grade 2        ... 78 times     GroupB, Grade 7        ... and so on      Unfortunately, we don't have cases of FORMAT C, but already aggregated data as in FORMAT F. I probably should add this: I implemented a solution, which uses (user entered) FORMAT F data using PHP/SQL: The PHP program constructs an SQL query, which computes the mean, the variance and the standard deviation of each Group. I even tested the solution using well-known sample data. The results are identical. Thus I'm sure the solution is correct. Since the computed statistical properties are sensible, I'd like to get a verification using a standard statistical package using the real FORMAT F data (not the test case data). My colleague is going to use SPSS.


title: Inference on linear terms in a GAM in the presence of smoothly modelled adjustment variables

text: In a multivariable generalized additive model with Gaussian family and identity link function, in which some explanatory variables are modelled linearly and some are modelled as smooth functions, to what extent is inference on the linear terms analogous to that on ordinary least squares coefficients? Based on my limited understanding of the GAM and of the linear model, and based on my experiences, it appears that they are practically the same.


title: Y v level for different additives; how to handle zero level

text: I want to model the effect of Level (continuous) of different Additives (categorical) on Y (continuous) by fitting a linear model, so my data is going to look like this:               Additive   Level   Y     A            1     2.3     A            2     4.5     A            3     6.7      <snip>     B            1     8.9     B            2     1.0      etc.      I could fit that with something like `Y ~ Level %in% Additive` But what about zero level? At zero level additive identity doesn't mean anything. I don't want to include another variable indicating whether _any_ additive is present because that would split off the zero level data. I want the zero level data to be used in estimating the regressions for _all_ the Additives, but I don't know how to enter the data, never mind how to specify the model! It feels as if this must be a recognised situation with "standard" approaches. I've hit it before but never resolved it to my own satisfaction, now it's arisen again. Can anyone point me in the right direction?


title: Calculating conditional probability

text: I've got the next question. Let's consider I have the pair of distributions: X|t ~ Binomial(n,t); t~Beta(a,b). Here n,a,b are known. I need to construct conditional probability to sample from it, t|X. In this example I can do this analytically, t|X ~ Beta(x+a,n-x+b) But how can I do this in general case, for any distributions? Numerically or somehow else.


title: Law of total expectation and how prove that two variables are independent

text: I am confused in understanding the formal definition for independence and how does this differ from the law of total expectations. I thought that $E[u|x]=E[u]$ formally means independence of $u$ and $x$. However, when I have $E[u|x]=a$ (const) and use the law of total expectations, $E[u]=E[E[u|x]]=a$ I obtain $E[u]=a$ Don't I have then $a=E[u|x]$ and $a=E[u] \implies a=E[u|x]=E[u] \implies \text{independence between x and u}$? Something does not add up, as this would mean that we can always establish independence between two variables. Where is my error in thinking?


title: Logistic regression with categorical data

text: I'm trying to apply logistic regression to the data with binary predictor. But some of my variables are numerical and some are categorical. If I just do this in R I get the model where for every categorical variable I have coefficients and p-values for for variable's possible values except for the first. How can I interprete such model? And what is the best way of finding best model for such problem? Regards, Oleg


title: How to score detector

text: I'm developing a classifier system to detect objects of interest in images. I want to report a score, but I'm a bit lost as to what the most fair and informative number is. Sensitivity and specificity seem to be commonly used, but while I have numbers for true positive, false positives and false negatives, I'm not sure how to properly calculate 'true negative'.


title: Looking for a proper minimization tool

text: I'm given a $M\times N$ ($M$ rows, $N$ cols) matrix $X$ and a $1\times M$ column vector $y$ ($M>N$). I need to find such a vector $x$ that minimizes the value of s, where $s= \sum_{i=1}^{M}\left [(y_i - \hat{y}_i)^2 \right] $ and where $\hat{y} = X \times x$ Basicaly, this is a set of linear constraints. There is also a constraint that for every $x_{i, j} \in X$ , the following is true: $0 \le x_{i,j} \lt 1$ It is known that $X$ is sparse (most of the elements are 0's). The typical number of rows in $X$ is 200 and the typical number of columns is 150 Is Simplex algorithm or gradient-based optimization methods well-suited to perform this job or should I try using global heuristic methods such as Genetic algorithm, taboo search etc? What Python tools would you suggest to attack this problem? Last, but not least: is there any other SE forum that is better suited for this questoin?


title: Analysing ratios of variables

text: Ratios (X/Y, e.g. body mass index) are variables with odd distributions. They have no means or moments (e.g. variances, skewness or kurtosis). Thus, my questions are:   1. How to compare between (or among) 2 or more groups with ratio variables?   2. How to use ratio variables as the independent variable or dependent variable in a regression analyses?   3. If:      * X/Y is an independent variable: is it better to use X and Y (or 1/Y) as independent variables separately?     * X/Y is a dependent variable: is it better to use X as a dependent variable and Y (or 1/Y) as an independent variable?   4. Are most studies on ratios statistically wrong?   5. Is the method described by a recent article Error propagation in calculated ratios. Clinical biochemistry, Vol. 40, No. 9-10. (June 2007), pp. 728-734) the way to a solution?


title: What distribution would lead to this highly peaked and skewed density plot?

text: The density of my data set was plotted in R as follows. ![enter image description here](http://i.stack.imgur.com/vUUFp.png) What kind of distribution would fit this data? As I am not experienced to tell by visualization, I can only guess it is not normally distributed. I shall test it with R: some guidance as a starting point is highly desirable.


title: Intuition behind why Stein's paradox only applies in dimensions $\ge 3$

text: Stein's Example shows that the maximum likelihood estimate of $n$ normally distributed variables with means $\mu_1,\ldots,\mu_n$ and variances $1$ is inadmissible (under a square loss function) iff $n\ge 3$. For a neat proof, see the first chapter of Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction by Bradley Effron. This was highly surprising to me at first, but there is some intuition behind why one might expect the standard estimate to be inadmissible (most notably, if $x\in N(\mu,1)$, then $\mathbb{E}||x||^2\approx ||\mu||^2+n$, as outlined in Stein's original paper, linked to below). My question is rather: What property of $n$-dimensional space (for $n\ge 3$) does $\mathbb{R}^2$ lack which facilitates Stein's example? Possible answers could be about the curvature of the $n$-sphere, or something completely different. In other words, why is the MLE admissible in $\mathbb{R}^2$? Edit 1: In response to mpiktas concern about 1.31 following from 1.30: (1) $E_\mu(||z-\hat{\mu}||^2)=E_\mu(S(\frac{N-2}{S})^2)=E_\mu(\frac{(N-2)^2}{S})$. $\hat{\mu_i} = (1-\frac{N-2}{S})z_i$ so $E_\mu(\frac{\partial\hat{\mu_i}}{\partial z_i} )=E_\mu( 1-\frac{N-2}{S}+2\frac{z_i^2}{S^2})$. Therefore we have: (2) Hence $2\sum_{i=1}^N E_\mu(\frac{\partial\hat{\mu_i}}{\partial z_i} )=2N-2E_\mu(\frac{N(N-2)}{S})+4E_\mu(\frac{(N-2)}{S})=2N-E_\mu\frac{2(N-2)^2}{S}$ Edit 2: In this paper, Stein proves that the MLE is admissible for $N=2$. (I added tags as well as I could, but surely someone could come up with better tags. If you can and have the permissions, please do add them)


title: Random forest highly imbalanced dataset: how to test _ create ROC curve?

text: I have a dataset containing 10,000 examples with 8 features. I would like to create a random forest to classify this dataset in two classes, more specific: substrate / no substrate. In this dataset there are approx 20 TRUE positives that I got from literature. I would like to predict other potential substrates. I created my negative data set (false) by resampling of the remaining examples/proteins. To overcome issues with imbalance, I used a down- sampling approach to create my random forest with the caret package using strata and class site parameters. Now I would like to compare this against random forest without down sampling by creating a ROC curve analysis. My question is now how to create my testing dataset. I could just take my TRUE positives and add 10,000 negatives (again by resampling) and run the random forest on it and create the ROC analysis. Would this not yield a highly optimistic result? I know it would be better to split my 20 true positive candidates, lets say in 15 to 5, and 5 would be for the testing set, but I'm afraid to lose sensitivity caused by so little true pos to generate random forest. On the other hand, I feel okay with my proposed testing set since both algorithms are just tested against each other. Thanks a lot.


title: Calculating the linear model with R

text: I need to calculate the linear model in R, i did the following: > summary(model) But what if I wanted to calculate only the first point? A bit stuck with this one... Many thanks! Here is the code used to create plot itself:               > elevation=c(12, 34, 32, 12, 11, 14, 56, 75, 43)      > snowfall=c(6, 52, 41, 25, 22, 9, 43, 67, 32)     > snowfallElevationPlot = data.frame(elevation=elevation, snowfall=snowfall)     > plot(elevation, snowfall)          > model <-lm(snowfall~elevation, data=snowfallElevationPlot)     > abline(model)


title: Correlation between two repeated measurements with several measurements at each repeat

text: It is fairly straightforward to calculate the variance of the paired t-test: $Var(\overline{D})=Var(\overline{Y_1})+Var(\overline{Y_2})-2*Covar[Var(\overline{Y_1})*Var(\overline{Y_2})]$ But how can you calculate this variance if you have at each of the two measurements a number of measurements per patient (to exclude measurement error). A real life example: a patient his fat is weighted by 4 doctors at the start of the trial, and by 4 different doctors at the ten of the trial. Thank you!


title: Fit poisson regression

text: I have a data set $Y$, $Y$ is assumed to follow a poisson distribution with mean = $\lambda$ . However, each element $y_i$ in $Y$ is accompanied with a covariate $x_i$. So here I want to use poisson regression to model this dataset. The link function is assumed to be log($\lambda(i)$)=$\beta_0$+$\beta_1x_i$. I then want to use glm(Y~X,family='poisson') to fit the model and get $\beta_0$ and $\beta_1$. After get $\beta_0$ and $\beta_1$, can I calculate the probability mass function given new $y$ and $x$. I want to know if all of my above thought make sense? Thanks


title: Time Series Forecasting with Daily Data: ARIMA with regressor

text: I'm using a daily time series of sales data that contains about 2 years of daily data points. Based on some of the online-tutorials / examples I tried to identify the seasonality in the data. It seems that there is a weekly, monthly and probably a yearly periodicity / seasonality. For example, there are paydays, particularly on 1st payday of the month effect that lasts for few days during the week. There are also some specific Holiday effects, clearly identifiable by noting the observations. Equipped with some of these observations, I tried the following:   1. ARIMA (with `Arima` and `auto.arima` from R-forecast package), using regressor (and other default values needed in the function). The regressor I created is basically a matrix of 0/1 values:     * 11 month (n-1) variables     * 12 holiday variables     * Could not figure out the payday part...since it's little more complicated effect than I thought. The payday effect works differently, depending on the weekday of the 1st of month. I used 7 (i.e., weekly frequency) to model the time series. I tried the test - forecasting 7 days at a time. The results are reasonable: average accuracy for a forecast of 11 weeks comes to weekly avg RMSE to 5%.   2. TBATS model (from R-forecast package) - using multiple seasonality (7, 30.4375, 365.25) and obviously no regressor. The accuracy is surprisingly better than the ARIMA model at weekly avg RMSE 3.5% . In this case, the model without ARMA errors perform slightly better. Now If I apply the coefficients for just the Holiday Effects from the ARIMA model described in #1, to the results of the TBATS model the weekly avg RMSE improves to 2.95% Now without having much background or knowledge on the underlying theories of these models, I'm in a dilemma whether this TBATS approach is even a valid one. Even though it's improving the RMSE significantly in the 11 weeks test, I'm wondering whether it can sustain this accuracy in the future. Or even if applying Holiday effects from ARIMA to the TBATS result is justifiable. Any thoughts from any / all the contributors will be highly appreciated. Link for Test Data Note: Do "Save Link As", to download the file.


title: Generation of random numbers from the multivariate normal distribution in spherical coordinates

text: Can any one point me to a research article or post the actual algorithm to generate random numbers from the multivariate normal distribution in spherical coordinates? I need both $N_p(0,I_p)$ , the standard multivariate normal case in spherical coordinates and $N_p(0,\Sigma)$ where $\Sigma \neq I_p$


title: How to measure weight similarity?

text: I'm doing some machine learning and get a set of optimum weights in the end. I'd like to verify that these weights are by and large the same no matter how many times I train on the data. I assume that my cost function has a global minimum. How can I measure the similarity between two set of weights that I got from training two times?


title: Luce choice axiom, question about conditional probability

text: I'm reading Luce (1959). Then I found this statement: > When a person chooses among alternatives, very often their responses appear > to be governed by probabilities that are conditioned on the choice set. But > ordinary probability theory with its standard definition of conditional > probability does not seem to be quite what is needed. An example illustrates > the difficulty. When deciding how to travel from home to another city, your > choice may be by airplane (a), bus (b), or car (c). Let A, B, C denote the > uncertain states of nature associated with form of travel. Note that if one > elects c all of the uncertainties of A and B remain because planes fly and > buses run whether or not you are on them. However, if you elect either a or > b , then your car remains in the garage and the set C is radically altered > from when the car is driven. So there really is no universal event > underlying the sources of uncertainty. > > The choice axiom of chapter 1 was introduced as a first attempt to construct > a probability-like theory of choice that by-passed the fixed, universal > sample space assumption. source: http://www.scholarpedia.org/article/Luce's_choice_axiom For me the probability measure is defined with the triplet $\Omega$, the sample space, a sigma-algebra $\mathcal{F}$ and finally a measure $P$. With respect to the foregoing example what seems to be the problem if I define: $\Omega = \{ \text{bus}, \text{car}, \text{airplane} \}$ One crucial assumption in common statistics is the ceteris paribus condition. Is this the reason we need to adjust basic probability theory in the context of choice behavior because the c.p. assumption is violated?


title: When is Gaussian mixture model more effective than Single Gaussian?

text: From my experiment, given training samples of around 600 each having vector length of 30, the single gaussian component performed much better unexpectedly. Having said that, I want to know how or when it is better to use mixture of Gaussian (any rule of thumb or whatever)


title: Motivations for Shi-Malik Algorithm

text: So I've been trying to make sense of the clustering algorithm on page 6 of this paper. ![Algorithm](http://i.stack.imgur.com/wE9PB.png) Are the "first" k eigenvalues they refer to the smallest eigenvalues? What are the $y_i$ exactly? I don't see the motivation for using them. If anyone could link other literature or explain how/why the first eigenvectors are being used, that would be great.


title: Rules to apply Monte Carlo simulation of p-values for chi-squared test

text: I'd like to understand the use of Monte Carlo simulation in the `chisq.test()` function in R. I have a qualitative variable which has 128 levels / classes. My sample size is 26 (I was not able to sample more "individuals"). So obviously, I will have some levels with 0 "individuals". But the fact is that I have only a very small number of classes represented out of the 127 possible. As I have heard that to apply chi-squared test we should have at least 5 individuals in each level (I do not completely understand the reason for that), I thought I had to use the `simulate.p.value` option to use Monte Carlo simulation to estimate the distribution and compute a p-value. Without Monte Carlo simulation, R gives me a p-value `< 1e-16`. With Monte Carlo simulation, it gives me a p-value at `4e-5`. I tried to compute the p-value with a vector of 26 ones and 101 zeros, and with Monte-Carlo simulation, I get a p-value at 1. Is it OK to state that, even if my sample size is small compared with the number of possible classes, the observed distribution is such that it is very unlikely that all possible classes exist at the same probability (1/127) in the real population?


title: Inference of the parameters of a linear congruential generator

text: Is there an R package or other open source software that solves the problem of inferring the parameters of a linear congruential generator? Such inference can be done as e.g. described here in a strict sense. On ther other hand, I would also be interested in statistical methods that can be applied in this context.


title: Fractional Polynomial interpreting interaction terms

text: **Background:** I have developed a logistic regression model where I am trying to analyze the effect of socio-economic data of a family on their probability of receiving a home loan. I have used fractional polynomials to achieve linearity in my continuous covariates. One such variable is "disposable income" which comes significant in my main effects model. On testing for interactions, the transformed disposable income (Idisp__1) 'interacted' with banksavings (not transformed through fractional polynomials) and interactions came to be significant. **Issue:** I am using a confidence interval graph as suggested by Hosmer and Lemeshow in _Applied Logistic Regression_ to analyze the odds ratio of the interaction together with the two variables (i.e., dispincome$~\times~$banksavings, dispincome and banksavings). When I keep disposable income constant and increase my banksavings by one unit I get significant confidence intervals ( i.e., 1 does not show up on the different confidence interval lines) and I get intuitive results. However when I keep my banksavings constant next time and now decide to increase my disposable income by 1 my confidence interval lines start touching (or crossing) 1, thereby, implying insignificance. **Question:** What does it mean when we have an interaction A$~\times~$B, and keeping B constant we get significant confidence intervals, but keeping A constant we get insignificant B confidence intervals? I am quite stumped and can't seem to find enough literature on the issue.


title: Building business intelligence tool -- where to start?

text: I'm currently a web developer with a backgound mainly in server side technologies however I've always had an interest in probability so recently I've been exploring building a dummy betting platform where people wager against insurance policies and essentially I'm aiming to use data from Google's index to set the odds -- essentially like the Webbot but not for finance. I know this is a massive endeavor so I don't expect to do it all but any suggestions on where to start are greatly appreciated. At the moment I'm brushing up on game theory, statistics, and some business intelligence tools however I wanted to know if anyone could recommend books or other resources which might be good starting points. Also any suggestions about business intelligence tools also are greatly appreciated.


title: Random parameter linear regression models

text: I am trying to estimate a random parameter linear regression models in LIMDEP (Hildreth, Houck and Swamy's Random Coeffficents Model). Since my data is x-sectional, therefore, I am not using Pds command (being x-sectional data), my commands are as follows:               regress;lhs=......     ;rhs= one,....................     ;Fcn= .......................     ;RPM     ;Pts=50     ;Halton$      LIMDEP works in two step: 1. In step-one LIMDEP use OLS and comes up with starting values for ramdom parameter model and outputs an initial model. My initial model(starting value) as R-Squred of 0.5. It also gives log-likelihood for starating model which in my case is -6600. 2. Instep-2 LIMDEP uses results from OLS model to estimate random parameter regression models. In this step LIMDEP uses maximum likelihood estimates, therefore it outputs log-likelihood at zero (-6900) and at convergence is (-6700).No more R-Squred values in step-2 as LIMDEP is using different technique. 3. I am having issue of the overall fit of the model. If I try to calculate Rouhe-Squred, it is less than 0.02 which is not very appealing. On the other hand I see that starting models used by LIMDEP to estimate random parameter linear regrssion models has an R-Squred of 0.5 which is good enough. In seceond step when LIMDEP gives final results it does not show R-Squred, instead it gives log-likelihood at zero (-6900) and at convergence (-6700).I have 7 total significant variable (four random parameter) all having T-value >2 and P-value<0.001). My question is how should I report my overall fit and is my model OK.


title: Alternative names for items in logistic regression

text: I am reading Applied Logistic Regression by Hosmer and Lemeshow. In chapter 1, they defined $D$ as $-2\ln(\text{likelihood ratio})$, which $D$ is also known as the _deviance_ , and $G$ as $$D(\text{model without the variable}) - D(\text{model with the variable}).$$ I wonder if $G$ also has an alternative name.   Thanks.


title: Listing of typical correlations across spectrum of topics and disciplines to assist with teaching meaning of correlation coefficient

text: ### Context When I teach students about the correlation coefficient, I like to give students a sense of how correlations map on to common associations that might be encountered in daily life and across various topics and disciplines. Some years back I read a psychology journal article (it might have been American Psychologist) that presented a table of meta analytic correlations across a wide range of topics and disciplines (e.g., 40 or 50 correlations from medical, psychological, economic, and other domains). The table aimed to highlight what are typical correlations in a given domain in order to provide context for the interpretation of a correlation coefficient. However, I can't remember the reference for the particular journal article. ### Question   * Does anyone know where a general listing of the magnitude of typical correlations across a wide range of disciplines and topics could be obtained?


I have been performing a relatively routine fit on data for some time. I fit the data to a fairly simple semi-empirical function. For each data point I perform repeated experiments to build a histogram and then fit using standard methods. However the situation is slightly more complicated than this, technically the variable x is a product of lots of values (~60-70) which I combine to create a single variable for the fit. There is in principle no reason that I could not fit to the full expression the only disadvantage being that I could no longer bin the data as each data point would be unique. I do know my experimental uncertainties quite well so I could quite easily fit with a Gaussian maximum likelihood. If I do this I increase my number of variables from one to 60-70. My question is this: are there any interesting effects/pitfalls that occur when fitting with a large number of variables? Although the number of free parameters does not change I feel like I am going from well sampling 1-D space to barely sampling a very high dimensional space, albeit it a very highly correlated one. The fit contains only two free parameters and I would have many hundred data points. If nothing else it would be interesting to see how the fit is biased by my combined variable.