title: Replicate simulation study from a paper and calculate the MSE in R

text: I have implemented a Gibbs Sampler for the **Bayesian Elastic Net** (BEN) according to this paper on Penalized Regression by Kyung et al.   In this paper, they execute a simulation study that has been used in other papers on Penalized Regression (LASSO, Bridge, Ridge) to compare the performance of the proposed models. Here are details of the simulation taken from the above mentioned paper: > We simulate data from the true model $$ y=X\beta+\sigma\epsilon > \quad\epsilon_i\,{\raise.17ex\hbox{$\scriptstyle\sim$}}\,\text{iid}\,N(0,1) > $$ We simulate data sets with $n=20$ to fit models and $n=200$ to compare > prediction errors of proposed models with eight predictors. We let > $\beta=(3,1.5,0,0,2,0,0,0)$ and $\sigma=3$. The pairwise correlation between > $x_i$ and $x_j$ was set to be $corr(i,j)=0.5^{|i-j|}$.   >  Later on they say, that for the prediction error, they calculate the > average mean squared error based on 50 replications. By average they mean > the median in this case. To simulate this data and calculate the MSE I've used following code in R:               # Number of observations     n.train <- 20     n.test  <- 200     # Error variance     sigma <- 3     # Pairwise correlation of X     cor <- 0.5     # Number of predictors     p <- 8     # Create training and test data set (package QRM and mvtnorm required)     Z <- equicorr(p, rho=cor)     X.train <- rmvnorm(n.train,sigma=Z)     X.test  <- rmvnorm(n.test,sigma=Z)     # Create error      error.train <- rnorm(n.train,mean=0,sd=1)     error.test  <- rnorm(n.test,mean=0,sd=1)     # Create beta     beta.true <- c(3,1.5,0,0,2,0,0,0)     # Create both responses     Y.train <- X.train %*% beta.true + sigma*error.train     Y.test  <- X.test %*% beta.true + sigma*error.test          # Fit the training data set with the BEN Gibbs Sampler     beta.ben <- BEN(X.train,Y.train, iter=11000, burn = 1000)     # Calculate the predicted response     Y.pred   <- X.test %*% beta.ben     # Calculate the mean squared error (MSE)     MSE      <- sum((Y.train - Y.pred)^2)/n.train      My problem is that my results are not even close to comparable to the ones in the paper which makes me doubt my simulation study "setup".   As one of the authors of the paper has uploaded the Gibbs Sampler code and I could check if I did something wrong, I know that the problem doesn't lie there. So my questions are:   1. Does anybody have experience with this kind of simulation study and can check if I did something wrong?   2. Is the MSE I calculate the same as the one used in the paper? In researching on this topic I found many different ways to calculate the MSE and it was also sometimes used but actually the mean squared prediction error was meant. For example the Wikipedia article on MSE alone lists three variations. I don't need help with coding, rather more information on how this simulation is typically excecuted so I can figure out what I'm doing wrong.


title: Difference between Variance Inflation Factor (VIF) and kappa in R?

text: I am running a regression analyis in r:               fit <- lm(Cost ~ Slope + YardDist, data = test)      I want to test the two independent variables for multicollinearity. I tested it with `vif()` (from the car package) and `kappa()`.               > vif(fit)        Slope YardDist      1.000121 1.000121      > kappa(fit)     [1] 11631.87      VIF tells me there is no multicollinearity and kappa tells me there is very high multicollinearity. What is the difference between both and which one is 'right'?


title: Obtaining covariance matrix from correlation matrix

text: I am trying to figure out how to convert a correlation matrix (R) to a covariance matrix (S) for input into a random number generator that only accepts S (`rmvnorm("mvtnorm")` in R)               library("mvtnorm")           TRUTH= 0.8 # target correlation value between X1 and X2     R <- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1)))     V <- diag(c(sqrt(1), sqrt(1))) # diagonal matrix of sqrt(variances)     S <- V %*% R %*% V     cor(rmvnorm(100, sigma=S) )          # repeat this to get an idea of the variance around Pearson's estimator      Instance where variances are not equal to 1               V <- diag(c(sqrt(3), sqrt(2)))      S <- V %*% R %*% V     cor(rmvnorm(100, sigma=S) )      This seems to be correct, but I would like expert criticism.


title: Illness-death multistate survival model

text: I am attempting to use and illness-death multi-state scenario to model a time- dependent covariate in a competing risks analysis following Beyersmann et al. 2012 (Competing Risks and Multi-state Models with R, Springer New York, 2012). In general all subjects would start in the "healty state" = 0 with a possibility to transition into the "sick state" = 1 (the time-dependent variable). Either of these states may then transition into one of two competing absorbing states "discharge" = 2, "death" = 3. My specific question is can a patient start in either state 0 or 1 or do all patients need to start in state 0 for this model to work. My data consists of patients starting in both the healthy and the sick states at the beginning of observation. Thanks.


title: Customer Lifetime Value based on RFM (Mixing logistic and linear regression)

text: I'm trying to build a simple model to compute Customer Lifetime Value (CLV) based on Recency, Frequency and Monetary (R,F,M). Browsing the web, I've found this question RFM & customer lifetime value modeling in R and after having looked at the answer given, some doubts have arisen:   1. One of the approaches consists of the following: > _Another popular approach is to build a slightly more complex model for > predicting monetary value based on two sub-models: one for probability of > response (e.g. using logistic regression as a function of RFM), and the > other for revenue conditional on response (again, it could be as simple as a > linear model of RFM). Expected monetary value is the product of the two > predictions._ > I already had my logistic regression output, but I don't know how to mix the response given by it (probability of response) in order to run, for instance, a linear regression for revenue. As you might realize, in the answer quoted above, it's said that it's done conditionally on response, but how can be this carried out? How can I mix those models?   1. Is it OK if I regress revenue (monetary) on the probability _P_ given by logistic regression? I'm afraid this variable is somehow included since I obtained it by running a model with _R, F_ and _M_ as predictor variables.   2. Would you please provide me with some online resources (papers) to get more of this topic?


title: Could I use Tukey's HSD test as a post-hoc test for a within-groups ANOVA test?

text: I understand Tukey's HSD is the most common post-hoc test for between-groups ANOVA. It appears some do not recommend or select it as a post-hoc test for within-groups ANOVA. Instead, they prefer Bonferroni test or Sidak test. But, both of those seem to just multiply whatever _p_ -value you would get from a paired t-test by the number of hypotheses you are testing. This approach seems less than robust. Instead, if I could I would rather use Tukey's HSD test in this situation [within-groups ANOVA]. Could I do that? Or are there specific reasons not to.


title: Experimental design: what to do when factors cannot be held at constant levels?

text: I am looking at some data with a 5 factors and a response variable. The experiment was designed as a factorial experiment, with observations at different levels of each factor. One thing that I have notices, though, is that there is some variability in the readings for the factor levels. For example, motor RPM is a variable, and readings were supposed to be recorded at 800 RPM, 1000 RPM and 1200 RPM. The experimenters were able to set the RPM at these levels for the experiment. There is an RPM monitor that keeps track of the RPM during the experiment, and is recorded along with the data. However, looking at the data, we see that there are slight variations. For example, when set at 800 RPM, we have readings that range between around 790 and 810 RPM. This occurs for all five of our factor variables. The RPM variability is negligible, but some factors are varying more than 15 percent of the desired factor level. This may be due to some error in the readings of the monitoring equipment. My question is, can I ignore the variability in the factor levels and assume that we read the observations at those factor levels, or can I take into account the factor level variability when doing the analysis?


title: Poisson likelihood and zero counts in expected value

text: The following lines are taken from a lecture found here. I have a very small question regarding this explanation of Poisson likelihood. _The Poisson likelihood can be used any time your data come in discrete intervals (which we'll call “counts”), and the counts are independent of each other. Schematically, we imagine dividing data space up into “bins”, which could be bins in energy channel of our detector, location on the sky, time of arrival, or any of a number of other things. Suppose that **in a particular model $m$, you expect there to be $m_i$ counts in bin $i$**. Then if the model is correct the likelihood of actually observing $d_i$ counts in bin $i$ of the data is, from the Poisson distribution,_ $L_i = m_i^{d_i}/d_i! \; exp(-m_i)$ _Note that $m_i$ can be any positive real number, whereas $d_i$ must be an integer. Note also that the sum of $L_i$ from $d_i = 0$ to ∞ is 1. The likelihood for the whole data set is the product of the likelihoods for each bin:_ $L = \Pi\; m_i^{d_i}/d_i! \; exp(-m_i)$ _This becomes better and better approximated by a Gaussian as $m_i$ increases._ (...) _If one uses Poisson likelihoods, small numbers are ﬁne. In fact, if you can manage, the best way to represent your data is to have bins so tiny that you expect either **0 or 1 count per bin**._ The question I have is: what happens if one of the **expected values** $m_i$ is zero (ie: the model predicts zero counts in that bin $i$)? Wouldn't that make the final $L$ value null? # Edit After reading @leonbloy and @Henry's answer I think I should explain a bit more my particular problem with this analysis and what it is that I'm clearly not understanding. In my case I have **two histograms** representing the _modeled_ data and the _observed_ data. So the $[d_i, \; i=1,n]$ histogram is made out of actual observations (ie: the counts in each bin $d_i$ are obtained through an experiment) and the histogram $[m_i, \; i=1,n]$ is produced by a model: this would mean that the $m_i$ counts in my case are **integers** that can actually **take the value 0**. Now, this is where I'm clearly getting something wrong because I read in that article a line that says " _Note that $m_i$ can be any positive real number_ " and @leonbloy's answer also points in that direction, when in my case the $m_i$ values are **all** integers and do **take the value 0** a lot of times. I must be skipping a crucial step concerning the calculation of the $m_i$ values. So how should I calculate the $m_i$ values if my model gives me **counts per bin**? @Henry's answer seems to contradict @leonbloy's answer since he says $m_i=0$ **is a possibility** whereas @leonbloy says **it is not**. In any case, Henry seems to be assuming the whenever $m_i=0$ then $d_i=0$ also so that one obtains $L_i=1$ because $0⁰=1$. This is not my case, since $m_i$ could be 0 but it doesn't mean $d_i$ is forced to be 0 also (my _observed_ histogram, $[d_i , \; i=1,n]$, and my _modeled_ histogram, $[m_i , \; i=1,n]$, are derived independently of each other) # Edit 2 I think I'm beginning to understand, but I'm still a bit lost. Please take a look at this excerpt taken from this paper: _(...) Therefore the selected test to perform on these histogram is a log likelihood test for Poisson statistics (Eidelman et al. 2004)_: $-2ln(\lambda(\theta))=2\, \sum\limits_{i=1}^N \, (\nu_i(\theta)-n_i+n_i\, ln\frac{n_i}{\nu_i(\theta)})$ (8) _In this formula $\theta$ is the set of unknown parameters one wants to derive, $n = (n_1 , n_2, . . . , n_N)$ is the data vector containing the **observations** with $N$ the number of bins in a histogram. $\nu$ are the **expected values** , which are derived from the histograms of the **modelled data** and are therefore dependent on $\theta$. When $n_i = 0$, the last term in Eq. (8) is set to zero (...)_ As you can see, the analysis is the same (the difference here is that the last equation shows the likelihood _ratio_ ), now the _expected values_ are $\nu_i(\theta)$. From what I've understood so far that parameter $\nu_i(\theta)$ should be **independent** of the number of counts in the bins and have **just one value** for each histogram. If this is so then why use the $i$ subscript? And how should I calculate that value for each histogram?


title: Thoughts on model self-penalization amidst difficult parameter estimation

text: It is well accepted that one should account for model complexity when performing model comparisons, and the general procedure is to penalize more complex models more strongly. While this makes sense when the parameters of a given model are easily estimated (i.e. analytically, as with the mean, variance, etc), it occurs to me that if parameter estimation is a more difficult endeavour then more complex models may to some degree self-penalize. That is, if parameter estimation requires search of a parameter space, presumably larger parameter spaces are more difficult to search and therefore any finite search algorithm is more likely (as the parameter space expands) to terminate prior to finding the point of global maximum likelihood. Has this idea been considered in the statistical literature at all?


title: Interpreting the difference between lognormal and power law distribution (network degree distribution)

text: First off, I'm not a statistician. However, I have been doing statistical network analysis for my PhD. As part of the network analysis, I plotted a Complementary Cumulative Distribution Function (CCDF) of network degrees. What I found was that, unlike conventional network distributions (e.g. WWW), the distribution is best fitted by a lognormal distribution. I did try to fit it against a power law and using Clauset et al's Matlab scripts, I found that the tail of the curve follows a power law with a cut-off. ![enter image description here](http://i.stack.imgur.com/TM2je.jpg) Dotted line represents power law fit. Purple line represents log-normal fit. Green line represents exponential fit. What I'm struggling to understand is what this all mean? I've read this paper by Newman which slightly touches on this topic: http://arxiv.org/abs/cond- mat/0412004 Below is my wild guess: If the degree distribution follows a power law distribution, I understand that it means there is linear preferential attachment in the distribution of links and network degree (rich gets richer effect or Yules process). Am I right in saying that with the lognormal distribution I'm witnessing, there is sublinear preferential attachment at the beginning of the curve and becomes more linear towards the tail where it can be fitted by a power law? Also, since a log-normal distribution occurs when the logarithm of the random variable (say X) is normally distributed, does this mean that in a log-normal distribution, there are more small values of X and less large values of X than a random variable that follows a power law distrbution would have? More importantly, with regards to network degree distribution, does a log- normal preferential attachment still suggest a scale-free network? My instinct tells me that since the tail of the curve can be fitted by a power law, the network can still be concluded as exhibiting scale-free characteristics. Your help is much appreciated. Best regards, Michael


title: Limiting joint distribution of estimators; Functional Statistics; Influence curves;

text: Let $X_1,...,X_n$ iid r.v. with distribution F, with mean $\mu$ and median $\theta$.Assume that $Var(X_i)=\sigma^2$ and $F'(\theta)>0$. If $\hat{\mu}_n$ is the sample mean, and $\hat{\theta}_n$ the sample median, then using the influence curve/functions, show that $\left( \begin{array}{ccc} \sqrt{n} (\hat{\mu}_n-\mu) \\\ \sqrt{n} (\hat{\theta}_n-\theta) \end{array} \right)\xrightarrow[d]{}N(\mathbf{0},C)$, with $C$ a variance-covariance matrix. Well, we know that $\sqrt{n} (\hat{\theta}_n-\theta)\xrightarrow[d]{}N(0,(2F'(\theta))^{-2})$ (Here I used some properties of influence curves of parameters that give 'linear approximation' similar to delta method) and also that $\sqrt{n} (\hat{\mu}_n-\mu)\xrightarrow[d]{}N(0,\sigma^2) $ So, I was thinking that I could immediately conclude what the exercise asks, remaining only the entries of C matrix to be explicitly determined. Can I do this? Thus, $\mathbf{C}=\left( \begin{array}{ccc} \sigma^2 & \lim Cov(\hat{\mu}_n,\hat{\theta}_n) \\\ \lim Cov(\hat{\theta}_n,\hat{\mu}_n) & (2F'(\theta))^{-2} \end{array} \right)$ Next, using properties/definition of influence curves (IC), and one other assumption, $\sqrt{n} (\hat{\mu}_n-\mu)=\frac{1}{\sqrt{n}}\sum IC_{\mu}(X_i,F)+R_1=\frac{1}{\sqrt{n}}\sum (X_i-\mu)+R_1$ $\sqrt{n} (\hat{\theta}_n-\theta)=\frac{1}{\sqrt{n}}\sum IC_{\theta}(X_i,F)+R_2=\frac{1}{\sqrt{n}}\sum \frac{sign(X_i-\theta)}{2F'(\theta)}+R_2$ With $R_1,R_2\rightarrow_p0$. Using the RHS expressions I've tried computing the Covariance in the matrix C, and obtained $\frac{\mu-2E(X_iI(X_i<\theta))}{2F'(\theta)}$ Is this deduction correct? Any help would be appreciated.


title: MCMC algorithm to estimate beta and variance

text: I'm looking for a generic MCMC algorithm for a linear model. I've been reading a lot of articles online and they are so confusing. I am hoping to understand how the method works theoretically **Can anybody help me with a step by step algorithm for implementing a MCMC method to estimate beta and variance ($\sigma^2$)?**


title: How to evaluate Likert scale data changes over multiple surveys of the same group?

text: I have five surveys of the same group of students over a semester. Each survey uses a 5-point Likert scale. The first and last survey contain some questions dealing with the beginning and end of the class (first impressions, final impressions), but most of the questions are identical for all four or five of the surveys. I want to evaluate the statistical significance of changes to students' responses over time. Unfortunately statistics is not my strong suit. I know of the t-test, but that seems to only be applicable to two groups of data (please correct me if I'm wrong). How should I go about evaluating this data? Is a repeated measures one-way ANOVA appropriate?


title: Having trouble with hidden layer neurons in neural network

text: I have some doubts related to my neural network implementation. I have 750 features and 98 outputs. I have number of samples = 5000. Now I used cross validation to select the number of neurons in the hidden layer(1 hidden layer). I used 10 fold cross validation. I got that number of neurons = 30 performed the best which is kind of strange. I couldn't explain what's going on here. Any ideas/ suggestions? How could just 30 neurons do better with such high dimensional data I used tansig activation for the hidden layer and logsig for the output. The outputs used for training were scaled to the range of [0 1]. The inputs where standardized to zero mean and unit variance.


title: Relation of slopes of predictors when they are correlated in linear regression

text: Suppose we have two independent variables and one response variable, we can fit three different models: \begin{align*} Y &= \beta_{00} + \beta_{10} X_1 + \epsilon \\\ Y &= \beta_{01} + \beta_{20} X_2 + \epsilon \\\ Y &= \beta_{02} + \beta_{11} X_1 + \beta_{21} X_2 + \epsilon \end{align*} When $X_1$ and $X_2$ are correlated, we know the estimates $b_{10}$ will not be equal to $b_{11}$, and similarly $b_{20}$ will not be equal to $b_{21}$. But what exactly are the relationships between these four estimates? What about the p-values? Particularly, if $X_1$ is highly correlated with $X_2$, by experience and intuition I tend to see $b_{10}$ and $b_{20}$ and their p-values quite close to each other. Is there a formula to support this observation?


title: Categorical independent variable and logistic regression

text: I am creating a logistic regression model as follows. The dependent variable is outcome of a game (Win/Lost) and the independent variable is the degree of MOON on the day of match. So when I take a match, I get the value of Moon in degrees (1 to 360). I feel that degree is just a measurement from a fixed line of reference point and hence I have to consider the position of Moon as a categorical variable by dividing the 360 degrees into groups and noting the group in which Moon existed during a match. So I divided the degrees into 12 groups of 30 degrees (based on sun sign) each and note the position of Moon. Am I right? Can I still divide the whole 360 degrees into 27 divisions (83 or 249 divisions) to get exact information? Can I have 27 or 83 possible outcome for a single categorical variable? Let me know if I need to explain my question further. As of now i have a sample of 900 entries.


title: Complete substantive examples of reproducible research using R

text: **The Question:** Are there any good examples of reproducible research using R that are freely available online? **Ideal Example:** Specifically, ideal examples would provide:   * The raw data (and ideally meta data explaining the data),   * All R code including data import, processing, analyses, and output generation,   * Sweave or some other approach for linking the final output to the final document,   * All in a format that is easily downloadable and compilable on a reader's computer. Ideally, the example would be a journal article or a thesis where the emphasis is on an actual applied topic as opposed to a statistical teaching example. **Reasons for interest:** I'm particularly interested in applied topics in journal articles and theses, because in these situations, several additional issues arise:   * Issues arise related to data cleaning and processing   * Issues arise related to managing metadata   * Journals and theses often have style guide expectations regarding the appearance and formatting of tables and figures   * Many journals and theses often have a wide range of analyses which raise issues regarding workflow (i.e., how to sequence analyses) and processing time (e.g., issues of caching analyses, etc.). Seeing complete working examples could provide good instructional material for researchers starting out with reproducible research.


title: What data qualities could one use to evaluate the likelihood of a paper's result being correct?

text: I'm in the initial stages of brainstorming a way to provide simplified, at-a- glance summaries of scientific papers encapsulated in a few statistics, and possibly a value on a scale based on weighted versions of those statistics in order to provide a quick, approximate visual for the relative likelihood that the data within is correct versus a contrasting paper. For example, if I tracked the replication status of a result on a scale (from lowest to highest) of "Did Not Replicate, No Replication Attempted, Successfully Replicated", then a person looking at two papers with opposing conclusions would assume that the paper with status "Did Not Replicate" is less likely to be correct than the paper with the status "Successfully Replicated" (or "No Replication Attempted"). The assumption is that, all other metrics being equal, being "higher" on the scale is corresponds to an increased likelihood. **What I'm curious about is, what are more metrics (in experimental design or otherwise) that can be used as a rule of thumb to "rank" the relative likelihood of a paper's results being correct versus another paper on the same subject? Is this even possible?** My proposed metrics so far:   * _**Replication Status_** : Did Not Replicate, No Replication Attempted, Successfully Replicated   * _**Type of Study_** : Retrospective, Prospective, Randomized Experiment   * _**(Experiments Only) Blindness_** : None, Blind, Double Blind   * _**Population Selection_** : Non-random, Random   * _**(Surveys Only) Survey Timeframe_** : Past, Immediate Past, Real-time   * _**Sample Size_**   * _**Confidence Interval_**


title: Average rate versus average period

text: I am working on a project where I need to do workload balancing, and to see how well my approach is working I want to compute both the average amount of time between requests for a workload and the average amount of time for a requested workload to be obtained. Here is what is bothering me. My data consists of a list of time intervals (e.g., 0.4 seconds, 1.3 seconds, etc.). From this data, I can compute either the average rate at which an event is taking place by taking the _reciprocal_ of each time interval to turn it into an instantaneous rate and averaging over that, or I can compute the average period by directly averaging over the time interval data. Doing the first is equivalent to computing the harmonic mean of the data, and doing the second involves computing the arithmetic mean of the data, so the two analyses get different results. What confuses me is that, although the two means get different results, a rate is equal to the inverse of a period, and so the two results seems like they should be related via reciprocation, even though they are not. Another option that I have is to use the geometric mean which will obtain consistent results for the rate and the period (as reciprocation commutes with taking the root), but my limited understanding says that one should only use the geometric mean in cases where you are taking the mean of something like growth rates. So in short, it is not clear to me which of these means --- harmonic, arithmetic, geometric, or something entirely different --- should be used in this situation. Any insight would be appreciated.


title: Use median for clustering, instead of average value

text: I have a List of Files on which I run some statistical methods to get a rating between 0 and 1. Now i need to cluster the files for similarities. The easiest way i thought that would be is to use the average value of the rating and cluster by n*sigma. So the biggest part of all files will be in the range of `[avg-sigma, avg+sigma]`, then the next part will be `(avg+sigma, avg + 2*sigma]` and so on... the problem is now if the values are very wide spreaded and have a distribution very close to one side of the range the clustering will be not very good. e.g. ten values of 0.01 and one value of 0.9 gives me avg of 0.1 but a median of 0.01... so a typical problem here... now i think it would be better to use the median as a starting point but can i use the standard deviation to calculate my clusters? Or is there even a better method then using StdDev to create clusters?


title: Probit analysis shows no significant variables

text: In my probit output all p values are insignificant. I want at least some of my variables to be significant. How I correct it? My variables are: crossbreed, age2, gender, education experience2, fullparttime, localbreed, farmsystem, ownland2, hirelabor, expensemoyr2, milkincome2, vetanaryattendence2, vetanaryservice, traninggot, furthertraning, health.


title: Spatial autocorrelation of non-homogeneous point data in irregularly shaped sample plots

text: Let's assume we registered the position of all individuals of different plant species within 3 irregularly shaped sample plots that are very close the each other. As an example, I put 2 pictures with 2 different species below. Each grid corresponds to 100m, and each circle represents one plant (the size of the circles is not important for this question). ![Species 1](http://i.stack.imgur.com/1xYpP.png) ![Species 2](http://i.stack.imgur.com/WZYsf.png) We want to analyse to what extent the occurrence of individuals depends on various ecological variables. In a number of similar studies, the authors provided a quantitative measure of the degree of intraspecies global spatial autocorrelation for each species, usually Moran's I. However, I'm not sure how we could do this in our case, because:   * A global analysis of spatial autocorrelation, such as Moran's I, assumes homogeneity within the study area, which is clearly not the case.    * The shape of the sample plots is quite irregular.   * Many species are completely absent in one of the 3 plots. We could try and find clusters at a local level using local spatial autocorrelation, such as LISA. If I understood LISA correctly, we would have to separate our study area in different homogeneous sub-units, which seems somehow arbitrarily. What Information about intraspecies spatial auto-correlation could we provide based on our data? Would this even make sense, from a statistical point of view, or should we better just describe the pattern qualitatively? By the way, we are using R for the analysis. Edit: Maybe some additional information about our study would help to answer the question. So far, we used present-only species distribution modelling techniques (like MaxEnt). Results are generally satisfying, with good model performance and predictions that make sense from an ecological point of view. However, some authors warn that spatial autocorrelation between occurrences could influence the model. Therefore, we would like to present a _simple_ measure (if possible) of spatial autocorrelation to describe the data. Edit2: Regarding the suggestion to use Ripley's K or L-function: As far as I know, those functions assume that underlying point processes are spatially homogeneous (stationary), which wouldn't be correct in our case. K and L-functions for inhomogeneous point pattern have been developed in recent years, which take non-constant intensity at each site in the study area into account (see for example here for an emulation in R). However, for our analysis, we are not really interested in spatial correlation _after_ accounting for different intensities, but rather if different intensities are spatially auto-correlated (I don't know if this makes completely sense).


title: Supervised learning algorithms for classification, that we should read first

text: What are the main supervised learning algorithms for classification (more than 2 classes), that we should learn first when we are beginners in that domain ? It is good if you can also give easy-to-understand tutorials for the algorithms that you propose here.


title: T-test when observations are years

text: I would like to know whether I can use the t-test or the non-parametric equivalent test when I have years as observations. Suppose I want to compare the profitability levels of two companies and I have profitability data over 10 years for each of the company. I thought I could use the company as if it was a treatment (i.e. company A or Not-Company A), and the years as cases. However, I fear this may pose problems of dependency, since the observations over time are likely to not be independent as would be, say, patients drawn randomly from a population. I would very much appreciate any help, and would be happy to provide additional explanation if needed.


title: quantile in scipy library

text: I am confused what is quantile in scipy.stats.norm.pdf function. Definition of quantile says that k-th of q-quantile is essentially value which divides population into k/q and (q-k)/q parts. But if I call               scipy.stats.norm.pdf([0,1,2,3,4], 2, 9)      then I will get vector `v` with 5 numbers. What `v[i]' mean in this case? probability of normal RV been equal`i`? Could you clarify it please?


title: Expected value of a transformed random variable

text: I would be interested in computing the expected value of the following random variable $Y$. Let $X$ be either a $\mathrm{Bin}(p,N)$ or a $\mathrm{Hyp}(n,m,N)$ with $m$ number of successes (I would be interested in both proofs). Is there a way to compute the **expected value** of $$ Y = X \log X $$ I tried to go through the moment generating function but I haven't had any result yet.


title: How can I tell which location is the best?

text: I have a list of the scores from various locations on different exams listed below. How can I compare the different test and scores to show which location is the best over all? It would really help if someone could inform which type of test to do and why. I first thought to use the Wilcoxon signed-rank test, but the data is not in pairs.                            English  AP CalcBC    AP PhyscsB   AP Chem     AP Econ      City      Rank  Score Rank Score  Rank  Score Rank  Score  Rank Score          (Beijing)   6   36.33   4   39.14   1   35.54   1   47.94   3   46.55     (Hefei)     5   38.95   11  26.31   9   18.61   8   27.83   8   30.52     (Hexi)     na      na   8   33.04   na     na   6   37.7    7   32.46     (Huzhou)    4   39.06   1   45.47   na     na   na     na   na     na     (Jiangyin)  9   29.9    9   30      5   24.4    7   31.75   6   36.7     (Jinan)     1   47.72   6   37.2    3   28.91   3   41.55   2   47.94     (Nantong)   3   39.67   5   38.45   2   32.59   5   38.9    4   45.05     (Shanghai)  na     na   3   39.44   6   23.18   2   44.4    5   43.86     (Suzhou)    8   30.19   10  27.5    7   21.75   9   26      9   19.87     (Yangzhou)  2   43.59   2   40.45   4   26.7    4   41.45   1   50.73     (Yixing)    7   33.26   7   35.91   8   21.48   10  24.7    na     na


title: Calculating Gini coefficient with unbound income brackets?

text: I need to compute the Gini coefficient on some population data arranged in income brackets: for example               $0->$1000 : 10000 people     $1000->$10000: 50000 people      My problem is that the last bracket is unbounded i.e it's in the form:               $1000000 < : 500 people      Is there any way to calculate the Gini coefficient given this data?


title: Is there a practical limit to the size of a piece of data with 0 statistical redundancy?

text: I was just wondering, just because I tend to wonder about these sort of things in my spare time, if there are any data pieces that are 1Gig+ that have 0 Statistical Redundancy? (i.e., uncompressible with lossless compression.)   Is it even possible for a file larger than a few bytes to have such a trait? Is that property possible at all, or merely theoretical? Anyway, please let me know if such a thing would be possible, and, if feasible, please link me to a place where I could view such a piece of data. Google gave me nothing. Might as well include the Wikipedia page about statistical redundancy: en.wikipedia.org/wiki/Redundancy (information theory)


title: Algorithm to produce autocorrelated uniformly distributed number

text: I would like to produce a time-series of autocorrelated probabilities (with a predefined mean level of autocorrelation). I've spotted this and this which I believe should give me what I'm looking for, but I'm unable to access them. If anyone could give me some hint on this 'sum of uniforms' method it would be appreciated. Many thanks.


title: How to generate an orthogonal vector

text: I have two vectors $a$, $b$ each containing $10$ random numbers from the standard normal distributions. I want to generate another vector $C$ of $10$ numbers from the standard distribution where $\mathbb{E}(a\cdot C)=\mathbb{E}(b\cdot C)=0$ (orthogonal), where $\mathbb{E}$ is the expectation.


title: What are good techniques and resources for teaching Bayes theorem?

text: My friend and I want to do a hands on tutorial on Bayes theorem for the Seattle LessWrong group. Neither of us have done this before, so we're searching for prior art; techniques that other people have tried before and descriptions of how they turned out. What are good techniques and resources for teaching bayes theorem? Reports of both successes and failures are useful, I'd like to know what _not_ to do in addition to what to do. The audience is a group of ~8 programmers and natural science students. They'll be smart and capable but not necessarily used to doing much math.


title: OLS Regression: Interaction effect is significant by itself, not when included in full model--ok to report?

text: I have a sample with 400 cases. When I run my full model, which includes 13 predictors, the interaction term is non-significant. However, when I run a model only including the three variables involved in the interaction (two main effects + interaction effect), the interaction effect is significant. I've tracked the problem down to one of the other IVs in the full model, which is extremely highly correlated with one of the items in the interaction term (r=.75), although collinearity tests indicate no problems with multicollinearity. Both the items included in the interaction term and the one which is highly correlated are composite scales, consisting of between four and nine items each. My question is: When reporting results from my analyses, I have a hypothesis about the interaction. Is it okay to talk about the significant interaction effects that exist in the absence of controlling for this other IV, which while not theoretically important is still important to include in the regressions. If so, how do I justify the inclusion?


title: How do I study the "correlation" between a continuous variable and a categorical variable?

text: What's a meaningful "correlation" measure to study the relation between the such two types of variables? In R, how to do it? Thanks a lot!


title: Example sets for Shlosser distinct value estimator?

text: I've implemented Shlosser's estimator for the number of distinct values in a population given a sample. Now I want to make sure my algorithm is correctly implemented. Are there any example sets out there that would allow me to verify my computation matches the canonical value?


title: How to conduct a meta-analysis using raw data?

text: I think meta-analyses are a great way of exploring a hypothesis when available evidence is heterogeneous. Usually however when conducting a meta-analysis one puts aggregated data in a model, possibly losing information. Therefore I am considering gathering data from published research from the authors in order to be able to conduct analysis on 'raw data' rather then aggregated data (my field is neuroscience/psychology/psychiatry). ### Questions   * Where can I learn more about conducting such analyses?   * Are researchers typically willing to share their data?   * Would ethical guidelines preclude such data sharing?


title: Correcting p values for multiple tests where tests are correlated (genetics)

text: I have p values from a lot of tests and would like to know whether there is actually something significant after correcting for multiple testing. The complication: my tests are not independent. The method I am thinking about (a variant of Fisher's Product Method, Zaykin et al., Genet Epidemiol, 2002) needs the correlation between the p values. In order to estimate this correlation, I am currently thinking about bootstrapping cases, running the analyses and correlating the resulting vectors of p values. Does anyone have a better idea? Or even a better idea for my original problem (correcting for multiple testing in correlated tests)? Background: I am logistically regressing whether or not my subjects are suffering from a particular disease on the interaction between their genotype (AA, Aa or aa) and a covariate. However, the genotype is actually a lot (30-250) of Single Nucleotide Polymorphisms (SNPs), which are certainly not independent but in Linkage Disequilibrium.


title: Why would predicted values be normally-distributed when the actual values are uniform?

text: I'm building a supervised learning model where the target variable is a uniformly-distributed continuous value ranging from 0-1 (originally a rank value from 1-38000, then scaled down to 0-1). The 20 predictor variables are a mix of continuous and categorical variables, with no more than 6 levels for the latter. But when I apply various modeling techniques to the data (including regression, random forest, MARS, neural nets) using R, the predicted values for every model are shaped like a bell curve. What could cause the fitted values (bell shape) to be so consistently different from the actual ones (uniform shape)? And is there any way to force these models to produce a more uniformly-shaped prediction?


title: How to choose the best transformation to achieve linearity?

text: I want to do multiple linear regression and then to predict new values with little extrapolation. I have my response variable in the range from -2 to +7, and three predictors (the ranges about +10 - +200). The distribution is nearly normal. But the relationship between the response and the predictors is not linear, I see curves on the plots. For example like this: http://cs10418.userapi.com/u17020874/153949434/x_9898cf38.jpg I would like to apply a transformation to achieve linearity. I tried to transform the response variable by checking different functions and looking at the resulting plots to see a linear relationship between the response and predictors. And I found that there are many functions which can give me visible linear relationship. For example, functions $t_1=\log(y+2.5)$ $t_2=\frac{1}{\log(y+5)}$ $t_3=\frac{1}{y+5}$ $t_4=\frac{1}{(y+10)^3}$ $t_5=\frac{1}{(y+3)^\frac{1}{3}}$ etc. give the similar results: http://cs10418.userapi.com/u17020874/153949434/x_06f13dbf.jpg After I am going to back-transform the predicted values (for $t=\frac{1}{(y+10)^3}$ as $y’=\frac{1}{t^\frac{1}{3}}-10$ and so on). The distributions are more or less similar to normal. How can I choose the best transformation for my data? Is there a quantitative (and not very complicated) way to evaluate linearity? To prove that selected transformation is the best or to find it automatically if possible. Or the only way is to do the non-linear multiple regression? Thank you for any suggestion!


title: Effect size and statistical significance

text: Using Cohen's $d$, I am getting small and medium effect sizes for results that are statistically non-significant ($p>.05$). Does this make sense?


title: Mahalanobis distance when n>p

text: I have an _n by p_ matrix, where _n_ is genes and _p_ is patients. Anyone whose worked with such data knows that _n_ is always larger than _p_. Using feature selection I have gotten _n_ down to a more reasonable number, however _n_ is still greater than _p_. I would like to compute the similarity of the patients based on their genetic profiles; I could use the euclidean distance, however Mahalanobis seems more approperiate as it accounts for the correlation among the variables. The problem (as noted in this post) is that Mahalanobis distance, specifically the covariance matrix, doesn't work when _n_ > _p_. When I run Mahalanobis distance in R, the error I get is: > Error in solve.default(cov, ...) : system is computationally singular: > reciprocal condition number = 2.81408e-21 So far to try solve this, I've used PCA and instead of using genes, I use components and this seems to allow me to compute the Mahalanobis distance; 5 components represent about 80% of the variance, so now _n_ < _p_. **My questions are:** Can I use PCA to meaningfully get the Mahalanobis distance between patients, or is it inappropriate? Are there alternative distance metrics that work when _n_ > _p_ and there is also much correlation among _n_?


title: Multinomial probit in R

text: I ran MNP official package in R, like below.               data(japan)     res <- mnp(cbind(LDP, NFP, SKG, JCP) ~ gender + education + age, data = japan,     n.draws = 10000, verbose = TRUE)     summary(res)      and the result was               Coefficients:                        mean    std.dev.  2.5%   97.5%     (Intercept):LDP 0.615184 0.517157 -0.386151 1.6         (Intercept):NFP 0.689753 0.568109 -0.419521 1.79     (Intercept):SKG 0.133961 0.455960 -0.758883 1.02     gendermale:LDP  0.099748 0.152323 -0.194786 0.40     gendermale:NFP  0.216824 0.166103 -0.102108 0.54     gendermale:SKG  0.132661 0.134605 -0.127145 0.40     education:LDP  -0.107038 0.074792 -0.253483 0.04     education:NFP  -0.107222 0.082324 -0.270127 0.05     education:SKG  -0.003728 0.066429 -0.132496 0.13     age:LDP         0.013518 0.006122  0.001492 0.03         age:NFP         0.006948 0.006783 -0.006572 0.02         age:SKG         0.009653 0.005431 -0.000812 0.02      From this result, how can I know the mean value between female and LDP/NFP/SKG respectively? There is another option for gender; Female, but it does not appear in the result with the mean value.


title: Annealing MCMC with constant temperature?

text: Let's say I'm trying to sample from a posterior distribution, $p(\theta|x)$, which I suspect is very "bumpy" (10s to 1000s of small modes). Moreover, the evaluation of the posterior kernel $K(\theta|x) \propto p(\theta|x)$ is very expensive to evaluate (and not conjugate) and the data set is large enough such that I would like to avoid any technique based on parallel chains. However, I would like to avoid being stuck in any one mode for too long during the sampling process, so I employ a standard annealed MCMC technique (such as Metroplis-Hastings with an annealed acceptance probability) but I never cool the temperature parameter. The distribution I'm sampling from is thus "flattened", and easier to sample from in the case of many modes, but the stationary distribution for the Markov chain has clearly changed from the one defined by $K(\theta|x)$. Does anyone have any references about what happens to the stationary distribution of the Markov chain in this case?


title: What does the change in intercept mean if you allow for clustering in your data?

text: When I estimate a logistic binomial regression with an allowed to vary effect on a second and a third level, the overall intercept represents the average log-odds for any level one unit in any higher level unit: $\text{logit}(\pi_{ijk}) = \beta_0 + \mu_{jk} + \mu_{k}$ Where $\beta_0$ is the fixed intercept and $\mu_{jk}$ and $\mu_{k}$ are the random cluster effects. $\beta_0$ is estimated as -3.525 in this case (or 2.9% probability of success: $\pi_{ijk} = \frac{e^{\beta_{0}}}{1+e^{\beta_{0}}}$). However, the model where no random effects are allowed (or the general mean): $\text{logit}(\pi_{ijk}) = \beta_0$ gives -2.932 log-odds or 5% probability of success. Does the decrease of overall intercept mean there were more clusters with low probability than with high probability? PS: I also switch from MQL to PQL when going to the random intercept model (since PQL is not possible for the fixed effect only model). It could have to do with that too.


title: References for spatial modeling with Bayesian belief networks in medical applications

text: I want to do research in spatial data mining where I want the concept of Bayesian belief networks applied on a medical domain like, for example, cancer. I have been searching for recent papers in various good journals but I am not able to find good papers where all these three concepts are together. Can anybody help me in providing links to these kinds of papers? Any sort of help would be appreciated.


title: Blocking and confounding in replicated $2^k$ factorial design

text: Consider the $2^6$ factorial in 8 blocks of eight each runs with ABCD, ACE and ABEF as the independent effects chosen to be confounded with blocks. Generate a design. Find the other effect of confounding. I can solve the second part of the question. I found the solution of the first part in a solutions manual (question 7-11) for Montgomery's textbook on the _Design and Analysis of Experiments_ (Wiley, 2001), but I could not understand this table. I want some one to explain this table. If we change the number of blocks, how do we construct the design?


title: What is correct terminology when describing the number of response levels of a Likert scale?

text: I know that the Likert scale is the sum (or some sort of combination) of all the responses of a survey, test etc. and the actual numeric responses are called response levels. Is this true of all psychometric numeric responses? ie are those responses also called response levels or are they some times called a scale (appropriately called that is)? If you have any literature to back up what you say it would be appreciated.


title: Difference in KS,Lift,Concordance in validating logistic regression

text: What is main difference in KS,Lift,Concordance in validating logistic regression? What are the advantages and disadvantages in using these measures to see performance of the model?


title: Does having self and other report measures as predictors require dyad analysis?

text: I have three different measures of child anxiety (IVs). For each measure, I have a self-rating by the child and an other-rating by a parent. I want to use these measures to predict a particular diagnosis (DV).   * How do I take into account the dyadic nature of the parent-child relationship?  I plan to use hierarchical regression or path analysis to analyse the predictive relationship between the IVs and DV.


title: Plot of quantity and price in purchases of the same product - Should it be a power-law distribution?

text: I work for the Brazilian Government and I realized recently that most of the purchase data I have tested is a power law. I am wondering what does it mean? At first, I thought that the relation between price and volume bought should be similar to a straight line without using the log function. Then, I thought about it and it kind of makes sense to think that I get twice the discount as I increase the volume ten-fold, for instance. However, I am having a hard time actually understanding what this means. Could someone help out? To make it more clear. One of the data I have is the quantity vs price of diesel purchases made by the Government in 2012. So, for every purchase made by different agencies, I have the price paid per liter and the quantity it bought. When I plot this two in a simple scatterplot, I pretty much get a power-law distribution. Well, first, there are a few problems with my data, since some quantities are totally wrong (a lot of people buying just 1 liter for as much as thousands of dollars for that single liter!). Besides that, on the other end of the distribution, I have those buying a whole bunch of liters (thousands) and paying 0 for each liter! Exactly, paying ZERO, nothing, nada. So, you can imagine we have some typing issues on the user end when filling out the quantities. One thing they never get wrong is the total price, since they need that to pay the seller. Anyway, even if I remove these outliers, I still end up with a power-law distribution (at least it seems so). In the end of the day, I was just trying to understand if there is a correlation between quantity bought and price paid. And here goes my question. Should the high volume/low price be represented as a power-law distribution? In other words, if I have a power- law, does it mean that I have a high correlation between price and quantity? If I compute the `cor(log(data$valuePerLiter), log(data$quantity))`, I get something around -.70, which seems to be a decent correlation. However, if I compute the `cor(data$valuePerLiter,data$quantity)` I get a value close to 0. So, do I have a strong correlation between price and quantity or not? Is a power-law distribution expected in this scenario? Thanks in advance! Rommel


title: Filling NAs in a dataset with column medians in R

text: I have a massive dataset with titles across the top for the row. All the values in the dataset are number minus about 50 of them which are NA. I need to find the median of each column whilst somehow not selecting the title of the column. Then once I have the median of that row I need to replace all potential 'NA's with the median of the corresponding column! Link to image of dataset.


title: Is this regression significant?

text: Hi I get this output from R summary of an lm:               lm(formula = weight.nz ~ dChgs.nz)          Residuals:          Min       1Q   Median       3Q      Max      -15373.7   -664.4    243.3   1104.2   9137.2           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 1.853e+00  2.141e+01   0.087    0.931           dChgs.nz  7.036e+07  5.841e+06  12.046   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1814 on 7464 degrees of freedom     Multiple R-squared: 0.01907,    Adjusted R-squared: 0.01894      F-statistic: 145.1 on 1 and 7464 DF,  p-value: < 2.2e-16      I think small p suggest significance BUT r2 is tiny - so not significant. I have read about problems of r2 on other questions on this subject that suggests I use stdError as a guide, but I am not sure how to interpret it. How do I use 1814 to tell me whether my result is significant or not?


title: How does SVM prediction work?

text: As far as I understand, the training phase usually uses the dual optimization formulation where we can implicitly calculate the weight vector which defines the discriminant function. How about the prediction phase, how do we use these weights and the kernel function when a new test sample arrives? edit: I should clarify, I am interested in the nonlinear SVM.


title: What kind of functions can have non whole degrees?

text: Thanks for the help in advance. I am reading a technical report on a regression algorithm that reports a pair of functions as having a total degree of freedom of 5.4. I believe that both of these functions are splines, but it isn't directly specified in the report (they are graphed but the explicit functions are not given \ they could also be polynomials). How might a function, or a set of functions, have a total degree of freedom that is not a whole number? What would be an example function exhibiting this behavior?


title: prerequisites for z-standardization

text: I tried to search for an answer to this but I was not successful in finding a good one. I wanted to compare data from two modalities together (an EMG reflex and self- rating). I found that in the literature there were studies used that Rhudy et al., 2008 ECON. The only thing I couldn't find an answer for is: Are there any prerequisites for z transformation? How I can check if I can do that or not?


title: When does replication reveal fraud?

text: Found this question here, and I couldn't answser it and thought it would be a wonderful question for this site. > Here's a little thought experiment for your weekend pleasure. Consider the > following: > > Joe Scientist decides to conduct a study (call it Study A) to test the > hypothesis that a parameter D > 0 vs. the null hypothesis that D = 0. He > designs a study, collects some data, conducts an appropriate statistical > analysis and concludes that D > 0. This result is published in the Journal > of Awesome Results along with all the details of how the study was done. > > Jane Scientist finds Joe's study very interesting and tries to replicate his > findings. She conducts a study (call it Study B) that is similar to Study A > but completely independent of it (and does not communicate with Joe). In her > analysis she does not find strong evidence that D > 0 and concludes that she > cannot rule out the possibility that D = 0. She publishes her findings in > the Journal of Null Results along with all the details. > > From these two studies, which of the following conclusions can we make? > >   1. Study A is obviously a fraud. If the truth were that D > 0, then Jane > should have concluded that D > 0 in her independent replication. > >     1. Study B is obviously a fraud. If Study A were conducted properly, > then Jane should have reached the same conclusion. > >     2. Neither Study A nor Study B was a fraud, but the result for Study A > was a Type I error, i.e. a false positive. > >   2. Neither Study A nor Study B was a fraud, but the result for Study B was > a Type II error, i.e a false negative. > > > > I realize that there are a number of subtle details concerning why things > might happen but I've purposely left them out. My question is, based on the > information that you actually have about the two studies, what would you > consider to be the most likely case? What further information would you like > to know beyond what was given here?


title: How do I analyse a two-alternative forced choice test?

text: I have recently as part of my bachelor project at my university conducted a two-alternative forced choice test. I have looked at ways to say something with significance about which images (the test participants were given a set of 4 pairs of images, for each pair they had to choose one). While testing is emphasized in my education, I do not have enough knowledge about analysis of the results I have now. How do I go about figuring something significant out about the independent variables that I am manipulating (in separate scenarios)? Is logistic regression something I could use or are there methods that are easier to understand for someone who has only had one course of statistics? **The data (and edits about the data)** : The data consists of rows of test participants and four columns named as image pairs that the participants were choosing between. The values are the test participants' choice (between the two images in the pair) I have four independent variables, three of them scales and one nominal that I am changing in four scenarios respectively. Each scenario has _paths_ that are chosen "at random" (depending on the last number on the test participant's digital watch e.g. if the time is 13:37, the test participant takes path 7). The reason for the paths is a detail of the implementation in Google Forms to allow all image pair combinations to be presented but not all to every test participant. I would like to find a model for achieving either high probability of one outcome of the dependent binary variable.   * IV: Amount, Spread, Diff objects, Placement model (nominal)   * DV: Image chosen or not. e.i. if I want a high probability of YES/1/ACCEPT/etc answers (to the binary question: "Does the virtual environment seem inhabited to a greater extend than the other image in this pair of images?"). to an image, what need the ratios between the independent variables be? Here is a link a copy of the spreadsheet containing the data. https://docs.google.com/spreadsheet/ccc?key=0Am7BMwlOEiUNdGJhNU1uQlg5QUtxYW1ZNHhfUHlqS1E&usp=sharing


title: Probability distributions

text: What is the probability that a random number [0,1) will lie in [0.5,0.6]? As far as I know, there are 3 probability distributions in physics we usually use, but I am not sure when we can use Gaussian, Exponential and uniform distribution. Any examples?


title: Marginal effect clogit not significant

text: I developed a conditional logit model in Stata. The model is good and the variables are highly significant. Then I do mfx predicted (PU0) to determine the marginal effects of variables $\frac{dy}{dx}$. The problem is when I get the table of the marginal effects of variables in Stata: the significance of the variables becomes very small and insignificant. Why and therefore I can comment or not the marginal effects? Example coefficients of smoker 0.6666 and the significance P value 0.000 but when I put mfx compute to have $\frac{dy}{dx}$ smoker marginal effect is 1.2 and significance is $p= 0.521$.


title: model specification for SARIMA order (2,1,2)x(0,1,1) period 12

text: Good morning scholars. Please am fitting a Seasonal Arima model of this form: (2,1,2)x(0,1,1) period 12 but I don't know how will look like. Can anybody help me wit the model specification? Thanks.


title: On the use of the periodogram of a time series

text: I have a time series of hourly measurements for a duration of one year:               time <- 1:(365*24)/24     set.seed(1)     x <- rnorm(length(time))      The measurements are taken at hourly intervals. The main periodicity for the time series is that of a diurnal cycle i.e. the temperature follows a daily pattern. However, I would like to show that the 'power' of this expected frequency varies in time, i.e. being greater during some days in comparison to others. To achieve this I was thinking of running a window over the time series, computing a periodogram for each day (i.e. each 24 rows) showing that the periodicity of 1/24 explains less of the variation in the data during certain times of the year. I am new to R so have a rather basic knowledge of how to produce a script for this. Any advice would be appreciated.


title: Lag length selection Granger causality test

text: Consider G-Causality on two stationary time series vectors (call these variables $X$ and $Y$), each with 100+ observations. It's daily financial market time series data. I have reason to believe that there's reverse causality between these two variables (i.e., $X$ causes $Y$ as well as $Y$ causes $X$). I want to know the current consensus (or nearing consensus) regarding the best method of selecting lag length. A Google search has only revealed some 1984 and 1985 papers. Scholarpedia's entry on Granger causality says that it's O.K. to minimise AIC or BIC, but no references were provided for this claim so I don't want to code that up until I get confirmation. Or is selection based on qualitative reasoning? **Disclaimer:** Cross posted on talkstats.com.


title: What are good values to assume for correlation and sphericity when estimating sample size for a within-subjects design

text: I am planning a within-subjects study with four different conditions and I am using G*Power to calculate my needed sample size. Not being intimately acquainted with the intricacies of statistics in general (and repeated measures ANOVA in particular) what are some good assumptions for the _correlation among repeated measures_ and _nonsphericity correction_ values? What can I base these values on rather than guesswork?


title: How do I know what the slope of the line is by being provided Theta(1) (aka ø(1)) ? (Cost function and linear regression)

text: Andrew Ng is explaining linear regression and the cost function, and why we want to minimize Theta(1). However, he says that when Theta(1) = 1 then our slope of the line is (x,y) = (1,1), (2,2), (3,3). I don't understand how he computes that slope? I don't understand how to figure this out? How do you know that the slope of the line is that? http://cl.ly/image/1J1n2E151k2n For example, what happens when Theta(1) = 2, how would i figure out the slope? ==== For more detail, you can see: https://class.coursera.org/ml-2012-002/lecture/index, at 2:30 he mentions that if Theta(1) = 1, then the slope of the line is: (1,1), (2,2), (3,3), or http://cl.ly/image/1J1n2E151k2n


title: Calculating the Fisher Information of bivariate normal

text: I'm lost. If I estimated the a Gaussian mixture model, with a shared diagonal covariance, will the Fisher information of the means be $\Sigma^{-1}$ ?


title: Sample size that will give statistically significant results

text: I have a computer program which I am trying to optimize. Suppose that I can collect random test cases which I will use to run both the original and the modified version to see if there are any speedups. How can I calculate the sample size (number of test cases) which makes sure that the average speed-up/down I am getting is a significant one ? Any ideas is appreciated.


title: Multivariate normal distribution of regression coefficient?

text: While reading a textbook on regression I encountered the following paragraph: > The least squares estimate of a vector of linear regression coefficients > ($\beta$) is > > $$ \hat{\beta} = (X^{t}X)^{-1}{X^t}y $$ > > which, when viewed as a function of data $y$ (considering the predictors $X$ > as constants), is a linear combination of the data. Using the Central Limit > Theorem, it can be shown that the distribution of $\beta$ will be > approximately multivariate normal if the sample size is large. I'm definitely missing something from the text, but I don't understand how can a single $\beta$ value have a distribution? How are the multiple $\beta$ values generated to obtain the distribution referred to in the text?


title: Bootstrap vs other simulated data methods

text: In the mixed effect model, many statisticians would like to simulate or bootstrap data to create empirical confidence regions for fixed effect parameters and random effect parameters. Resampling (ie bootstapping) seems intuitive for me because it makes few assumptions about the nature of the data. As an alternative, some identify the multivariate distribution of a set of variables and draw at random from that distribution. My question is: Is there a principle where one would decide between one of these approaches? Is one of them always better?


title: how to estimate a generalized difference in difference model in panel data with stata

text: Please i need to know how to estimate the model bellow with stata: Log Xcit = α0 + α1 log PIBct + β0 Log liberalct + β1 log Liberalct×DepFini + β2 logLiberalct×Tangi + λc + λi + λt + εcit Liberalct is a binary variable equal to 1 in the year of and all years after an equity market liberalization event and 0 otherwise knowing that i have 4 countries that liberalised at different dates and 4 countries represent the control group.


title: How to interpret model diagnostics when doing linear regression in R?

text: I ran `lm()` on my data with models selected by individual `lm`'s of each characteristic and then combined the top $R^2$ based on $p$-value. For instance, the first few characteristics are taken, then the rest are evaluated if they have $p<.005$. My characteristics contain some duplication: for instance, I have a characteristic and its normalized variant in test P. My $p$-values are all very small but my diagrams do not look correct for R and T. (Referring to this blog post: Evaluating Linear Regression Model in R.) In test P (and T) there is one outlier according to Cooks Distance. How do I find and eliminate that instance? According to this tutorial on Using R for Linear Regression, > The plot in the upper left shows the residual errors plotted versus their > fitted values. The residuals should be randomly distributed around the > horizontal line representing a residual error of zero; that is, there should > not be a distinct trend in the distribution of points. Test P looks ok in the residual error but test R and T have a grouping what does that mean and how do I account for it? > The plot in the lower left is a standard Q-Q plot, which should suggest that > the residual errors are normally distributed. The scale-location plot in the > upper right shows the square root of the standardized residuals (sort of a > square root of relative error) as a function of the fitted values. Again, > there should be no obvious trend in this plot. Again Test P looks ok in the standard Q-Q plot but test R and T have a grouping what does that mean and how do I account for it? Also what is the coefficients on the output. I notice it lists the characteristics and a p value but i don't understand what it means. And finally how do I make predictions using the model I created? **Test P** F-statistic: 2.684 on 280 and 2221 DF, p-value: < 2.2e-16 ![enter image description here](http://i.stack.imgur.com/10A8r.png) **Test R** F-statistic: 3.691 on 258 and 2243 DF, p-value: < 2.2e-16 ![enter image description here](http://i.stack.imgur.com/jy6IR.png) **Test T** F-statistic: 4.029 on 268 and 2233 DF, p-value: < 2.2e-16 ![enter image description here](http://i.stack.imgur.com/bs69P.png) edit after running gls my p looks like this ![](http://i.stack.imgur.com/cUBGB.png)


title: When no model comparison, should I use REML vs ML?

text: I'm running LMM, and I will make no comparison of models. Could I ask which one should I use between REML and ML?


title: Condition number of covariance matrix

text: I am interested in generating a covariance matrix of dimension say 100. I managed to get a correlation matrix with finite condition number. To construct a covariance matrix I need to have standard deviations. I think for my case the most suitable one is to generate standard deviations from gamma distribution. So, it gives me small standard deviations as well as large standard deviations. As a result of that, the resulting covariance matrix has a very high condition number. I want to know whether the condition number can be affected by the scale of the variables and if I want to incorporate different scales in the covariance matrix how can I get a covariance matrix with a reasonable condition number? Any help or insight regarding this is highly appreciated.


title: Poisson Regression Nested Design

text: At the moment, I'm trying to do a 'Poisson Regression'. I have several years (2004 to 2012) with four seasons in each year. My dataset consists of the amount of carcasses found every 10 days, and the total carcasses for every season in every year. My goal is to see if there are any differences in carcasses between the seasons and between the years. Seasons is nested in years and I heard that Poisson regression is able to deal with nested datasets by including random effects. My question is how do include `Year` as a random effect in the R code. I used this code so far for the Poisson Regression:               mymodel = glm(counts~Season+Year, quasipoisson, mydataset)      I hope some could shed some light on this matter.


title: How to cross validate with optimal sample size selection?

text: My problem is as follows: I got a training set composed of texts, where 10% of observations are 1, and 90% are 0. To make it simpler, lets say that 1000 is 1, and 9000 is 0. I implement Naive Bayes and observe the results. Even before the Cross- Validation results depends strongly on the size of used observations - if I take 1000 of 1s, and about 3000 of 0, the results are the best. If I take all observations I have - 10000, the results are really poor. Now the question is: How to cross-validate with different amounts of 0 taken? If I take 1000 of 1s and 3000 of 0s, should I train it on 900 of 1s and 2700 of 0s and then test it on the rest 400 observations or rather test it on the whole of the rest - 100 1s and 6300 0s?


title: How genetic beta-diversity reflects traditional beta-diversity?

text: Im trying to figure out how betadisper {vegan, R}, and more specifically how the 'average distance to median' relates to betadiversity in a more traditional sense. betadisper calculates the variance found within communities/sites. I use distance matrices from sequence data where the distance corresponds to the raw base-pair difference. I have calculated betadiversity using the beta.multi (total=b.JAC or b.SOR) {betapart, R}. For calculating the traditional betadiversity I have clustered sequences into operational taxonomic units "OTU" (microbial substitute for species) and constructed a normal presence-absence matrix for sites-hosts. So, my first thought was that, a high value in betadiversity would correspond to very different averages from betadisper, but so is not necessarily the case. The average distances between sites for a host can show similarly low values although having a high betadiversity, and vice versa. I guess, similar averages of sites (for a host) as reported from betadisper doesn't mean the sites have to be compositionally similar (as betadisper don't look at the position of the centroid). Since, the base-pair distance Im using is directly related to which sequences cluster into which OTUs, I would imagine this being reflected in the "average distance to median" from betadisper. Also, could one aggregate the "average distance to medan" for the sites in order to derive at a single value (for each host) reflecting "traditional" betadiversity..? Not sure this makes sense, but any thoughts on this would be appreciated, thanks


title: Logistic regression

text: I want to categorize 20 responses from the BBC survey into average ability (chance level) and greater ability. I read somewhere that 25% above chance is no longer chance. Is it appropriate to use 12.5% above 50% (62.5%) as the category "greater ability" for logistic regression? Have I calculated this correctly? I am using the BBC survey "Spot the fake Smile" as the responses for ability. http://www.bbc.co.uk/science/humanbody/mind/index_surveys.shtml I have 20 predictor variables to create a composite variate of personality.


title: Why do we take the absolute value in a hypothesis test?

text: Recently I studied the definition of _p_ -value as > The _p_ -value associated with a test is the probability that we obtain the > observed value of the test statistic or a value that is more extreme in the > direction given by the alternative hypothesis, when $H_0$ (null hypothesis) > is true. For a two-sided alternative, the _p_ -value =$P_{H_0}[|T|\geq|t_0|]$ where $T$ is the test statistic and $t_0$ is the observed value of the test statistic. Why do we take the absolute values of $T$ and $t_0$?


title: Fisher's Exact: why no test statistic for 2x2 tables in SPSS

text: Whenever I have used the Fisher's Exact test on tables larger than 2x2 contingency tables, SPSS has produced a test statistic that I can quote in my work. However, when I do the same with 2x2 tables, no test statistic is produced. I want the way in which I report my results to appear standard throughout my thesis. Therefore do I (a) omit the test statistic that is provided for the larger tables or (b) find some way of calculating the test statistic? A similar issue was noted in a question entitled "Fisher's Exact test value" but the answers didn't address this issue specifically.


title: Graphic representation of data

text: I conducted exploratory analysis of my data and got Mean and standard deviation values from the original data and plotted the data using boxplots. I conducted regression models (GLM and GLMM) on the data (which allow to estimate marginal means) and I would like to know what I have to report, if the "raw" mean and SD or the values estimated by the model. Or if I should report both, for example the "raw" in a table and the estimated in a boxplot or vice-versa. Thank you!


title: Equivalence between Elastic Net formulations

text: According to Hastie's paper, the elastic net has two equivalent formulations: $\hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\\{ \sum_{i=1}^N(y_i-\sum_{j=1}^p x_{ij} \beta_j)^2 + \lambda_2 \sum_{j=1}^p \beta_j^2 + \lambda_1 \sum_{j=1}^p |\beta_j| \right\\}$ and $\hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\\{ \sum_{i=1}^N(y_i - \sum_{j=1}^p x_{ij} \beta_j)^2\right\\} \text{ s.t. } (1-\alpha)\sum_{j=1}^p |\beta_j| + \alpha\sum_{j=1}^p \beta_j^2$ where $\alpha = \lambda_2/(\lambda_1 + \lambda_2)$ My question is how to prove this equivalence formally. Ridge regression and the lasso also have these two possible formulations, but I could not find any reference where this equivalence is proven. A similar question I found in CrossValidated is this one Lagrangian relaxation in the context of ridge regression but I'm unable to understand Tristan's explanation. I have some understanding of Lagrange optimization theory, and I guess the answer is around those lines, but since all the papers treat the equivalence as obvious I would like to find a proper reference where this is explicitly demonstrated.


title: Is the notion of 'no-effect' flawed?

text: I am a student of mathematics, but my knowledge of probability and statistics is rather limited. So please forgive if this is a bad question. T-tests are a popular way to perform statistical inference. However the null- hypothesis of 'no-effect' seems flawed. Take any two variables: suppose they're 'amount of time person spent on cross-validated in 2012' and 'extent to which person likes ice-cream in 2013' and you'll most likely see an effect, as long as your sample size is enormous. At worst, you can regard the whole population as a single sample, in which case there is no uncertainty:, thus if you see a positive relationship you can reject the 'no effect' hypothesis and if you see a negative relationship you can also reject the 'no effect' hypothesis. So either way, the null-hypothesis is rejected. _Therefore, we can reject the null-hypothesis without collecting any data._ Is the notion of 'no effect' fundamentally flawed, or have I made a logical error? Of course you may say, "that's why we also report effect sizes," but that's kind of missing the point. If the concept of 'no effect' is flawed this badly, then we should forget this 'no-effect' business and use more sensible hypotheses.


title: Analysing the correlation of two ranked scales

text: So I am researching if self-esteem has an effect on hearing aid satisfaction in patients. The patients have filled out two questionnaires, the first measuring self esteem, the second measuring their hearing aid satisfaction. My hypothesis is that a higher self esteem score will correlate to a higher hearing aid satisfaction. The self esteem questionnaire is a ten part likert scale which asks the patients to ring whether they strongly agree, agree, disagree or strongly disagree with a statement. The hearing aid satisfaction questionnaire is similar but uses tick boxes for the patients to choose which option best describes their use. So 20 people have done the questionnaires and so I have 20 scores for self esteem and hearing aid satisfaction for each patient. I want to see if the two scores are correlated. What statistical test do I use? (I literally have no experience in stats at all)


title: How can I handle the non-response in this simple random sampling?

text: I am planning to conduct a simple random sampling for a population of $N$ patients with index $1,...N$. I have enough money to interview $n$ patients (sample size=$n$). In order to determine the sample patients, I first generate $n$ random integers(with replacement) and ask the patients with such index whether they would like to join the interview. However, due to some reasons, several patients (non-responder) were not able to join. My question is: what shall I do if in this case I have patients who can not join, but I also want to maximize the use of $n$ samples? I can think of: If 3 patients can not join, re-generate 3 random integers and assign three new patients. Repeat this process until all patients are willing to join. However, I will still record information from the 3 non-responders, and the total sample size is $n+3$. The non-responder's information will be analyzed later to see if there is any bias. Is this method reasonable?


title: Interpolation of time series data

text: I am working on using cubic spline interpolation in time series data. I used Galdolfo and Prachowney algorithms. Now how do I obtain estimates of models of cubic spline?


title: Finding a statistically different proportion

text: I have the following problem: given a proportion $p_1=\frac{X_1}{n}$, where $X_1$ is the number of succesful outcomes in $n$ tries ($X_1$ and $n$ are both > 100), I have to find the smallest proportion $p_2$ that is greater than and statistically different from $p_1$ at a given confidence level (e.g. 95%). An application of this would be to find out a success rate (i.e., $X_2$ successes in $n$ tries) that would ensure that a given enhancement to a system would actually be beneficial. I could find the confidence interval for $p_1$ but merely having $p_2$ outside this interval wouldn't be enough to assure these proportions would be different. Doubling the interval didn't seem right either. I also tried using the formulas for the confidence interval for the difference of proportions $c = p\mp z\cdot s_p$ where $p = p_1-p_2$ is the difference between the proportions and $s_p$ is given by $s_p = \sqrt{\frac{p_1(1-p_1)+p_2(1-p_2)}{n}}$ but I couldn't figure out how to get $p_2$ from them. Any hints on how I might proceed? For what it's worth, this is not homework, but a slightly reworded problem from a book I'm using for self-study. EDIT: Glen's answer put me in the right direction. It also occurred to me to brute force a solution, posted below.


title: Obtaining and sampling from the posterior predictive of a naive Bayes classifier

text: I have trained a naive Bayes classifier with on a dataset with a dichotomous outcome and multinomial attributes (predictors). I managed to get a Maximum _a posteriori_ (MAP) estimate which is good enough to determine class membership of a given instance of attributes. However, to calculate things like the Highest Density Interval (HDI) I would need to get the posterior distribution rather than a point estimate. How can I do that? And more specifically, how could I sample from that distribution? (via BUGS for instance). Any help is appreciated


title: Is conditional logit a specific form of GLM? And what are its specificities?

text: Background: For a project, I am fitting a conditional logit model where I have 5 control cases for every realized case. To do that I use the `clogit()` function in the package `survival`. I wanted to graph interactions with the `effects` package by John Fox et al. It turns out that this package can't handle `clogit` objects (output of `clogit()`). As I believed I remembered that conditional logit were a special case of GLM, I thought the clever/lazy way to get my interaction plots would be to refit the model using a fixed effects glm and then use `effect()`. The documentation of `clogit` seemed to confirm my intuition: > It turns out that the logliklihood for a conditional logistic regresson > model = loglik from a Cox model with a particular data structure. [...] When > a well tested Cox model routine is available many packages use this ‘trick’ > rather than writing a new software routine from scratch, and this is what > the clogit routine does. > > In detail, a stratified Cox model with each case/control group assigned to > its own stratum, time set to a constant, status of 1=case 0=control, and > using the exact partial likelihood has the same likelihood formula as a > conditional logistic regression. The clogit routine creates the necessary > dummy variable of times (all 1) and the strata, then calls coxph. Based on this description, it seems that I should be able to reproduce the stratification achieved through `strata()` by using a random intercept for each case/control group with `1|group` in `lmer()`. However, when I try, the results of `clogit` and `lmer` differ. One thing is that I probably have the wrong likelihood function. I don't really know how to specify this in `lmer` but more important, I am wondering what else I am missing. I wonder whether I am completely wrong or somewhat on the right track but missing some pieces? What I would like is to understand what are the difference in terms of how the model is fitted between a conditional logit and a regular one (I understand that might be quite a long answer, so a book reference would be a great start). The my usual references for regression (Gelman and Hill, 2007; Mills 2011) are somewhat silent on the subject.


title: How I calculate the confidence interval for nonparametric count data? (Is it possible in GraphPad Prism?)

text: hope you can help me with my problem. My results are counts. What i did is identifying if embryos had (yes = 1) or not (no = 0) malformations or delay in their development, within samples of 15 embryos. In one type of experiment, I repeated this 6 times (6 independent replicates); in another type, I repeated this 9 times (9 independent replicates), so my n is quite low. I observed a lot of features (between 6 and 12) in every embryo, every 24 hours for 5 days (24 hpf, 48 hpf, 72 hpf and 96 hpf) so I have a lot of data. My first idea was to do a "repeated measures ANOVA", but unfortunately, my data are absolutely NOT normal distributed; because they're **"binary" counts** (with minimal value 0 and maximum value 15), and I have **lots of ties**. I tried to use Friedman's test, but it really takes too long! (At least 30 minutes for every feature). I tried a lot of ways to normalize the data (even one I had never heard before: using the arcsin of the percentage square root...!), but nothing works. Then a lab mate told me that I could calculate the confidence interval for my controls, and take as "different" from that any value that is outside the CI. I have the GraphPad Prism program, and it says that for calculating CI, you assume independet and gaussian distributed values... So I got stuck again at the beginning. How can I calculate a non-parametric CI, using Prism? On the other hand, I was thinking about using bootstrap, which theory I understand very basically, but I have never done it, and don't know if I can do it with this program. Any help would be extremely appreciated!


title: If I'm not interested in the interaction, is there any reason to run a Two-Way ANOVA instead of two One-Way ANOVAs?

text: I mean any reason aside from the convenience of being able to complete the analysis within a single procedure.


title: Existence of co-variables

text: In a study, there are 3 independent variables (A, B, C). I want to manipulate only (A, B). The idea of not manipulating variable C is for it not to interfere with the other two manipulated variables A, B. In this way, the researcher is confident that only A or/and B are causing an effect on the dependent variable. (i) I am not sure how to describe the above in statistical language. How should I best put it? (ii) In this case, variable C is not regarded as a co-variate, as it is kept constant. Please advise. Thanks


title: Confusion on Maximum Likelihood's starting setup equation

text: This is my first post on CV. I am wrapping my head around Maximum Likelihood Estimation (MLE), and ironically, I understand everything about the algorithm, except the starting point. What confuses me is the following: All the examples I have seen of course, start with some vector of observations, x1, x[2], ... x[N]. The observations are iid. The objective is to estimate parameters (theta) of the PDF - ONE PDF - from which sired such observations. This part I get. We are allowed to estimate **ONE PDF** because X ('big X') is a random variable from which all those x's (small x's) are sired. As far as I know, big X's have PDFs. Small x's dont. Anyway, so the PDF whos theta we want to estimate is written like this: p(x1, x[2], x[3]... x[n] | theta) = p(x1|theta) * p(x[2]|theta) * ... p(x[N]|theta). Like I said if I accept this I know the rest of the algorithm, but this starting point is where I dont understand. Here all the sudden I now have **N PDFs**. This confuses me a lot: Like I said before, as far as I know, observations (like x[5]) dont have PDFs. Random variables do. So in this case the observations are 'small x', where as a random variable is 'big X'. We are trying to estimate the pdf of 'big X' - X. How is it that I now have an individual PDFfor each and every observation 'small x'? Every lecture video/article seem to use contradictory terms among themselves of 'random variable' and 'observation' while simultaneously taking the effort to point out that they are not the same and shouldnt be confused - so which one is it?? Thanks in advance! **EDIT for whuber:** For additional context, here is my current state of understanding: By the way, I am going to assume a Gaussian process for X from here on out for this example.   * A random variable, usually denoted 'big something' - X in this case, can have a PDF. A 'small something' is the one outcome - one observation of big X, we call this 'small x'.    * So therefore, we can ask "What is the probability of X = 4?". Small x = 4 here, and this has meaning, and we can write "P(X=4) = ?". Big X's PDF evaluated at a chosen small x's value, gives us the probability X come out as x. In other words, P(X=x).   * To me this means only X can have a PDF. (X represents a name for a process that can generate many possible values - he can take many values), while x, cannot have a PDF. (What is the PDF of the number 7?). This is how I see it.    * If I have one observation, (x = 7), then I understand that the jist of the MLE is to say, which guassian, with what mu and sigma, maximize the probability of this one number being generated from it? If I had all the time in the universe, I would: 1) Take a guassian PDF with some mu and some sigma, and then 2) plug x = 7 into it and get a probability of it being 7. Eventually I would realize that the best mu for my one x = 7 is in fact when mu = 7.    * If I have two observations, (x1 = 7, x2 = 9), then my question in english is, "What possible combo of mu and sigma did my X have, such that I could have gotten x1 = 7 and x2 = 9 with the highest odds?". To me this means it is nonsensical to assume that x1 and x2 were aberrations and on the tail end of a gaussian. So not knowing anything else, we assume that well, if they are here, they were the most likely. So what values of mu and sigma most likely led to those being generated?    * This now means ok, what is the probability of x1=7 AND x2=9? So P(x1 = 7 AND x2 = 9) = P(x1=7) * P(x2=9). We can do this because they are IID. **EDIT 2** I have attached an image of what I believe to me a Likelihood function for ONE observation, (x = 3), assuming it came from a gaussian X, based on my new understanding from the convos here. Interesting...looks like the light pattern out of the headlight of a car... ![enter image description here](http://i.stack.imgur.com/WiA9Y.jpg)


title: Singularity issues in multinomial logit model with differing choice sets

text: I am estimating a discrete choice model in which individuals choose which schools to attend. I have a large amount of data on individuals and schools. However, each particular school only appears in my data set one or two times since individuals face different sets of school choices based on where they live, etc. For example, individual A may choose between schools 1, 2, and 3 while individual B chooses between 4, 5, and 6. These choices should provide information on the importance of school attributes, but I cannot determine how to specify this model correctly. I have followed the advice in this post and ensured that all individual- specific characteristics appear after the |. There is no missing data in my sample. Unfortunately I cannot post my data directly given its somewhat sensitive nature. I appreciate any advice!


title: How do I create a random variance matrix (x'x)/N with confidence that it will be invertible?

text: I see some similar threads, but no answer that I know how to apply to this question. I wish to create a variance matrix for a (pseudo)random matrix $X$ with dimensions $K \times N$, where each row $k$ represents $N$ draws of a random variable with mean $\mu_k$, $\sigma^2_k=10$, and 10 was arbitrarily selected. The goal is to create $X$ such that $(X'X)/N$ is computationally invertible. I'm using R... so that defines the bounds of "computationally invertible". So far, I've done the following, which I assumed would be sufficient. The last line points out to me that I am probably confused about how easy it is to acquire linear dependence.               N<-50     K<-round(runif(1,1,10))     ex<-rnorm(1,0,10)     x<-matrix(rnorm(N,ex,10)) #one row is already created.     tick<-0 #cheap trick to make sure your loop repeated successfully.     for(i in 2:K){     ex<-rnorm(1,0,10) #the means should each be different, should help?     holding<-as.matrix(rnorm(N,ex,1)) #it's an Nx1, k=1     x<-cbind(x,holding)     tick<-tick+1     }     print(tick)     x<-t(x) #x is now K by N     det(cov(x)) #It's computationally 0 with great consistency.      Primary question: How do I make one such random matrix $X$ such that it does have an invertible variance matrix? Aside: Could someone explain why this is so frequently singular, since I thought by LLN each column/row should converge to their uncorrelated state of linear independence? Is the "computationally" singular restriction so powerful? I've noticed that it's actually working sometimes if I set $N$ as arbitrarily small, such as 5.


title: Why will a statistic be significant with sufficiently large samples unless the population effect is exactly zero?

text: From Wikipedia > **Given a sufficiently large sample size, a statistical comparison will > always show a significant difference unless the population effect size is > exactly zero.** > > For example, a sample Pearson correlation coefficient of 0.1 is strongly > statistically significant if the sample size is 1000. Reporting only the > significant p-value from this analysis could be misleading if a correlation > of 0.1 is too small to be of interest in a particular application. I was wondering why "given a sufficiently large sample size, a statistical comparison will always show a significant difference unless the population effect size is exactly zero"? Thanks and regards!


title: Items from two constructs load on one factor in factor analysis

text: I am trying to determine the factor structure of a set of 84 items. Exploratory factor analysis using varimax rotation was conducted to estimate the underlying factor structure for the sample data. Two rounds of exploratory factor analysis were conducted. Items with loading below the generally accepted 0.60 threshold and those that highly cross-loaded on more than one factor were rejected. However, the overall results of the exploratory factor analysis shows that items for two different constructs load on a same factor. ### Questions   * **How should I interpret the fact that items for two constructs load on the one factor?**   * **Is this a problem?**   * **If yes, what could I do?**


title: What does the Scale parameter mean in linear regression?

text: I have started to use the `GENLIN` procedure in SPSS more than any of the specific dialogues, but I don't understand the Scale parameter or why it has the effects it does on the regression results. Here's an SPSS example: _code block 1_ Normally, if I wanted to see if _set_ has an effect on the linear parameters $b_0$ and $b_1$, I'd could just do linear regression with the interaction term _set1X_ : _code block 2_ The results are obvious: _set_ affects slope but not the intercept of my model (as I'd expect since I generated these sample data that way). But if I use `GENLIN` instead of `REGRESSION`, my sig. values and 95%CIs for the parameters are different: _code block 3_ This Scale parameter (regressed by `GENLIN` as 1.127 with a SE of 0.3563) seems to make the difference. If I change the scale to `PEARSON` (for Pearson's chi sq.), `DEVIANCE` (?), or a constant like 1, I get all different answers. **Summary** : what is Scale in Generalized Linear Models (`GENLIN` in SPSS) and how should I handle it? Why doesn't OLS Linear regression use such a parameter? How do I know how to assign scale? ## Code block 1               data list list /X set Y.     begin data.     1   1   7.0     2   1   7.8     3   1   12.4     4   1   14.8     5   1   19.0     6   1   22.7     7   1   24.4     8   1   25.5     9   1   29.5     10  1   31.0     1   2   7.9     2   2   12.7     3   2   14.3     4   2   20.1     5   2   20.8     6   2   26.5     7   2   30.9     8   2   35.8     9   2   38.0     10  2   43.7     end data.     dataset name exampleData WINDOW=front.     variable level X (scale) Y (scale) set (nominal).          compute set1 = (set=1).     compute set1X = set1*X.     EXECUTE.      ## Code block 2               REGRESSION       /MISSING LISTWISE       /STATISTICS COEFF OUTS CI(95) R ANOVA CHANGE       /CRITERIA=PIN(.05) POUT(.10)       /NOORIGIN        /DEPENDENT Y       /METHOD=FORWARD X set1 set1X.      ## Code block 3               * Generalized Linear Models.     GENLIN Y BY set (ORDER=ASCENDING) WITH X       /MODEL X set*X INTERCEPT=YES      DISTRIBUTION=NORMAL LINK=IDENTITY       /CRITERIA SCALE=MLE COVB=MODEL PCONVERGE=1E-006(ABSOLUTE) SINGULAR=1E-012 ANALYSISTYPE=3(WALD)          CILEVEL=95 CITYPE=WALD LIKELIHOOD=FULL       /MISSING CLASSMISSING=EXCLUDE       /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION.


title: Assessing or comparing correlations at multiple locations

text: I have 10 independent locations (separated by at least 5 km) where I have collected abundance data for 17 species of birds. In particular, I want to look for correlations between one particular species and the 16 other species. Each site was sampled multiple times (separated by 2-8 weeks), some only 17 times while others were sampled 50+ times. I know I can perform Spearman correlations for each location (most data are not normally distributed) and calculate r for each pair. However, as might be expected, not every location shares the same significant relationships and some are even conflicting (some + and some -). Some significant results seem rather spurious (significant based on only 2 observations of one species). Is there some way I can compile all these different correlation results into a single test to determine whether there is a general trend for the one species to correlate with others? Am I doing this all wrong? I'm floundering and would welcome any advice.


title: Monte Carlo method

text: I need to write a report on what is Monte Carlo method suitable for and for what it is not suitable for, however, I can't find a lot of material regarding this information. What is it suitable for apart for integration, physics and travelling salesman problem? And for what it is not?